File content as of revision 1:5ec4a7674069:

#include "mbed.h"
#define ITMAX 100
#define CGOLD 0.3819660
#define GOLD 1.618034
#define GLIMIT 100.0
#define TINY 1.0e-20
#define SHIFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
#define ZEPS 1.0e-10
Serial pc(USBTX,USBRX);
double f(double x);
double brent(double min,double mid,double max,double tol);
double SIGN(double x,double y);
void mnbrak();//引数は全てグローバルで指定するとする
double FMAX(double x,double y);

double ax=0.0,bx=1.0,cx=0.0;
double fa,fb,fc;

void mnbrak()
{
    double ulim,u,r,q,fu,dum,fa,fb,fc;
    
    fa=f(ax);
    fb=f(bx);
    pc.printf("fa=%lf,fb=%lf\r\n",fa,fb);
    if(fb>fa)
    {
        SHIFT(dum,ax,bx,dum);
        SHIFT(dum,fb,fa,dum);
        }
    cx=bx+GOLD*(bx-ax);
    fc=f(cx);
    while (fb>fc)
    {
        r=(bx-ax)*(fb-fc);
        q=(bx-cx)*(fb-fa);
        u=bx-((bx-cx)*q-(bx-cx)*r)/
        (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));
        ulim=bx+GLIMIT*(cx-bx);
        
        if((bx-u)*(u-cx)>0.0)
        {
            fu=f(u);
            if(fu<fc)
            {
                ax=bx;
                bx=u;
                fa=fb;
                fb=fu;
                return;
                }
            else if(fu>fb)
            {
                cx=u;
                fc=fu;
                return;
                }
             u=cx*+GOLD*(cx-bx);
             fu=f(u);       
            }
        else if((cx-u)*(u-ulim)>0.0)
        {
            fu=f(u);
            if(fu<fc)
            {
                SHIFT(bx,cx,u,cx+GOLD*(cx-bx));
                SHIFT(fb,fc,fu,f(u));
                }
            
            }
        else if((u-ulim)*(ulim-cx)>=0.0)
        {
            u=ulim;
            fu=f(u);
            }
        else
        {
            u=cx+GOLD*(cx-bx);
            fu=f(u);
            }
        SHIFT(ax,bx,cx,u);
        SHIFT(fa,fb,fc,fu);
        }
    }
double brent(double ax,double bx,double cx,double tol)
{
    pc.printf("bernt start");
    int iter;
    double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm,xmin;
    double e=0.0;
    a=(ax <cx ? ax : cx);
    b=(ax >cx ? ax : cx);
    x=w=v=bx;
    fw=fv=fx=f(x);
    for(iter=1;iter<=ITMAX;iter++)
    {
        //pc.printf("bernt loop");
        xm=0.5*(a+b);
        tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
        pc.printf("x =%lf,w = %lf,u = %lf\n\r",x,w,u);
        if(fabs(x-xm)<=(tol2-0.5*(b-a)))
        {
            //pc.printf("bernt out");
            xmin=x;
             pc.printf("xmin=%lf\r\n",x);
             pc.printf("fx=%lf\r\n",fx);
            return xmin;
        }
        if(fabs(e)>tol1)
        {
            r=(x-w)*(fx-fv);
            q=(x-v)*(fx-fw);
            p=(x-v)*q-(x-w)*r;
            q=2.0*(q-r);
            if(q>0.0)p=-p;
            q=fabs(q);
            etemp=e;
            e=d;
            if(fabs(p)>=fabs(0.5*q*etemp)||p<=q*(a-x)||p>=q*(b-x))
            {  d=CGOLD*(e= (x>=xm ? a-x : b-x));}
            else
            {
               d=p/q;
               u=x+d;
               if(u-a < tol2 || b-u < tol2)
               {d=SIGN(tol1,xm-x);}
            }
        }
        else
        {
            d=CGOLD*(e= (x>=xm ? a-x : b-x));
            }    
        u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
        fu=f(u);///////注意
         if(fu <= fx)
         {
            if(u >= x)a=x; else b=x;
            SHIFT(v,w,x,u);
            SHIFT(fv,fw,fx,fu);
          }
         else{
          if(u < x){a=u;} 
          else {b=u;}
          if(fu <= fw || w==x)
          {
              v=w;
              w=u;
              fv=fw;
              fw=fu;
              }    
          else if (fu <= fv || v==x || v==w)
          {
            v=u;
            fv=fu;
          }
        }
        
        }
        pc.printf("xmin=%lf\r\n",x);
        pc.printf("fx=%lf\r\n",fx);
    return xmin;
    }
//極小値を求めたい関数を定義    
double f(double x){
    double x_return;
    x_return=x*x*x-x;
    return x_return;
    }   
double SIGN(double x,double y)
{
    double x_return;
    x_return=abs(x);
    if(y<0.0)x_return=-x_return;
    printf("f");
    return x_return; 
    } 
double FMAX(double x, double y)
{
    if(x>y){return x;}
    if(y>x){return y;}
    }        
int main() {
   pc.baud(921600);
   pc.printf("loop start\r\n");
   mnbrak();
   pc.printf("(a,b,c)=(%lf,%lf,%lf)\r\n",ax,bx,cx);
   pc.printf("(fa,fb,fc)=(%lf,%lf,%lf)\r\n",fa,fb,fc);
   brent(ax,bx,cx,0.000001);
   return 0;
}