!******************************************************************************
!SINGLE PRECISION MATHS LIBRARY AFTER CODY AND WAITE,
!SOFTWARE MANUAL FOR THE ELEMENTARY FUNCTIONS, PRENTICE HALL 1980.
!WRITTEN BY D.B.TAYLOR, JANUARY 1986. CODE FOR IEEE SINGLE PRECISION


!THE FLOATING POINT RECIPEES HAVE BEEN SELECTED AND INDEPENDANT ROUTINES 
!HAVE BEEN WRITTEN FOR EACH ENTRY WITH THE EXCEPTION OF ALOG10 WHICH CALLS ALOG,
! ATAN2 WHICH CALLS  ATAN,  ACOS AND  ACOS WHICH CALL SQRT, AND  SINH,
! COSH AND  TANH WHICH CALL  EXP.

!THE ERROR EXITS TO THE DEFINED ERROR EXCEPTION NUMBERS ARE IN LINE WITH
!FORTRAN STANDARDS EXCEPT THAT THE  EXP FUNCTION DOES NOT FLAG UNDERFLOW
!AND ZERO**POSITIVE IS ACCEPTED BY POWRR. 

!THE LIMITS FOR THE ARGUEMENTS ARE THOSE DESCRIBED BY CODY AND WAITE,
!SUITABLY INTERPRETED FOR BYTE ARITHEMETIC, WITH THE EXCEPTION THAT THE RANGES
!FOR  SIN/ COS AND  TAN/ COTAN ARE EXTENDED TO (-5.274E7 TO 5.274E7)
!AND (-2.633E7 TO 2.633E7) RESPECTIVELY.
!AT THE EXTREMAL VALUES ONLY ONE HEX DIGIT IS ACCURATE. ONLY WITHIN
!THE RESTRICTED RANGES, (-12867 TO 12867 ) FOR  SIN/ COS AND 
!(-6433 TO 6433) FOR  TAN/ COTAN 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 DOUBLE PRECISION (LONG) ARITHMETIC IN ALOG, EXP, POWRR,
!POWRI, SIN, COS, TAN AND COTAN. ON AMDAHL THE VALUE OF THE INTEGER
!CONSTANT 'IMPROUNDS' IS SET TO 1 TO IMPLY THAT THE IMP CONVERSION FROM
!LONGREAL TO REAL 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 LONGREAL TO REAL.
!******************************************************************************
!THE FOLLOWING FUNCTIONS ARE PROVIDED ON THIS FILE. THE NUMBERED GROUPS
!ARE TESTED BY THE TEST SUBROUTINE STEST'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 %EXTERNALREALFN  SQRT  %ALIAS "F_SQRT" (%REALNAME ARG)   
! 2 %EXTERNALREALFN ALOG   %ALIAS "F_LOG"   (%REALNAME ARG)
! 2 %EXTERNALREALFN ALOG10 %ALIAS "F_LOG10" (%REALNAME ARG)
! 3 %EXTERNALREALFN  EXP   %ALIAS "F_EXP"  (%REALNAME ARG)
! 4 %EXTERNALREALFN POWRR  %ALIAS "F_POWRR" (%REAL ARG1,ARG2)     
!   %EXTERNALREALFN POWRI  %ALIAS "F_POWRI" %C
!                                          (%REAL A,%INTEGER N)
! 5 %EXTERNALREALFN  SIN   %ALIAS "F_SIN"  (%REALNAME ARG)
! 5 %EXTERNALREALFN  COS   %ALIAS "F_COS"  (%REALNAME ARG)
! 6 %EXTERNALREALFN  TAN   %ALIAS "F_TAN"  (%REALNAME ARG)
! 6 %EXTERNALREALFN  COTAN %ALIAS "F_COT"(%REALNAME ARG)
! 7 %EXTERNALREALFN  ASIN  %ALIAS "F_ASIN" (%REALNAME ARG)
! 7 %EXTERNALREALFN  ACOS  %ALIAS "F_ACOS" (%REALNAME ARG)
! 8 %EXTERNALREALFN  ATAN  %ALIAS "F_ATAN" (%REALNAME ARG)
! 8 %EXTERNALREALFN  ATAN2 %ALIAS "F_ATAN2"(%REALNAME ARGY,ARGX)
! 9 %EXTERNALREALFN  SINH  %ALIAS "F_SINH" (%REALNAME ARG)
! 9 %EXTERNALREALFN  COSH  %ALIAS "F_COSH" (%REALNAME ARG)
!10 %EXTERNALREALFN  TANH  %ALIAS "F_TANH" (%REALNAME ARG)
!11 %EXTERNALREALFN  GAMMA %ALIAS "F_GAMMA"(%REALNAME  FX)
!   %EXTERNALREALFN ALGAMA %ALIAS "F_LGAMMA"(%REALNAME ARG)
!   %EXTERNALREALFN  ERF   %ALIAS "F_ERF"  (%REALNAME ARG)
!   %EXTERNALREALFN  ERFC  %ALIAS "F_ERFC" (%REALNAME ARG)
!******************************************************************************
%CONSTANTINTEGER TEST= 0, IMPROUNDS= 1
!SINGLE PRECISION ERROR NUMBERS
%CONSTANTINTEGER %C
 ALOGSMALL   = 21, { ALOG ARG NEGATIVE OR ZERO}
 SQRTNEG     = 18, { SQRT ARG NEGATIVE}
 EXPLARGE    = 19, { EXP ARG GREATER THAN 174.6731}
 EXPSMALL    = 20, { EXP ARG LESS THAN -180.2182, NOT USED}
 SINLARGE    = 26, { SIN ARG MAGNITUDE GREATER THAN 5.274E7}
 ASINLARGE   = 29, { ASIN ARG MAGNITUDE GREATER THAN 1.0}
 COSLARGE    = 27, { COS ARG MAGNITUDE GREATER THAN 5.274E7}
 ACOSLARGE   = 30, { ACOS ARG MAGNITUDE GREATER THAN 1.0}
 TANLARGE    = 28, { TAN/ COTAN ARG MAGNITUDE GREATER THAN 2.633E7}
 COTANSMALL  = 41, { COTAN ARG SMALL IN MAGNITUDE}
 ARCTAN2ZERO = 31, { ATAN2 ARGS ZERO}
 SINHLARGE   = 32, { SINH ARG MAGNITUDE GREATER THAN 175.3662 }
 COSHLARGE   = 33, { COSH ARG MAGNITUDE GREATER THAN 175.3662 }
 POWERNEG    = 23, { NEGATIVE S.P. VALUE RAISED TO A NON-INTEGER POWER}
 POWERZERO   = 22, { S.P. ZERO RAISED TO NON-POSITIVE POWER   }
 POWERBIG    = 24, { S.P. VALUES RAISED TO S.P. POWER TOO BIG}
 POWERSMALL  = 25, { S.P. VALUE TO S.P. POWER UNDERFLOWS. NOT USED}
 GAMLARGE    = 42, { GAMMA ARG MAGNITUDE GREATER THAN 57.0}
 GAMNEGINT   = 43, { GAMMA ARG NEAR ZERO OR NEGATIVE INTEGER}
 LGAMNEG     = 44, { ALGAMA ARG IS NEGATIVE}
 LGAMLARGE   = 45  { ALGAMA ARG IS GREATER THAN 4.29E73}

