PROGRAM badluk(input,output);
LABEL 1,2;
CONST
   zon=-5.0;
   iybeg=1900;
   iyend=2000;
VAR
   timzon,frac: real;
   ic,icon,idwk,im: integer;
   iyyy,jd,jday,n: integer;
(*$I MODFILE.PAS *)
(*$I JULDAY.PAS *)
(*$I FLMOON.PAS *)
BEGIN
   timzon := zon/24.0;
   writeln('Full moons on Friday the 13th from',iybeg:5,' to',iyend:5);
   FOR iyyy := iybeg to iyend DO BEGIN
      FOR im := 1 to 12 DO BEGIN
         jday := julday(im,13,iyyy);
         idwk := (jday+1) MOD 7;
         IF (idwk = 5)  THEN BEGIN
            n := trunc(12.37*(iyyy-1900+
               trunc((im-0.5)/12.0)));
            icon := 0;
1:            flmoon(n,2,jd,frac);
            frac := 24.0*(frac+timzon);
            IF (frac < 0.0) THEN BEGIN
               jd := jd-1;
               frac := frac+24.0
            END;
            IF (frac > 12) THEN BEGIN
               jd := jd+1;
               frac := frac-12.0
            END ELSE BEGIN
               frac := frac+12.0
            END;
            IF (jd = jday) THEN BEGIN
               writeln;
               writeln(im:2,'/',13:2,'/',iyyy:4);
               writeln('Full moon ',frac:5:1,
                  ' hrs after midnight (EST).');
               GOTO 2 END
            ELSE BEGIN
               IF (jday >= jd) THEN ic := 1 ELSE ic := -1;
               IF (ic = -icon) THEN GOTO 2;
               icon := ic;
               n := n+ic
            END;
            GOTO 1;
2:         END
      END
   END
END.
PROCEDURE balanc(VAR a: glnpnp; n: integer);
(* Programs using routine BALANC should define the type
TYPE
   glnpnp = ARRAY [1..np,1..np] OF real;
where 'np by np' is the physical dimension of the array to be analyzed. *)
CONST
   radix=2.0;
VAR
   last,j,i: integer;
   s,r,g,f,c,sqrdx: real;
BEGIN
   sqrdx := sqr(radix);
   REPEAT
      last := 1;
      FOR i := 1 to n DO BEGIN
         c := 0.0;
         r := 0.0;
         FOR j := 1 to n DO
            IF (j <> i) THEN BEGIN
               c := c+abs(a[j,i]);
               r := r+abs(a[i,j])
            END;
         IF ((c <> 0.0) AND (r <> 0.0)) THEN BEGIN
            g := r/radix;
            f := 1.0;
            s := c+r;
            WHILE (c < g) DO BEGIN
               f := f*radix;
               c := c*sqrdx
            END;
            g := r*radix;
            WHILE (c > g) DO BEGIN
               f := f/radix;
               c := c/sqrdx
            END;
            IF ((c+r)/f < 0.95*s) THEN BEGIN
               last := 0;
               g := 1.0/f;
               FOR j := 1 to n DO a[i,j] := a[i,j]*g;
               FOR j := 1 to n DO a[j,i] := a[j,i]*f
            END
         END
      END;
   UNTIL (last <> 0)
END;
PROCEDURE bcucof(y,y1,y2,y12: gl4array; d1,d2: real; VAR c: gl4by4);
(* Programs using routine BCUCOF must define the types
TYPE
   gl4array = ARRAY [1..4] OF real;
   gl4by4 = ARRAY [1..4,1..4] OF real;
in the main routine. They must also declare the variables
VAR
   glflag: boolean;
   wt: ARRAY [1..16,1..16] OF real;
and initialize glflag to true. The values of wt are read from the file
bcucof.dat whose contents are listed at the end of this routine. The procedure
GLOPEN assigns bcucof.dat to infile and opens the file for reading. *)
VAR
   l,k,j,i: integer;
   xx,d1d2: real;
   cl,x: ARRAY[1..16] OF real;
   infile: text;
