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;
PROCEDURE ddpoly(c: glcarray; nc: integer; x: real;
       VAR pd: glpdarray; nd: integer);
(* Programs using routine DDPOLY must define the types
TYPE
   glcarray = ARRAY [1..nc] OF integer;
   glpdarray = ARRAY [1..nd] OF integer;
in the main routine. *)
VAR
   nnd,j,i: integer;
   cnst: real;
BEGIN
   pd[1] := c[nc];
   FOR j := 2 to nd DO BEGIN
      pd[j] := 0.0
   END;
   FOR i := nc-1 DOWNTO 1 DO BEGIN
      IF (nd < (nc+1-i)) THEN nnd := nd ELSE nnd := nc+1-i;
      FOR j := nnd DOWNTO 2 DO BEGIN
         pd[j] := pd[j]*x+pd[j-1]
      END;
      pd[1] := pd[1]*x+c[i]
   END;
   cnst := 2.0;
   FOR i := 3 to nd DO BEGIN
      pd[i] := cnst*pd[i];
      cnst := cnst*i
   END
END;
PROCEDURE des(input,key: gl64array; VAR newkey: integer;
         isw: integer; VAR jotput: gl64array);
(* Programs using routine DES must define the type
TYPE
   gl32array = ARRAY [1..32] OF integer;
   gl64array = ARRAY [1..64] OF integer;
   gl48array = ARRAY [1..48] OF real;
in the main routine. They must also declare the variables
VAR
   ip,ipm: gl64array;
   desflg: boolean;
and initialize desflg to true. The values of ip and ipm are read from
the file desinp.dat whose contents are listed at the end of this routine.
The procedure GLOPEN is used to assign desinp.dat to infile and open the file
for reading. *)
VAR
   j,ii,ic,i: integer;
   titmp,icf: gl32array;
   itmp: gl64array;
   kns: ARRAY[1..48,1..16] OF integer;
   tkns: gl48array;
   infile: text;
BEGIN
   IF desflg THEN BEGIN
      desflg := false;
      glopen(infile,'desinp.dat');
      FOR i := 1 to 64 DO read(infile,ip[i]);
      FOR i := 1 to 64 DO read(infile,ipm[i]);
      close(infile)
   END;
   IF (newkey <> 0) THEN BEGIN
      newkey := 0;
      FOR i := 1 to 16 DO BEGIN
         ks(key,i,tkns);
         FOR j := 1 to 48 DO kns[j,i] := tkns[j]
      END
   END;
   FOR j := 1 to 64 DO itmp[j] := input[ip[j]];
   FOR i := 1 to 16 DO BEGIN
      ii := i;
      IF (isw = 1) THEN ii := 17-i;
      FOR j := 1 to 48 DO tkns[j] := kns[j,ii];
      FOR j := 1 to 32 DO titmp[j] := itmp[32+j];
      cyfun(titmp,tkns,icf);
      FOR j := 1 to 32 DO BEGIN
         ic := icf[j]+itmp[j];
         itmp[j] := itmp[j+32];
         itmp[j+32] := abs(ic MOD 2)
      END
   END;
   FOR j := 1 to 32 DO BEGIN
      ic := itmp[j];
      itmp[j] := itmp[j+32];
      itmp[j+32] := ic
   END;
   FOR j := 1 to 64 DO jotput[j] := itmp[ipm[j]]
END;
(* Contents of the file desinp.dat
58 50 42 34 26 18 10 2 60 52 44 36 28 20 12 4 62 54 46
38 30 22 14 6 64 56 48 40 32 24 16 8 57 49 41 33 25 17 9 1 59 51
43 35 27 19 11 3 61 53 45 37 29 21 13 5 63 55 47 39 31 23 15 7
40 8 48 16 56 24 64 32 39 7 47 15 55 23 63 31 38 6 46 14
54 22 62 30 37 5 45 13 53 21 61 29 36 4 44 12 52 20 60 28 35 3
43 11 51 19 59 27 34 2 42 10 50 18 58 26 33 1 41 9 49 17 57 25 *)
PROCEDURE ks(key: gl64array; n: integer; VAR kn: gl48array);
(* Programs using routine KS must define the types
TYPE
   gl48array = ARRAY [1..48] OF integer;
   gl56array = ARRAY [1..56] OF integer;
   gl64array = ARRAY [1..64] OF integer;
in the main routine. They must also declare the variables
VAR
   glicd,ipc1: gl56array;
   ipc2: gl48array;
   ksflg: boolean;
and initialize ksflg to true. The values of ipc1 and ipc2 are read from
the file ksinpu.dat whose contents are given at the end of this routine. The
procedure GLOPEN assigns ksinpu.dat to infile and opens the file for reading. *)
VAR
   j,it,id,ic,i: integer;
   infile: text;
