!****************************************************************************** !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) !(38 and 39 have been suppressed for conformity) !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 IS PUT INTO TEST MODE BY SETTING THE INTEGER 'TEST' TO 1. !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 TEST= 0, 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, not used} IPOWERZERO = 39, { INTEGER ZERO RAISED TO NON-POSITIVE POWER, not used} 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= -128, 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 TERROR (%INTEGERNAME TYPE) %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 %IF TEST=1 %THEN TERROR(DSQRTNEG) %ELSE 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 %IF TEST=1 %THEN TERROR( DLOGSMALL ) %ELSE ERROR( DLOGSMALL ) %IF 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))<>(32-Q)) INTEGER(ADDR(IN)+4)=(INTEGER(ADDR(IN)+4)<174.6731 %START %IF IN<-180.2182 %THENRESULT= ZERO %IF IN>174.6731 %START %IF TEST=1 %THEN TERROR(DEXPLARGE) %ELSE ERROR( DEXPLARGE ) %RESULT= GREATEST %FINISH %FINISH %IF BYTEINTEGER(ADDR(W))=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 %IF TEST=1 %THEN TERROR( DPOWERNEG ) %ELSE ERROR( DPOWERNEG ) X= -X %FINISHELSESTART %IF Y<=ZERO %START %IF Y=ZERO %START %IF TEST=1 %THEN TERROR( DPOWERZERO ) %C %ELSE ERROR( DPOWERZERO ) %FINISHELSESTART %IF TEST=1 %THEN TERROR( DPOWERNEG ) %ELSE 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))<>(32-Q)) INTEGER(ADDR(G)+4)=(INTEGER(ADDR(G)+4)<>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 SAFETYLAB2: !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 SAFETYLAB3: IW1= INTPT(W1) !ALLOW OVERFLOW TO CAUSE TERMINATION HERE. !%IF IW1>IW1MAX %START !PRINTSTRING(" IW1 IS >4030.") %AND NEWLINES(2) !%IF IW1ZERO %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 %IF TEST=1 %THEN TERROR( DPOWERLARGE ) %ELSE ERROR( DPOWERLARGE ) ! NOTE THAT THERE SHOULD REALLY BE A EXCEPTION FOR DPOWEROVER. %RESULT= GREATEST %FINISH %IF N<0 %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 %IF TEST=1 %THEN TERROR( DPOWERZERO ) %ELSE ERROR( DPOWERZERO ) %RESULT= ZERO %FINISH %FINISH %IF ARG=ZERO %START %IF NARG>0 %THEN %RESULT= ZERO %ELSESTART %IF TEST=1 %THEN TERROR( DPOWERZERO ) %ELSE 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 %IF TEST=1 %THEN TERROR( DSINLARGE ) %ELSE 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))SINLIMIT %START %IF TEST=1 %THEN TERROR( DCOSLARGE ) %ELSE 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 %IF TEST=1 %THEN TERROR( DTANLARGE ) %ELSE ERROR( DTANLARGE ) %END {OF TANARGBIG, LOCAL} %ROUTINE COTANARGSMALL %IF TEST=1 %THEN TERROR( DCOTANSMALL ) %ELSE 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) 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 XSETGREATEST %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)HALF %START FLAG= 1 - FLAG %IF Y>ONE %START %IF TEST=1 %THEN TERROR( DASINLARGE ) %ELSE 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 YHALF %START FLAG= 1 - FLAG %IF Y>ONE %START %IF TEST=1 %THEN TERROR( DACOSLARGE ) %ELSE 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 YSETQUADRANT %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 ARGONE %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)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 %IF TEST=1 %THEN TERROR( DARCTAN2ZERO )%ELSE 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 %IF TEST=1 %THEN TERROR( DSINHLARGE ) %ELSE 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 ARGONE %START W= Y - LNV %IF W>MAXEXPARG %START %IF TEST=1 %THEN TERROR( DCOSHLARGE ) %ELSE 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 FFXMAX %START !" ABSOLUTE VALUE OF DGAMMA ARG > 57." %IF TEST=1 %THEN TERROR( DGAMLARGE ) %ELSE ERROR( DGAMLARGE ) %IF X>ZERO %THEN %RESULT= GAMMAOFFXMAX %ELSE %RESULT= ZERO %FINISH RX= ONE %IF TWO<= X<= THREE %THEN ->G2 %IF XG3 G1: X= X - ONE RX= RX*X %IF X<3 %THEN ->G2 ->G1 G3: %IF MOD(X)XSMALL %START %IF TEST=1 %THEN TERROR(DGAMNEGINT) %ELSE ERROR( DGAMNEGINT ) %RESULT= GAMMAOFFXMAX %FINISHELSESTART %IF Y=ZERO %THEN %RESULT= ONE/XVSMALL %ELSE %RESULT= -ONE/XVSMALL %FINISHELSE %RESULT= ONE/X %FINISH %FINISH RX= RX/X X= X + ONE %IF XG3 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 %IF TEST=1 %THEN TERROR( DLGAMNEG ) %ELSE 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 %IF TEST=1 %THEN TERROR( DLGAMLARGE ) %ELSE 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 XVTWO %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>1 %IF N#0 %THEN X= X*X %ELSE %EXIT %REPEAT %IF Y=0 %START %IF TEST=1 %THEN TERROR( IPOWLARGE ) %ELSE ERROR( IPOWLARGE ) %RESULT= 0 %FINISHELSE %RESULT= Y %END {OF POWII} !****************************************************************************** %ENDOFFILE