BEGIN
   IF glflag THEN BEGIN
      glflag := FALSE;
      glopen(infile,'bcucof.dat');
      FOR i := 1 to 16 DO FOR k := 1 to 16 DO read(infile,wt[k,i]);
      close(infile)
   END;
   d1d2 := d1*d2;
   FOR i := 1 to 4 DO BEGIN
      x[i] := y[i];
      x[i+4] := y1[i]*d1;
      x[i+8] := y2[i]*d2;
      x[i+12] := y12[i]*d1d2
   END;
   FOR i := 1 to 16 DO BEGIN
      xx := 0.0;
      FOR k := 1 to 16 DO xx := xx+wt[i,k]*x[k];
      cl[i] := xx
   END;
   l := 0;
   FOR i := 1 to 4 DO
      FOR j := 1 to 4 DO BEGIN
          l := l+1;
          c[i,j] := cl[l]
      END
END;
(* Contents of the file bcucof.dat
1 0 -3 2 0 0 0 0 -3 0 9 -6 2 0 -6 4 0 0 0 0 0 0 0 0 3 0 -9 6 -2 0 6 -4
0 0 0 0 0 0 0 0 0 0 9 -6 0 0 -6 4 0 0 3 -2 0 0 0 0 0 0 -9 6 0 0 6 -4
0 0 0 0 1 0 -3 2 -2 0 6 -4 1 0 -3 2 0 0 0 0 0 0 0 0 -1 0 3 -2 1 0 -3 2
0 0 0 0 0 0 0 0 0 0 -3 2 0 0 3 -2 0 0 0 0 0 0 3 -2 0 0 -6 4 0 0 3 -2
0 1 -2 1 0 0 0 0 0 -3 6 -3 0 2 -4 2 0 0 0 0 0 0 0 0 0 3 -6 3 0 -2 4 -2
0 0 0 0 0 0 0 0 0 0 -3 3 0 0 2 -2 0 0 -1 1 0 0 0 0 0 0 3 -3 0 0 -2 2
0 0 0 0 0 1 -2 1 0 -2 4 -2 0 1 -2 1 0 0 0 0 0 0 0 0 0 -1 2 -1 0 1 -2 1
0 0 0 0 0 0 0 0 0 0 1 -1 0 0 -1 1 0 0 0 0 0 0 -1 1 0 0 2 -2 0 0 -1 1 *)
PROCEDURE bcuint(y,y1,y2,y12: gl4array; x1l,x1u,x2l,x2u,x1,x2: real;
        VAR ansy,ansy1,ansy2: real);
(* Programs using procedure BCUINT must define the types
TYPE
   gl4array = ARRAY [1..4] OF real;
   gl4by4 = ARRAY [1..4,1..4] OF real;
in the main routine. *)
VAR
   i: integer;
   t,u,d1,d2: real;
   c: gl4by4;
BEGIN
   d1 := x1u-x1l;
   d2 := x2u-x2l;
   bcucof(y,y1,y2,y12,d1,d2,c);
   IF ((x1u = x1l) OR (x2u = x2l)) THEN BEGIN
      writeln('pause in routine BCUINT - bad input'); readln
   END;
   t := (x1-x1l)/d1;
   u := (x2-x2l)/d2;
   ansy := 0.0;
   ansy2 := 0.0;
   ansy1 := 0.0;
   FOR i := 4 DOWNTO 1 DO BEGIN
      ansy := t*ansy+((c[i,4]*u+c[i,3])*u+c[i,2])*u+c[i,1];
      ansy2 := t*ansy2+(3.0*c[i,4]*u+2.0*c[i,3])*u+c[i,2];
      ansy1 := u*ansy1+(3.0*c[4,i]*t+2.0*c[3,i])*t+c[2,i]
   END;
   ansy1 := ansy1/d1;
   ansy2 := ansy2/d2
