!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'



!
!
!
!*
!     ERROR ROUTINE
!*
!
!
!






%SYSTEMLONGREALFN IARCTAN(%LONGREAL X2, X1)
      %LONGREAL XSQ, CONSTANT
      %LONGREAL Y, YSQ, X
      %INTEGER SIGN, INV, 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/(YSQ+S4))))))) %C
           )+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 %THEN %RESULT=-RES %ELSE %RESULT=RES
%FINISH
      %IF X>=20.101 %START
%IF MINUS=1 %THEN %RESULT=-1.0 %ELSE %RESULT=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) %ELSE %RESULT=ICOS(X)/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 %C
Y=ISQRT(1.0-X*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 %THEN %RESULT=-ARKSN %ELSE %RESULT=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 %AND %C
%RESULT=X*(1.0+XPOW*(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 
      %EXTERNAL %LONG %REAL %FN GAMMAFN(%LONG %REAL X)
         %LONG %REAL          RX,G,RT
         %INTEGER             K
         %CONST %LONG %REAL %ARRAY   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
%CONST%LONGREAL %ARRAY 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 %THEN %START
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
%EXTERNALLONGREALFN ERFN(%LONGREAL X)
      %LONGREAL A0, A1, B0, B1, Z, SUM, T, SA, SB, TEMP
      %INTEGER S, N, SC


      %CONSTLONGREAL EPS=0.0000000000000001;!@-16
      %CONSTLONGREAL SQPI=1.772453850905516
      %IF X=0 %THEN %RESULT=0 
      %IF X>0 %THEN S=1 %ELSE S=-1 %AND X=MOD(X)
      Z=X*X
      %IF X<1 %THENSTART 
         SUM=X
         T=X
         N=0
E1:      N=N+1
         T=-T*Z*(2*N-1)/(N*(2*N+1))
         SUM=SUM+T
         %IF MOD(T)>EPS %THEN ->E1
         SUM=2*SUM/SQPI
      %FINISHELSESTART 
         A0=0
         A1=1
         B0=1
         B1=X
         SUM=A1/B1
         N=1
E2:      N=N+1
SB=1
E3:
TEMP=X*A1+(N-1)*A0/2
%IF MOD(TEMP)>1000 %THEN SA=1000 %ELSE SA=1
A0=A1/SA
A1=TEMP/SA
TEMP=X*B1+(N-1)*B0/2
%IF MOD(TEMP)>1000 %THEN SB=1000 %ELSE SB=0
B0=B1/SB
B0=TEMP/SB
N=N+SC
%IF SC=1 %THEN SC=0 %AND ->E3
         T=(A1/B1)*(SA/SB)
         %IF MOD(SUM-T)>EPS %THENSTART 
            SUM=T
            ->E2
         %FINISH 
         SUM=1-IEXP(-Z)*T/SQPI
      %FINISH 
%RESULT=S*SUM
%END 
%EXTERNALLONGREALFN ERFNC(%LONGREAL X)
      %RESULT =1-ERFN(X)
%END 



%EXTERNALINTEGERFN BITS(%INTEGER X)
!     %INTEGER I, J
!     J=0
!      %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 %C
  , 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, %C
  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 SOLVE LN EQ(%LONGREALARRAYNAME A, B, %INTEGER N, %C
  %LONGREALNAME DET)
      %LONGREAL AMAX, CH
      %INTEGER I, J, J MAX, S
      ->L3 %IF N>0
      PRINTSTRING('
      SOLVE LN EQ DATA FAULT: N=')
      WRITE(N, 2);  NEWLINE;  %STOP 
L3:   ->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
         A MAX=0;  J MAX=0
         %CYCLE J=I, 1, N
            ->L4 %IF MOD(A(J, I))<=AMAX
            AMAX=MOD(A(J, I));  JMAX=J
L4:      %REPEAT 
         ->L5 %IF J MAX=I
         DET=-DET
         ->L6 %IF J MAX#0
         DET=0;  ->L2
L6:      %CYCLE J=I, 1, N
            CH=A(I, J)
            A(I, J)=A(J MAX, J)
            A(J MAX, J)=CH
         %REPEAT 
         CH=B(I)
         B(I)=B(J MAX)
         B(J MAX)=CH
L5:      CH=A(I, I)
         DET=DET*CH
         %CYCLE J=I+1, 1, N
            A MAX=A(J, I)/CH
            %CYCLE S=I+1, 1, N
               A(J, S)=A(J, S)-A(I, S)*A MAX
            %REPEAT 
            B(J)=B(J)-B(I)*A MAX
         %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 DIV MATRIX(%LONGREALARRAYNAME A, B, %INTEGER N, M, %C
  %LONGREALNAME DET)
      %COMMENT A=INV(B)A:BNXN, ANXM
      %LONGREAL AMAX, CH
      %INTEGER I, J, JMAX, S, K
      ->L3 %IF N>0
      PRINTSTRING('
      DIV MATRIX DATA FAULT N=')
      WRITE(N, 2)
      NEWLINE;  %STOP 
L3:   ->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 J MAX=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
      ->L10 %IF N>0
      PRINTSTRING('
MATRIX BOUND ZERO OR NEGATIVE')
      NEWLINE
      %MONITOR 
      %STOP 
L10:  
      %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, %C
  %LONGREALNAME DET)
      %COMMENT A=INV B USING DIV MATRIX
      ->L3 %IF N>0
      PRINTSTRING('
      INVERT DATA FAULT N=')
      WRITE(N, 2);  NEWLINE;  %STOP 
L3:   UNIT(A, N)
      DIV MATRIX(A, B, N, N, DET)
%END 


%EXTERNALLONGREALFN DET(%LONGREALARRAYNAME A, %INTEGER N)
      %LONGREALARRAY B(1:N);  %LONGREAL DET
      %INTEGER I
      ->L10 %IF N>0
      PRINTSTRING('
MATRIX BOUND ZERO OR NEGATIVE')
      NEWLINE
      %MONITOR 
      %STOP 
L10:  
      %CYCLE I=1, 1, N
         B(I)=0
      %REPEAT 
      SOLVE LN EQ(A, B, N, DET)
      %RESULT =DET
%END 


%EXTERNALROUTINE NULL(%LONGREALARRAYNAME A, %INTEGER N, M)
      %INTEGER I, J
      ->L10 %IF N>0 %AND M>0
      PRINTSTRING('
MATRIX BOUND ZERO OR NEGATIVE')
      NEWLINE
      %MONITOR 
      %STOP 
L10:  
      %CYCLE I=1, 1, N
         %CYCLE J=1, 1, M
            A(I, J)=0
         %REPEAT 
      %REPEAT 
%END 


%EXTERNALROUTINE ADD MATRIX(%LONGREALARRAYNAME A, B, C, %INTEGER N, M)
      %INTEGER I, J
      ->L10 %IF N>0 %AND M>0
      PRINTSTRING('
MATRIX BOUND ZERO OR NEGATIVE')
      NEWLINE
      %MONITOR 
      %STOP 
L10:  
      %CYCLE I=1, 1, N
         %CYCLE J=1, 1, M
            A(I, J)=B(I, J)+C(I, J)
         %REPEAT 
      %REPEAT 
%END 


%EXTERNALROUTINE SUB MATRIX(%LONGREALARRAYNAME A, B, C, %INTEGER N, M)
      %INTEGER I, J
      ->L10 %IF N>0 %AND M>0
      PRINTSTRING('
MATRIX BOUND ZERO OR NEGATIVE')
      NEWLINE
      %MONITOR 
      %STOP 
L10:  
      %CYCLE I=1, 1, N
         %CYCLE J=1, 1, M
            A(I, J)=B(I, J)-C(I, J)
         %REPEAT 
      %REPEAT 
%END 


%EXTERNALROUTINE COPY MATRIX(%LONGREALARRAYNAME A, B, %INTEGER N, M)
      %INTEGER I, J
      ->L10 %IF N>0 %AND M>0
      PRINTSTRING('
MATRIX BOUND ZERO OR NEGATIVE')
      NEWLINE
      %MONITOR 
      %STOP 
L10:  
      %CYCLE I=1, 1, N
         %CYCLE J=1, 1, M
            A(I, J)=B(I, J)
         %REPEAT 
      %REPEAT 
%END 


%EXTERNALROUTINE MULT MATRIX(%LONGREALARRAYNAME A, B, C, %INTEGER N, %C
  P, M)
      %COMMENT A=B*C, A IS N X M
      %INTEGER I, J, K
      %LONGREAL R
      ->L10 %IF N>0 %AND M>0 %AND P>0
      PRINTSTRING('
MATRIX BOUND ZERO OR NEGATIVE')
      NEWLINE
      %MONITOR 
      %STOP 
L10:  
      %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 MULT TR MATRIX(%LONGREALARRAYNAME A, B, C, %INTEGER %C
  N, P, M)
      %LONGREAL R
      %COMMENT A=B*C, A IS N X M
      %INTEGER I, J, K
      ->L10 %IF N>0 %AND M>0 %AND P>0
      PRINTSTRING('
MATRIX BOUND ZERO OR NEGATIVE')
      NEWLINE
      %MONITOR 
      %STOP 
L10:  
      %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 TRANS MATRIX(%LONGREALARRAYNAME A, B, %INTEGER N, M)
      %COMMENT AN X M, B M X N
      %INTEGER I, J
      ->L10 %IF N>0 %AND M>0
      PRINTSTRING('
MATRIX BOUND ZERO OR NEGATIVE')
      NEWLINE
      %MONITOR 
      %STOP 
L10:  
      %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 
!*
%SYSTEMROUTINESPEC SSERR(%INTEGER ERROR)
!
!*
!*
!*
!
! 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='"'
%SIGNAL %EVENT 4,1;     !SYMBOL IN DATA
L2: S = 0; DELIM=I
L3: -> L4 %IF NEXT SYMBOL = DELIM
%IF LEN=MAXLEN %THEN %SIGNAL%EVENT 6,1
L5: READ ITEM(T)
   S = S.T
LEN=LEN+1
  -> L3
L4: SKIP SYMBOL
   -> L5 %IF NEXT SYMBOL = DELIM
%END

%ENDOFFILE