! 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