END;
FUNCTION bessi(n: integer; x: real): real;
CONST
   iacc=40;
   bigno=1.0e10;
   bigni=1.0e-10;
VAR
   bi,bim,bip,tox,ans: real;
   j,m: integer;
BEGIN
   IF  (n < 2) THEN BEGIN
      writeln('pause in routine BESSI');
      writeln('index n is less than 2'); readln
   END;
   tox := 2.0/x;
   bip := 0.0;
   bi := 1.0;
   m := 2*(n+trunc(sqrt(iacc*n)));
   FOR j := m DOWNTO 1 DO BEGIN
      bim := bip+j*tox*bi;
      bip := bi;
      bi := bim;
      IF (abs(bi) > bigno) THEN BEGIN
         ans := ans*bigni;
         bi := bi*bigni;
         bip := bip*bigni
      END;
      IF (j=n) THEN ans := bip
   END;
   bessi := ans*bessi0(x)/bi
END;
FUNCTION bessi0(x: real): real;
VAR
   ax: real;  y,ans: double;
BEGIN
   IF (abs(x) < 3.75) THEN BEGIN
      y := sqr(x/3.75);
      ans := 1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492+y*
         (0.2659732+y*(0.360768e-1+y*0.45813e-2)))))  END
   ELSE BEGIN
      ax := abs(x);  y := 3.75/ax;
      ans := (exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
         +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
         +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
         +y*0.392377e-2))))))))  END;
   bessi0 := sngl(ans)
END;
FUNCTION bessi1(x: real): real;
VAR
   ax: real;  y,ans: double;
BEGIN
   IF (abs(x) < 3.75) THEN BEGIN
      y := sqr(x/3.75);
      ans := x*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
         +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3)))))) END
   ELSE BEGIN
      ax := abs(x); y := 3.75/ax;
      ans := 0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1-y*0.420059e-2));
      ans := 0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
         +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
      ans := (exp(ax)/sqrt(ax))*ans  END;
   bessi1 := sngl(ans)
END;
FUNCTION bessj(n: integer; x: real): real;
CONST
   iacc=40;
   bigno=1.0e10;
   bigni=1.0e-10;
VAR
   bj,bjm,bjp,sum,tox,ans: real;
   j,jsum,m: integer;
BEGIN
   IF (n < 2) THEN BEGIN
      writeln('pause in BESSJ'); readln
   END;
   tox := 2.0/x;
   IF (x > 1.0*n) THEN BEGIN
      bjm := bessj0(x);
      bj := bessj1(x);
      FOR j := 1 to n-1 DO BEGIN
         bjp := j*tox*bj-bjm;
         bjm := bj;
         bj := bjp
      END;
      ans := bj
   END ELSE BEGIN
      m := 2*((n+trunc(sqrt(1.0*(iacc*n)))) DIV 2);
      jsum := 0;
      sum := 0.0;
      bjp := 0.0;
      bj := 1.0;
      FOR j := m DOWNTO 1 DO BEGIN
         bjm := j*tox*bj-bjp;
         bjp := bj;
         bj := bjm;
         IF (abs(bj) > bigno) THEN BEGIN
            bj := bj*bigni;
            bjp := bjp*bigni;
            ans := ans*bigni;
            sum := sum*bigni
         END;
         IF (jsum <> 0) THEN sum := sum+bj;
         jsum := 1-jsum;
         IF (j = n) THEN ans := bjp
      END;
      sum := 2.0*sum-bj;
      ans := ans/sum
   END;
   bessj := ans
END;
FUNCTION bessj0(x: real): real;
VAR
   ax,xx,z: real;  y,ans,ans1,ans2: double;
