! 03/06/82 IROUT31S HAS ERFN & ERFC REIMPLEMENTED USING ICL METHOD
!02/06/82 CONVERTED TO IMP90 BY PDS
!17/11/80 - IROUT29S - CAPACITY EXCEEDED CHECKING ADDED TO READSTRING
!COPIED FROM ERCS08.B44SCE_IROUT27A ON 26.10.78
!LRANDOM ADDED 24.1.80
!UNUSED SPECS AND UNUSED ROUTINE FILL REMOVED 10/8/79
!SOME MATHS ROUTINE FAULT NOS ALTERED 19/12/78
!IEXPTEN RENAMED AS EXPTEN
!ILOGTEN RENAMED AS LOGTEN
!ICOT RENAMED AS COT
!IHYPSIN RENAMED AS HYPSIN
!IHYPCOS RENAMED AS HYPCOS
!IHYPTAN RENAMED AS HYPTAN
!THE FOLLOWING ROUTINES HAVE BEEN REMOVED:
!DA AND SQ ROUTINES,FROMSTRING,CLOSESTREAM,IFDBINARY
!RFDBINARY,CPUTIME,TIME,DATE,SETMARGINS
!SPECS ADDED FOR CPUTIME, TIME , AND DATE
!* MODIFIED 10/08/78 AGK   (ADDED MLIBERR ROUTINE)
!* MODIFIED 20/02/78 ERROR MESSAGE VALUES
!* MODIFIED 29/09/78 LCG TO REMOVE MLIBERR
!* MODIFIED 17/2/77  LAB   (APPENDED SQ & DA FILE ROUTINES AND
!                           DATE AND TIME ROUTINES FOR IMP AND FORTE)
!* MODIFIED  16/12/76  GEM (RANDOM)
%SYSTEMROUTINESPEC MOVE( %INTEGER L,F,T)
%EXTERNALLONGREALFNSPEC CPUTIME
%EXTERNALSTRINGFNSPEC TIME
%EXTERNALSTRINGFNSPEC DATE
%SYSTEMROUTINESPEC MLIBERR( %INTEGER ERROR)
%SYSTEMROUTINESPEC IOCP( %INTEGER EP,N)
!******** MODIFIED 02:07:76 15.15  LCG    (ALGLRTS,MATHFNS CONCATONATED
!                                         ,COMPLEX ROUTINES
!                                         & DUPLICATES FOR FORTRAN)
%SYSTEMLONGREALFNSPEC ISIN( %LONGREAL X)
%SYSTEMLONGREALFNSPEC ICOS( %LONGREAL X)
%SYSTEMLONGREALFNSPEC IEXP( %LONGREAL X)
%SYSTEMLONGREALFNSPEC ISQRT( %LONGREAL X)
%SYSTEMLONGREALFNSPEC ILOG( %LONGREAL X)
%CONSTLONGREAL LOG10A=2.3025850929940456840179914546843642076011
%OWNLONGREAL PI=3.141592653589793238462643
%CONSTLONGREAL R1=R'41C98867F42983DF'
%CONSTLONGREAL R2=R'C2562FB2813C6014'
%CONSTLONGREAL R3=R'C1146D547FED8A3D'
%CONSTLONGREAL R4=R'C0157BD961F06C89'
%CONSTLONGREAL S1=R'421B189E39236635'
%CONSTLONGREAL S2=R'4168EE1BDE0C3700'
%CONSTLONGREAL S3=R'41224E7F3CBDFE41'
%CONSTLONGREAL S4=R'41144831DAFBF542'
%CONSTLONGREAL RT3=R'411BB67AE8584CAA'
%CONSTLONGREAL PIBY6=R'40860A91C16B9B2C'
%CONSTLONGREAL PIBY2M1=R'40921FB54442D184'
%CONSTLONGREAL RT3M1=R'40BB67AE8584CAA7'
%CONSTLONGREAL TANPIBY12=R'404498517A7B3559'
%SYSTEMLONGREALFN IARCTAN( %LONGREAL X2,X1)
%LONGREAL CONSTANT
%LONGREAL Y,YSQ,X
%INTEGER SIGN,SIGN1,SIGN2
      SIGN = 1; SIGN1 = 1; SIGN2 = 1
      CONSTANT = 0
      %IF X2 < 0 %THEN SIGN2 =  - 1 %AND X2 =  - X2
      %IF X1 < 0 %THEN SIGN1 =  - 1 %AND X1 =  - X1
      %IF X1 = 0 %AND X2 = 0 %START
         MLIBERR(71);                    ! ERROR(1, 38, 1, 0, 1, 0)
         X = 0
      %FINISHELSESTART
         %IF X1 > X2 %START
            SIGN =  - 1
            Y = X2 / X1
         %FINISHELSE Y = X1 / X2
         %IF Y > TANPIBY12 %THEN Y = (RT3M1 * Y - 1.0 + Y) / (Y + RT3) %AND %C
            CONSTANT = PIBY6
         YSQ = Y * Y
         X = Y * (R1 / (YSQ + S1 + (R2 / (YSQ + S2 + (R3 / (YSQ + S3 + (R4 / %C
            (YSQ + S4)))))))) + CONSTANT
         %IF SIGN =  - 1 %THEN X = 1.0 - X + PIBY2M1
         %IF SIGN1 = 1 %START
            %IF SIGN2 =  - 1 %THEN X = PI - X
         %FINISHELSESTART
            %IF SIGN2 =  - 1 %THEN X = X - PI %ELSE X =  - X
         %FINISH
      %FINISH
      %RESULT = X
