!******************************************************************************
!DOUBLE PRECISION MATHS LIBRARY AFTER CODY AND WAITE,
!SOFTWARE MANUAL FOR THE ELEMENTARY FUNCTIONS, PRENTICE HALL 1980.
!WRITTEN BY D.B.TAYLOR, FEBRUARY 1986 FOR USE WITH SEPARATE SINGLE
!PRECISON ROUTINES.
!THE ENTRY F_POWII IS ALSO PRROVIDED.
!(IT NEEDS EXCEPTION NUMBERS FOR IPOWEXPNEG , IPOWERZERO AND IPOWLARGE
!ASSUMED TO BE 38, 39 AND 40)
!THE FLOATING POINT RECIPEES HAVE BEEN SELECTED AND INDEPENDANT ROUTINES
!HAVE BEEN WRITTEN FOR EACH ENTRY WITH THE EXCEPTION OF DLOG10 WHICH CALLS DLOG,
!DATAN2 WHICH CALLS DATAN, DACOS AND DACOS WHICH CALL DSQRT, AND DSINH,
!DCOSH AND DTANH WHICH CALL DEXP.
!THE ERROR EXITS TO THE DEFINED ERROR EXCEPTION NUMBERS ARE IN LINE WITH
!FORTRAN STANDARDS EXCEPT THAT THE DEXP FUNCTION DOES NOT FLAG UNDERFLOW
!AND ZERO**POSITIVE IS ACCEPTED BY POWDD.
!THE LIMITS FOR THE ARGUEMENTS ARE THOSE DESCRIBED BY CODY AND WAITE,
!SUITABLY INTERPRETED FOR BYTE ARITHEMETIC, WITH THE EXCEPTION THAT THE RANGES
!FOR DSIN/DCOS AND DTAN/DCOTAN ARE EXTENDED TO (-7.205D16 TO 7.205D16)
!AND (-3.521D15 TO 3.521D15) RESPECTIVELY.
!AT THE EXTREMAL VALUES ONLY ONE HEX DIGIT IS ACCURATE. ONLY WITHIN
!THE RESTRICTED RANGES, (-12.108D8 TO 12.108D8 ) FOR DSIN/DCOS AND
!(-4.2165D8 TO 4.2166D8) FOR DTAN/DCOTAN ARE ERRORS CONFINED TO THE BOTTOM
!HEX DIGIT (MORE OR LESS).
!the software test can be run if a version of f_mle is provided which
!does not call diagnostics and termination.
!THE SOFTWARE SATISFIES THE TESTS GIVEN BY CODY AND WAITE
!BUT IS SENSITIVE TO COMPILER CHANGES AND THE TESTS SHOULD BE RERUN IF
!THE SOURCE IS RECOMPILED AND THE SUMMARY TEST RESULTS COMPARED. THE TESTS
!ARE PROVIDED ON A SEPARATE FILE.
!THIS CODE USES EXTENDED PRECISION (LONGLONG) ARITHMETIC IN POWDD AND
!POWDI. ON AMDAHL THE VALUE OF THE INTEGER CONSTANT 'IMPROUNDS' IS SET
!TO 1 TO IMPLY THAT THE IMP CONVERSION FROM LONGLONG TO LONG EMBODIES
!ROUNDING. ON BYTE MACHINES, SUCH AS ICL 2900, ON WHICH IMP DOES NOT
!ROUND HERE, 'IMPROUNDS' SHOULD BE SET TO 0 TO AWAKEN IN-LINE CODE
!TO ROUND LONGLONG TO LONG.
!******************************************************************************
!THE FOLLOWING FUNCTIONS ARE PROVIDED ON THIS FILE. THE NUMBERED GROUPS
!ARE TESTED BY THE TEST SUBROUTINE DTEST'N', WHERE 'N' APPEARS IN THE LEFT
!MARGIN. FOR N= 1 TO 10, THE CORRESPONDING CHAPTERS OF 'CODY AND WAITE' ARE
!NUMBERED N+3, I.E. CHAPS 4 TO 13.
! 1 %EXTERNALLONGREALFN DSQRT %ALIAS "F_DSQRT" (%LONGREALNAME ARG)
! 2 %EXTERNALLONGREALFN DLOG %ALIAS "F_DLOG" (%LONGREALNAME ARG)
! 2 %EXTERNALLONGREALFN DLOG10 %ALIAS "F_DLOG10"(%LONGREALNAME ARG)
! 3 %EXTERNALLONGREALFN DEXP %ALIAS "F_DEXP" (%LONGREALNAME ARG)
! 4 %EXTERNALLONGREALFN POWDD %ALIAS "F_POWDD" (%LONGREAL ARG1,ARG2)
! %EXTERNALLONGREALFN POWDI %ALIAS "F_POWDI" %C
! (%LONGREAL A,%INTEGER N)
! %EXTERNALINTEGERFN POWII %ALIAS "F_POWII" (%INTEGER IARG,NARG)
! 5 %EXTERNALLONGREALFN DSIN %ALIAS "F_DSIN" (%LONGREALNAME ARG)
! 5 %EXTERNALLONGREALFN DCOS %ALIAS "F_DCOS" (%LONGREALNAME ARG)
! 6 %EXTERNALLONGREALFN DTAN %ALIAS "F_DTAN" (%LONGREALNAME ARG)
! 6 %EXTERNALLONGREALFN DCOTAN %ALIAS "F_DCOT"(%LONGREALNAME ARG)
! 7 %EXTERNALLONGREALFN DASIN %ALIAS "F_DASIN" (%LONGREALNAME ARG)
! 7 %EXTERNALLONGREALFN DACOS %ALIAS "F_DACOS" (%LONGREALNAME ARG)
! 8 %EXTERNALLONGREALFN DATAN %ALIAS "F_DATAN" (%LONGREALNAME ARG)
! 8 %EXTERNALLONGREALFN DATAN2 %ALIAS "F_DATAN2"(%LONGREALNAME ARGY,ARGX)
! 9 %EXTERNALLONGREALFN DSINH %ALIAS "F_DSINH" (%LONGREALNAME ARG)
! 9 %EXTERNALLONGREALFN DCOSH %ALIAS "F_DCOSH" (%LONGREALNAME ARG)
!10 %EXTERNALLONGREALFN DTANH %ALIAS "F_DTANH" (%LONGREALNAME ARG)
!11 %EXTERNALLONGREALFN DGAMMA %ALIAS "F_DGAMMA"(%LONGREALNAME FX)
! %EXTERNALLONGREALFN DLGAMA %ALIAS "F_DLGAMMA"(%LONGREALNAME FX)
! %EXTERNALLONGREALFN DERF %ALIAS "F_DERF" (%LONGREALNAME ARG)
! %EXTERNALLONGREALFN DERFC %ALIAS "F_DERFC" (%LONGREALNAME ARG)
!******************************************************************************
%CONSTANTINTEGER IMPROUNDS= 1
!DOUBLE PRECISION ERROR NUMBERS
%CONSTANTINTEGER %C
DLOGSMALL = 1, { DLOG ARG NEGATIVE OR ZERO.}
DSQRTNEG = 2, { DSQRT ARG NEGATIVE}
DEXPLARGE = 3, { DEXP ARG GREATER THAN 174.6731}
DEXPSMALL = 4, { DEXP ARG LESS THAN -180.2182, NOT USED}
DSINLARGE = 5, { DSIN ARG MAGNITUDE GREATER THAN 7.205D16}
DASINLARGE = 6, { DASIN ARG MAGNITUDE GREATER THAN 1.0}
DCOSLARGE = 7, { DCOS ARG MAGNITUDE GREATER THAN 7.205D16 }
DACOSLARGE = 8, { DACOS ARG MAGNITUDE GREATER THAN 1.0}
DTANLARGE = 9, { DTAN/DCOTAN ARG MAGNITUDE GREATER THAN 3.521D15}
DCOTANSMALL = 10, { DCOTAN ARG SMALL IN MAGNITUDE}
DARCSINLARGE = 11, { NOT USED}
DARCCOSLARGE = 12, { NOT USED}
DARCTAN2ZERO = 13, { DATAN2 ARGS ZERO.}
DSINHLARGE = 14, { DSINH ARG MAGNITUDE GREATER THAN 175.3662 }
DCOSHLARGE = 15, { DCOSH ARG MAGNITUDE GREATER THAN 175.3662 }
DPOWERNEG = 16, { NEGATIVE D.P. VALUE RAISED TO A NON-INTEGER POWER}
DPOWERZERO = 17, { D.P. ZERO RAISED TO NON-POSITIVE POWER }
DGAMLARGE = 34, { DGAMMA ARG MAGNITUDE GREATER THAN 57.0}
DGAMNEGINT = 35, { DGAMMA ARG NEAR ZERO OR NEGATIVE INTEGER}
DLGAMNEG = 36, { DLGAMA ARG IS NEGATIVE}
DLGAMLARGE = 37, { DLGAMA ARG IS GREATER THAN 4.29D73}
IPOWEXPNEG = 38, { INTEGER RAISED TO NEGATIVE INTEGER POWER}
IPOWERZERO = 39, { INTEGER ZERO RAISED TO NON-POSITIVE POWER}
IPOWLARGE = 40, { INTEGER TO INTEGER POWER OVERFLOWS}
DPOWERLARGE = 46 { D.P. VALUE RAISED TO TOO LARGE A D.P. POWER }
%CONSTANTINTEGER EXPEXCESS= 64, MAXEXP= 127, MINEXP= 0, EPSEXP= 49
%CONSTANTLONGREAL GREATEST = R'7FFFFFFFFFFFFFFF',
SMALLEST = R'0010000000000000',
MAXEXPARG= 174.6731,
MINEXPARG= -180.2182,
SINLIMIT = 7.205@16,
TANLIMIT = 3.521@15,
ONE = 1.0,
TWO = 2.0,
THREE = 3.0,
TWENTY = 20.0,
HALF = R'4080000000000000',
QUART = R'4040000000000000',
ZERO = 0.0,
SQRTHALF = R'40B504F333F9DE65',
RTEPS = R'3A10000000000000',
TRICK1 = R'4F00000000000000',
TRICK2 = R'4E10000000000000',
TRICK3 = R'4E00000000000000',
PI = R'413243F6A8885A31',
PIBY2 = R'411921FB54442D18',
PIBY3 = R'4110C152382D7366',
PIBY4 = R'40C90FDAA22168C2',
PIBY6 = R'40860A91C16B9B2C'
!GREATEST= 7.2370055773322608@ 75,
!SQRTHALF= 7.0710678118654753@ -1,
!RTEPS= 1.49 @-8,
!PI = 3.14159 26535 89793 23846 26433 8328,
!PIBY2= 1.57079 63267 94896 61923,
!PIBY3= 1.04719 75511 96597 74615,
!PIBY4= 0.78539 81633 97448 30962; !PI/4= 0.62207 73250 42055 06043 (OCTAL)
!PIBY6= 0.52359 87755 98298 87308
%EXTERNALROUTINESPEC ERROR %ALIAS "F_MLE"(%INTEGER TYPE)
!******************************************************************************
%EXTERNALLONGREALFN DSQRT %ALIAS "F_DSQRT" (%LONGREALNAME ARG)
%CONSTANTREAL C1= 1.161322, C2= 0.172924, C3= 0.175241
%LONGREAL IN,Y
%INTEGER E
IN= ARG
%IF IN<=ZERO %THENSTART
IN= -IN
%IF IN=ZERO %THENRESULT= ZERO
ERROR( DSQRTNEG )
%RESULT= ZERO
%FINISH
! OBTAIN HEX EXPONENT AND FRACTION IN .0625 TO 1.0
E= BYTEINTEGER(ADDR(IN)) - EXPEXCESS
BYTEINTEGER(ADDR(IN))= EXPEXCESS
! OBTAIN AN REAL*4 APPROX TO 2*SQRT AND PERFORM 3 NEWTON-RAPHSON STEPS,
! THE LAST MODIFIED TO IMPROVE ROUNDING.
INTEGER(ADDR(Y)+4)= 0
REAL(ADDR(Y))= C1 + REAL(ADDR(IN)) - C2/(REAL(ADDR(IN)) + C3)
Y= QUART*Y + IN/Y
Y= HALF*(Y + IN/Y)
Y= Y + HALF*(IN/Y - Y)
! CORRECT THE EXPONENT IF NECESSARY
%IF E&1#0 %THEN Y= 4.0*Y %AND E= E-1
%IF E#0 %THEN BYTEINTEGER(ADDR(Y))= BYTEINTEGER(ADDR(Y)) + E//2
%RESULT= Y
%END {OF DSQRT}
!******************************************************************************
%EXTERNALLONGREALFN DLOG %ALIAS "F_DLOG"(%LONGREALNAME ARG)
!LOGE2= R'40B17217F7D1CF7A'
%CONSTANTLONGREAL C1= R'40B1800000000000',
C2= R'BDDE8082E3086543',
A0= R'C2401FFC4ACED669',
A1= R'4210624A2016AFED',
A2= R'C0CA20AD9AB5E947',
B0= R'C33017FD381B20EF',
B1= R'43138083FA15267E',
B2= R'C223AB0096CF9C1E'
!C1= 0.69335 93750,
!C2= -0.21219 44400 54690 58277 @-3,
!A0= -0.64124 94342 37455 81147 @ 2,
!A1= 0.16383 94356 30215 34222 @ 2,
!A2= -0.78956 11288 74912 57267,
!B0= -0.76949 93210 84948 79777 @ 3,
!B1= 0.31203 22209 19245 32844 @ 3,
!B2= -0.35667 97773 90346 46171 @ 2
%CONSTANTINTEGERARRAY QVAL(1:15)= 3,2,2,1,1,1,1,0,0,0,0,0,0,0,0
%INTEGER P,Q
%LONGREAL IN,D,Z,W,PRESULT
IN= ARG
%IF IN<= ZERO %THENSTART
ERROR( DLOGSMALL )
%IF IN<ZERO %THEN IN= -IN %ELSE %RESULT= -GREATEST
%FINISH
P= BYTEINTEGER(ADDR(IN)) - EXPEXCESS
BYTEINTEGER(ADDR(IN))= 0
Q= QVAL(INTEGER(ADDR(IN))>>20)
!%IF SHORTF>=X'00400000' %START
! %IF SHORTF>=X'00800000' %THEN Q= 0 %ELSE Q= 1
!%FINISHELSESTART
! %IF SHORTF>=X'00200000' %THEN Q= 2 %ELSE Q= 3
!%FINISH
INTEGER(ADDR(IN))=(INTEGER(ADDR(IN))<<Q)!(INTEGER(ADDR(IN)+4)>>(32-Q))
INTEGER(ADDR(IN)+4)=(INTEGER(ADDR(IN)+4)<<Q)
BYTEINTEGER(ADDR(IN))= X'40'
%IF IN<SQRTHALF %START
Z= IN - HALF
D= Z*HALF + HALF
PRESULT= FLOAT(4*P-Q-1)
%FINISHELSESTART
Z= (IN - HALF) - HALF
D= IN*HALF + HALF
PRESULT= FLOAT(4*P-Q)
%FINISH
IN= Z/D
W= IN*IN
Z= ((A2*W + A1)*W + A0)
D= ((W + B2)*W + B1)*W + B0
Z= (W*Z)/D
Z= IN + (IN*Z)
%RESULT= ((PRESULT*C2 +Z) + PRESULT*C1)
%END {OF DLOG}
!******************************************************************************
%EXTERNALLONGREALFN DLOG10 %ALIAS "F_DLOG10"(%LONGREALNAME ARG)
%CONSTANTLONGREAL C3= R'406F2DEC549B9439'
!C3= 0.43429 44819 03251 82765
%RESULT= C3*DLOG(ARG)
%END {OF DLOG10}
!******************************************************************************
%EXTERNALLONGREALFN DEXP %ALIAS "F_DEXP"(%LONGREALNAME ARG)
!EXPE2= R'40B17217F7D1CF7A'
%CONSTANTLONGREAL %C
C1= R'40B1800000000000',
C2= R'BDDE8082E3086543',
C3= R'41171547652B82FE',
P0= R'403FFFFFFFFFFFFF',
P1= R'3F1C70E46FB3F6E0',
P2= R'3D1152A46F58DC1C',
Q0= HALF,
Q1= R'3FE38C738A128D95',
Q2= R'3E207F32DFBC7747'
!C1= 0.69335 93750,
!C2= -0.21219 44400 54690 58277 @-3
!C3= 0.14426 95040 88896 341 @1,
!P0 = 0.24999 99999 99999 993 @ 0,
!P1 = 0.69436 00015 11792 853 @-2,
!P2 = 0.16520 33002 68279 130 @-4,
!Q0= HALF,
!Q1 = 0.55553 86669 69001 188 @-1,
!Q2 = 0.49586 28849 05541 294 @-3
!P0= 0.25,
!P1= 0.75753 18015 94227 76666 @-2,
!P2= 0.31555 19276 56846 46356 @-4,
!Q0= HALF,
!Q1= 0.56817 30269 85512 21787 @-1,
!Q2= 0.63121 89437 43985 03557 @-3,
!Q3= 0.75104 02839 98700 46114 @-6
%CONSTANTLONGREALARRAY POW2(1:3)= 2.0,4.0,8.0
%INTEGER P,Q,N
%LONGREAL IN,D,Z,W,XN,G
IN= ARG
W= MOD(IN)
%IF W>174.6731 %START
%IF IN<-180.2182 %THENRESULT= ZERO
%IF IN>174.6731 %START
ERROR( DEXPLARGE )
%RESULT= GREATEST
%FINISH
%FINISH
%IF BYTEINTEGER(ADDR(W))<EPSEXP %THEN %RESULT= 1.0
N= INT(IN*C3)
XN= FLOAT(N)
G= (IN - XN*C1) - XN*C2
Z= G*G
W= ((P2*Z + P1)*Z + P0)*G
D= (Q2*Z + Q1)*Z + Q0
W= HALF + (W/(D - W))
N= N + 1
%IF N>=0 %THEN P= N>>2 %ELSE P= X'C0000000' ! (N>>2)
Q= N & 3
!%IF N>=0 %THEN P= N//4 %ELSE P= (N-3)//4
!Q= N - 4*P
LABELFORCOMPILER:
%IF Q>0 %THEN W= W*POW2(Q)
BYTEINTEGER(ADDR(W))= BYTEINTEGER(ADDR(W)) + P
%RESULT= W
%END {OF DEXP}
!******************************************************************************
%EXTERNALLONGREALFN POWDD %ALIAS "F_POWDD"(%LONGREAL ARG1,ARG2)
! ARRAY OF 2**((1-J)/16), J= 1,2,...,17 ROUNDED TO WORKING PRECISION
%CONSTANTLONGREALARRAY A1(1:17)= %C
R'4110000000000000',
R'40F5257D152486CC',
R'40EAC0C6E7DD2439',
R'40E0CCDEEC2A94E1',
R'40D744FCCAD69D6B',
R'40CE248C151F8481',
R'40C5672A115506DB',
R'40BD08A39F580C37',
R'40B504F333F9DE65',
R'40AD583EEA42A14B',
R'40A5FED6A9B15139',
R'409EF5326091A112',
R'409837F0518DB8A9',
R'4091C3D373AB11C3',
R'408B95C1E3EA8BD7',
R'4085AAC367CC487B',
R'4080000000000000'
! ARRAY 2**((1-2*J)/16) - A1(2*J), J= 1,2,...,8. I.E. CORRECTION TO A1.
%CONSTANTLONGREALARRAY A2(1:8)= %C
R'322C7B9D0C7AED98',
R'3211065895048DD3',
R'B21C1DCA7C706A0D',
R'B241577EE04992F1',
R'B239B67F57370A66',
R'B2525F6EE0F61447',
R'32360FD6D8E0AE5A',
R'3214C5C95B8C2154'
%CONSTANTLONGREAL K= R'4071547652B82FE1'
!K= LOG2 - 1 =0.44269504088896341
%CONSTANTLONGREAL %C
P1= R'401555555555554D',
P2= R'3F333333333C101C',
P3= R'3E924921713C6D5C',
P4= R'3E1C78FDDB4AFC28',
Q1= R'40B17217F7D1CF79',
Q2= R'403D7F7BFF05899C',
Q3= R'3FE35846B8181368',
Q4= R'3F276556DC263B30',
Q5= R'3E5761F8635F367C',
Q6= R'3DA17BD701C263A0',
Q7= R'3CFA76EF1C9663FD'
!P1 = 0.83333 33333 33332 11405 @-1,
!P2 = 0.12500 00000 05037 99174 @-1,
!P3 = 0.22321 42128 59242 58967 @-2,
!P4 = 0.43445 77567 21631 19635 @-3,
!Q1 = 0.69314 71805 59945 29629 @ 0,
!Q2 = 0.24022 65069 59095 37056 @ 0,
!Q3 = 0.55504 10866 40855 95326 @-1,
!Q4 = 0.96181 29059 51724 16964 @-2,
!Q5 = 0.13333 54131 35857 84703 @-2,
!Q6 = 0.15400 29044 09897 64601 @-3,
!Q7 = 0.14928 85268 05956 08186 @-4
%CONSTANTLONGREALARRAY INVPOW2(1:4)= %C
R'4080000000000000',
R'4040000000000000',
R'4020000000000000',
R'4010000000000000'
%CONSTANTINTEGERARRAY QVAL(1:15)= 3,2,2,1,1,1,1,0,0,0,0,0,0,0,0
%CONSTANTINTEGER IW1MAX= 4030, IW1MIN= -4159
%INTEGER PP,P,Q,N,M,IW1
%LONGREAL D,Z,R,V,G,X,Y,U1,U2,W,W1,W2,Y1,Y2
%LONGLONGREAL UU1,UU2,WW
X= ARG1; Y= ARG2
%IF X<=ZERO %START
%IF X<>ZERO %START
ERROR( DPOWERNEG )
X= -X
%FINISHELSESTART
%IF Y<=ZERO %START
%IF Y=ZERO %START
ERROR( DPOWERZERO )
%FINISHELSESTART
ERROR( DPOWERNEG)
%FINISH
%RESULT= GREATEST
%FINISHELSERESULT= ZERO
%FINISH
%FINISH
M= BYTEINTEGER(ADDR(X)) - EXPEXCESS
G= X
BYTEINTEGER(ADDR(G))= 0
Q= QVAL(INTEGER(ADDR(G))>>20)
INTEGER(ADDR(G))=(INTEGER(ADDR(G))<<Q)!(INTEGER(ADDR(G)+4)>>(32-Q))
INTEGER(ADDR(G)+4)=(INTEGER(ADDR(G)+4)<<Q)
BYTEINTEGER(ADDR(G))= X'40'
P= 1
%IF G<=A1(9) %THEN P= 9
%IF G<=A1(P+4) %THEN P= P + 4
%IF G<=A1(P+2) %THEN P= P + 2
PP= P + 1
Z= (G - A1(PP)) - A2(PP>>1)
D= HALF*(G + A1(PP))
Z= Z/D
V= Z*Z
R= (((P4*V + P3)*V + P2)*V + P1)*V*Z
R= R + K*R
U2= (R + Z*K) + Z
U1= FLOAT((4*M - Q)*16 - P)
%IF BYTEINTEGER(ADDR(U1))#0 %THEN %C
BYTEINTEGER(ADDR(U1))= BYTEINTEGER(ADDR(U1)) - 1
!Y1= Y + TRICK3
!Y2= Y - Y1
!W= U2*Y + U1*Y2
!W1= W + TRICK3
!W2= W - W1
!W= W1 + U1*Y1
!W1= W + TRICK3
!W2= W2 + (W - W1)
!W= W2 + TRICK3
!IW1= INTPT(16.0*(W1 + W))
!W2= W2 - W
UU1= U1
UU2= U2
WW= Y*(UU1 + UU2)
W= WW
%IF IMPROUNDS=0 %START
WW= WW + (WW - W)
W= WW
%FINISH
W1= W + TRICK3
W2= WW - W1
!IW1= INTPT(16.0*W1)
BYTEINTEGER(ADDR(W1))= BYTEINTEGER(ADDR(W1)) + 1
IW1= INTPT(W1)
!ALLOW OVERFLOW TO CAUSE TERMINATION HERE.
!%IF IW1>IW1MAX %START
!PRINTSTRING(" IW1 IS >4030.") %AND NEWLINES(2)
!%IF IW1<IW1MIN %START
!PRINTSTRING(" IW1 IS <-4159.") %AND NEWLINES(2)
!%IF INTEGER(ADDR(W2))#0 %AND %C
! BYTEINTEGER(ADDR(W2))&X'80'=0 %START
%IF W2>ZERO %START
W2= W2 - INVPOW2(4)
IW1= IW1 + 1
%FINISH
!%IF IW1<0 %THEN I= 0
!N= IW1/16 + I
!P= N*16 - IW1 + 1
!M= N/4 + I
!Q= M*4 - N
%IF IW1<0 %START
N= ((-IW1)>>4)
M= -(N>>2)
P= -(N<<4) - IW1 + 1
Q= -((-M)<<2) + N
%FINISHELSESTART
N= IW1>>4 + 1
M= N>>2 + 1
P= N<<4 - IW1 + 1
!M= (IW1 + X'00000050')>>6
!P= 17 - (IW1 & X'0000000F')
Q= M<<2 - N
%FINISH
Z= ((((((Q7*W2 + Q6)*W2 + Q5)*W2 + Q4)*W2 + Q3)*W2 + Q2)*W2 + Q1)*W2
Z= A1(P) + A1(P)*Z
%IF Q>0 %THEN Z= Z*INVPOW2(Q)
N= BYTEINTEGER(ADDR(Z)) + M
%IF N>MAXEXP %START
ERROR( dpowerLARGE )
%RESULT= GREATEST
%FINISH
%IF N<minexp %START
!PRINTSTRING(" EXPONENTIATION UNDERFLOW."); NEWLINES(2)
%RESULT= ZERO
%FINISH
BYTEINTEGER(ADDR(Z))= N
%RESULT= Z
%END {OF POWDD}
!******************************************************************************
%EXTERNALLONGREALFN POWDI %ALIAS "F_POWDI" %C
(%LONGREAL ARG, %INTEGER NARG)
%CONSTANTLONGLONGREAL LONE= 1.0
!
! CALCULATES ARG**NARG
!
%LONGLONGREAL XX,YY
%LONGREAL Y
%INTEGER N
%IF NARG<0 %THEN N= -NARG %ELSE N= NARG
%IF N=0 %START
%IF ARG#ZERO %THEN %RESULT= ONE %ELSESTART
ERROR( DPOWERZERO )
%RESULT= ZERO
%FINISH
%FINISH
%IF ARG=ZERO %START
%IF NARG>0 %THEN %RESULT= ZERO %ELSESTART
ERROR( DPOWERZERO )
%RESULT= ZERO
%FINISH
%FINISH
%IF ARG=ONE %THEN %RESULT= ONE
%IF ARG=-ONE %START
%IF N&1=1 %THEN %RESULT= -ONE %ELSE %RESULT= ONE
%FINISH
!
! USER POWER FUNCTION FOR REAL EXPONENT IF MAGNITUDE OF EXPONENT > 64
!
%IF N>64 %START
%IF ARG>=ZERO %THEN %RESULT= POWDD( ARG,FLOAT(NARG)) %ELSESTART
%IF N&1#1 %THEN %RESULT= POWDD(-ARG,FLOAT(NARG)) %C
%ELSE %RESULT= -POWDD(-ARG,FLOAT(NARG))
%FINISH
%FINISHELSESTART
!
! USE EXTENDED PRECISION FOR LOWER POWERS (<=64)
! ALLOW OVERFLOW TO SIGNAL FAILURE ON DIVISION OR MULT
!
XX= ARG
YY= LONE
%CYCLE
%IF N&1#0 %THEN YY= YY*XX
N = N>>1
%IF N#0 %THEN XX= XX*XX %ELSE %EXIT
%REPEAT
Y= YY
%IF IMPROUNDS=0 %START
YY= YY + (YY - Y)
Y= YY
%FINISH
%IF NARG<0 %THEN %RESULT= ONE/Y %ELSE %RESULT= Y
%FINISH
%END {OF POWDI}
!******************************************************************************
! CONSTANTS FOR DSIN AND DCOS. NOTE THAT RANGE IS EXTENDED INTO 'GRAINY'
! AREA. SEE NOTE IN INTRODUCTION.
! ADDING TRICK1 TO A LONGREAL GIVES FLOAT(INT(..))
! ADDING TRICK2 FORCES LOWEST INTEGER HEX DIGIT INTO BOTTOM OF 2ND WORD.
! THESE ARE NEEDED FOR LARGE INTEGERS WHICH CANNOT FIT IN 4 BYTES.
!CODY AND WAITE USE SINLIMIT = 2.10828771@8
%CONSTANTLONGREAL %C
DC1 = 3.14160 15625,
DC2 = -0.89089 10206 76153 73566 @-5,
PIINV = 0.31830 98861 83790 67154,
D1T4 = R'C0AAAAAAAAAAAAA9',
D2T4 = R'3F88888888888583',
D3T4 = R'BE34034034027C34',
D4T4 = R'3CB8EF1D29278319',
D5T4 = R'BB1AE6454B5DC0B5',
D6T4 = R'392C2478D0D5AD50',
D7T4 = R'B735C841B811F653',
D8T4 = R'353101FED3502BA1'
! THE D..T4 ARE ESSENTIALLY 4 TIMES THE D.., WITH A LITTLE MASSAGING.
! D1 = -0.16666 66666 66666 65052,
! D2 = 0.83333 33333 33316 50314 @-2,
! D3 = -0.19841 26984 12018 40457 @-3,
! D4 = 0.27557 31921 01527 56119 @-5,
! D5 = -0.25052 10679 82745 84544 @-7,
! D6 = 0.16058 93649 03715 89114 @-9,
! D7 = -0.76429 17806 89104 67734 @-12,
! D8 = 0.27204 79095 78888 46175 @-14
%EXTERNALLONGREALFN DSIN %ALIAS "F_DSIN" (%LONGREALNAME ARG)
%LONGREAL XN,F,G
%LONGREAL X1,X2,NN,FF
%INTEGERNAME N
N== INTEGER(ADDR(NN)+4)
F= MOD(ARG)
%IF F>SINLIMIT %START
ERROR ( DSINLARGE )
%RESULT= ZERO
%FINISH
!N= INT(F*PIINV)
!XN= FLOAT(N)
XN= (F*PIINV) + HALF
XN= XN + TRICK1
NN= XN + TRICK2
X1= F + TRICK1
X2= F - X1
G= X1 + X2
FF= ((X1 - XN*DC1) + X2) - XN*DC2
%IF N&1=0 %THEN F= FF %ELSE F= -FF
%IF MOD(F)>=RTEPS %THENSTART
G= F*F
F= F + QUART*F*G*(((((((D8T4*G+D7T4)*G+D6T4)*G+D5T4)*G+D4T4) %C
*G+D3T4)*G+D2T4)*G+D1T4)
%FINISH
%IF BYTEINTEGER(ADDR(ARG))<X'80' %THENRESULT= F %ELSERESULT= -F
%END {OF DSIN}
!******************************************************************************
%EXTERNALLONGREALFN DCOS %ALIAS "F_DCOS" (%LONGREALNAME ARG)
%LONGREAL XN,F,G,Y
%LONGREAL X1,X2,NN
%INTEGERNAME N
N== INTEGER(ADDR(NN)+4)
F= MOD(ARG)
Y= F + PIBY2
%IF Y>SINLIMIT %START
ERROR( DCOSLARGE )
%RESULT= ZERO
%FINISH
!N= INT(Y*PIINV)
!XN= FLOAT(N) - HALF
XN= (Y*PIINV) + HALF
XN= XN + TRICK1
NN= XN + TRICK2
XN= XN - HALF
! ON MACHINES WITH SLOW DOUBLE PRECISION USE TWO MULTS,
! OTHERWISE USE DOUBLE PRECISION WITH ROUNDING.
!X1= FLOAT(INT(F))
X1= F + TRICK1
X2= F - X1
F= ((X1 - XN*DC1) + X2) - XN*DC2
%IF MOD(F)>=RTEPS %THENSTART
G= F*F
F= F + QUART*F*G*(((((((D8T4*G+D7T4)*G+D6T4)*G+D5T4)*G+D4T4) %C
*G+D3T4)*G+D2T4)*G+D1T4)
%FINISH
%IF N&1=0 %THENRESULT= F %ELSERESULT= -F
%END {OF DCOS}
!******************************************************************************
! CONSTANTS FOR DTAN AND DCOTAN
%CONSTANTLONGREAL %C
TOVPI = R'40A2F9836E4E4415',
TCC1 = R'4119220000000000',
TCC2 = R'BC4ABBBD2E7B9676',
TCP1 = R'C022256BCA9A11FF',
TCP2 = R'3EE0741531DD56F6',
TCP3 = R'BD12BAB72EA2C724',
TCQ1 = R'C0777AC11FEF6755',
TCQ2 = R'3F691E7A85F88565',
TCQ3 = R'BE146F6499094841',
TCQ4 = R'3B85BBA783B3C748'
!TOVPI= 0.63661 97723 67581 34308,
!TCC1 = 1.57080 07812 50000,
!TCC2 = -0.44544 55103 38076 8678 @-5,
!TCP1 = -0.13338 35000 64219 60681 @ 0,
!TCP2 = 0.34248 87823 58905 89960 @-2,
!TCP3 = -0.17861 70734 22544 26711 @-4,
!TCQ1 = -0.46671 68333 97552 94240 @ 0,
!TCQ2 = 0.25663 83228 94401 12864 @-1,
!TCQ3 = -0.31181 53190 70100 27307 @-3,
!TCQ4 = 0.49819 43399 37865 12270 @-6
!******************************************************************************
%ROUTINE TANARGBIG
ERROR( DTANLARGE )
%END {OF TANARGBIG, LOCAL}
%ROUTINE COTANARGSMALL
ERROR( DCOTANSMALL )
%END {OF COTANARGSMALL, LOCAL}
!******************************************************************************
%EXTERNALLONGREALFN DTAN %ALIAS "F_DTAN" (%LONGREALNAME ARG)
!TANLIMIT = 3.521 @ 15
! CODY AND WAITE USE TANLIMIT = 4.216574280 @ 8
%LONGREAL X,X1,X2,F,G,W,Z,XN,NN
%INTEGERNAME N
N== INTEGER(ADDR(NN) + 4)
X= ARG
%IF BYTEINTEGER(ADDR(X))&X'80'=0 %START
%IF X>TANLIMIT %THEN TANARGBIG %AND ->SETGREATEST
XN= (X*TOVPI) + HALF
XN= XN + TRICK1
NN= XN + TRICK2
%FINISHELSESTART
%IF X<-TANLIMIT %THEN TANARGBIG %AND ->SETGREATEST
XN= X*TOVPI - HALF
XN= XN + TRICK1
NN= XN - TRICK2
%FINISH
X1= X + TRICK1
X2= X - X1
F= ((X1 - XN*TCC1) + X2) - XN*TCC2
%IF MOD(F)<RTEPS %THEN W= F %AND Z= ONE %C
%ELSESTART
G= F*F
W= F + F*G*((TCP3*G + TCP2)*G + TCP1)
Z= (((TCQ4*G + TCQ3)*G + TCQ2)*G + TCQ1)*G + ONE
%FINISH
%IF N&1#0 %THEN %RESULT= -Z/W %ELSE %RESULT= �W/Z
SETGREATEST: %RESULT= GREATEST
%END {OF DTAN}
!******************************************************************************
%EXTERNALLONGREALFN DCOTAN %ALIAS "F_DCOT"(%LONGREALNAME ARG)
%CONSTANTLONGREAL EPS1= R'0210000000000001'
!EPS1= 0.1381786...= R'0210000000000001', JUST > 1/GREATEST
!TANLIMIT = 3.521 @ 15
! CODY AND WAITE USE TANLIMIT = 4.216574280 @ 8
%LONGREAL X,X1,X2,F,G,W,Z,XN,NN
%INTEGERNAME N
N== INTEGER(ADDR(NN) + 4)
X= ARG
%IF BYTEINTEGER(ADDR(X))&X'80'=0 %START
%IF X<EPS1 %THEN COTANARGSMALL %AND ->SETGREATEST
%IF X>TANLIMIT %THEN TANARGBIG %AND ->SETGREATEST
XN= (X*TOVPI) + HALF
XN= XN + TRICK1
NN= XN + TRICK2
%FINISHELSESTART
%IF X>-EPS1 %THEN COTANARGSMALL %AND ->SETGREATEST
%IF X<-TANLIMIT %THEN TANARGBIG %AND ->SETGREATEST
XN= X*TOVPI - HALF
XN= XN + TRICK1
NN= XN - TRICK2
%FINISH
X1= X + TRICK1
X2= X - X1
F= ((X1 - XN*TCC1) + X2) - XN*TCC2
%IF MOD(F)<RTEPS %THEN W= F %AND Z= ONE %C
%ELSESTART
G= F*F
W= F + F*G*((TCP3*G + TCP2)*G + TCP1)
Z= (((TCQ4*G + TCQ3)*G + TCQ2)*G + TCQ1)*G + ONE
%FINISH
%IF N&1#0 %THEN %RESULT= -W/Z %ELSE %RESULT= Z/W
SETGREATEST: %RESULT= GREATEST
%END {OF DCOTAN}
!******************************************************************************
! CONSTANTS FOR DASIN AND DACOS
%CONSTANTLONGREAL %C
ASCP1 = R'C21B5E55A83A0A62',
ASCP2 = R'4239354E6C15A914',
ASCP3 = R'C227B059534DB53E',
ASCP4 = R'41A270BB2761C939',
ASCP5 = R'C0B25DEDAF30F326',
ASCQ0 = R'C2A43601F15C3E62',
ASCQ1 = R'431A124F101EB843',
ASCQ2 = R'C317DDCEFC56A848',
ASCQ3 = R'4296F3E4B2C8E37D',
ASCQ4 = R'C217D2E86EF9861F'
!ASCP1 = -0.27368 49452 41642 55994 @ 2,
!ASCP2 = 0.57208 22787 78917 31407 @ 2,
!ASCP3 = -0.39688 86299 75048 77339 @ 2,
!ASCP4 = 0.10152 52223 38064 63645 @ 2,
!ASCP5 = -0.69674 57344 73506 48411,
!ASCQ0 = -0.16421 09671 44985 60795 @ 3,
!ASCQ1 = 0.41714 43024 82604 12556 @ 3,
!ASCQ2 = -0.38186 30336 17501 49284 @ 3,
!ASCQ3 = 0.15095 27084 10306 04719 @ 3,
!ASCQ4 = -0.23823 85915 36702 38830 @ 2
!******************************************************************************
%EXTERNALLONGREALFN DASIN %ALIAS "F_DASIN"(%LONGREALNAME ARG)
%CONSTANTLONGREAL A1= PIBY4
%LONGREAL G,Y,W
%INTEGER FLAG
FLAG= 0
Y= MOD(ARG)
%IF Y>HALF %START
FLAG= 1 - FLAG
%IF Y>ONE %START
ERROR(DASINLARGE)
%RESULT= GREATEST
%FINISH
G= (HALF - Y) + HALF
G= G*HALF
Y= DSQRT(G)
Y= Y + Y; { MULTIPLY BY 2 }
Y= -Y; { Y MUST HAVE BEEN +VE SINCE SQRT IS +VE }
->RATFUN
%FINISH
%IF Y<RTEPS %THEN %RESULT= ARG %ELSE G= Y*Y
RATFUN:
W= ((((ASCP5*G + ASCP4)*G + ASCP3)*G + ASCP2)*G + ASCP1)*G
G= ((((G + ASCQ4)*G + ASCQ3)*G + ASCQ2)*G + ASCQ1)*G + ASCQ0
W= Y + (Y*W/G)
%IF FLAG= 1 %THEN W= (A1 + W) + A1
%IF ARG<ZERO %THEN W= -W
%RESULT= W
%END {OF DASIN}
!******************************************************************************
%EXTERNALLONGREALFN DACOS %ALIAS "F_DACOS"(%LONGREALNAME ARG)
%CONSTLONGREAL A1= 0.78539 81633 97448 30962; !P1/4= 0.62207 73250 42055 06043 (OCTAL)
%CONSTLONGREAL B0= PIBY2
%LONGREAL G,Y,W
%INTEGER FLAG
FLAG= 1
Y= MOD(ARG)
%IF Y>HALF %START
FLAG= 1 - FLAG
%IF Y>ONE %START
ERROR (DACOSLARGE)
%RESULT= GREATEST
%FINISH
G= (HALF - Y) + HALF
G= G*HALF
Y= DSQRT(G)
Y= Y + Y; { MULTIPLY BY 2 }
Y= -Y; { Y MUST HAVE BEEN +VE SINCE SQRT IS +VE }
->RATFUN
%FINISH
%IF Y<RTEPS %THEN W= Y %AND ->SETQUADRANT %ELSE G= Y*Y
RATFUN:
W= ((((ASCP5*G + ASCP4)*G + ASCP3)*G + ASCP2)*G + ASCP1)*G
G= ((((G + ASCQ4)*G + ASCQ3)*G + ASCQ2)*G + ASCQ1)*G + ASCQ0
W= Y + (Y*W/G)
SETQUADRANT:
%IF ARG<ZERO %START
%IF FLAG= 0 %THEN W= (B0 + W) + B0 %ELSE W= (A1 + W) + A1
%FINISH %ELSE %IF FLAG= 1 %THEN W= (A1 - W) + A1 %ELSE W= -W
%RESULT= W
%END {OF DACOS}
!******************************************************************************
%EXTERNALLONGREALFN DATAN %ALIAS "F_DATAN"(%LONGREALNAME ARG)
%CONSTANTLONGREAL %C
ROOT3= R'411BB67AE8584CAA',
RT3M1= R'40BB67AE8584CAA7',
TMRT3= R'404498517A7B3559',
P0 = R'C1DB05328830E70F',
P1 = R'C214817FB9E2BCCB',
P2 = R'C187E9FAE46B531A',
P3 = R'C0D66BD6CD8C3DE9',
Q0 = R'422910F979892B53',
Q1 = R'42562848102DB692',
Q2 = R'423B9414641B47AD',
Q3 = R'41F0624F0A563883'
!ROOT3= 1.73205 08075 68877 29353,
!RT3M1= 0.73205 08075 68877 29353,
!TMRT3= 0.26794 91924 31122 70647,
!P0 = -0.13688 76889 41919 26929 @2,
!P1 = -0.20505 85519 58616 51981 @2,
!P2 = -0.84946 24035 13206 83534 @1,
!P3 = -0.83758 29936 81500 59274 @0,
!Q0 = 0.41066 30668 25757 81263 @2,
!Q1 = 0.86157 34959 71302 42515 @2,
!Q2 = 0.59578 43614 25973 44465 @2,
!Q3 = 0.15024 00116 00285 76121 @2
%CONSTANTLONGREALARRAY A(0:3)= ZERO,PIBY6,PIBY2,PIBY3
%LONGREAL F,G,W
%INTEGER N,S
F= ARG
S= BYTEINTEGER(ADDR(F))&X'80'
%IF S#0 %THEN BYTEINTEGER(ADDR(F))= BYTEINTEGER(ADDR(F))&X'7F'
%IF F>ONE %THEN F= ONE/F %AND N= 2 %ELSE N= 0
%IF F>TMRT3 %START
F= (((RT3M1*F - HALF) - HALF) + F)/(ROOT3 + F)
N= N + 1
%FINISH
%IF MOD(F)<RTEPS %THEN W= F %ELSESTART
G= F*F
W= (((P3*G + P2)*G + P1)*G + P0)*G
W= F + F*W/((((G + Q3)*G + Q2)*G + Q1)*G + Q0)
%FINISH
%IF N>1 %THEN W= -W
%IF N>0 %THEN W= W + A(N)
%IF S#0 %THEN W= -W
%RESULT= W
%END {OF DATAN}
!******************************************************************************
%EXTERNALLONGREALFN DATAN2 %ALIAS "F_DATAN2"(%LONGREALNAME ARGY,ARGX)
%LONGREAL W,INY,INX
%INTEGER EXPDIF,EXPYP64,EXPXP64
INY= ARGY
INX= ARGX
%IF INTEGER(ADDR(INY))=0 %OR INTEGER(ADDR(INX))=0 %START
%IF INTEGER(ADDR(INY))=0 %AND INTEGER(ADDR(INX))=0 %START
ERROR( DARCTAN2ZERO )
%RESULT= ZERO
%FINISH
%IF INTEGER(ADDR(INY))=0 %THEN W= ZERO
%IF INTEGER(ADDR(INX))=0 %THEN W= PIBY2 %AND ->OVERFLOW
%FINISHELSESTART
EXPYP64= BYTEINTEGER(ADDR(INY))&X'7F'
EXPXP64= BYTEINTEGER(ADDR(INX))&X'7F'
EXPDIF= EXPYP64 - EXPXP64
%IF EXPDIF> 60 %THEN W= PIBY2 %AND ->OVERFLOW %ELSE %C
%IF EXPDIF<-61 %THEN W= ZERO %ELSE W= MOD(INY/INX) %AND W=DATAN(W)
%FINISH
%IF BYTEINTEGER(ADDR(INX))&X'80'#0 %THEN W= PI - W
OVERFLOW:
%IF BYTEINTEGER(ADDR(INY))&X'80'#0 %THEN W= -W
%RESULT= W
%END {OF DATAN2}
!******************************************************************************
%EXTERNALLONGREALFN DSINH %ALIAS "F_DSINH"(%LONGREALNAME ARG)
%CONSTANTLONGREAL WBAR= 19.95, { 0.35*(56 + 1) }
LNV = R'40B1730000000000',
VOV2M1= R'3CE80897581016B3',
VSQINV= R'403FFF8BFC520EE7',
P0 = R'C555E44D594CD063',
P1 = R'C42D2B856D28291F',
P2 = R'C2A3C20B1C2DF253',
P3 = R'C0CA273DC37E40D7',
Q0 = R'C620359D017CCE25',
Q1 = R'448D42B91DB2F638',
Q2 = R'C3115BC381C97FF2'
!LNV= 0.69316 10107 42187 50000, {0.542714 (OCTAL)}
!VOV2M1= 0.13830 27787 96019 02638 @-4, { V/2 - 1}
!VSQINV= 0.24999 30850 04514 99336 @ 0,
!P0= -0.35181 28343 01771 17881 @ 6,
!P1= -0.11563 52119 68517 68270 @ 5,
!P2= -0.16375 79820 26307 51372 @ 3,
!P3= -0.78966 12741 73570 99479,
!Q0= -0.21108 77005 81062 71242 @ 7,
!Q1= 0.36162 72310 94218 36460 @ 5,
!Q2= -0.27773 52311 96507 01667 @ 3
%LONGREAL Y,W,Z,F,RF
Y= MOD(ARG)
%IF Y>ONE %START
W= Y - LNV
%IF W>MAXEXPARG %START
ERROR ( DSINHLARGE )
%IF ARG>ZERO %THEN %RESULT= GREATEST %ELSE %RESULT= -GREATEST
%FINISHELSESTART
Z= DEXP(W)
%IF W<=WBAR %THEN Z= Z - VSQINV/Z
Z= Z + VOV2M1*Z
%IF ARG<ZERO %THEN %RESULT= -Z %ELSE %RESULT= Z
%FINISH
%FINISHELSE %IF Y<RTEPS %THEN %RESULT= ARG %ELSESTART
F= ARG*ARG
RF= F*(((P3*F + P2)*F + P1)*F + P0)/(((F + Q2)*F + Q1)*F + Q0)
%RESULT= ARG + (ARG*RF)
%FINISH
%END {OF DSINH}
!******************************************************************************
%EXTERNALLONGREALFN DCOSH %ALIAS "F_DCOSH"(%LONGREALNAME ARG)
%CONSTANTLONGREAL WBAR= 19.95, { 0.35*(56 + 1) }
LNV = R'40B1730000000000',
VOV2M1= R'3CE80897581016B3',
VSQINV= R'403FFF8BFC520EE7'
!LNV= 0.69316 10107 42187 50000, {0.542714 (OCTAL)}
!VOV2M1= 0.13830 27787 96019 02638 @-4,
!VSQINV= 0.24999 30850 04514 99336 @ 0
%LONGREAL Y,W,Z
Y= MOD(ARG)
%IF Y>ONE %START
W= Y - LNV
%IF W>MAXEXPARG %START
ERROR ( DCOSHLARGE )
%RESULT= GREATEST
%FINISH
Z= DEXP(W)
%IF W<=WBAR %THEN Z= Z + VSQINV/Z
%RESULT= Z + VOV2M1*Z
%FINISHELSESTART
Z= DEXP(Y)
%RESULT= HALF*Z + HALF/Z
%FINISH
%END {OF DCOSH}
!******************************************************************************
%EXTERNALLONGREALFN DTANH %ALIAS "F_DTANH"(%LONGREALNAME ARG)
%CONSTLONGREAL XBIG = 21.16, {= (LN(2) + (14+1)*LN(16))/2}
LN3BY2 = R'408C9F53D5681855',
P0 = R'C364D69726F87862',
P1 = R'C26339D686E97332',
P2 = R'C0F6E14677DD3BF4',
Q0 = R'4412E83C574E9692',
Q1 = R'438B9C5A680F8796',
Q2 = R'4270BEA787AFDFEA'
! LN3BY2 = 0.54930 61443 34054 84570, {LN(3)/2 }
! P0 = -0.16134 11902 39962 28053 @ 4,
! P1 = -0.99225 92967 22360 83313 @ 2,
! P2 = -0.96437 49277 72254 69787,
! Q0 = 0.48402 35707 19886 88686 @ 4,
! Q1 = 0.22337 72071 89623 12926 @ 4,
! Q2 = 0.11274 47438 05349 49335 @ 3
%LONGREAL F,ANS,G,R
F= MOD(ARG)
%IF F>XBIG %THEN ANS= ONE %ELSE %IF F>LN3BY2 %START
F= F + F
ANS= HALF - (ONE / (DEXP(F) + ONE))
ANS= ANS + ANS
%FINISH %ELSE %IF F<RTEPS %THEN ANS= F %ELSE %START
G= F*F
R= (((P2*G + P1)*G + P0)*G)/(((G + Q2)*G + Q1)*G + Q0)
ANS= F + (F*R)
%FINISH
%IF ARG<ZERO %THEN %RESULT= -ANS %ELSE %RESULT= ANS
%END {OF DTANH}
!******************************************************************************
%EXTERNALLONGREALFN DGAMMA %ALIAS "F_DGAMMA"(%LONGREALNAME FX)
%CONSTANTLONGREAL FXMAX= 57.0,
GAMMAOFFXMAX= 7.1@74,
XSMALL= 1.0@-17,
XVSMALL= 1.4@-76
%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
%LONGREAL X,Y,RX,G,Z
%INTEGER K
X= FX
Y= MOD(X)
%IF Y>FXMAX %START
!" ABSOLUTE VALUE OF DGAMMA ARG > 57."
ERROR( DGAMLARGE )
%IF X>ZERO %THEN %RESULT= GAMMAOFFXMAX %ELSE %RESULT= ZERO
%FINISH
RX= ONE
%IF TWO<= X<= THREE %THEN ->G2
%IF X<TWO %THEN ->G3
G1: X= X - ONE
RX= RX*X
%IF X<3 %THEN ->G2
->G1
G3: %IF MOD(X)<XSMALL %START
%IF Y>XSMALL %START
ERROR( DGAMNEGINT )
%RESULT= GAMMAOFFXMAX
%FINISHELSESTART
%IF Y<XVSMALL %START
!" ** DGAMMA ARG TOO NEAR ZERO OR NEG INTEGER."
ERROR( DGAMNEGINT )
%IF X>=ZERO %THEN %RESULT= ONE/XVSMALL %ELSE %RESULT= -ONE/XVSMALL
%FINISHELSE %RESULT= ONE/X
%FINISH
%FINISH
RX= RX/X
X= X + ONE
%IF X<TWO %THEN ->G3
G2: G= -0.00000051132627267
Z= X - TWO
! %CYCLE K= 15,-1,1
! G= G*RT+A(K)
! %REPEAT
G= ((((((((((((((G*Z + A(15))*Z + A(14))*Z + A(13))*Z + A(12)) %C
*Z + A(11))*Z + A(10))*Z + A(9))*Z + A(8))*Z + A(7))*Z + A(6))%C
*Z + A(5))*Z + A(4))*Z + A(3))*Z + A(2))*Z + A(1)
%RESULT= RX*G
%END {OF DGAMMA}
!*
!******************************************************************************
%EXTERNALLONGREALFN DLGAMA %ALIAS "F_DLGAMMA"(%LONGREALNAME ARG)
%LONGREAL G,T,X,Y
%INTEGER I,M
%CONSTANTLONGREAL %C
XSMALL= 1.0@-17, XBIG = 7.7@7, LNR2PI= 9.18938533204672742@-1
! RANGE DEPENDENT CONSTANTS
! FOR IBM 360/370 AND SIMILAR MACHINES
%CONSTANTLONGREAL %C
XVBIG= 4.29@+73, GBIG = 7.231@+75
! FOR DEC10, HONEYWELL, UNIVAC 1100 (S.P.)
!R2 %CONSTANTREAL %C
!R2 XVBIG= 2.05@36, GBIG= 1.69@38
! FOR CDC 7600/CYBER
!R4 %CONSTANTLONGREAL %C
!R4 XVBIG= 1.72@+319, GBIG= 1.26@+322
! FOR UNIVAC 1100 (D.P.)
!R5 %CONSTANTLONGREAL %C
!R5 XVBIG= 1.28@305, GBIG 8.98@+307
!
X= ARG
%IF X<=XSMALL %START {VERY SMALL RANGE}
%IF X>ZERO %THEN %RESULT= -DLOG(X) %ELSESTART
ERROR( DLGAMNEG )
%RESULT= ZERO
%FINISH
%FINISH
%IF X<=15.0 %START {MAIN SMALL X RANGE}
M = INTPT(X) {X IS POSITIVE}
T= X - FLOAT(M)
M= M - 1
G= ONE
%IF M#0 %START
%IF M<0 %THEN G= G/X %ELSESTART
%CYCLE I= 1,1,M
G= G*(X - FLOAT(I))
%REPEAT
%FINISH
%FINISH
T= TWO*T - ONE
! EXPANSION (0026) EVALUATED AS Y(T) --PRECISION 17E.18
Y= (((((((-1.46381209600000000@-11*T %C
+4.26560716800000000@-11)*T -4.01499750400000000@-11)*T %C
+1.27679856640000000@-10)*T -6.13513953280000000@-10)*T %C
+1.82243164160000000@-9)*T -5.11961333760000000@-9)*T %C
+1.53835215257600000@-8)
Y= ((((((( Y*T -4.64774927155200000@-8)*T %C
+1.39383522590720000@-7)*T -4.17808776355840000@-7)*T %C
+1.25281466396672000@-6)*T -3.75499034136576000@-6)*T %C
+1.12524642975590400@-5)*T -3.36375833240268800@-5)
Y= (((((((( Y*T %C
+1.00928148823365120@-4)*T -2.96890121633200000@-4)*T %C
+9.15785997288933120@-4)*T -2.42259538436268176@-3)*T %C
+9.04033494028101968@-3)*T -1.34118505705967765@-2)*T %C
+1.03703363422075456@-1)*T +1.61691987244425092@-2)*T %C
+8.86226925452758013@-1
Y= Y*G
%RESULT= DLOG(Y)
%FINISHELSESTART
%IF X<=XBIG %START {MAIN LARGE X RANGE}
T= 450.0/(X*X) - ONE
! EXPANSION (0059) EVALUATED AS Y(T) --PRECISION 08E.09
! EXPANSION (0059) EVALUATED AS Y(T) --PRECISION 17E.18
Y= (((( -9.94561064728159347@-17*T +2.00201927337982364@-14)*T %C
-6.45101975779653651@-12)*T +3.89978899876484712@-9)*T %C
-6.16502049453716986@-6)*T +8.33271644065786580@-2
%RESULT= (X-HALF)*DLOG(X) - X + LNR2PI + Y/X
%FINISHELSESTART
%IF X<=XVBIG %THEN {ASYMPTOTIC LARGE X RANGE} %C
%RESULT= (X-HALF)*DLOG(X) - X + LNR2PI %ELSESTART
ERROR( DLGAMLARGE )
%RESULT= GBIG
%FINISH
%FINISH
%FINISH
%END {OF DLGAMA}
!******************************************************************************
%EXTERNALLONGREALFN DERF %ALIAS "F_DERF"(%LONGREALNAME ARG)
! ERF(X)
%INTEGER J
%LONGREAL XV,X,X2,BJP2,BJP1,BJ
!
! *** IMPLEMENTATION DEPENDENT FEATURES ***
!
!
%CONSTANTLONGREAL XUP= 6.25@0, SQRTPI= 1.7724538509055160@0
%CONSTANTLONGREALARRAY C(1:18)= %C
1.9449071068178803@0,4.20186582324414@-2,-1.86866103976769@-2,
5.1281061839107@-3,-1.0683107461726@-3,1.744737872522@-4,
-2.15642065714@-5,1.7282657974@-6,-2.00479241@-8,
-1.64782105@-8,2.0008475@-9,2.57716@-11,-3.06343@-11,
1.9158@-12,3.703@-13,-5.43@-14,-4.0@-15,1.2@-15
%CONSTANTLONGREALARRAY D(1:17)= %C
1.4831105640848036@0,-3.010710733865950@-1,6.89948306898316@-2,
-1.39162712647222@-2,2.4207995224335@-3,-3.658639685849@-4,
4.86209844323@-5,-5.7492565580@-6,6.113243578@-7,
-5.89910153@-8,5.2070091@-9,-4.232976@-10,3.18811@-11,
-2.2361@-12,1.467@-13,-9.0@-15,5.0@-16
!
!
! NO FAILURE EXITS
X= ARG
XV= MOD(X)
%IF XV<XUP %START
!ELSE GO TO 120
%IF XV>TWO %START
! ELSE GO TO 60
X2= TWO - TWENTY/(XV+THREE)
!
! SUMMATION
BJ = X2*C(18) + C(17)
BJP2= X2*BJ - C(18) + C(16)
BJP1= X2*BJP2 - BJ + C(15)
BJ = X2*BJP1 - BJP2 + C(14)
BJP2= X2*BJ - BJP1 + C(13)
BJP1= X2*BJP2 - BJ + C(12)
BJ = X2*BJP1 - BJP2 + C(11)
BJP2= X2*BJ - BJP1 + C(10)
BJP1= X2*BJP2 - BJ + C(9)
BJ = X2*BJP1 - BJP2 + C(8)
BJP2= X2*BJ - BJP1 + C(7)
BJP1= X2*BJP2 - BJ + C(6)
BJ = X2*BJP1 - BJP2 + C(5)
BJP2= X2*BJ - BJP1 + C(4)
BJP1= X2*BJP2 - BJ + C(3)
BJ = X2*BJP1 - BJP2 + C(2)
BJP2= X2*BJ - BJP1 + C(1)
X2= HALF*(BJP2-BJP1)/XV*DEXP(-X*X)/SQRTPI
%IF X>ZERO %THEN %RESULT= ONE - X2 %ELSE %RESULT= X2 - ONE
!
%FINISHELSESTART
X2= X*X - TWO
! SUMMATION
BJP2= X2*D(17) + D(16)
BJP1= X2*BJP2 - D(17) + D(15)
BJ = X2*BJP1 - BJP2 + D(14)
BJP2= X2*BJ - BJP1 + D(13)
BJP1= X2*BJP2 - BJ + D(12)
BJ = X2*BJP1 - BJP2 + D(11)
BJP2= X2*BJ - BJP1 + D(10)
BJP1= X2*BJP2 - BJ + D(9)
BJ = X2*BJP1 - BJP2 + D(8)
BJP2= X2*BJ - BJP1 + D(7)
BJP1= X2*BJP2 - BJ + D(6)
BJ = X2*BJP1 - BJP2 + D(5)
BJP2= X2*BJ - BJP1 + D(4)
BJP1= X2*BJP2 - BJ + D(3)
BJ = X2*BJP1 - BJP2 + D(2)
BJP2= X2*BJ - BJP1 + D(1)
%RESULT= HALF*(BJP2-BJP1)*X
!
%FINISH
%FINISHELSESTART
%IF X>ZERO %THEN %RESULT= ONE %ELSE %RESULT= -ONE
%FINISH
%END {OF DERF}
!******************************************************************************
%EXTERNALLONGREALFN DERFC %ALIAS "F_DERFC"(%LONGREALNAME ARG)
!
! COMPLEMENT OF ERROR FUNCTION ERFC(X)
!
! *** IMPLEMENTATION DEPENDENT FEATURES ***
%LONGREAL X,Y,T
!
! PRECISION DEPENDENT CONSTANTS
%CONSTANTLONGREAL XHI= 13.0@0, XLO= -6.25@0
!
! TEST EXTREME EXITS
X= ARG
%IF X>=XHI %THEN %RESULT= ZERO
%IF X<=XLO %THEN %RESULT= TWO
!
! EXPANSION ARGUMENT
T= ONE - 7.5@0/(MOD(X)+3.75@0)
!
!
!
! EXPANSION (0021) EVALUATED AS Y(T) --PRECISION 16E
Y= +1.455897212750385@-1+T*( -2.734219314954260@-1+ %C
T*( +2.260080669166197@-1+T*( -1.635718955239687@-1+ %C
T*( +1.026043120322792@-1+T*( -5.480232669380236@-2+ %C
T*( +2.414322397093253@-2+T*( -8.220621168415435@-3+ %C
T*( +1.802962431316418@-3+T*( -2.553523453642242@-5+ %C
T*( -1.524627476123466@-4+T*( +4.789838226695987@-5+ %C
T*( +3.584014089915968@-6+T*( -6.182369348098529@-6+ %C
T*( +7.478317101785790@-7+T*( +6.575825478226343@-7+ %C
T*( -1.822565715362025@-7+T*( -7.043998994397452@-8+ %C
T*( +3.026547320064576@-8+T*( +7.532536116142436@-9+ %C
T*( -4.066088879757269@-9+T*( -5.718639670776992@-10+ %C
T*( +3.328130055126039@-10))))))))))))))))))))))
!
Y= DEXP(-X*X)*Y
%IF X<ZERO %THEN %RESULT= TWO - Y %ELSE %RESULT= Y
%END {OF DERFC}
!******************************************************************************
!******************************************************************************
%EXTERNALINTEGERFN POWII %ALIAS "F_POWII"(%INTEGER IARG,NARG)
!
! CALCULATES IARG**NARG
!
%INTEGER X,Y
%INTEGER N
X= IARG
Y= 1
N= NARG
%IF N=0 %START
%IF X#0 %THEN %RESULT= 1 %ELSESTART
ERROR( IPOWERZERO )
%RESULT= 0
%FINISH
%FINISH
%IF X=1 %THEN %RESULT= 1
%IF X=-1 %START
%IF N<0 %THEN N= -N
%IF N&1=1 %THEN %RESULT= -1 %ELSE %RESULT= 1
%FINISH
%IF N<0 %START
ERROR( IPOWEXPNEG )
%RESULT= 0
%FINISH
%IF X=0 %THEN %RESULT= 0
%CYCLE
%IF N&1#0 %THEN Y= Y*X
N = N>>1
%IF N#0 %THEN X= X*X %ELSE %EXIT
%REPEAT
%IF Y=0 %START
ERROR( IPOWLARGE )
%RESULT= 0
%FINISHELSE %RESULT= Y
%END {OF POWII}
!******************************************************************************
%ENDOFFILE