BEGIN
   IF (abs(x) < 8.0) THEN BEGIN
      y := sqr(x);
      ans1 := 57568490574.0+y*(-13362590354.0+y*(651619640.7
         +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
      ans2 := 57568490411.0+y*(1029532985.0+y*(9494680.718
         +y*(59272.64853+y*(267.8532712+y*1.0))));
      bessj0 := sngl(ans1/ans2)  END
   ELSE BEGIN
      ax := abs(x); z := 8.0/ax; y := sqr(z); xx := ax-0.785398164;
      ans1 := 1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
         +y*(-0.2073370639e-5+y*0.2093887211e-6)));
      ans2 := -0.1562499995e-1+y*(0.1430488765e-3
         +y*(-0.6911147651e-5+y*(0.7621095161e-6
         -y*0.934945152e-7)));
      ans := sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
      bessj0 := sngl(ans)  END
END;
FUNCTION bessj1(x: real): real;
VAR
   ax,xx,z: real;  y,ans,ans1,ans2: double;
FUNCTION sign(x: real): real;
   BEGIN
      IF x >= 0.0 THEN sign := 1.0
      ELSE sign := -1.0;
   END;
BEGIN
   IF (abs(x) < 8.0) THEN BEGIN
      y := sqr(x);
      ans1 := x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
         +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
      ans2 := 144725228442.0+y*(2300535178.0+y*(18583304.74
         +y*(99447.43394+y*(376.9991397+y*1.0))));
      bessj1 := sngl(ans1/ans2)  END
   ELSE BEGIN
      ax := abs(x); z := 8.0/ax; y := sqr(z); xx := ax-2.356194491;
      ans1 := 1.0+y*(0.183105e-2+y*(-0.3516396496e-4
         +y*(0.2457520174e-5+y*(-0.240337019e-6))));
      ans2 := 0.04687499995+y*(-0.2002690873e-3
         +y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6)));
      ans := sqrt(0.636619772/ax)*(cos(xx)*ans1
         -z*sin(xx)*ans2)*sign(x);
      bessj1 := sngl(ans)  END
END;
FUNCTION bessk(n: integer; x: real): real;
VAR
   tox,bkp,bkm,bk: real;
   j: integer;
BEGIN
   IF (n < 2) THEN BEGIN
      writeln('pause in routine BESSK');
      writeln('index n less than 2'); readln
   END;
   tox := 2.0/x;
   bkm := bessk0(x);
   bk := bessk1(x);
   FOR j := 1 to n-1 DO BEGIN
      bkp := bkm+j*tox*bk;
      bkm := bk;
      bk := bkp
   END;
   bessk := bk
END;
FUNCTION bessk0(x: real): real;
VAR
   y,ans: double;
BEGIN
   IF (x <= 2.0) THEN BEGIN
      y := x*x/4.0;
      ans := (-ln(x/2.0)*bessi0(x))+(-0.57721566+y*(0.42278420
         +y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2
         +y*(0.10750e-3+y*0.74e-5))))))  END
   ELSE BEGIN
      y := (2.0/x);
      ans := (exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1
         +y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
         +y*(-0.251540e-2+y*0.53208e-3))))))  END;
   bessk0 := sngl(ans)
END;
FUNCTION bessk1(x: real): real;
VAR
   y,ans: double;
BEGIN
   IF  (x <= 2.0)  THEN BEGIN
      y := x*x/4.0;
      ans := (ln(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144
         +y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1
         +y*(-0.110404e-2+y*(-0.4686e-4)))))))  END
   ELSE BEGIN
      y := 2.0/x;
      ans := (exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619
         +y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2
         +y*(0.325614e-2+y*(-0.68245e-3)))))))  END;
   bessk1 := sngl(ans)
END;
FUNCTION bessy(n: integer; x: real): real;
VAR
   by,bym,byp,tox: real;
   j: integer;
BEGIN
   IF (n < 2) THEN BEGIN
      writeln('pause in BESSY - index n less than 2'); readln
   END;
   tox := 2.0/x;
   by := bessy1(x);
   bym := bessy0(x);
   FOR j := 1 to n-1 DO BEGIN
      byp := j*tox*by-bym;
      bym := by;
      by := byp
   END;
   bessy := by
END;
FUNCTION bessy0(x: real): real;
VAR
   xx,z: real;  y,ans,ans1,ans2: double;