%CONSTANTINTEGER  EXPEXCESS= 64, MAXEXP= 127, MINEXP= -128, EPSEXP= -49

%CONSTANTREAL     GREATEST = R'7FFFFFFF',
                  SMALLEST = R'00100000',
                  MAXEXPARG= 174.6731,
                  MINEXPARG= -180.2182,
                  SINLIMIT = 5.274@7,
                  TANLIMIT = 2.633@7,
                  ONE      =  1.0,
                  TWO      =  2.0,
                  THREE    =  3.0,
                  TWENTY   =  20.0,
                  HALF     = R'40800000',
                  QUART    = R'40400000',
                  ZERO     =  0.0,
                  SQRTHALF = R'40B504F3',
                  RTEPS    = R'3E100000',
                  TRICK1   = R'47000000',
                  TRICK2   = R'46100000',
                  TRICK3   = R'46000000',
                  PI       = R'413243F7',
                  PIBY2    = R'411921FB',
                  PIBY3    = R'4110C152',
                  PIBY4    = R'40C90FDB',
                  PIBY6    = R'40860A92'

%CONSTANTLONGREAL %C
LZERO= 0.0,
LONE= 1.0,
LTWO= 2.0,
LTHREE= 3.0,
DPI= 3.141592653589793238462643,
LPIBY2= R'411921FB54442D18',
LLOGE2= R'40B17217F7D1CF7A'

!GREATEST= 7.2370055773322608@ 75,
!SQRTHALF= 7.0710678118654753@ -1,
!RTEPS= 2.44140625 @-4,
!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 SERROR (%INTEGERNAME TYPE)
%EXTERNALROUTINESPEC ERROR %ALIAS "F_MLE"(%INTEGER TYPE)

!******************************************************************************
%EXTERNALREALFN  SQRT %ALIAS "F_SQRT" (%REALNAME ARG)     
%CONSTANTREAL %C
C1= 1.161322, C2= 0.172924, C3= 0.175241

%REAL IN,Y
%INTEGER E
!
IN= ARG
   %IF IN<=0.0 %THENSTART
      IN= -IN
      %IF IN=ZERO %THENRESULT= ZERO
      %IF TEST=1 %THEN SERROR( SQRTNEG)  %ELSE ERROR( SQRTNEG )
      %RESULT= ZERO
   %FINISH

!     OBTAIN HEX EXPONENT AND FRACTION IN .0625 TO 1.0

E= BYTEINTEGER(ADDR(IN)) - EXPEXCESS
BYTEINTEGER(ADDR(IN))= EXPEXCESS