%END
%EXTERNALLONGREALFN HYPTAN( %LONGREAL X)
%CONSTLONGREAL A1=676440.765
%CONSTLONGREAL A2=45092.124
%CONSTLONGREAL A3=594.459
%CONSTLONGREAL B1=2029322.295
%CONSTLONGREAL B2=947005.29
%CONSTLONGREAL B3=52028.55
%CONSTLONGREAL B4=630.476
%LONGREAL XSQ,NUM,DENOM,RES
%INTEGER MINUS
      MINUS = 0
      %IF X <=  - 0.54931 %THEN MINUS = 1 %AND X =  - X %ELSE MINUS = 0
      %IF(X >= 0.54931 %AND X < 20.101) %START
         RES = 1 - (2.0 / (IEXP(2.0 * X) + 1.0))
         %IF MINUS = 1 %THENRESULT =  - RES %ELSERESULT = RES
      %FINISH
      %IF X >= 20.101 %START
         %IF MINUS = 1 %THENRESULT =  - 1.0 %ELSERESULT = 1.0
      %FINISH
      XSQ = X * X
      NUM = XSQ * (A1 + XSQ * (A2 + XSQ * (A3 + XSQ)))
      DENOM = B1 + XSQ * (B2 + XSQ * (B3 + XSQ * (B4 + XSQ)))
      %RESULT = (1 - NUM / DENOM) * X
%END
%EXTERNALLONGREALFN COT( %LONGREAL X)
%LONGREAL DENOM
%CONSTLONGREAL MAX=1.0@15
%CONSTLONGREAL NZERO=0.00000000000001;   !@-14
      %IF X > MAX %THEN MLIBERR(66) %ELSESTART
         DENOM = ISIN(X)
         %IF MOD(DENOM) <= NZERO %THEN MLIBERR(67) %ELSERESULT = ICOS(X) / %C
            DENOM
      %FINISH
      %STOP
%END
%SYSTEMLONGREALFN IRADIUS( %LONGREAL X,Y)
!      %IF X*X>7@75 %OR X*X+Y*Y>7@75 %THEN MLIBERR(70)
                                         ! IMERR(25, 10)
      %RESULT = ISQRT(X * X + Y * Y)
%END
%SYSTEMLONGREALFN IARCSIN( %LONGREAL X)
%LONGREAL ARKSN,G,F,Q,Y,Y2
%LONGREALARRAY DT(1:42)
%INTEGER I,J,M,M2
%CONSTLONGREAL DROOT=1.4142135623730950488
%CONSTLONGREALARRAY DA(0:21)=0.0,
                                   1.48666649328710346896,
                                   0.03885303371652290716,
                                   0.00288544142208447113,
                                   0.00028842183344755366,
                                   0.00003322367192785279,
                                   0.00000415847787805283,
                                   0.00000054965045259742,
                                   0.00000007550078449372,
                                   0.00000001067193805630,
                                   0.00000000154218037928,
                                   0.00000000022681145985,
                                   0.00000000003383885639,
                                   0.00000000000510893752,
                                   0.00000000000077911392,
                                   0.00000000000011983786,
                                   .00000000000001856973,
                                   .00000000000000289619,
                                   .00000000000000045428,
                                   .00000000000000007162,
                                   .00000000000000001134,
                                   .00000000000000000180
      ARKSN = 0.0
      G = 0.0
      F = 0.0
      Q = 0.0
      I = 0
      J = 1
      DT(1) = 1.0
      %IF MOD(X) > 1.0 %THEN MLIBERR(72)
                                         ! IMERR(25, 1)
      %IF(MOD(X) >= 1. / DROOT %AND MOD(X) <= 1.0) %THEN Y = ISQRT(1.0 - X * %C
         X) * DROOT %ELSE Y = X * DROOT
      DT(2) = Y
      ARKSN = ARKSN + 0.74333324664350173448
      Y2 = Y * 2
      %CYCLE M=2,1,21
         G = ARKSN
         M2 = M << 1
         DT(M2 - 1) = Y2 * DT(M2 - 2) - DT(M2 - 3)
         DT(M2) = Y2 * DT(M2 - 1) - DT(M2 - 2)
         ARKSN = ARKSN + DA(M) * DT(M2 - 1)
         %IF MOD(G - ARKSN) < 0.0000000000000000005 %THEN ->RESULT
      %REPEAT
