PROCEDURE indexx(n: integer; arrin: glsarray; VAR indx: gliarray);
(* Programs using routine INDEXX must define the types
TYPE
   glsarray = ARRAY [1..np] OF real;
   gliarray = ARRAY [1..np] OF integer;
in the main routine, with np >= n.   *)
LABEL 99;
VAR
   l,j,ir,indxt,i: integer;
   q: real;
BEGIN
   FOR j := 1 to n DO BEGIN
      indx[j] := j
   END;
   l := (n DIV 2) + 1;
   ir := n;
   WHILE true DO BEGIN
      IF (l > 1) THEN BEGIN
            l := l-1;
            indxt := indx[l];
            q := arrin[indxt]
      END ELSE BEGIN
         indxt := indx[ir];
         q := arrin[indxt];
         indx[ir] := indx[1];
         ir := ir-1;
         IF (ir = 1) THEN BEGIN
            indx[1] := indxt;
            GOTO 99
         END
      END;
      i := l;
      j := l+l;
      WHILE (j <= ir) DO BEGIN
         IF (j < ir) THEN BEGIN
             IF (arrin[indx[j]] < arrin[indx[j+1]]) THEN j := j+1
         END;
         IF (q < arrin[indx[j]]) THEN BEGIN
            indx[i] := indx[j];
            i := j;
            j := j+j
         END ELSE
            j := ir+1
      END;
      indx[i] := indxt
   END;
99:   END;
FUNCTION irbit1(VAR iseed: integer): integer;
(* This routine runs much faster if you can perform bitwise logical operations
on integers.  For example, here is a TURBO Pascal version:
CONST
   ib1=1; ib3=4; ib5=16; ib14=8192;
VAR
   newbit: boolean;
BEGIN
   newbit := (iseed AND ib14) <> 0;
   IF ((iseed AND ib5) <> 0) THEN newbit :=  NOT newbit;
   IF ((iseed AND ib3) <> 0) THEN newbit :=  NOT newbit;
   IF ((iseed AND ib1) <> 0) THEN newbit :=  NOT newbit;
   iseed := iseed SHL 1;
   IF (newbit) THEN BEGIN
      irbit1 := 1; iseed := iseed OR ib1
   END ELSE BEGIN
      irbit1 := 0; iseed := iseed AND (NOT ib1)
   END
END; *)
(* Here is the slower version for other Pascal systems: *)
CONST
   ib1=1;
   ib3=4;
   ib5=16;
   ib14=8192;    (* Values chosen not to overflow 2-byte integers *)
VAR
   mask: integer;
   newbit: boolean;
FUNCTION iand(i1,i2: integer): integer;
   VAR
      i: integer;
   BEGIN
      IF ((i1=0) OR (i2=0)) THEN iand := 0
      ELSE BEGIN
         i := ord(odd(i1) AND odd(i2));
         i1 := i1 DIV 2; i2 := i2 DIV 2;
         iand := 2*iand(i1,i2) + i
      END
   END;
FUNCTION inot(ib: integer): integer;
   BEGIN inot := maxint-ib END;
FUNCTION ior(i1,i2: integer): integer;
   VAR
      i: integer;
   BEGIN
      IF ((i1=0) AND (i2=0)) THEN ior := 0
      ELSE BEGIN
         i := ord(odd(i1) OR odd(i2));
         i1 := i1 DIV 2; i2 := i2 DIV 2;
         ior := 2*ior(i1,i2) + i
      END
   END;
BEGIN
   mask := maxint DIV 2;
   newbit := iand(iseed,ib14) <> 0;
   IF (iand(iseed,ib5) <> 0) THEN newbit :=  NOT newbit;
   IF (iand(iseed,ib3) <> 0) THEN newbit :=  NOT newbit;
   IF (iand(iseed,ib1) <> 0) THEN newbit :=  NOT newbit;
   irbit1 := 0;
   iseed := iand(2*iand(mask,iseed),inot(ib1));
   IF (newbit) THEN BEGIN
      irbit1 := 1;
      iseed := ior(iseed,ib1);
   END
END;
FUNCTION irbit2(VAR iseed: integer): integer;
(* This routine runs much faster if you can perform bitwise logical operations
on integers.  For example, here is a TURBO Pascal version:
CONST
   ib1=1; ib3=4; ib5=16; ib14=8192; mask = 21;
BEGIN
   IF (iseed AND ib14) <> 0 THEN BEGIN
      iseed := ((iseed XOR mask) SHL 1) OR ib1;
      irbit2 := 1
   END ELSE BEGIN
      iseed := (iseed SHL 1) AND (NOT ib1);
      irbit2 := 0
   END
END; *)
(* Here is the slower version for other Pascal systems: *)
CONST
   ib1=1;
   ib3=4;
   ib5=16;
   ib14=8192;   (* Values chosen not to overflow 2-byte integers *)
   mask = 21;   (* ib1+ib3+ib5 *)
VAR
   mask1: integer;
FUNCTION iand(i1,i2: integer): integer;
   VAR
      i: integer;
   BEGIN
      IF ((i1=0) OR (i2=0)) THEN iand := 0
      ELSE BEGIN
         i := ord(odd(i1) AND odd(i2));
         i1 := i1 DIV 2; i2 := i2 DIV 2;
         iand := 2*iand(i1,i2) + i
      END
   END;
FUNCTION inot(i: integer): integer;
   BEGIN inot := maxint-i END;
FUNCTION ior(i1,i2: integer): integer;
   VAR
      i: integer;
   BEGIN
      IF ((i1=0) AND (i2=0)) THEN ior := 0
      ELSE BEGIN
         i := ord(odd(i1) OR odd(i2));
         i1 := i1 DIV 2; i2 := i2 DIV 2;
         ior := 2*ior(i1,i2) + i
      END
   END;
FUNCTION ieor(i1,i2: integer): integer;
   VAR
      i: integer;
   BEGIN
      IF ((i1=0) AND (i2=0)) THEN ieor := 0
      ELSE BEGIN
         i := ord((odd(i1) AND NOT odd(i2))
               OR (odd(i2) AND NOT odd(i1)));
         i1 := i1 DIV 2; i2 := i2 DIV 2;
         ieor := 2*ieor(i1,i2) + i
      END
   END;
BEGIN
   mask1 := maxint DIV 2;
   IF (iand(iseed,ib14) <> 0) THEN BEGIN
      iseed := ior(2*iand(ieor(iseed,mask),mask1),ib1);
      irbit2 := 1
   END ELSE BEGIN
      iseed := iand(2*iand(iseed,mask1),inot(ib1));
      irbit2 := 0
   END
END;