!     OBTAIN AN APPROX TO 2*SQRT AND PERFORM 2 NEWTON-RAPHSON STEPS,
!     THE LAST MODIFIED TO IMPROVE ROUNDING.


Y= C1 + IN - C2/(IN + C3)
Y= QUART*Y + IN/Y
Y= Y + HALF*(IN/Y -Y)

!     RESET EXPONENT OF Y TO CORRECT VALUE 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  SQRT}
!******************************************************************************
%EXTERNALREALFN ALOG %ALIAS "F_LOG"(%REALNAME ARG)

%CONSTANTREAL C1= R'40B18000',
C2= R'BDDE8083',
A0= R'C08D7E3D',
B0= R'C16A1F9D'

!C1=  0.69335 93750,
!C2= -0.21219 44400 54690 58277 @-3,
!A0= -0.55270 74855 @ 0,
!B0= -0.66327 18214 @ 1
!B1=  0.10000 00000 @ 1
!
%CONSTANTINTEGERARRAY QVAL(1:15)= 3,2,2,1,1,1,1,0,0,0,0,0,0,0,0         
!
%INTEGER P,Q
%REAL IN,W,Z
%LONGREAL ZZ,PRESULT
!
IN= ARG
%IF IN<= 0.0 %THENSTART
  %IF TEST=1 %THEN SERROR( ALOGSMALL ) %ELSE ERROR( ALOGSMALL )
  !PRINTSTRING(" ** LOG ARG NEG.") %AND NEWLINE
  %IF IN<0.0 %THEN IN= -IN %ELSE %RESULT= -GREATEST
%FINISH
!
P= BYTEINTEGER(ADDR(IN)) - EXPEXCESS
BYTEINTEGER(ADDR(IN))= 0
Q= QVAL(INTEGER(ADDR(IN))>>20)
INTEGER(ADDR(IN))=(INTEGER(ADDR(IN))<<Q)
BYTEINTEGER(ADDR(IN))= X'40'
%IF IN<SQRTHALF %START
  Z= IN - HALF
  W= Z*HALF + HALF
  PRESULT= FLOAT(4*P-Q-1)
%FINISHELSESTART
  Z= (IN - HALF) - HALF
  W= IN*HALF + HALF
  PRESULT= FLOAT(4*P-Q)
%FINISH

IN= Z/W
W= IN*IN
Z= (W*A0)/(W + B0)
Z= IN + (IN*Z)
!%RESULT= ((PRESULT*C2 + Z) + PRESULT*C1)
ZZ= PRESULT*LLOGE2 + Z
%IF IMPROUNDS=0 %START
  Z= ZZ
  ZZ= ZZ + (ZZ - Z)