BEGIN
  IF ksflg THEN BEGIN
      ksflg := false;
      glopen(infile,'ksinpu.dat');
      FOR i := 1 to 56 DO read(infile,ipc1[i]);
      FOR i := 1 to 48 DO read(infile,ipc2[i]);
      close(infile)
   END;
   IF (n = 1) THEN BEGIN
      FOR j := 1 to 56 DO BEGIN
         glicd[j] := key[ipc1[j]]
      END
   END;
   it := 2;
   IF ((n = 1) OR (n = 2) OR (n = 9) OR (n = 16)) THEN it := 1;
   FOR i := 1 to it DO BEGIN
      ic := glicd[1];
      id := glicd[29];
      FOR j := 1 to 27 DO BEGIN
         glicd[j] := glicd[j+1];
         glicd[j+28] := glicd[j+29]
      END;
      glicd[28] := ic;
      glicd[56] := id
   END;
   FOR j := 1 to 48 DO kn[j] := glicd[ipc2[j]]
END;
(* Contents of the file ksinpu.dat
57 49 41 33 25 17 9 1 58 50 42 34 26 18 10 2 59 51 43 35 27 19 11 3 60 52 44 36
63 55 47 39 31 23 15 7 62 54 46 38 30 22 14 6 61 53 45 37 29 21 13 5 28 20 12 4
14 17 11 24 1 5 3 28 15 6 21 10 23 19 12 4 26 8 16 7 27 20 13 2
41 52 31 37 47 55 30 40 51 45 33 48 44 49 39 56 34 53 46 42 50 36 29 32 *)
PROCEDURE cyfun(ir: gl32array; k: gl48array; VAR iout: gl32array);
(* Programs using routine CYFUN must define the types
TYPE
   gl32array = ARRAY [1..32] OF integer;
   gl48array = ARRAY [1..48] OF integer;
in the main routine. They must also declare the variables
VAR
   iet: gl48array;
   ipp: gl32array;
   is: ARRAY[1..16,1..4,1..8] OF integer;
   ibin: ARRAY[1..4,1..16] OF integer;
   cyflg: boolean;
and initialize cyflg to true. The values of iet, ipp, is and ibin are read
from the file cyfuni.dat. The procedure GLOPEN assigns cyfuni.dat to
infile and opens the file for reading. *)
VAR
   jj,irow,icol,kk,iss,ki,j: integer;
   ie: gl48array;
   itmp: gl32array;
   infile: text;
BEGIN
   IF cyflg THEN BEGIN
      cyflg := false;
      glopen(infile,'cyfuni.dat');
      FOR j := 1 to 48 DO read(infile,iet[j]);
      FOR j := 1 to 32 DO read(infile,ipp[j]);
      FOR jj := 1 to 8 DO
         FOR ki := 1 to 4 DO
            FOR j := 1 to 16 DO read(infile,is[j,ki,jj]);
      FOR j := 1 to 16 DO
         FOR ki := 1 to 4 DO read(infile,ibin[ki,j]);
      close(infile)
   END;
   FOR j := 1 to 48 DO ie[j] := (((ir[iet[j]]+k[j]) MOD 2)+2) MOD 2;
   FOR jj := 1 to 8 DO BEGIN
      j := 6*jj-5;
      irow := 2*ie[j]+ie[j+5];
      icol := 8*ie[j+1]+4*ie[j+2]+2*ie[j+3]+ie[j+4];
      iss := is[icol+1,irow+1,jj];
      kk := 4*(jj-1);
      FOR ki := 1 to 4 DO itmp[kk+ki] := ibin[ki,iss+1]
   END;
   FOR j := 1 to 32 DO iout[j] := itmp[ipp[j]]
