C----------------------------------------------------------------------
C FORTRAN 77 program to test RJBESL
C
C Method:
C
C Two different accuracy tests are used. In the first interval,
C function values are compared against values generated with the
C multiplication formula, where the Bessel values used in the
C multiplication formula are obtained from the function program.
C In the remaining intervals, function values are compared
C against values generated with a local Taylor series expansion.
C Derivatives in the expansion are expressed in terms of the
C first two Bessel functions, which are in turn obtained from
C the function program.
C
C Data required
C
C None
C
C Subprograms required from this package
C
C MACHAR - An environmental inquiry program providing
C information on the floating-point arithmetic
C system. Note that the call to MACHAR can
C be deleted provided the following five
C parameters are assigned the values indicated
C
C IBETA - the radix of the floating-point system
C IT - the number of base-IBETA digits in the
C significant of a floating-point number
C EPS - the smallest positive floating-point
C number such that 1.0+EPS .NE. 1.0
C XMAX - the largest finite floating-point number
C
C REN(K) - a function subprogram returning random real
C numbers uniformly distributed over (0,1)
C
C
C Intrinsic functions required are:
C
C ABS, DBLE, INT, LOG, MAX, REAL, SQRT
C
C Reference: "Performance evaluation of programs for certain
C Bessel functions", W. J. Cody and L. Stoltz,
C ACM Trans. on Math. Software, Vol. 15, 1989,
C pp 41-48.
C
C "Use of Taylor series to test accuracy of function
C programs," W. J. Cody and L. Stoltz, submitted for
C publication.
C
C Latest modification: September 6, 1989
C
C Author: W. J. Cody
C Mathematics and Computer Science Division
C Argonne National Laboratory
C Argonne, IL 60439
C
C Acknowledgement: this program is a minor modification of the test
C driver for RIBESL whose primary author was Laura Stoltz.
C
C----------------------------------------------------------------------
INTEGER I,IBETA,IEXP,II,III,IOUT,IRND,IT,J,JT,J1,J2,K,KK,K1,
1 K2,K3,LAST,M,MACHEP,MAXEXP,MB,MINEXP,MVR,N,NCALC,NDUM,
2 NDX,NDX2,NEGEP,NGRD,NK,NO1,NUM
CS REAL
CD DOUBLE PRECISION
1 A,AIT,ALBETA,ALPHA,ALPHSQ,A1,AR1,AR2,B,BETA,C,CONV,D,DEL,
2 DELTA,DERIV,EPS,EPSNEG,F,FXMX,G,HALF,HUND,ONE,REN,RTPI,R6,
3 R7,SIX,SIXTEN,SUM,TEN,TWO,T1,T2,U,U2,W,X,XL,XLAM,XLARGE,
4 XMAX,XMIN,XN,X1,X99,Y,YSQ,Z,ZERO,ZZ
DIMENSION AR1(11,6),AR2(13,9),G(5),NDX(24),NDX2(8),U(560),U2(560)
CS DATA ZERO,HALF,ONE,TWO,SIX/0.0E0,0.5E0,1.0E0,2.0E0,6.0E0/,
CS 1 TEN,SIXTEN,HUND,X99/10.0E0,1.6E1,1.0E2,-999.0E0/,
CS 2 XLAM,XLARGE/1.03125E0,1.0E4/,
CS 3 C,RTPI/0.9189385332E0,0.6366E0/
CD DATA ZERO,HALF,ONE,TWO,SIX/0.0D0,0.5D0,1.0D0,2.0D0,6.0D0/,
CD 1 TEN,SIXTEN,HUND,X99/10.0D0,1.6D1,1.0D2,-999.0D0/,
CD 2 XLAM,XLARGE/1.03125D0,1.0D4/,
CD 3 C,RTPI/0.9189385332D0,0.6366D0/
C----------------------------------------------------------------------
C Arrays related to expansion of the derivatives in terms
C of the first two Bessel functions.
C----------------------------------------------------------------------
DATA NDX/9,7,5,3,1,8,6,4,2,7,5,3,1,6,4,2,5,3,1,4,2,3,1,2/
DATA NDX2/5,9,13,16,19,21,23,24/
CS DATA AR1/0.0E0,-1.0E0,0.0E0,1.0E0,0.0E0,1.0E0,-3.0E0,0.0E0,-2.0E0,
CS 1 1.2E1,0.0E0,1.0E0,0.0E0,-1.0E0,-1.0E0,2.0E0,0.0E0,
CS 2 2.0E0,-6.0E0,1.0E0,-7.0E0,2.4E1,0.0E0,0.0E0,1.0E0,0.0E0,
CS 3 -3.0E0,0.0E0,-2.0E0,1.1E1,0.0E0,1.2E1,-5.0E1,0.0E0,
CS 4 0.0E0,0.0E0,0.0E0,1.0E0,0.0E0,0.0E0,-6.0E0,0.0E0,-2.0E0,
CS 5 3.5E1,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,1.0E0,
CS 6 0.0E0,0.0E0,-1.0E1,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,
CS 7 0.0E0,0.0E0,0.0E0,0.0E0,1.0E0/
CS DATA AR2/-1.0E0,9.0E0,-6.0E1,0.0E0,3.0E0,-5.1E1,3.6E2,0.0E0,
CS 1 1.0E0,-1.8E1,3.45E2,-2.52E3,0.0E0,0.0E0,-3.0E0,3.3E1,
CS 2 -1.2E2,-1.0E0,1.5E1,-1.92E2,7.2E2,0.0E0,4.0E0,-9.6E1,
CS 3 1.32E3,-5.04E3,0.0E0,3.0E0,-7.8E1,2.74E2,0.0E0,-2.7E1,
CS 4 5.7E2,-1.764E3,0.0E0,-4.0E0,2.46E2,-4.666E3,1.3068E4,
CS 5 0.0E0,0.0E0,1.8E1,-2.25E2,0.0E0,3.0E0,-1.5E2,1.624E3,
CS 6 0.0E0,0.0E0,-3.6E1,1.32E3,-1.3132E4,0.0E0,0.0E0,-3.0E0,
CS 7 8.5E1,0.0E0,0.0E0,4.5E1,-7.35E2,0.0E0,0.0E0,6.0E0,
CS 8 -5.5E2,6.769E3,0.0E0,0.0E0,0.0E0,-1.5E1,0.0E0,0.0E0,
CS 9 -3.0E0,1.75E2,0.0E0,0.0E0,0.0E0,6.0E1,-1.96E3,0.0E0,
CS a 0.0E0,0.0E0,1.0E0,0.0E0,0.0E0,0.0E0,-2.1E1,0.0E0,0.0E0,
CS b 0.0E0,-4.0E0,3.22E2,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,
CS c 0.0E0,1.0E0,0.0E0,0.0E0,0.0E0,0.0E0,-2.8E1,0.0E0,0.0E0,
CS d 0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,
CS e 0.0E0,1.0E0/
CD DATA AR1/0.0D0,-1.0D0,0.0D0,1.0D0,0.0D0,1.0D0,-3.0D0,0.0D0,-2.0D0,
CD 1 1.2D1,0.0D0,1.0D0,0.0D0,-1.0D0,-1.0D0,2.0D0,0.0D0,
CD 2 2.0D0,-6.0D0,1.0D0,-7.0D0,2.4D1,0.0D0,0.0D0,1.0D0,0.0D0,
CD 3 -3.0D0,0.0D0,-2.0D0,1.1D1,0.0D0,1.2D1,-5.0D1,0.0D0,
CD 4 0.0D0,0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,-6.0D0,0.0D0,-2.0D0,
CD 5 3.5D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,1.0D0,
CD 6 0.0D0,0.0D0,-1.0D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,
CD 7 0.0D0,0.0D0,0.0D0,0.0D0,1.0D0/
CD DATA AR2/-1.0D0,9.0D0,-6.0D1,0.0D0,3.0D0,-5.1D1,3.6D2,0.0D0,
CD 1 1.0D0,-1.8D1,3.45D2,-2.52D3,0.0D0,0.0D0,-3.0D0,3.3D1,
CD 2 -1.2D2,-1.0D0,1.5D1,-1.92D2,7.2D2,0.0D0,4.0D0,-9.6D1,
CD 3 1.32D3,-5.04D3,0.0D0,3.0D0,-7.8D1,2.74D2,0.0D0,-2.7D1,
CD 4 5.7D2,-1.764D3,0.0D0,-4.0D0,2.46D2,-4.666D3,1.3068D4,
CD 5 0.0D0,0.0D0,1.8D1,-2.25D2,0.0D0,3.0D0,-1.5D2,1.624D3,
CD 6 0.0D0,0.0D0,-3.6D1,1.32D3,-1.3132D4,0.0D0,0.0D0,-3.0D0,
CD 7 8.5D1,0.0D0,0.0D0,4.5D1,-7.35D2,0.0D0,0.0D0,6.0D0,
CD 8 -5.5D2,6.769D3,0.0D0,0.0D0,0.0D0,-1.5D1,0.0D0,0.0D0,
CD 9 -3.0D0,1.75D2,0.0D0,0.0D0,0.0D0,6.0D1,-1.96D3,0.0D0,
CD a 0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,-2.1D1,0.0D0,0.0D0,
CD b 0.0D0,-4.0D0,3.22D2,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,
CD c 0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,0.0D0,-2.8D1,0.0D0,0.0D0,
CD d 0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,
CD e 0.0D0,1.0D0/
DATA IOUT/6/
C----------------------------------------------------------------------
C Statement function for integer to float conversion
C----------------------------------------------------------------------
CS CONV(NDUM) = REAL(NDUM)
CD CONV(NDUM) = DBLE(NDUM)
C----------------------------------------------------------------------
C Determine machine parameters and set constants
C----------------------------------------------------------------------
CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
1 MAXEXP,EPS,EPSNEG,XMIN,XMAX)
BETA = CONV(IBETA)
AIT = CONV(IT)
ALBETA = LOG(BETA)
A = ZERO
B = TWO
JT = 0
DELTA = XLAM - ONE
F = (DELTA) * (XLAM+ONE) * HALF
C----------------------------------------------------------------------
C Random argument accuracy tests
C----------------------------------------------------------------------
DO 300 J = 1, 4
N = 2000
XN = CONV(N)
K1 = 0
K2 = 0
K3 = 0
X1 = ZERO
A1 = ZERO
R6 = ZERO
R7 = ZERO
DEL = (B - A) / XN
XL = A
DO 200 I = 1, N
X = DEL * REN(JT) + XL
110 ALPHA = REN(JT)
C----------------------------------------------------------------------
C Carefully purify arguments
C----------------------------------------------------------------------
IF (J .EQ. 1) THEN
Y = X/XLAM
ELSE
Y = X - DELTA
END IF
W = SIXTEN * Y
T1 = W + Y
T1 = W + T1
Y = T1 - W
Y = Y - W
IF (J .EQ. 1) THEN
X = Y * XLAM
MB = 11
ELSE
X = Y + DELTA
MB = 2
END IF
CALL RJBESL(Y,ALPHA,MB,U2,NCALC)
CALL RJBESL(X,ALPHA,MB,U,NCALC)
IF (J .EQ. 1) THEN
C----------------------------------------------------------------------
C Accuracy test is based on the multiplication theorem
C----------------------------------------------------------------------
D = -F*Y
MB = NCALC - 2
XMB = CONV(MB)
SUM = U2(MB+1)
IND = MB
DO 125 II = 2, MB
SUM = SUM * D / XMB + U2(IND)
IND = IND - 1
XMB = XMB - ONE
125 CONTINUE
ZZ = SUM * D + U2(IND)
ZZ = ZZ * XLAM ** ALPHA
ELSE
C----------------------------------------------------------------------
C Accuracy test is based on local Taylor's series expansion
C----------------------------------------------------------------------
IF (ABS(U(1)) .LT. ABS(U2(1))) THEN
Z = X
X = Y
Y = Z
DELTA = X - Y
DO 130 II = 1, MB
Z = U(II)
U(II) = U2(II)
U2(II) = Z
130 CONTINUE
END IF
C----------------------------------------------------------------------
C Filter out cases where function values or derivatives are small
C----------------------------------------------------------------------
W = SQRT(RTPI/X)/SIXTEN
Z = MIN(ABS(U2(1)),ABS(U2(2)))
IF (Z .LT. W) GO TO 110
YSQ = Y * Y
ALPHSQ = ALPHA * ALPHA
MB = 8
J1 = MB
XJ1 = CONV(J1+1)
IEXP = 0
NK = 13
NUM = 2
DO 180 II = 1, MB
IF (NK .EQ. 0) THEN
NK = 11
NUM = 1
END IF
K = 9 - J1
IF (K .GT. 1) THEN
NO1 = NDX2(K-1) + 1
ELSE
NO1 = 1
END IF
MVR = NO1
LAST = NDX2(K)
K = LAST - NO1 + 1
C----------------------------------------------------------------------
C Group J(ALPHA) terms in the derivative
C----------------------------------------------------------------------
DO 160 III = 1, K
J2 = NDX(MVR)
IF (NUM .EQ. 1) THEN
G(III) = AR1(NK,J2)
ELSE
G(III) = AR2(NK,J2)
END IF
IF (J2 .GT. 1) THEN
157 J2 = J2 - 1
IF (NUM .EQ. 1) THEN
G(III) = G(III) * ALPHA + AR1(NK,J2)
ELSE
G(III) = G(III) * ALPHA + AR2(NK,J2)
END IF
IF (J2 .GT. 1) GO TO 157
END IF
MVR = MVR + 1
NK = NK - 1
160 CONTINUE
T1 = G(1)
DO 162 III = 2, K
T1 = T1 / YSQ + G(III)
162 CONTINUE
IF (IEXP .EQ. 1) T1 = T1 / Y
C----------------------------------------------------------------------
C Group J(ALPHA+1) terms in the derivative
C----------------------------------------------------------------------
IEXP = 1 - IEXP
NK = NK + K
MVR = NO1
KK = K
DO 165 III = 1, K
J2 = NDX(MVR)
M = MOD(J2,2)
IF (M .EQ. 1) J2 = J2 - 1
IF (J2 .GE. 2) THEN
IF (NUM .EQ. 1) THEN
G(III) = AR1(NK,J2)
ELSE
G(III) = AR2(NK,J2)
END IF
163 J2 = J2 - 2
IF (J2 .GE. 2) THEN
IF (NUM .EQ. 1) THEN
G(III) = G(III) * ALPHSQ +
1 AR1(NK,J2)
ELSE
G(III) = G(III) * ALPHSQ +
1 AR2(NK,J2)
END IF
GO TO 163
END IF
ELSE
KK = III - 1
END IF
MVR = MVR + 1
NK = NK - 1
165 CONTINUE
T2 = G(1)
DO 167 III = 2, KK
T2 = T2 / YSQ + G(III)
167 CONTINUE
IF (IEXP .EQ. 1) T2 = T2 / Y
DERIV = U2(1) * T1 - U2(2) * T2
IF (J1 .EQ. 8) THEN
SUM = DERIV
ELSE
SUM = SUM * DELTA / XJ1 + DERIV
END IF
J1 = J1 - 1
XJ1 = XJ1 - ONE
180 CONTINUE
ZZ = SUM * DELTA + U2(1)
END IF
MB = 2
Z = U(1)
C----------------------------------------------------------------------
C Accumulate Results
C----------------------------------------------------------------------
W = (Z - ZZ) / Z
IF (W .GT. ZERO) THEN
K1 = K1 + 1
ELSE IF (W .LT. ZERO) THEN
K3 = K3 + 1
END IF
W = ABS(W)
IF (W .GT. R6) THEN
R6 = W
X1 = X
A1 = ALPHA
FXMX = Z
END IF
R7 = R7 + W * W
XL = XL + DEL
200 CONTINUE
C----------------------------------------------------------------------
C Gather and print statistics for test
C----------------------------------------------------------------------
K2 = N - K1 - K3
R7 = SQRT(R7/XN)
IF (J .EQ. 1) THEN
WRITE (IOUT,1000)
ELSE
WRITE (IOUT,1001)
END IF
WRITE (IOUT,1010) N,A,B
WRITE (IOUT,1011) K1,K2,K3
WRITE (IOUT,1020) IT,IBETA
IF (R6 .NE. ZERO) THEN
W = LOG(R6)/ALBETA
ELSE
W = X99
END IF
WRITE (IOUT,1021) R6,IBETA,W,X1,A1
WRITE (IOUT,1024) FXMX
W = MAX(AIT+W,ZERO)
WRITE (IOUT,1022) IBETA,W
IF (R7 .NE. ZERO) THEN
W = LOG(R7)/ALBETA
ELSE
W = X99
END IF
WRITE (IOUT,1023) R7,IBETA,W
W = MAX(AIT+W,ZERO)
WRITE (IOUT,1022) IBETA,W
C----------------------------------------------------------------------
C Initialize for next test
C----------------------------------------------------------------------
A = B
IF (J .EQ. 1) THEN
B = TEN
ELSE IF (J .EQ. 2) THEN
B = B + B
ELSE IF (J .EQ. 3) THEN
A = A + TEN
B = A + TEN
END IF
300 CONTINUE
C----------------------------------------------------------------------
C Test of error returns
C
C First, test with bad parameters
C----------------------------------------------------------------------
WRITE (IOUT, 2006)
X = ONE
ALPHA = ONE + HALF
MB = 5
CALL RJBESL(X,ALPHA,MB,U,NCALC)
WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC
ALPHA = HALF
MB = -MB
CALL RJBESL(X,ALPHA,MB,U,NCALC)
C----------------------------------------------------------------------
C Last tests are with extreme parameters
C----------------------------------------------------------------------
X = ZERO
ALPHA = ONE
MB = 2
CALL RJBESL(X,ALPHA,MB,U,NCALC)
WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC
X = -ONE
ALPHA = HALF
MB = 5
CALL RJBESL(X,ALPHA,MB,U,NCALC)
WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC
C----------------------------------------------------------------------
C Determine largest safe argument
C----------------------------------------------------------------------
WRITE (IOUT, 2015)
X = XLARGE * (ONE - SQRT(SQRT(EPS)))
MB = 2
CALL RJBESL(X,ALPHA,MB,U,NCALC)
WRITE (IOUT, 2012) X
WRITE (IOUT, 2014) NCALC,U(1)
X = XLARGE * (ONE + SQRT(SQRT(EPS)))
MB = 2
CALL RJBESL(X,ALPHA,MB,U,NCALC)
WRITE (IOUT, 2012) X
WRITE (IOUT, 2013)
WRITE (IOUT, 2014) NCALC,U(1)
WRITE (IOUT, 2020)
STOP
C----------------------------------------------------------------------
1000 FORMAT('1Test of J(X,ALPHA) vs Multiplication Theorem'//)
1001 FORMAT('1Test of J(X,ALPHA) vs Taylor series'//)
1010 FORMAT(I7,' Random arguments were tested from the interval ',
1 '(',F5.2,',',F5.2,')'//)
1011 FORMAT(' J(X,ALPHA) was larger',I6,' times,'/
1 15X,' agreed',I6,' times, and'/
1 11X,'was smaller',I6,' times.'//)
1020 FORMAT(' There are',I4,' base',I4,
1 ' significant digits in a floating-point number'//)
1021 FORMAT(' The maximum relative error of',E15.4,' = ',I4,' **',
1 F7.2/4X,'occurred for X =',E13.6,' and NU =',E13.6)
1022 FORMAT(' The estimated loss of base',I4,
1 ' significant digits is',F7.2//)
1023 FORMAT(' The root mean square relative error was',E15.4,
1 ' = ',I4,' **',F7.2)
1024 FORMAT(4x,'with J(X,ALPHA) = ',E13.6)
2006 FORMAT('1Check of Error Returns'///
1 ' The following summarizes calls with indicated parameters'//
2 ' NCALC different from MB indicates some form of error'//
3 ' See documentation for RJBESL for details'//
4 7X,'ARG',12X,'ALPHA',6X,'MB',6X,'B(1)',6X,'NCALC'//)
2011 FORMAT(2E15.7,I5,E15.7,I5//)
2012 FORMAT(' RJBESL will be called with the argument',E13.6)
2013 FORMAT(' This should trigger an error message.')
2014 FORMAT(' NCALC returned the value',I5/
1 ' and RJBESL returned U(1) = ',E13.6/)
2015 FORMAT(' Tests near the largest acceptable argument for RJBESL'/)
2020 FORMAT(' This concludes the tests.')
C ---------- Last line of RJBESL test program ----------
END
SUBROUTINE MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
1 MAXEXP,EPS,EPSNEG,XMIN,XMAX)
C----------------------------------------------------------------------
C This Fortran 77 subroutine is intended to determine the parameters
C of the floating-point arithmetic system specified below. The
C determination of the first three uses an extension of an algorithm
C due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some,
C but not all, of the improvements suggested by M. Gentleman and S.
C Marovich, CACM 17 (1974), pp. 276-277. An earlier version of this
C program was published in the book Software Manual for the
C Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall,
C Englewood Cliffs, NJ, 1980.
C
C The program as given here must be modified before compiling. If
C a single (double) precision version is desired, change all
C occurrences of CS (CD) in columns 1 and 2 to blanks.
C
C Parameter values reported are as follows:
C
C IBETA - the radix for the floating-point representation
C IT - the number of base IBETA digits in the floating-point
C significand
C IRND - 0 if floating-point addition chops
C 1 if floating-point addition rounds, but not in the
C IEEE style
C 2 if floating-point addition rounds in the IEEE style
C 3 if floating-point addition chops, and there is
C partial underflow
C 4 if floating-point addition rounds, but not in the
C IEEE style, and there is partial underflow
C 5 if floating-point addition rounds in the IEEE style,
C and there is partial underflow
C NGRD - the number of guard digits for multiplication with
C truncating arithmetic. It is
C 0 if floating-point arithmetic rounds, or if it
C truncates and only IT base IBETA digits
C participate in the post-normalization shift of the
C floating-point significand in multiplication;
C 1 if floating-point arithmetic truncates and more
C than IT base IBETA digits participate in the
C post-normalization shift of the floating-point
C significand in multiplication.
C MACHEP - the largest negative integer such that
C 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that
C MACHEP is bounded below by -(IT+3)
C NEGEPS - the largest negative integer such that
C 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that
C NEGEPS is bounded below by -(IT+3)
C IEXP - the number of bits (decimal places if IBETA = 10)
C reserved for the representation of the exponent
C (including the bias or sign) of a floating-point
C number
C MINEXP - the largest in magnitude negative integer such that
C FLOAT(IBETA)**MINEXP is positive and normalized
C MAXEXP - the smallest positive power of BETA that overflows
C EPS - FLOAT(IBETA)**MACHEP.
C EPSNEG - FLOAT(IBETA)**NEGEPS.
C XMIN - the smallest non-vanishing normalized floating-point
C power of the radix, i.e., XMIN = FLOAT(IBETA)**MINEXP
C XMAX - the largest finite floating-point number. In
C particular XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
C Note - on some machines XMAX will be only the
C second, or perhaps third, largest number, being
C too small by 1 or 2 units in the last digit of
C the significand.
C
C Latest modification: May 30, 1989
C
C Author: W. J. Cody
C Mathematics and Computer Science Division
C Argonne National Laboratory
C Argonne, IL 60439
C
C----------------------------------------------------------------------
INTEGER I,IBETA,IEXP,IRND,IT,ITEMP,IZ,J,K,MACHEP,MAXEXP,
1 MINEXP,MX,NEGEP,NGRD,NXRES
CS REAL
CD DOUBLE PRECISION
1 A,B,BETA,BETAIN,BETAH,CONV,EPS,EPSNEG,ONE,T,TEMP,TEMPA,
2 TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO
C----------------------------------------------------------------------
CS CONV(I) = REAL(I)
CD CONV(I) = DBLE(I)
ONE = CONV(1)
TWO = ONE + ONE
ZERO = ONE - ONE
C----------------------------------------------------------------------
C Determine IBETA, BETA ala Malcolm.
C----------------------------------------------------------------------
A = ONE
10 A = A + A
TEMP = A+ONE
TEMP1 = TEMP-A
IF (TEMP1-ONE .EQ. ZERO) GO TO 10
B = ONE
20 B = B + B
TEMP = A+B
ITEMP = INT(TEMP-A)
IF (ITEMP .EQ. 0) GO TO 20
IBETA = ITEMP
BETA = CONV(IBETA)
C----------------------------------------------------------------------
C Determine IT, IRND.
C----------------------------------------------------------------------
IT = 0
B = ONE
100 IT = IT + 1
B = B * BETA
TEMP = B+ONE
TEMP1 = TEMP-B
IF (TEMP1-ONE .EQ. ZERO) GO TO 100
IRND = 0
BETAH = BETA / TWO
TEMP = A+BETAH
IF (TEMP-A .NE. ZERO) IRND = 1
TEMPA = A + BETA
TEMP = TEMPA+BETAH
IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2
C----------------------------------------------------------------------
C Determine NEGEP, EPSNEG.
C----------------------------------------------------------------------
NEGEP = IT + 3
BETAIN = ONE / BETA
A = ONE
DO 200 I = 1, NEGEP
A = A * BETAIN
200 CONTINUE
B = A
210 TEMP = ONE-A
IF (TEMP-ONE .NE. ZERO) GO TO 220
A = A * BETA
NEGEP = NEGEP - 1
GO TO 210
220 NEGEP = -NEGEP
EPSNEG = A
C----------------------------------------------------------------------
C Determine MACHEP, EPS.
C----------------------------------------------------------------------
MACHEP = -IT - 3
A = B
300 TEMP = ONE+A
IF (TEMP-ONE .NE. ZERO) GO TO 320
A = A * BETA
MACHEP = MACHEP + 1
GO TO 300
320 EPS = A
C----------------------------------------------------------------------
C Determine NGRD.
C----------------------------------------------------------------------
NGRD = 0
TEMP = ONE+EPS
IF ((IRND .EQ. 0) .AND. (TEMP*ONE-ONE .NE. ZERO)) NGRD = 1
C----------------------------------------------------------------------
C Determine IEXP, MINEXP, XMIN.
C
C Loop to determine largest I and K = 2**I such that
C (1/BETA) ** (2**(I))
C does not underflow.
C Exit from loop is signaled by an underflow.
C----------------------------------------------------------------------
I = 0
K = 1
Z = BETAIN
T = ONE + EPS
NXRES = 0
400 Y = Z
Z = Y * Y
C----------------------------------------------------------------------
C Check for underflow here.
C----------------------------------------------------------------------
A = Z * ONE
TEMP = Z * T
IF ((A+A .EQ. ZERO) .OR. (ABS(Z) .GE. Y)) GO TO 410
TEMP1 = TEMP * BETAIN
IF (TEMP1*BETA .EQ. Z) GO TO 410
I = I + 1
K = K + K
GO TO 400
410 IF (IBETA .EQ. 10) GO TO 420
IEXP = I + 1
MX = K + K
GO TO 450
C----------------------------------------------------------------------
C This segment is for decimal machines only.
C----------------------------------------------------------------------
420 IEXP = 2
IZ = IBETA
430 IF (K .LT. IZ) GO TO 440
IZ = IZ * IBETA
IEXP = IEXP + 1
GO TO 430
440 MX = IZ + IZ - 1
C----------------------------------------------------------------------
C Loop to determine MINEXP, XMIN.
C Exit from loop is signaled by an underflow.
C----------------------------------------------------------------------
450 XMIN = Y
Y = Y * BETAIN
C----------------------------------------------------------------------
C Check for underflow here.
C----------------------------------------------------------------------
A = Y * ONE
TEMP = Y * T
IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460
K = K + 1
TEMP1 = TEMP * BETAIN
IF ((TEMP1*BETA .NE. Y) .OR. (TEMP .EQ. Y)) THEN
GO TO 450
ELSE
NXRES = 3
XMIN = Y
END IF
460 MINEXP = -K
C----------------------------------------------------------------------
C Determine MAXEXP, XMAX.
C----------------------------------------------------------------------
IF ((MX .GT. K+K-3) .OR. (IBETA .EQ. 10)) GO TO 500
MX = MX + MX
IEXP = IEXP + 1
500 MAXEXP = MX + MINEXP
C----------------------------------------------------------------------
C Adjust IRND to reflect partial underflow.
C----------------------------------------------------------------------
IRND = IRND + NXRES
C----------------------------------------------------------------------
C Adjust for IEEE-style machines.
C----------------------------------------------------------------------
IF (IRND .GE. 2) MAXEXP = MAXEXP - 2
C----------------------------------------------------------------------
C Adjust for machines with implicit leading bit in binary
C significand, and machines with radix point at extreme
C right of significand.
C----------------------------------------------------------------------
I = MAXEXP + MINEXP
IF ((IBETA .EQ. 2) .AND. (I .EQ. 0)) MAXEXP = MAXEXP - 1
IF (I .GT. 20) MAXEXP = MAXEXP - 1
IF (A .NE. Y) MAXEXP = MAXEXP - 2
XMAX = ONE - EPSNEG
IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG
XMAX = XMAX / (BETA * BETA * BETA * XMIN)
I = MAXEXP + MINEXP + 3
IF (I .LE. 0) GO TO 520
DO 510 J = 1, I
IF (IBETA .EQ. 2) XMAX = XMAX + XMAX
IF (IBETA .NE. 2) XMAX = XMAX * BETA
510 CONTINUE
520 RETURN
C---------- Last line of MACHAR ----------
END
FUNCTION REN(K)
C---------------------------------------------------------------------
C Random number generator - based on Algorithm 266 by Pike and
C Hill (modified by Hansson), Communications of the ACM,
C Vol. 8, No. 10, October 1965.
C
C This subprogram is intended for use on computers with
C fixed point wordlength of at least 29 bits. It is
C best if the floating-point significand has at most
C 29 bits.
C
C Latest modification: May 30, 1989
C
C Author: W. J. Cody
C Mathematics and Computer Science Division
C Argonne National Laboratory
C Argonne, IL 60439
C
C---------------------------------------------------------------------
INTEGER IY,J,K
CS REAL CONV,C1,C2,C3,ONE,REN
CD DOUBLE PRECISION CONV,C1,C2,C3,ONE,REN
DATA IY/100001/
CS DATA ONE,C1,C2,C3/1.0E0,2796203.0E0,1.0E-6,1.0E-12/
CD DATA ONE,C1,C2,C3/1.0D0,2796203.0D0,1.0D-6,1.0D-12/
C---------------------------------------------------------------------
C Statement functions for conversion between integer and float
C---------------------------------------------------------------------
CS CONV(J) = REAL(J)
CD CONV(J) = DBLE(J)
C---------------------------------------------------------------------
J = K
IY = IY * 125
IY = IY - (IY/2796203) * 2796203
REN = CONV(IY) / C1 * (ONE + C2 + C3)
RETURN
C---------- Last card of REN ----------
END