RESULT:%IF 1.0 / DROOT <= MOD(X) <= 1.0 %START
         ARKSN = PI / 2.0 - ARKSN * Y
         %IF X < 0.0 %THENRESULT =  - ARKSN %ELSERESULT = ARKSN
      %FINISHELSERESULT = Y * ARKSN
%END
%SYSTEMLONGREALFN IARCCOS( %LONGREAL X)
      %IF MOD(X) > 1.0 %THEN MLIBERR(73)
                                         ! IMERR(25, 2)
      %RESULT = PI / 2.0 - IARCSIN(X)
%END
%EXTERNALLONGREALFN HYPSIN( %LONGREAL X)
%LONGREAL Y,Z,XPOW
      %IF MOD(X) > 172.694 %THEN MLIBERR(74)
                                         ! IMERR(25, 4)
      %IF MOD(X) < 0.3465736 %THEN XPOW = X * X %ANDRESULT = X * (1.0 + XPOW %C
         * (0.16666505 + 0.00836915 * XPOW))
      Y = 0.5 * IEXP(X)
      Z = 0.5 * IEXP( - X)
      %RESULT = Y - Z
%END
%EXTERNALLONGREALFN HYPCOS( %LONGREAL X)
%LONGREAL Y,Z
      %IF MOD(X) > 172.694 %THEN MLIBERR(75)
                                         ! IMERR(25, 5)
      Y = 0.5 * IEXP(X)
      Z = 0.5 * IEXP( - X)
      %RESULT = Y + Z
%END
%EXTERNALLONGREALFN EXPTEN( %LONGREAL X)
%INTEGER I
      I = INTPT(X)
      %IF X - I = 0.0 %THENRESULT = 10.0 ** I
                                         !10@XZ=E@(X*LN10)
      %RESULT = IEXP(LOG10A * X)
%END
%EXTERNALLONGREALFN LOGTEN( %LONGREAL X)
      %IF X = 1.0 %THENRESULT = 0.0
      %RESULT = 0.4342944819032518276511289 * ILOG(X)
%END
%EXTERNALLONGREALFN GAMMAFN( %LONGREAL X)
%LONGREAL RX,G,RT
%INTEGER K
%CONSTLONGREALARRAY A(1:15)=0.99999999999999990,
                            0.42278433509851812,
                            0.41184033042198148,
                            0.08157691940138868,
                            0.07424900794340127,
                           -0.00026695102875553,
                            0.01115381967190670,
                           -0.00285150124303465,
                            0.00209975903507706,
                           -0.00090834655742005,
                            0.00046776781149650,
                           -0.00020644763191593,
                            0.00008155304980664,
                           -0.00002484100538487,
                            0.00000510635920726
      %IF X > 57 %THEN MLIBERR(65)
      RX = 1
      %IF 2 <= X <= 3 %THEN ->G2
      %IF X < 2 %THEN ->G3
G1:   X = X - 1
      RX = RX * X
      %IF X < 3 %THEN ->G2
      ->G1
G3:   %IF MOD(X) < 1@-11 %THEN MLIBERR(65)
      RX = RX / X
      X = X + 1
      %IF X < 2 %THEN ->G3
G2:   G =  - 0.00000051132627267
      RT = X - 2
      %CYCLE K=15, - 1,1
         G = G * RT + A(K)
      %REPEAT
      %RESULT = RX * G
%END
%EXTERNALLONGREALFN LOGGAMMA( %LONGREAL X)
%LONGREAL U,YT,YR,FX,LX,MODD,SX
%INTEGER IND,K
%CONSTLONGREALARRAY A(1:7)=0.083333333333333,
                          -0.002777777777778,
                           0.000793650793651,
                          -0.000595238095238,
                           0.000841750841751,
                          -0.001917520917527,
                           0.006410256410256
      YT = 0.0
      YR = 4.2913@73
      %IF X <= YT %THEN MLIBERR(63)
      %IF X >= YR %THEN MLIBERR(64)
      U = X
      FX = 0
      %IF X < 0 %THEN IND = 1 %ELSE IND =  - 1
      X = MOD(X)
      %IF X >= 13 %THEN ->G1
      %IF X < 2@-10 %THEN ->G6