BEGIN
   IF (x < 8.0) THEN BEGIN
      y := sqr(x);
      ans1 := -2957821389.0+y*(7062834065.0+y*(-512359803.6
         +y*(10879881.29+y*(-86327.92757+y*228.4622733))));
      ans2 := 40076544269.0+y*(745249964.8+y*(7189466.438
         +y*(47447.26470+y*(226.1030244+y*1.0))));
      ans := (ans1/ans2)+0.636619772*bessj0(x)*ln(x);
      bessy0 := sngl(ans)  END
   ELSE BEGIN
      z := 8.0/x; y := sqr(z); xx := x-0.785398164;
      ans1 := 1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
         +y*(-0.2073370639e-5+y*0.2093887211e-6)));
      ans2 := -0.1562499995e-1+y*(0.1430488765e-3
         +y*(-0.6911147651e-5+y*(0.7621095161e-6+y*(-0.934945152e-7))));
      ans := sin(xx)*ans1+z*cos(xx)*ans2;
      ans := sqrt(0.636619772/x)*ans;
      bessy0 := sngl(ans)  END
END;
FUNCTION bessy1(x: real): real;
VAR
   xx,z: real;  y,ans,ans1,ans2: double;
BEGIN
   IF (x < 8.0) THEN BEGIN
      y := sqr(x);
      ans1 := x*(-0.4900604943e13+y*(0.1275274390e13
         +y*(-0.5153438139e11+y*(0.7349264551e9
         +y*(-0.4237922726e7+y*0.8511937935e4)))));
      ans2 := 0.2499580570e14+y*(0.4244419664e12
         +y*(0.3733650367e10+y*(0.2245904002e8
         +y*(0.1020426050e6+y*(0.3549632885e3+y*1.0)))));
      ans := (ans1/ans2)+0.636619772*(bessj1(x)*ln(x)-1.0/x);
      bessy1 := sngl(ans) END
   ELSE BEGIN
      z := 8.0/x; y := sqr(z); xx := x-2.356194491;
      ans1 := 1.0+y*(0.183105e-2+y*(-0.3516396496e-4
         +y*(0.2457520174e-5+y*(-0.240337019e-6))));
      ans2 := 0.04687499995+y*(-0.2002690873e-3
         +y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6)));
      ans := sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
      bessy1 := sngl(ans) END
END;
FUNCTION beta(VAR z,w: real): real;
BEGIN
   beta := exp(gammln(z)+gammln(w)-gammln(z+w))
END;
FUNCTION betacf(a,b,x: real): real;
LABEL 1;
CONST
   itmax=100;
   eps=3.0e-7;
VAR
   tem,qap,qam,qab,em,d: real;
   bz,bpp,bp,bm,az,app: real;
   am,aold,ap: real;
   m: integer;
BEGIN
   am := 1.0;
   bm := 1.0;
   az := 1.0;
   qab := a+b;
   qap := a+1.0;
   qam := a-1.0;
   bz := 1.0-qab*x/qap;
   FOR m := 1 to itmax DO BEGIN
      em := m;
      tem := em+em;
      d := em*(b-m)*x/((qam+tem)*(a+tem));
      ap := az+d*am;
      bp := bz+d*bm;
      d := -(a+em)*(qab+em)*x/((a+tem)*(qap+tem));
      app := ap+d*az;
      bpp := bp+d*bz;
      aold := az;
      am := ap/bpp;
      bm := bp/bpp;
      az := app/bpp;
      bz := 1.0;
      IF ((abs(az-aold)) < (eps*abs(az))) THEN GOTO 1
   END;
   writeln('pause in BETACF');
   writeln('a or b too big, or itmax too small'); readln;
1:   betacf := az
END;
FUNCTION betai(a,b,x: real): real;
VAR
   bt: real;