END;
(* Contents of the file cyfuni.dat
32 1 2 3 4 5 4 5 6 7 8 9 8 9 10 11 12 13 12 13 14 15 16 17
16 17 18 19 20 21 20 21 22 23 24 25 24 25 26 27 28 29 28 29 30 31 32 1
16 7 20 21 29 12 28 17 1 15 23 26 5 18 31 10
2 8 24 14 32 27 3 9 19 13 30 6 22 11 4 25
14 4 13 1 2 15 11 8 3 10 6 12 5 9 0 7 0 15 7 4 14 2 13 1 10 6 12 11 9 5 3 8
4 1 14 8 13 6 2 11 15 12 9 7 3 10 5 0 15 12 8 2 4 9 1 7 5 11 3 14 10 0 6 13
15 1 8 14 6 11 3 4 9 7 2 13 12 0 5 10 3 13 4 7 15 2 8 14 12 0 1 10 6 9 11 5
0 14 7 11 10 4 13 1 5 8 12 6 9 3 2 15 13 8 10 1 3 15 4 2 11 6 7 12 0 5 14 9
10 0 9 14 6 3 15 5 1 13 12 7 11 4 2 8 13 7 0 9 3 4 6 10 2 8 5 14 12 11 15 1
13 6 4 9 8 15 3 0 11 1 2 12 5 10 14 7 1 10 13 0 6 9 8 7 4 15 14 3 11 5 2 12
7 13 14 3 0 6 9 10 1 2 8 5 11 12 4 15 13 8 11 5 6 15 0 3 4 7 2 12 1 10 14 9
10 6 9 0 12 11 7 13 15 1 3 14 5 2 8 4 3 15 0 6 10 1 13 8 9 4 5 11 12 7 2 14
2 12 4 1 7 10 11 6 8 5 3 15 13 0 14 9 14 11 2 12 4 7 13 1 5 0 15 10 3 9 8 6
4 2 1 11 10 13 7 8 15 9 12 5 6 3 0 14 11 8 12 7 1 14 2 13 6 15 0 9 10 4 5 3
12 1 10 15 9 2 6 8 0 13 3 4 14 7 5 11 10 15 4 2 7 12 9 5 6 1 13 14 0 11 3 8
9 14 15 5 2 8 12 3 7 0 4 10 1 13 11 6 4 3 2 12 9 5 15 10 11 14 1 7 6 0 8 13
4 11 2 14 15 0 8 13 3 12 9 7 5 10 6 1 13 0 11 7 4 9 1 10 14 3 5 12 2 15 8 6
1 4 11 13 12 3 7 14 10 15 6 8 0 5 9 2 6 11 13 8 1 4 10 7 9 5 0 15 14 2 3 12
13 2 8 4 6 15 11 1 10 9 3 14 5 0 12 7 1 15 13 8 10 3 7 4 12 5 6 11 0 14 9 2
7 11 4 1 9 12 14 2 0 6 10 13 15 3 5 8 2 1 14 7 4 10 8 13 15 12 9 0 3 5 6 11
0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 1 0 1 0 0 0 1 0 1 0 1 1 0 0 1 1 1
1 0 0 0 1 0 0 1 1 0 1 0 1 0 1 1 1 1 0 0 1 1 0 1 1 1 1 0 1 1 1 1 *)
FUNCTION df1dim(x: real): real;
(* Programs using routine DF1DIM must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
where n is the dimension of vector x. They must also
define the variables
VAR
   ncom: integer;
   pcom,xicom: glnarray
in the main routine, and externally assign them values. *)
VAR
   df1: real;
   j: integer;
   xt,df: glnarray;
BEGIN
   FOR j := 1 to ncom DO BEGIN
      xt[j] := pcom[j]+x*xicom[j]
   END;
   dfunc(xt,df);
   df1 := 0.0;
   FOR j := 1 to ncom DO BEGIN
      df1 := df1+df[j]*xicom[j]
   END;
   df1dim := df1
END;
PROCEDURE dfpmin(VAR p: glnarray; n: integer; ftol: real;
       VAR iter: integer; VAR fret: real);
(* Programs using routine DFPMIN must supply a
FUNCTION fnc(p: glnarray):real; and a
PROCEDURE dfnc(p: glnarray; VAR g: glnarray);
which evaluate a function and its gradient. They must
also define the types
TYPE
   glnarray = ARRAY [1..n] OF real;
   glnbyn = ARRAY [1..n,1..n] OF real;
in the main routine. *)
LABEL 99;
CONST
   itmax=200;
   eps=1.0e-10;
VAR
   j,i,its: integer;
   fp,fae,fad,fac: real;
   xi,g,dg: glnarray;
   hdg: glnarray;
   hessin: glnbyn;
BEGIN
   fp := fnc(p);
   dfnc(p,g);
   FOR i := 1 to n DO BEGIN
      FOR j := 1 to n DO BEGIN
         hessin[i,j] := 0.0
      END;
      hessin[i,i] := 1.0;
      xi[i] := -g[i]
   END;
   FOR its := 1 to itmax DO BEGIN
      iter := its;
      linmin(p,xi,n,fret);
      IF ((2.0*abs(fret-fp)) <= (ftol*(abs(fret)+abs(fp)+eps)))
         THEN GOTO 99;
      fp := fret;
      FOR i := 1 to n DO BEGIN
         dg[i] := g[i]
      END;
      fret := fnc(p);
      dfnc(p,g);
      FOR i := 1 to n DO BEGIN
         dg[i] := g[i]-dg[i]
      END;
      FOR i := 1 to n DO BEGIN
         hdg[i] := 0.0;
         FOR j := 1 to n DO BEGIN
            hdg[i] := hdg[i]+hessin[i,j]*dg[j]
         END
      END;
      fac := 0.0;
      fae := 0.0;
      FOR i := 1 to n DO BEGIN
         fac := fac+dg[i]*xi[i];
         fae := fae+dg[i]*hdg[i]
      END;
      fac := 1.0/fac;
      fad := 1.0/fae;
      FOR i := 1 to n DO BEGIN
         dg[i] := fac*xi[i]-fad*hdg[i]
      END;
      FOR i := 1 to n DO BEGIN
         FOR j := 1 to n DO BEGIN
            hessin[i,j] := hessin[i,j]+fac*xi[i]*xi[j]
            -fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
         END
      END;
      FOR i := 1 to n DO BEGIN
         xi[i] := 0.0;
         FOR j := 1 to n DO BEGIN
            xi[i] := xi[i]-hessin[i,j]*g[j]
         END
      END
   END;
   writeln('pause in routine DFPMIN');
   writeln('too many iterations'); readln;