%FINISH
%RESULT= ZZ
%END {OF ALOG}
!******************************************************************************
%EXTERNALREALFN ALOG10 %ALIAS "F_LOG10"(%REALNAME ARG)
%CONSTANTREAL C3= R'406F2DEC' {549B9439'}
!C3= 0.43429 44819 03251 82765
%RESULT= C3*ALOG(ARG)
%END {OF ALOG10}
!******************************************************************************
%EXTERNALREALFN  EXP %ALIAS "F_EXP"(%REALNAME ARG)

!
%CONSTANTREAL %C
C3= R'41171547',
P0= R'40400000',
P1= R'3F110A61',
Q0= R'407FFFFE',
Q1= R'3FCCBF5B'

%CONSTANTREAL C1= R'40B00000', C2= R'3F17217F'
!
!C3= 0.14426 95040 88896 341 @1,
!P0 = 0.24999 99995 0 @ 0,
!P1 = 0.41602 88626 8 @-2,
!Q1 = 0.49987 17877 8 @-1

%CONSTANTREALARRAY POW2(1:3)= 2.0,4.0,8.0

%INTEGER P,Q,N
%LONGREAL GG,XXN
%REAL IN,D,Z,W,XN,G

IN= ARG
%IF IN<-180.2182 %THENRESULT= 0.0
%IF IN>174.6731 %THENSTART
  !PRINTSTRING(" ** EXP ARG > 174.6731") %AND NEWLINE
  %IF TEST=1 %THEN SERROR( EXPLARGE) %ELSE ERROR(  EXPLARGE )
  %RESULT= GREATEST
%FINISH

P= (X'7F' & BYTEINTEGER(ADDR(IN)))
%IF P<49 %THEN %RESULT= 1.0
N= INT(IN*C3)
XXN= FLOAT(N)
N= N + 1
%IF N>=0 %THEN P= N>>2 %ELSE P= (N>>2) ! X'C0000000'
!Q= N - 4*P
Q= N & X'00000003'
!%IF N>=0 %THEN P= N//4 %ELSE P= (N-3)//4
!Q= N - 4*P
!
! SINGLE PREC OPTION
!G= (IN - XN*C1) - XN*C2
!
!DOUBLE PREC OPTION
GG= IN
GG= GG - lLOGE2*XXN
G= GG
! ROUNDING CODE OMITTED ON AMDAHL
%IF IMPROUNDS=0 %START
  GG= GG + (GG - G)
  G= GG
%FINISH
!
Z= G*G
W= (P1*Z + P0)*G
D= Q1*Z + HALF
W= HALF + (W/(D - W))
%IF Q>0 %THEN W= W*POW2(Q)
! LABEL NEEDED TO GET CORRECT RESULTS ??
LABELFORCOMPILER:
BYTEINTEGER(ADDR(W))= BYTEINTEGER(ADDR(W)) + P
%RESULT= W
%END {OF  EXP}

!******************************************************************************
%EXTERNALREALFN POWRR %ALIAS "F_POWRR"(%REAL ARG1,ARG2)

!  ARRAY OF 2**((1-J)/16), J= 1,2,...,17 ROUNDED TO WORKING PRECISION

%CONSTANTREALARRAY A1(1:17)= %C
R'41100000',
R'40F5257D',
R'40EAC0C7',
R'40E0CCDF',
R'40D744FD',
R'40CE248C',
R'40C5672A',
R'40BD08A4',
R'40B504F3',
R'40AD583F',
R'40A5FED7',
R'409EF532',
R'409837F0',
R'4091C3D3',
R'408B95C2',
R'4085AAC3',
R'40800000'

!  ARRAY 2**((1-2*J)/16) - A1(2*J), J= 1,2,...,8. I.E. CORRECTION TO A1.

%CONSTANTREALARRAY A2(1:8)= %C
R'3A152486',
R'BA13D56B',
R'3A151F84',
R'BA60A7F3',
R'BA15BD5E',
R'3A6091A1',
R'3A73AB11',
R'3A67CC48'


%CONSTANTREAL K= R'40715476'
!K= LOG2 - 1 =0.44269504088896341

%CONSTANTREAL  %C
P1= 0.83357 541 @-1,
Q1= 0.69314 675 @ 0,
Q2= 0.24018 510 @ 0,
Q3= 0.54360 383 @-1

%CONSTANTREALARRAY 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
%REAL D,Z,R,V,G,X,Y,U1,U2,W,W1,W2,Y1,Y2
%LONGREAL UU1,UU2,WW

X= ARG1; Y= ARG2
%IF X<=ZERO %START
   %IF X<>ZERO %START
      %IF TEST=1 %THEN SERROR( POWERNEG ) %ELSE ERROR( POWERNEG )
      X= -X
   %FINISHELSESTART
      %IF Y<=ZERO %START
         %IF Y=ZERO %START
            %IF TEST=1 %THEN SERROR( POWERZERO ) %C
                       %ELSE ERROR( POWERZERO )
         %FINISHELSESTART
            %IF TEST=1 %THEN SERROR( POWERNEG ) %ELSE ERROR( POWERNEG)
         %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)
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= 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= ((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 SERROR(  POWERBIG ) %ELSE ERROR(  POWERBIG )
%RESULT= GREATEST
%FINISH
%IF N<0 %START
!PRINTSTRING(" EXPONENTIATION UNDERFLOW."); NEWLINES(2)
!%IF TEST=1 %THEN SERROR( POWERSMALL ) %ELSE ERROR ( POWERSMALL )       
%RESULT= ZERO
%FINISH
BYTEINTEGER(ADDR(Z))= N
%RESULT= Z
%END {OF POWRR}

!******************************************************************************
%EXTERNALREALFN POWRI %ALIAS "F_POWRI" %C
         (%REAL ARG, %INTEGER NARG)
%CONSTANTLONGREAL LONE= 1.0
!
!     CALCULATES ARG**NARG
!
%LONGREAL XX,YY
%REAL 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 SERROR( POWERZERO ) %ELSE ERROR( POWERZERO )
    %RESULT= ZERO
  %FINISH
%FINISH

%IF ARG=ZERO %START
  %IF NARG>0 %THEN %RESULT= ZERO %ELSESTART
    %IF TEST=1 %THEN SERROR( POWERZERO ) %ELSE ERROR( POWERZERO )
    %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 FUNCION FOR REAL EXPONENT IF MAGNITUDE OF EXPONENT > 64
!
%IF N>64 %START
  %IF ARG>=ZERO %THEN %RESULT=  POWRR( ARG,FLOAT(NARG)) %ELSESTART
    %IF N&1#1   %THEN %RESULT=  POWRR(-ARG,FLOAT(NARG)) %C
                %ELSE %RESULT= -POWRR(-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 POWRI}
!******************************************************************************
!  CONSTANTS FOR  SIN AND  COS. NOTE THAT RANGE IS EXTENDED INTO 'GRAINY'
!  AREA. SEE NOTE IN INTRODUCTION.


!CODY AND WAITE USE SINLIMIT =  12867 = pi*16^3

%CONSTANTREAL %C
              SC1   =  3.14062500,
              SC2   =  9.67653701@-4,
              PIINV =  0.318309903,
              R1T4  =  R'C0AAAAB0',
              R2T4  =  R'3F888798',
              R3T4  =  R'BE33ECE8',
              R4T4  =  R'3CAE92B4'
!
!     THESE ARE ESSENTIALLY R1,R2,R3 AND R4 TIMES 4 WITH MASSAGING
!     OF R1T4 AND R2T4 TO MINIMIZE RMS RELATIVE ERROR.
!
!     CODY AND WAITE USE
!      -0.166666567, 0.833302457@-2,-0.198074180@-3, 0.260190313@-5
!     THESE ARE ROUNDED IN HEX. REMEZ ALG FOR MIN REL ERROR AND
!     NO ERROR AT 0 AND PI/2 GIVES
!      -0.166666588, 0.833304665@-2,-0.198080850@-3, 0.26022725@-5
!
!
%EXTERNALREALFUNCTION SIN %ALIAS "F_SIN" (%REALNAME ARG)
%REAL XN,F,G
%LONGREAL FF
%INTEGER N

F= MOD(ARG)
%IF F>SINLIMIT %START
  %IF TEST=1 %THEN SERROR(  SINLARGE ) %ELSE ERROR (  SINLARGE )        
%RESULT= ZERO
%FINISH
!
N= INT(F*PIINV)
XN= FLOAT(N)
FF= F
FF= FF - DPI*XN
%IF IMPROUNDS=0 %START
  F= FF
  FF= FF + (FF - F)
%FINISH
%IF N&1=0 %THEN F= FF %ELSE F= -FF
!     NOTE AMDAHL IMP COMPILER ROUNDS HERE.
%IF MOD(F)>=RTEPS %THENSTART
G= F*F
F= F + QUART*F*G*(((R4T4*G + R3T4)*G + R2T4)*G + R1T4)               
%FINISH
%IF BYTEINTEGER(ADDR(ARG))<X'80' %THENRESULT= F %ELSERESULT= -F
%END {OF SIN}
!
%EXTERNALREALFUNCTION COS %ALIAS "F_COS" (%REALNAME ARG)

%REAL XN,F,G,Y
%LONGREAL FF
%INTEGER N

F= MOD(ARG)
Y= F + PIBY2
%IF Y>SINLIMIT %START
 %IF TEST=1 %THEN SERROR(  COSLARGE ) %ELSE ERROR(  COSLARGE )
%RESULT= ZERO
%FINISH
!
N= INT(Y*PIINV)
XN= FLOAT(N) - HALF
!
!
FF= F
FF= FF - DPI*XN
F= FF
!     NOTE AMDAHL IMP COMPILER ROUNDS HERE.
%IF IMPROUNDS=0 %START
  FF= FF + (FF - F)
  F= FF
%FINISH
%IF MOD(F)>=RTEPS %THENSTART
G= F*F
F= F + QUART*F*G*(((R4T4*G + R3T4)*G + R2T4)*G + R1T4)               
%FINISH
%IF N&1=0 %THENRESULT= F %ELSERESULT= -F
%END {OF  COS}
!******************************************************************************
! CONSTANTS FOR  TAN AND  COTAN

%CONSTANTREAL %C
TOVPI  = R'40A2F983',
TCC1   = R'41192000',
TCC2   = R'3E1FB544'

!TOVPI=  0.63661 97723 67581 34308,
!TCC1   =  1.57031 25,
!TCC2   =  4.83826 79489 7 @-4

%CONSTANTREAL %C
TCP1   = -0.95801 7723 @ -1,
TCQ1   = -0.42913 5777 @  0,
TCQ2   =  0.97168 5835 @ -2


!******************************************************************************
%ROUTINE TANARGBIG
%IF TEST=1 %THEN SERROR(  TANLARGE ) %ELSE ERROR(  TANLARGE )
%END {OF TANARGBIG, LOCAL}
%ROUTINE COTANARGSMALL
%IF TEST=1 %THEN SERROR(  COTANSMALL ) %ELSE ERROR(  COTANSMALL )
%END {OF COTANARGSMALL, LOCAL}

!******************************************************************************
%EXTERNALREALFN  TAN   %ALIAS "F_TAN"  (%REALNAME ARG)

!   CODY AND WAITE USE TANLIMIT =  6433 = (pi/2)*16^3

%LONGREAL XX
%REAL X,F,G,W,Z,XN
%INTEGER N

X= ARG
%IF MOD(X)>TANLIMIT %THEN TANARGBIG %AND ->SETGREATEST
N= INT(X*TOVPI)
XN= FLOAT(N)

XX= X
XX= XX - XN*LPIBY2
F= XX
%IF IMPROUNDS=0 %START
  XX= XX + (XX - F)
  F = XX
%FINISH
%IF MOD(F)<RTEPS %THEN W= F %AND Z= ONE %C
%ELSESTART
   G= F*F
   W= F + F*G*TCP1
   Z= (TCQ2*G + TCQ1)*G + ONE
%FINISH

%IF N&1#0 %THEN %RESULT= -Z/W %ELSE %RESULT= W/Z

SETGREATEST: %RESULT= GREATEST
%END {OF  TAN}
!******************************************************************************

%EXTERNALREALFN  COTAN %ALIAS "F_COT"(%REALNAME ARG)

%CONSTANTREAL EPS1= R'02100001'
!EPS1= 0.1381786...= R'02100001', JUST > 1/GREATEST

!   CODY AND WAITE USE TANLIMIT =  4.216574280 @ 8

%LONGREAL XX
%REAL X,F,G,W,Z,XN
%INTEGER N

X= ARG
W= MOD(X)
%IF W<EPS1 %THEN COTANARGSMALL %AND ->SETGREATEST
%IF W>TANLIMIT %THEN TANARGBIG %AND ->SETGREATEST
N= INT(X*TOVPI)
XN= FLOAT(N)

XX= X
XX= XX - XN*LPIBY2
F= XX
%IF IMPROUNDS=0 %START
  XX= XX + (XX - F)
  F = XX
%FINISH
%IF MOD(F)<RTEPS %THEN W= F %AND Z= ONE %C
%ELSESTART
   G= F*F
   W= F + F*G*TCP1
   Z= (TCQ2*G + TCQ1)*G + ONE
%FINISH

%IF N&1#0 %THEN %RESULT= -W/Z %ELSE %RESULT= Z/W
SETGREATEST: %RESULT= GREATEST
%END {OF  COTAN}
!******************************************************************************
!   CONSTANTS FOR  ASIN AND  ACOS

%CONSTANTREAL %C
ASCP1=  0.93393 5835 @ 0,
ASCP2= -0.50440 0557 @ 0,
ASCQ0=  0.56036 3004 @+1,
ASCQ1= -0.55484 6723 @+1

!******************************************************************************
%EXTERNALREALFN  ASIN %ALIAS "F_ASIN"(%REALNAME ARG)

%CONSTANTREAL A1= PIBY4
%REAL G,Y,W
%INTEGER FLAG

FLAG= 0
Y= MOD (ARG)
%IF Y>HALF %START
  FLAG= 1 - FLAG
  %IF Y>ONE %START
    %IF TEST=1 %THEN SERROR(  ASINLARGE ) %ELSE ERROR( ASINLARGE)
    %RESULT= GREATEST
  %FINISH
  G= (HALF - Y) + HALF
  G= G*HALF
  Y=  SQRT(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= (ASCP2*G + ASCP1)*G
G= (G + ASCQ1)*G + ASCQ0
W= Y + (Y*W/G)
%IF FLAG= 1 %THEN W= (A1 + W) + A1
%IF ARG<0.0 %THEN W= -W
%RESULT= W
%END {OF ASIN}

!******************************************************************************
%EXTERNALREALFN  ACOS %ALIAS "F_ACOS"(%REALNAME ARG)
    
%CONSTREAL A1= PIBY4, B0= PIBY2

%REAL G,Y,W
%INTEGER FLAG

FLAG= 1
Y= MOD(ARG)
%IF Y>HALF %START
   FLAG= 1 - FLAG
   %IF Y>ONE %START
       %IF TEST=1 %THEN SERROR(  ACOSLARGE ) %ELSE ERROR ( ACOSLARGE)
       %RESULT= GREATEST
   %FINISH
   G= (HALF - Y) + HALF
   G= G*HALF
   Y=  SQRT(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= (ASCP2*G + ASCP1)*G
G= (G + ASCQ1)*G + ASCQ0
W= Y + (Y*W/G)

SETQUADRANT:
%IF ARG<0.0 %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  ACOS}
!******************************************************************************
%EXTERNALREALFN  ATAN %ALIAS "F_ATAN"(%REALNAME ARG)
%CONSTANTREAL %C
ROOT3= R'411BB67B',
RT3M1= R'40BB67AF',
TMRT3= R'40449851'

!ROOT3=  1.73205 08075 68877 29353,
!RT3M1=  0.73205 08075 68877 29353,
!TMRT3=  0.26794 91924 31122 70647

%CONSTANTREAL %C
                  P0 = -0.47083 25141 @+0, 
                  P1 = -0.50909 58253 @-1, 
                  Q0 =  0.14125 00740 @+1 

%CONSTANTREALARRAY A(0:3)= ZERO,PIBY6,PIBY2,PIBY3

%REAL 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= F + F*(P1*G + P0)*G/(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  ATAN}
!******************************************************************************
%EXTERNALREALFN  ATAN2 %ALIAS "F_ATAN2"(%REALNAME ARGY,ARGX)
%REAL 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 SERROR( ARCTAN2ZERO )%ELSE ERROR( ARCTAN2ZERO )      
%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= ATAN(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  ATAN2}
!******************************************************************************
%EXTERNALREALFN  SINH  %ALIAS "F_SINH"(%REALNAME ARG)

%CONSTANTREAL WBAR= 8.75,  { 0.35*(24 + 1) }
LNV   = R'40B17300',
VOV2M1= R'3CE80897',
VSQINV= R'403FFF8C'

!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

%CONSTANTREAL %C
     P0 = -0.71379 3159 @+1, 
     P1 = -0.19033 3399 @+0, 
     Q0 = -0.42827 7109 @+2 

%REAL Y,W,Z,F,RF

Y= MOD(ARG)
%IF Y>ONE %START 
  W= Y - LNV
  %IF W>MAXEXPARG %START
    %IF TEST=1 %THEN SERROR(  SINHLARGE ) %ELSE ERROR (  SINHLARGE )    
    %IF ARG>ZERO %THEN %RESULT= GREATEST %ELSE %RESULT= -GREATEST   
  %FINISHELSESTART
    Z=  EXP(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*(P1*F + P0)/(F + Q0)  
  %RESULT= ARG + (ARG*RF)
%FINISH
%END {OF  SINH}

!******************************************************************************
%EXTERNALREALFN  COSH %ALIAS "F_COSH"(%REALNAME ARG)

%CONSTANTREAL WBAR= 8.75,  { 0.35*(24 + 1) }
LNV   = R'40B17300',
VOV2M1= R'3CE80897',
VSQINV= R'403FFF8C'

!LNV=  0.69316 10107 42187 50000, {0.542714 (OCTAL)}
!VOV2M1=  0.13830 27787 96019 02638 @-4,
!VSQINV=  0.24999 30850 04514 99336 @ 0
 
%REAL Y,W,Z
  
Y= MOD(ARG)
%IF Y>ONE %START
  W= Y - LNV
  %IF W>MAXEXPARG %START
    %IF TEST=1 %THEN SERROR(  COSHLARGE ) %ELSE ERROR (  COSHLARGE )    
    %RESULT= GREATEST
  %FINISH
  Z=  EXP(W)
  %IF W<=WBAR %THEN Z= Z + VSQINV/Z
  %RESULT= Z + VOV2M1*Z
%FINISHELSESTART
  Z=  EXP(Y)
  %RESULT= HALF*Z + HALF/Z
%FINISH
%END {OF  COSH}
!******************************************************************************
%EXTERNALREALFN  TANH %ALIAS "F_TANH"(%REALNAME ARG)


%CONSTREAL    XBIG = 10.05, {= (LN(2) + (6+1)*LN(16))/2}
LN3BY2 = R'408C9F53'

!                LN3BY2 =  0.54930 61443 34054 84570    {LN(3)/2 }   

%CONSTANTREAL %C
   P0    = -0.82377 28127 @+0, 
   P1    = -0.38310 10665 @-2, 
   Q0    =  0.24713 19654 @+1 

%REAL F,ANS,G,R
 
F= MOD(ARG)
%IF F>XBIG %THEN ANS= ONE %ELSE %IF F>LN3BY2 %START
    F= F + F
    ANS= HALF - (ONE / ( EXP(F) + ONE))
    ANS= ANS + ANS
%FINISH %ELSE %IF F<RTEPS %THEN ANS= F %ELSE %START
    G= F*F
    R= (P1*G + P0)*G/(G + Q0)
    ANS= F + (F*R)
%FINISH
%IF ARG<0.0 %THEN %RESULT= -ANS %ELSE %RESULT= ANS
%END {OF  TANH}
!******************************************************************************
%EXTERNALREALFN  GAMMA %ALIAS "F_GAMMA"(%REALNAME ARG)          
%CONSTANTREAL ARGMAX= 57.0,
GAMMAOFARGMAX= 7.1@74,
XSMALL= 1.0@-7,
XVSMALL= 1.4@-76

%CONSTREALARRAY A(1:8)=%C
R'41100000',
R'406C3BC1',
R'40696ACA',
R'4015008C',
R'401283C4',
R'3F10EB02',
R'3F161BC6',
R'3E64B7BB'

!0.99999998,
!0.42278677,
!0.41178574,
!0.82039600@-01,
!0.72323079@-01,
!0.41303702@-02,
!0.53975815@-02,
!0.15368304@-02 

%REAL X,Y,RX,G,Z

      X= ARG
      Y= MOD(X)
      %IF Y>ARGMAX %START
!" ABSOLUTE VALUE OF  GAMMA ARG > 57."      
      %IF TEST=1 %THEN SERROR( GAMLARGE ) %ELSE ERROR( GAMLARGE )
%IF X>ZERO %THEN %RESULT=  GAMMAOFARGMAX %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
%IF TEST=1 %THEN SERROR(GAMNEGINT) %ELSE ERROR( GAMNEGINT )
%RESULT=  GAMMAOFARGMAX
%FINISHELSESTART
%IF Y<XVSMALL %START
!" **  GAMMA ARG TOO NEAR ZERO OR NEG INTEGER." 
%IF TEST=1 %THEN SERROR( GAMNEGINT ) %ELSE ERROR( GAMNEGINT )
%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:   
      Z= X - TWO
      G= ((((((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  GAMMA}
!*
!******************************************************************************
%EXTERNALREALFN ALGAMA %ALIAS "F_LGAMMA"(%REALNAME ARG)
%REAL  G,T,X,Y
%INTEGER I,M
%CONSTANTREAL %C
XSMALL= 1.0@-8, XBIG= 1.2@+3, LNR2PI= 9.18938533@-1

!     RANGE DEPENDENT CONSTANTS
!     FOR IBM 360/370 AND SIMILAR MACHINES
%CONSTANTREAL %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   %CONSTANTREAL %C
!R4   XVBIG= 1.72@+319, GBIG= 1.26@+322
!     FOR UNIVAC 1100 (S.P.)
!R5   %CONSTANTREAL %C
!R5   XVBIG= 1.28@305, GBIG 8.98@+307
!
X= ARG
%IF X<=XSMALL %START {VERY SMALL RANGE}
  %IF X>ZERO %THEN %RESULT= -ALOG(X) %ELSESTART
    %IF TEST=1 %THEN SERROR( LGAMNEG ) %ELSE ERROR( LGAMNEG )
    %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 08E.09
  Y= (((((((((((  +1.88278283@-6*T  -5.48272091@-6)*T %C
 +1.03144033@-5)*T  -3.13088821@-5)*T  +1.01593694@-4)*T %C
 -2.98340924@-4)*T  +9.15547391@-4)*T  -2.42216251@-3)*T %C
 +9.04037536@-3)*T  -1.34119055@-2)*T  +1.03703361@-1)*T %C
 +1.61692007@-2)*T  +8.86226925@-1
  
  Y= Y*G
  %RESULT= ALOG(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

  Y= (  +3.89980902@-9*T  -6.16502533@-6)*T  +8.33271644@-2

  %RESULT= (X-HALF)*ALOG(X) - X + LNR2PI + Y/X

%FINISHELSESTART
  %IF X<=XVBIG %THEN {ASYMPTOTIC LARGE X RANGE} %C
    %RESULT= (X-HALF)*ALOG(X) - X + LNR2PI %ELSESTART
      %IF TEST=1 %THEN SERROR( LGAMLARGE ) %ELSE ERROR( LGAMLARGE )
      %RESULT= GBIG
    %FINISH
%FINISH
%FINISH
%END {OF ALGAMA}
!******************************************************************************
%EXTERNALREALFN  ERF  %ALIAS "F_ERF"(%REALNAME ARG)
!     ERF(X)
%REAL XV,X,X2,BJP2,BJP1,BJ
!
!     *** IMPLEMENTATION DEPENDENT FEATURES ***
!***** DELETE THIS LINE WHEN CORRECT CODE HAS BEEN ACTIVATED *****
!
%CONSTANTREAL XUP= 4.0@0, SQRTPI= 1.772454@0
%CONSTANTREALARRAY C(1:8)= %C
1.944907@0,4.2019@-2,-1.8687@-2,5.129@-3,-1.068@-3,
1.74@-4,-2.1@-5,2.0@-6
%CONSTANTREALARRAY D(1:8)= %C
1.483110@0,-3.01071@-1,6.8995@-2,-1.3916@-2,2.421@-3,
-3.66@-4,4.9@-5,-6.0@-6
!
!
!     NO FAILURE EXITS
X= ARG
XV= MOD(X)
%IF XV<XUP %START
 %IF XV>TWO %START
    X2= TWO - TWENTY/(XV+THREE)
!
!     SUMMATION
    BJP2= X2*C(8)        + C(7)
    BJP1= X2*BJP2 - C(8) + 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* EXP(-X*X)/SQRTPI
    %IF X>ZERO %THEN %RESULT= ONE - X2 %ELSE %RESULT= X2 - ONE
!
  %FINISHELSESTART
    X2= X*X - TWO
!     SUMMATION
    BJP2= X2*D(8)        + D(7)
    BJP1= X2*BJP2 - D(8) + 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 ERF}
!******************************************************************************
%EXTERNALREALFN  ERFC  %ALIAS "F_ERFC"(%REALNAME ARG)
!
!     COMPLEMENT OF ERROR FUNCTION ERFC(X)
!
!     *** IMPLEMENTATION DEPENDENT FEATURES ***
%REAL X,Y,T
!
!     PRECISION DEPENDENT CONSTANTS
%CONSTANTREAL XHI= 9.5@0, XLO= -4.5@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 08E
Y=    +1.4558972@-1+T*(  -2.7342192@-1+T*(  +2.2600824@-1+ %C
 T*(  -1.6357229@-1+T*(  +1.0260225@-1+T*(  -5.4799165@-2+ %C
 T*(  +2.4151896@-2+T*(  -8.2316935@-3+T*(  +1.7866301@-3+ %C
 T*(  -6.4127909@-6+T*(  -1.3874589@-4+ %C
 T*(  +3.1475326@-5)))))))))))
!
!
Y=  EXP(-X*X)*Y
%IF X<ZERO %THEN %RESULT= TWO - Y %ELSE %RESULT= Y
%END {OF ERFC}

!******************************************************************************
%ENDOFFILE