BEGIN
   IF ((x < 0.0) OR (x > 1.0)) THEN BEGIN
      writeln('pause in routine BETAI'); readln
   END;
   IF ((x = 0.0) OR (x = 1.0)) THEN bt := 0.0
   ELSE bt := exp(gammln(a+b)-gammln(a)-gammln(b)
           +a*ln(x)+b*ln(1.0-x));
   IF (x < ((a+1.0)/(a+b+2.0))) THEN
      betai := bt*betacf(a,b,x)/a
   ELSE betai := 1.0-bt*betacf(b,a,1.0-x)/b
END;
FUNCTION bico(n,k: integer): real;
BEGIN
   bico := round(exp(factln(n)-factln(k)-factln(n-k)));
END;
PROCEDURE bksub(ne,nb,jf,k1,k2: integer; VAR c: glcarray;
      nci,ncj,nck: integer);
(* Programs using routine BKSUB must define the type
TYPE
   glcarray = ARRAY [1..nci,1..ncj,1..nck] OF real;
in the main routine. *)
VAR
   nbf,kp,k,j,i: integer;
   xx: real;
BEGIN
   nbf := ne-nb;
   FOR k := k2 DOWNTO k1 DO BEGIN
      kp := k+1;
      FOR j := 1 to nbf DO BEGIN
         xx := c[j,jf,kp];
         FOR i := 1 to ne DO BEGIN
            c[i,jf,k] := c[i,jf,k]-c[i,j,k]*xx
         END
      END
   END;
   FOR k := k1 to k2 DO BEGIN
      kp := k+1;
      FOR i := 1 to nb DO BEGIN
         c[i,1,k] := c[i+nbf,jf,k]
      END;
      FOR i := 1 to nbf DO BEGIN
         c[i+nb,1,k] := c[i,jf,kp]
      END
   END
END;
FUNCTION bnldev(pp: real; n: integer; VAR idum: integer): real;
LABEL 1;
CONST
   pi=3.141592654;
VAR
   am,em,en,g,angle: real;
   oldg,p,pc,bnl: real;
   pclog,plog,pold,sq,t,y: real;
   j,nold: integer;
BEGIN
   nold := -1;
   pold := -1.0;
   IF (pp <= 0.5) THEN p := pp ELSE p := 1.0-pp;
   am := n*p;
   IF (n < 25) THEN BEGIN
      bnl := 0.0;
      FOR j := 1 to n DO BEGIN
         IF (ran3(idum) < p) THEN bnl := bnl+1.0
      END
   END ELSE IF (am < 1.0) THEN BEGIN
      g := exp(-am);
      t := 1.0;
      FOR j := 0 to n DO BEGIN
         t := t*ran3(idum);
         IF (t < g) THEN GOTO 1
      END;
      j := n;
1:      bnl := j
   END ELSE BEGIN
      IF (n <> nold) THEN BEGIN
         en := n;
         oldg := gammln(en+1.0);
         nold := n
      END;
      IF (p <> pold) THEN BEGIN
         pc := 1.0-p;
         plog := ln(p);
         pclog := ln(pc);
         pold := p
      END;
      sq := sqrt(2.0*am*pc);
      REPEAT
         REPEAT
            angle := pi*ran3(idum);
            y := sin(angle)/cos(angle);
            em := sq*y+am;
         UNTIL ((em >= 0.0) AND (em < en+1.0));
         em := trunc(em);
         t := 1.2*sq*(1.0+sqr(y))*exp(oldg-gammln(em+1.0)
            -gammln(en-em+1.0)+em*plog+(en-em)*pclog);
      UNTIL (ran3(idum) <= t);
      bnl := em
   END;
   IF (p <> pp) THEN bnl := n-bnl;
   bnldev := bnl
END;
FUNCTION brent(ax,bx,cx,tol: real; VAR xmin: real): real;
(* Programs using routine BRENT must supply an external function
func(x:real):real whose minimum is to be found. *)
LABEL 1,2,3;
CONST
   itmax=100;
   cgold=0.3819660;
   zeps=1.0e-10;
