dbrent.pas

来自「Delphi Pascal 数据挖掘领域算法包 数值算法大全」· PAS 代码 · 共 112 行

FUNCTION dbrent(ax,bx,cx,tol: real; VAR xmin: real): real;
(* Programs using routine DBRENT must define the external
functions func(x:real):real and dfunc(x:real):real which are,
respectively, the function whose minimum is to be found,
and its derivative. *)
LABEL 1,2,3;
CONST
   itmax=100;
   zeps=1.0e-10;
VAR
   a,b,d,d1,d2: real;
   du,dv,dw,dx: real;
   e,fu,fv,fw,fx: real;
   iter: integer;
   olde,tol1,tol2: real;
   u,u1,u2,v,w,x,xm: real;
   ok1,ok2: boolean;
FUNCTION sign(a,b: real): real;
   BEGIN
      IF (b > 0.0) THEN sign := abs(a) ELSE sign := -abs(a)
   END;
BEGIN
   IF ax < cx THEN a := ax ELSE a := cx;
   IF ax > cx THEN b := ax ELSE b := cx;
   v := bx;
   w := v;
   x := v;
   e := 0.0;
   fx := func(x);
   fv := fx;
   fw := fx;
   dx := dfunc(x);
   dv := dx;
   dw := dx;
   FOR iter := 1 TO itmax DO BEGIN
      xm := 0.5*(a+b);
      tol1 := tol*abs(x)+zeps;
      tol2 := 2.0*tol1;
      IF (abs(x-xm) <= (tol2-0.5*(b-a))) THEN GOTO 3;
      IF (abs(e) > tol1) THEN BEGIN
         d1 := 2.0*(b-a);
         d2 := d1;
         IF (dw <> dx) THEN  d1 := (w-x)*dx/(dx-dw);
         IF (dv <> dx) THEN  d2 := (v-x)*dx/(dx-dv);
         u1 := x+d1;
         u2 := x+d2;
         ok1 := ((a-u1)*(u1-b) > 0.0) AND (dx*d1 <= 0.0);
         ok2 := ((a-u2)*(u2-b) > 0.0) AND (dx*d2 <= 0.0);
         olde := e;
         e := d;
         IF (NOT (ok1 OR ok2)) THEN GOTO 1
         ELSE IF (ok1 AND ok2) THEN BEGIN
            IF (abs(d1) < abs(d2)) THEN BEGIN
               d := d1
            END ELSE BEGIN
               d := d2
            END
         END ELSE IF (ok1) THEN BEGIN
            d := d1
         END ELSE BEGIN
            d := d2
         END;
         IF (abs(d) > abs(0.5*olde)) THEN GOTO 1;
         u := x+d;
         IF (((u-a) < tol2) OR ((b-u) < tol2)) THEN BEGIN
            d := sign(tol1,xm-x)
         END;
         GOTO 2
      END;
1:      IF (dx >= 0.0) THEN e := a-x ELSE e := b-x;
      d := 0.5*e;
2:      IF (abs(d) >= tol1) THEN BEGIN
         u := x+d;
         fu := func(u)
      END ELSE BEGIN
         u := x+sign(tol1,d);
         fu := func(u);
         IF (fu > fx) THEN GOTO 3
      END;
      du := dfunc(u);
      IF (fu <= fx) THEN BEGIN
         IF (u >= x) THEN a := x ELSE b := x;
         v := w;
         fv := fw;
         dv := dw;
         w := x;
         fw := fx;
         dw := dx;
         x := u;
         fx := fu;
         dx := du
      END ELSE BEGIN
         IF (u < x) THEN a := u ELSE b := u;
         IF ((fu <= fw) OR (w = x)) THEN BEGIN
            v := w;
            fv := fw;
            dv := dw;
            w := u;
            fw := fu;
            dw := du
         END ELSE IF ((fu < fv) OR (v = x) OR (v = w)) THEN BEGIN
            v := u;
            fv := fu;
            dv := du
         END
      END
   END;
   writeln('pause in routine DBRENT - too many iterations');
3:   xmin := x;
   dbrent := fx
END;