G2:   LX = LOG(X * X) / 2
      FX = FX + LX
      X = X + 1
      %IF X < 13 %THEN ->G2
G1:   %IF IND =  - 1 %THEN ->G5
      K = INTPT(X + 5@-8)
      %IF MOD(X - K) < 2@-10 %THEN ->G6
G5:   MODD = X * X
      LX =  - 2.955065359477124@-2
      SX = LX * X / MODD
      %CYCLE K=7, - 1,1
         LX = (SX * X) / MODD + A(K)
         SX = (LX * X) / MODD
      %REPEAT
      SX = SX - X - FX + 0.918938533204673
      LX = LOG(X * X) / 2
      SX = SX + (X - 0.5) * LX
      %IF IND = 1 %THENSTART
         FX = SIN(PI * U)
         X = U * FX
         MODD = PI / (X * X)
         LX = LOG(X * MODD * X * MODD) / 2
         SX = LX - SX
      %FINISH
      %RESULT = SX
      ->G7
G6:   %RESULT = 1@17
G7:%END
%LONGREALFNSPEC ERFNC( %LONGREAL INARG)
%EXTERNALLONGREALFN ERFN( %LONGREAL INARG)
!***********************************************************************
!*    EVALUATES THE RRORFN USING ICLS METHOD                           *
!***********************************************************************
%INTEGER SIGN
%LONGREAL X,RESULT
%CONSTLONGREAL L1=0.47
%CONSTLONGREAL L2=6.0924
%CONSTLONGREAL TBYRTPI=R'41120DD750420B69'{2//SQRT(PI)}
%CONSTLONGREAL LAST=R'40F4DE9BD37A6F4D'{22/23}
      *MPSR_X'4080';                     ! MASK OUT UNDERFLOW
                                         ! NO NEEDED IF CALLED FROM IMP
      %IF INARG <= 0 %START
         %IF INARG = 0 %THENRESULT = 0 %ELSE SIGN = 1 %AND X =  - INARG
      %FINISHELSE SIGN = 0 %AND X = INARG
      %IF X < L1 %START
         RESULT =  - (X * X)
         RESULT = TBYRTPI * EXP(RESULT) * (X / (1 + (2 * RESULT / (3 - (4 * %C
            RESULT / (5 + (6 * RESULT / (7 - (8 * RESULT / (9 + (10 * RESULT %C
            / (11 - (12 * RESULT / (13 * (14 * RESULT / (15 - (16 * RESULT / %C
            (17 + (18 * RESULT / (19 - (20 * RESULT / (21 + (LAST * %C
            RESULT)))))))))))))))))))))))
      %FINISHELSEIF X < L2 %THEN RESULT = 1 - ERFNC(X) %ELSE RESULT = 1
      %IF SIGN = 1 %THEN RESULT =  - RESULT
      %RESULT = RESULT
%END
%EXTERNALLONGREALFN ERFNC( %LONGREAL INARG)
!***********************************************************************
!*    THE COMPLIMENTARY ERROR FN SUCH THAT ERFNC(X)+ERFN(X)=1          *
!***********************************************************************
%INTEGER SIGN;                           ! OR FOR GT>0 1 OTHERWISE
%LONGREAL X,RESULT
%CONSTLONGREAL ZERO=0
%CONSTLONGREAL ONE=1
%CONSTLONGREAL TWO=2
%CONSTLONGREAL L1=0.47{LIMITS FOR 3 RANGES}
%CONSTLONGREAL L2=8.0
%CONSTLONGREAL L3=13.306
%CONSTLONGREAL PL0=0.182633488422951126@4
%CONSTLONGREAL PL1=0.289802932921676556@4
%CONSTLONGREAL PL2=0.232043959025163525@4
%CONSTLONGREAL PL3=0.114326207070388617@4
%CONSTLONGREAL PL4=0.368519615471001064@3
%CONSTLONGREAL PL5=0.770816173036842861@2
%CONSTLONGREAL PL6=0.967580788298726540@1
%CONSTLONGREAL PL7=0.564187782550739741
%CONSTLONGREAL QL0=0.182633488422951126@4
%CONSTLONGREAL QL1=0.495882756472114071@4
%CONSTLONGREAL QL2=0.608954242327244355@4
%CONSTLONGREAL QL3=0.442961280388368273@4
%CONSTLONGREAL QL4=0.209438436778953959@4
%CONSTLONGREAL QL5=0.661736120710765347@3
%CONSTLONGREAL QL6=0.137125596050062220@3
%CONSTLONGREAL QL7=0.171498094362760785@2
%CONSTLONGREAL PH0=0.297886562639399289@1
%CONSTLONGREAL PH1=0.740974060596474179@1
%CONSTLONGREAL PH2=0.616020985310963054@1
%CONSTLONGREAL PH3=0.501904972678426746@1
%CONSTLONGREAL PH4=0.127536664472996595@1
%CONSTLONGREAL PH5=0.564189583547755074
%CONSTLONGREAL QH0=0.336907520698275277@1
%CONSTLONGREAL QH1=0.960896532719278787@1
%CONSTLONGREAL QH2=0.170814407474660043@2
%CONSTLONGREAL QH3=0.120489519278551290@2
%CONSTLONGREAL QH4=0.939603401623505415@1
%CONSTLONGREAL QH5=0.226052852076732697@1
      *MPSR_X'4080';                     ! MASK OUT UNDERFLOW
                                         ! NOT NEEDED IF CALLED FROM IMP
      %IF INARG <= 0 %THENSTART
         %IF INARG = 0 %THENRESULT = ONE %ELSE SIGN = 1 %AND X = INARG
      %FINISHELSESTART
         SIGN = 0; X = INARG
      %FINISH
      %IF X < L1 %THENRESULT = ONE - ERFN(INARG)
      %IF X < L2 %START
         RESULT =  - (X * X)
         RESULT = (((((((PL7 * X + PL6) * X + PL5) * X + PL4) * X + PL3) * X %C
            + PL2) * X + PL1) * X + PL0) / ((((((((X + QL7) * X + QL6) * X + %C
            QL5) * X + QL4) * X + QL3) * X + QL2) * X + QL1) * X + QL0) * %C
            EXP(RESULT)
         %IF SIGN = 1 %THEN RESULT = TWO - RESULT
      %FINISHELSEIF X < L3 %START
         RESULT =  - (X * X)
         RESULT = (((((PH5 * X + PH4) * X + PH3) * X + PH2) * X + PH1) * X + %C
            PH0) / ((((((X + QH5) * X + QH4) * X + QH3) * X + QH2) * X + %C
            QH1) * X + QH0) * EXP(RESULT)
         %IF SIGN = 1 %THEN RESULT = TWO - RESULT
      %FINISHELSEIF SIGN = 1 %THEN RESULT = TWO %ELSE RESULT = 0
      %RESULT = RESULT
%END
%EXTERNALINTEGERFN BITS( %INTEGER X)
!     %INTEGER I, J
!     WHILE X#0 %THEN J=J+X&1 %AND X=X>>1
!     %RESULT =J
      *LSS_X
      *LB_0
      *JAT_4,<OUT>
AGN:  *SHZ_X
      *USH_1
      *ADB_1
      *JAF_4,<AGN>
OUT:  *LSS_ %B
       *EXIT_-64
%END;                                    ! of BITS.
%EXTERNALINTEGERFN PARITY( %INTEGER I)
      %RESULT = 1 - (I & 1) * 2
%END
%EXTERNALINTEGERFN SHIFTC( %INTEGER X,Y)
%INTEGER Z
      %IF - 32 <= Y <= 32 %THEN ->L3
      PRINTSTRING("
      ILLEGAL SHIFTC")
      %MONITOR
      %STOP
L3:   ->L2 %IF Y <= 0
      *LSS_X
      *ROT_Y
      *EXIT_-64
L2:   %RESULT = X %IF Y = 0
      *LSS_X
      *LUH_X
      *USH_Y
      *STUH_ %B
       *EXIT_-64
%END
%EXTERNALROUTINE ISOCARD( %BYTEINTEGERARRAYNAME CARD)
      IOCP(10,ADDR(CARD(1)))
%END
%EXTERNALLONGREALFN RFDISO( %BYTEINTEGERARRAYNAME CARD, %INTEGER COL1,COL2,
    %INTEGERNAME ERROR)
%INTEGER I,J,K,L,M,SCALE
%INTEGER SIGN
%REAL RES
      SIGN = 0
      ->L1 %IF COL1 > COL2
      M = 1; ERROR = 0; RES = 0; SCALE = 0; J = 1
      %CYCLE I=COL1,1,COL2
         K = CARD(I)
         ->L5 %IF K = ' '
         ->L7 %IF K = '+'
         ->L4 %IF K = '-'
         ->L2 %UNLESS '0' <= K <= '9'
         RES = RES * 10 + (K & 15)
         SCALE = 10 * SCALE
         M = 0
         ->L3
L5:      ->L3 %IF M = 1
         ERROR = 1
         RES = 10 * RES
         SCALE = 10 * SCALE
         ->L3
L4:      J =  - 1
L7:      ->L8 %IF SIGN = 1 %OR M = 0; SIGN = 1
L3:   %REPEAT
      ERROR = M ! ERROR
      %RESULT = J * RES %UNLESS SCALE # 0
      %RESULT = J * RES / SCALE
L1:   ERROR = 3
      %RESULT = 0
L2:   ->L7 %IF CARD(I) = '&'
      ->L6 %IF CARD(I) = '.'
L8:   ERROR = 2
      %RESULT = I
L6:   ->L8 %IF SCALE # 0; SCALE = 1
      ->L3
%END
%EXTERNALINTEGERFN IFDISO( %BYTEINTEGERARRAYNAME CARD, %INTEGER COL1,COL2,
    %INTEGERNAME ERROR)
%INTEGER I,J,K,RES
%INTEGER SIGN
%INTEGER M
      ->L1 %IF COL1 > COL2
      RES = 0
      ERROR = 0
      SIGN = 0
      M = 1
      J = 1
      %CYCLE I=COL1,1,COL2
         K = CARD(I)
         ->L5 %IF K = ' '
         ->L6 %IF K = '+'
         ->L4 %IF K = '-'
         ->L2 %UNLESS '0' <= K <= '9'
         RES = 10 * RES + K & 15
         M = 0
         ->L3
L5:      ->L3 %IF M = 1
         ERROR = 1
         RES = 10 * RES
         ->L3
L4:      J =  - 1
L6:      ->L7 %IF SIGN = 1 %OR M = 0
         SIGN = 1
L3:   %REPEAT
      ERROR = M ! ERROR
      %RESULT = J * RES
L1:   ERROR = 3
      %RESULT = 0
L2:   ->L6 %IF K = '&'
L7:   ERROR = 2
      %RESULT = I
%END
%EXTERNALROUTINE SOLVELNEQ( %LONGREALARRAYNAME A,B, %INTEGER N,
    %LONGREALNAME DET)
%LONGREAL AMAX,CH
%INTEGER I,J,JMAX,S
      MLIBERR(X'A07') %UNLESS N > 0
      ->L1 %IF N > 1
      DET = A(1,1)
      ->L2 %IF DET = 0
      B(1) = B(1) / DET
      ->L2
L1:   DET = 1
      %CYCLE I=1,1,N - 1
         AMAX = 0; JMAX = 0
         %CYCLE J=I,1,N
            ->L4 %IF MOD(A(J,I)) <= AMAX
            AMAX = MOD(A(J,I)); JMAX = J
L4:      %REPEAT
          ->L5 %IF JMAX = I
         DET =  - DET
         ->L6 %IF JMAX # 0
         DET = 0; ->L2
L6:      %CYCLE J=I,1,N
            CH = A(I,J)
            A(I,J) = A(JMAX,J)
            A(JMAX,J) = CH
         %REPEAT
         CH = B(I)
         B(I) = B(JMAX)
         B(JMAX) = CH
L5:      CH = A(I,I)
         DET = DET * CH
         %CYCLE J=I + 1,1,N
            AMAX = A(J,I) / CH
            %CYCLE S=I + 1,1,N
               A(J,S) = A(J,S) - A(I,S) * AMAX
            %REPEAT
            B(J) = B(J) - B(I) * AMAX
         %REPEAT
      %REPEAT
      CH = A(N,N)
      DET = DET * CH
      ->L2 %IF DET = 0
      B(N) = B(N) / CH
      %CYCLE I=N - 1, - 1,1
         CH = B(I)
         %CYCLE J=I + 1,1,N
            CH = CH - A(I,J) * B(J)
         %REPEAT
         B(I) = CH / A(I,I)
      %REPEAT
L2:%END
%EXTERNALROUTINE DIVMATRIX( %LONGREALARRAYNAME A,B, %INTEGER N,M,
    %LONGREALNAME DET)
                                         ! A=INV(B)A:BNXN, ANXM
%LONGREAL AMAX,CH
%INTEGER I,J,JMAX,S,K
      MLIBERR(X'A07') %UNLESS N > 0
      ->L1 %IF N > 1
      DET = B(1,1)
      ->L2 %IF DET = 0
      %CYCLE I=1,1,M
         A(1,I) = A(1,I) / DET
      %REPEAT
       ->L2
L1:   DET = 1
      %CYCLE I=1,1,N - 1
         AMAX = 0; JMAX = 0
         %CYCLE J=I,1,N
            ->L4 %IF MOD(B(J,I)) <= AMAX
            AMAX = MOD(B(J,I)); JMAX = J
L4:      %REPEAT
          ->L5 %IF JMAX = I
         DET =  - DET
         ->L6 %IF JMAX # 0
         DET = 0; ->L2
L6:      %CYCLE J=I,1,N
            CH = B(I,J)
            B(I,J) = B(JMAX,J)
            B(JMAX,J) = CH
         %REPEAT
         %CYCLE K=1,1,M
            CH = A(I,K)
            A(I,K) = A(JMAX,K)
            A(JMAX,K) = CH
         %REPEAT
L5:      CH = B(I,I)
         DET = DET * CH
         %CYCLE J=I + 1,1,N
            AMAX = B(J,I) / CH
            %CYCLE S=I + 1,1,N
               B(J,S) = B(J,S) - B(I,S) * AMAX
            %REPEAT
            %CYCLE K=1,1,M
               A(J,K) = A(J,K) - A(I,K) * AMAX
            %REPEAT
         %REPEAT
      %REPEAT
      CH = B(N,N)
      DET = DET * CH
      ->L2 %IF DET = 0
      %CYCLE K=1,1,M
         A(N,K) = A(N,K) / CH
      %REPEAT
      %CYCLE I=N - 1, - 1,1
         AMAX = B(I,I)
         %CYCLE K=1,1,M
            CH = A(I,K)
            %CYCLE J=I + 1,1,N
               CH = CH - B(I,J) * A(J,K)
            %REPEAT
            A(I,K) = CH / AMAX
         %REPEAT
      %REPEAT
L2:%END
%EXTERNALROUTINE UNIT( %LONGREALARRAYNAME A, %INTEGER N)
%INTEGER I,J
      MLIBERR(X'A07') %UNLESS N > 0
      %CYCLE I=1,1,N
         %CYCLE J=1,1,N
            A(I,J) = 0
         %REPEAT
         A(I,I) = 1
      %REPEAT
%END
%EXTERNALROUTINE INVERT( %LONGREALARRAYNAME A,B, %INTEGER N, %LONGREALNAME DET)
                                         ! A=INV B USING DIV MATRIX
      MLIBERR(X'A07') %UNLESS N > 0
      UNIT(A,N)
      DIVMATRIX(A,B,N,N,DET)
%END
%EXTERNALLONGREALFN DET( %LONGREALARRAYNAME A, %INTEGER N)
%LONGREALARRAY B(1:N); %LONGREAL DET
%INTEGER I
      MLIBERR(X'A07') %UNLESS N > 0
      %CYCLE I=1,1,N
         B(I) = 0
      %REPEAT
      SOLVELNEQ(A,B,N,DET)
      %RESULT = DET
%END
%EXTERNALROUTINE NULL( %LONGREALARRAYNAME A, %INTEGER N,M)
%INTEGER I,J
      MLIBERR(X'A07') %UNLESS N > 0 %AND M > 0
      %CYCLE I=1,1,N
         %CYCLE J=1,1,M
            A(I,J) = 0
         %REPEAT
      %REPEAT
%END
%EXTERNALROUTINE ADDMATRIX( %LONGREALARRAYNAME A,B,C, %INTEGER N,M)
%INTEGER I,J
      MLIBERR(X'A07') %UNLESS N > 0 %AND M > 0
      %CYCLE I=1,1,N
         %CYCLE J=1,1,M
            A(I,J) = B(I,J) + C(I,J)
         %REPEAT
      %REPEAT
%END
%EXTERNALROUTINE SUBMATRIX( %LONGREALARRAYNAME A,B,C, %INTEGER N,M)
%INTEGER I,J
      MLIBERR(X'A07') %UNLESS N > 0 %AND M > 0
      %CYCLE I=1,1,N
         %CYCLE J=1,1,M
            A(I,J) = B(I,J) - C(I,J)
         %REPEAT
      %REPEAT
%END
%EXTERNALROUTINE COPYMATRIX( %LONGREALARRAYNAME A,B, %INTEGER N,M)
%INTEGER I,J
      MLIBERR(X'A07') %UNLESS N > 0 %AND M > 0
      %CYCLE I=1,1,N
         %CYCLE J=1,1,M
            A(I,J) = B(I,J)
         %REPEAT
      %REPEAT
%END
%EXTERNALROUTINE MULTMATRIX( %LONGREALARRAYNAME A,B,C, %INTEGER N,P,M)
                                         ! A=B*C, A IS N X M
%INTEGER I,J,K
%LONGREAL R
      MLIBERR(X'A07') %UNLESS N > 0 %AND M > 0 %AND P > 0
      %CYCLE I=1,1,N
         %CYCLE J=1,1,M
            R = 0
            %CYCLE K=1,1,P
               R = R + B(I,K) * C(K,J)
            %REPEAT
            A(I,J) = R
         %REPEAT
      %REPEAT
%END
%EXTERNALROUTINE MULTTRMATRIX( %LONGREALARRAYNAME A,B,C, %INTEGER N,P,M)
%LONGREAL R
                                         ! A=B*C, A IS N X M
%INTEGER I,J,K
      MLIBERR(X'A07') %UNLESS N > 0 %AND M > 0 %AND P > 0
      %CYCLE I=1,1,N
         %CYCLE J=1,1,M
            R = 0
            %CYCLE K=1,1,P
               R = R + B(I,K) * C(J,K)
            %REPEAT
            A(I,J) = R
         %REPEAT
      %REPEAT
%END
%EXTERNALROUTINE TRANSMATRIX( %LONGREALARRAYNAME A,B, %INTEGER N,M)
                                         ! AN X M, B M X N
%INTEGER I,J
      MLIBERR(X'A07') %UNLESS N > 0 %AND M > 0
      %CYCLE I=1,1,N
         %CYCLE J=1,1,M
            A(I,J) = B(J,I)
         %REPEAT
      %REPEAT
%END
%EXTERNALREALFN RANDOM( %INTEGERNAME I, %INTEGER N)
%REAL X
%INTEGER J
      ->L2 %IF N < 0; ->L1 %IF N > 1
!      %CONTROL 0;                       !  I=(65539*I)&X'7FFFFFFF'
!  %CONTROL X'1111'
      J = I
      *LSS_65539
      *IMYD_J
      *AND_X'000000007FFFFFFF'
      *STUH_ %B
       *ST_J
      I = J
      %RESULT = 0.0000000004656613 * I
L1:   X = 0
      %CYCLE N=1,1,N; X = X + RANDOM(I,1); %REPEAT; %RESULT = X
L2:   PRINTSTRING("NEGATIVE ARGUMENT IN RANDOM"); %MONITOR
      %STOP
%END
!*
%EXTERNALLONGREALFN LRANDOM( %INTEGERNAME I, %INTEGER N)
%LONGREAL X
%INTEGER J
      ->L2 %IF N < 0; ->L1 %IF N > 1
!      %CONTROL 0;                       !  I=(65539*I)&X'7FFFFFFF'
!  %CONTROL X'1111'
      J = I
      *LSS_65539
      *IMYD_J
      *AND_X'000000007FFFFFFF'
      *STUH_ %B
       *ST_J
      I = J
      %RESULT = 0.0000000004656613 * I
L1:   X = 0
      %CYCLE N=1,1,N; X = X + LRANDOM(I,1); %REPEAT; %RESULT = X
L2:   PRINTSTRING("NEGATIVE ARGUMENT IN LRANDOM"); %MONITOR
      %STOP
%END;                                    !OF LRANDOM
!*
!
! FORTRAN ROUTINES FOR DATE AND TIME CALLABLE ONLY FROM FORTE
!*
%EXTERNALROUTINE CPUTIM( %LONGREALNAME X)
      X = CPUTIME
%END
%EXTERNALROUTINE CTIME( %LONGREALNAME X)
%STRING(10)T
      T = TIME
      MOVE(8,ADDR(T) + 1,ADDR(X))
      BYTEINTEGER(ADDR(X) + 2) = ':'
      BYTEINTEGER(ADDR(X) + 5) = ':'
%END
%EXTERNALROUTINE HDATE( %LONGREALNAME X)
%STRING(10)D
      D = DATE
      MOVE(8,ADDR(D) + 1,ADDR(X))
%END
!*
%EXTERNALROUTINE READSTRING( %INTEGER DR0,DR1)
!WHEN CALLED THIS ROUTINE TAKES A %STRINGNAME PARAMETER
%INTEGER I,DELIM,LEN,MAXLEN; %STRING(1)T
%STRINGNAME S
      LEN = 0;                           !CURRENT LENGTH
      S == STRING(DR1)
      MAXLEN = (DR0 & X'FF') - 1;        !BOUND LESS 1 BYTE FOR LENGTH BYTE
L1:   READSYMBOL(I)
      ->L1 %IF I = 10 %OR I = 32 %OR I = 25
      ->L2 %IF I = '''' %OR I = '"'
      %SIGNALEVENT 4,2;                  !SYMBOL INSTEAD OF STRING
L2:   S = ""; DELIM = I
L3:   ->L4 %IF NEXTSYMBOL = DELIM
      %IF LEN = MAXLEN %THENSIGNALEVENT 6,1
L5:   READITEM(T)
      S = S.T
      LEN = LEN + 1
      ->L3
L4:   SKIPSYMBOL
      ->L5 %IF NEXTSYMBOL = DELIM
%END
%ENDOFFILE