VAR
   a,b,d,e,etemp: real;
   fu,fv,fw,fx: real;
   iter: integer;
   p,q,r,tol1,tol2: real;
   u,v,w,x,xm: real;
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;
   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
         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) THEN  p := -p;
         q := abs(q);
         etemp := e;
         e := d;
         IF((abs(p) >= abs(0.5*q*etemp)) OR (p <= q*(a-x))
            OR (p >= q*(b-x))) THEN GOTO 1;
         d := p/q;
         u := x+d;
         IF (((u-a) < tol2) OR ((b-u) < tol2)) THEN d := sign(tol1,xm-x);
         GOTO 2
      END;
1:      IF (x >= xm)  THEN e := a-x ELSE e := b-x;
      d := cgold*e;
2:      IF (abs(d) >= tol1)  THEN u := x+d ELSE u := x+sign(tol1,d);
      fu := func(u);
      IF (fu <= fx)  THEN BEGIN
         IF (u >= x)  THEN a := x ELSE b := x;
         v := w;
         fv := fw;
         w := x;
         fw := fx;
         x := u;
         fx := fu
      END ELSE BEGIN
         IF (u < x)  THEN a := u ELSE b := u;
         IF ((fu <= fw) OR (w = x))  THEN BEGIN
            v := w;
            fv := fw;
            w := u;
            fw := fu
         END ELSE IF ((fu <= fv) OR (v = x) OR (v = 2)) THEN BEGIN
            v := u;
            fv := fu
         END
      END
   END;
   writeln('pause in routine BRENT - too many iterations');
3:   xmin := x;
   brent := fx
END;
PROCEDURE bsstep(VAR y: glyarray; dydx: glyarray; nv: integer; VAR x: real;
      htry,eps: real; yscal: glyarray; VAR hdid,hnext: real);
(* Programs using routine BSSTEP must define the type
TYPE
   glyarray = ARRAY [1..nv] OF real;
in the main routine. *)
LABEL 99;
CONST
   imax=11;
   nuse=7;
   one=1.0e0;
   shrink=0.95e0;
   grow=1.2e0;
VAR
   j,i: integer;
   xsav,xest,h,errmax: real;
   ysav,dysav,yseq,yerr: glyarray;
   nseq: ARRAY [1..imax] OF integer;
BEGIN
   nseq[1] := 2; nseq[2] := 4; nseq[3] := 6;
   nseq[4] := 8; nseq[5] := 12; nseq[6] := 16;
   nseq[7] := 24; nseq[8] := 32; nseq[9] := 48;
   nseq[10] := 64; nseq[11] := 96;
   h := htry;
   xsav := x;
   FOR i := 1 to nv DO BEGIN
      ysav[i] := y[i];
      dysav[i] := dydx[i]
   END;
   WHILE true DO BEGIN
      FOR i := 1 to imax DO BEGIN
         mmid(ysav,dysav,nv,xsav,h,nseq[i],yseq);
         xest := sqr(h/nseq[i]);
         rzextr(i,xest,yseq,y,yerr,nv,nuse);
         errmax := 0.0;
         FOR j := 1 to nv DO BEGIN
            IF (errmax < abs(yerr[j]/yscal[j])) THEN
               errmax := abs(yerr[j]/yscal[j])
         END;
         errmax := errmax/eps;
         IF (errmax < one) THEN BEGIN
            x := x+h;
            hdid := h;
            IF (i = nuse) THEN hnext := h*shrink
            ELSE IF (i = nuse-1) THEN hnext := h*grow
            ELSE hnext := (h*nseq[nuse-1])/nseq[i];
            GOTO 99
         END
      END;
      h := 0.25*h;
      IF (((imax-nuse) DIV 2) > 0) THEN BEGIN
         FOR i := 1 to ((imax-nuse) DIV 2) DO h := h/2
      END;
      IF ((x+h) = x) THEN BEGIN
         writeln('pause in routine BSSTEP');
         writeln('step size underflow'); readln
      END
   END;
99:   END;