99:   END;
PROCEDURE difeq(k,k1,k2,jsf,is1,isf: integer; indexv: glindex;
      ne: integer; VAR s: glsarray; nsi,nsj: integer;
      y: glyarray; nyj,nyk: integer);
(* Programs using routine DIFEQ must define the types
TYPE
   glindex = ARRAY [1..nyj] OF integer;
   glsarray = ARRAY [1..nsi,1..nsj] OF real;
   glyarray = ARRAY [1..nyj,1..nyk] OF real;
and also declare the global quantities
CONST
   m=41;
VAR
   h,c2,anorm: real;
   mm,n: integer;
   x: ARRAY [1..m] OF real;
in the main routine.   *)
VAR
   temp2,temp: real;
BEGIN
   IF (k = k1) THEN BEGIN
      IF (((n+mm) MOD 2) = 1) THEN BEGIN
         s[3,3+indexv[1]] := 1.0;
         s[3,3+indexv[2]] := 0.0;
         s[3,3+indexv[3]] := 0.0;
         s[3,jsf] := y[1,1]
      END ELSE BEGIN
         s[3,3+indexv[1]] := 0.0;
         s[3,3+indexv[2]] := 1.0;
         s[3,3+indexv[3]] := 0.0;
         s[3,jsf] := y[2,1]
      END
   END ELSE IF (k > k2) THEN BEGIN
      s[1,3+indexv[1]] := -(y[3,m]-c2)/(2.0*(mm+1.0));
      s[1,3+indexv[2]] := 1.0;
      s[1,3+indexv[3]] := -y[1,m]/(2.0*(mm+1.0));
      s[1,jsf] := y[2,m]-(y[3,m]-c2)*y[1,m]/(2.0*(mm+1.0));
      s[2,3+indexv[1]] := 1.0;
      s[2,3+indexv[2]] := 0.0;
      s[2,3+indexv[3]] := 0.0;
      s[2,jsf] := y[1,m]-anorm
   END ELSE BEGIN
      s[1,indexv[1]] := -1.0;
      s[1,indexv[2]] := -0.5*h;
      s[1,indexv[3]] := 0.0;
      s[1,3+indexv[1]] := 1.0;
      s[1,3+indexv[2]] := -0.5*h;
      s[1,3+indexv[3]] := 0.0;
      temp := h/(1.0-sqr(x[k]+x[k-1])*0.25);
      temp2 := 0.5*(y[3,k]+y[3,k-1])-c2*0.25*sqr(x[k]+x[k-1]);
      s[2,indexv[1]] := temp*temp2*0.5;
      s[2,indexv[2]] := -1.0-0.5*temp*(mm+1.0)*(x[k]+x[k-1]);
      s[2,indexv[3]] := 0.25*temp*(y[1,k]+y[1,k-1]);
      s[2,3+indexv[1]] := s[2,indexv[1]];
      s[2,3+indexv[2]] := 2.0+s[2,indexv[2]];
      s[2,3+indexv[3]] := s[2,indexv[3]];
      s[3,indexv[1]] := 0.0;
      s[3,indexv[2]] := 0.0;
      s[3,indexv[3]] := -1.0;
      s[3,3+indexv[1]] := 0.0;
      s[3,3+indexv[2]] := 0.0;
      s[3,3+indexv[3]] := 1.0;
      s[1,jsf] := y[1,k]-y[1,k-1]-0.5*h*(y[2,k]+y[2,k-1]);
      s[2,jsf] := y[2,k]-y[2,k-1]-temp*((x[k]+x[k-1])
         *0.5*(mm+1.0)*(y[2,k]+y[2,k-1])-temp2
         *0.5*(y[1,k]+y[1,k-1]));
      s[3,jsf] := y[3,k]-y[3,k-1]
   END
END;
