program main c*********************************************************************72 c cc toms597_test() tests toms597(). c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 16 June 2021 c implicit none call timestamp ( ) write ( *, '(a)' ) '' write ( *, '(a)' ) 'toms597_test():' write ( *, '(a)' ) ' Fortran77 version' write ( *, '(a)' ) ' Test toms597():' call ren_test ( ) call ribesl_test ( ) c c Terminate. c write ( *, '(a)' ) '' write ( *, '(a)' ) 'toms597_test():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) '' call timestamp ( ) stop 0 end subroutine ren_test ( ) c*********************************************************************72 c cc ren_test() tests ren(). c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 16 June 2021 c c Author: c c John Burkardt c implicit none integer i integer k double precision ren double precision u write ( *, '(a)' ) '' write ( *, '(a)' ) 'ren_test():' write ( *, '(a)' ) ' Fortran77 version' write ( *, '(a)' ) ' ren() returns a random value in [0,1].' write ( *, '(a)' ) '' k = 0 do i = 1, 10 u = ren ( k ) write ( *, '(2x,i2,2x,g14.6)' ) i, u end do return end subroutine ribesl_test ( ) c*********************************************************************72 c cc ribesl_test() tests ribesl(). c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 January 2016 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 SIGNIFICAND OF A FLOATING-POINT NUMBER C MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE C INTEGER SUCH THAT FLOAT(IBETA)**MINEXP C IS A POSITIVE FLOATING-POINT NUMBER C EPS - THE SMALLEST POSITIVE FLOATING-POINT C NUMBER SUCH THAT 1.0+EPS .NE. 1.0 C EPSNEG - THE SMALLEST POSITIVE FLOATING-POINT C NUMBER SUCH THAT 1.0-EPSNEG .NE. 1.0 C C REN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD Fortran SUBPROGRAMS REQUIRED FOR SINGLE PRECISION C C ABS, ALOG, AMAX1, FLOAT, IFIX, SQRT C C STANDARD Fortran SUBPROGRAMS REQUIRED FOR DOUBLE PRECISION C C DABS, DLOG, DMAX1, DBLE, FLOAT, IFIX, SNGL, DSQRT C C C LATEST REVISION - MAY 18, 1982 C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C------------------------------------------------------------------ implicit none integer i, ibeta, iexp, ii, iout, irnd, it, ize, j, jt, k1, k2, * k3, machep, maxexp, mb, minexp, n, nb, nbm1, ncalc, negep, ngrd double precision a, ait, albeta, aleps, alpha, a1, b, beta, c, * c1, c2, del, den, den1, dgamma, ren, eps, epsneg, half, one, r6, * r7, sum, ten, two, u, w, x, xl, xmax, xmin, xn, xnb, x1, x99, y, * yy, z, zero, zz dimension u(160) data zero, half, one, two, ten /0.0d0,0.5d0,1.0d0,2.0d0,10.0d0/ data x99, c1, c2 /-999.0d0,0.796875d0,1.0095608028653558798921d-3/ data iout /6/ c------------------------------------------------------------------ c c determine machine parameters and set constants c c------------------------------------------------------------------ call machar(ibeta, it, irnd, ngrd, machep, negep, iexp, minexp, * maxexp, eps, epsneg, xmin, xmax) beta = dble(float(ibeta)) albeta = dlog(beta) aleps = dlog(eps) ait = dble(float(it)) a = zero b = two n = 2000 xn = dble(float(n)) jt = 0 ize = 1 mb = 2 c----------------------------------------------------------------- c random argument accuracy tests c----------------------------------------------------------------- do 60 j=1,4 k1 = 0 k3 = 0 x1 = zero a1 = zero r6 = zero r7 = zero del = (b-a)/xn xl = a if (j.eq.1) nb = ifix(5.53e0-0.21e0*sngl(aleps)) if (j.eq.2) nb = ifix(7.26e0-0.27e0*sngl(aleps)) if (j.eq.3) nb = ifix(12.93e0-0.33e0*sngl(aleps)) xnb = dble(float(nb)) nbm1 = nb - 1 c do 50 i=1,n x = del*ren(jt) + xl if (j.eq.4) go to 20 c------------------------------------------------------------------ c first three tests are based on maclaurin series c------------------------------------------------------------------ alpha = ren(jt) y = x/two yy = y*y call ribesl(x, alpha, mb, ize, u, ncalc) z = u(1) den = xnb den1 = den + alpha sum = one/den/den1 do ii=1,nbm1 den = den - one den1 = den1 - one sum = (yy*sum+one)/den/den1 end do y = y**alpha zz = (sum*yy*y+y)/dgamma(one+alpha) go to 30 c------------------------------------------------------------------ c last test is based on asymptotic form for alpha = 0.5 c------------------------------------------------------------------ 20 call ribesl(x, alpha, mb, ize, u, ncalc) z = u(1) zz = c/dsqrt(x) 30 w = (z-zz)/z if (w.gt.zero) k1 = k1 + 1 if (w.lt.zero) k3 = k3 + 1 w = dabs(w) if (w.le.r6) go to 40 r6 = w x1 = x a1 = alpha 40 r7 = r7 + w*w xl = xl + del 50 continue c------------------------------------------------------------------ c gather and print statistics for test c------------------------------------------------------------------ k2 = n - k3 - k1 r7 = dsqrt(r7/xn) write ( iout, '(a)' ) '' if (j.ne.4) write (iout,99999) if (j.eq.4) write (iout,99998) write ( iout, '(a)' ) '' write (iout,99997) n, a, b write (iout,99996) k1, k2, k3 write (iout,99995) it, ibeta w = x99 if (r6.ne.zero) w = dlog(dabs(r6))/albeta write (iout,99994) r6, ibeta, w, x1, a1 w = dmax1(ait+w,zero) write (iout,99993) ibeta, w w = x99 if (r7.ne.zero) w = dlog(dabs(r7))/albeta write (iout,99992) r7, ibeta, w w = dmax1(ait+w,zero) write (iout,99993) ibeta, w c------------------------------------------------------------------ c initialize for next test c------------------------------------------------------------------ a = b b = b + b if (j.eq.2) b = ten if (j.ne.3) go to 60 c = (c1+c2)*half alpha = half a = -aleps*half + two b = a + ten ize = 2 60 continue c----------------------------------------------------------------- c test of error returns c c first, test with bad parameters c----------------------------------------------------------------- write (iout,99991) x = one alpha = half nb = 5 ize = 1 u(1) = zero call ribesl(x, alpha, nb, ize, u, ncalc) write (iout,99990) x, alpha, nb, ize, u(1), ncalc x = -x u(1) = zero call ribesl(x, alpha, nb, ize, u, ncalc) write (iout,99990) x, alpha, nb, ize, u(1), ncalc x = -x alpha = one + half call ribesl(x, alpha, nb, ize, u, ncalc) write (iout,99990) x, alpha, nb, ize, u(1), ncalc alpha = half nb = -nb call ribesl(x, alpha, nb, ize, u, ncalc) write (iout,99990) x, alpha, nb, ize, u(1), ncalc nb = -nb ize = 5 call ribesl(x, alpha, nb, ize, u, ncalc) write (iout,99990) x, alpha, nb, ize, u(1), ncalc c----------------------------------------------------------------- c last tests are with extreme parameters c----------------------------------------------------------------- x = xmin alpha = zero nb = 150 ize = 1 u(1) = zero call ribesl(x, alpha, nb, ize, u, ncalc) write (iout,99990) x, alpha, nb, ize, u(1), ncalc x = ten + ten + ten u(1) = zero call ribesl(x, alpha, nb, ize, u, ncalc) write (iout,99990) x, alpha, nb, ize, u(1), ncalc x = xmax u(1) = zero call ribesl(x, alpha, nb, ize, u, ncalc) write (iout,99990) x, alpha, nb, ize, u(1), ncalc return c----------------------------------------------------------------- 99999 format ( 'test of i(x,alpha) vs series expansion' ) 99998 format ( 'test of exp(-x)*i(x,0.5) vs 1/sqrt(2*x*pi)' ) 99997 format (i7, 48h random arguments were tested from the interval , * 1h(, f5.1, 1h,, f5.1, 1h)//) 99996 format (22h i(x,alpha) was larger, i6, 7h times,/15x, 7h agreed, * i6, 11h times, and/11x, 11hwas smaller, i6, 7h times.//) 99995 format (10h there are, i4, 5h base, i4, 22h significant digits in, * 24h a floating-point number//) 99994 format (30h the maximum relative error of, e15.4, 3h = , i4, * 3h **, f7.2/4x, 16hoccurred for x =, e13.6, 10h and nu =, e13.6) 99993 format (27h the estimated loss of base, i4, 18h significant digit, * 4hs is, f7.2//) 99992 format (40h the root mean square relative error was, e15.4, 3h = , * i4, 3h **, f7.2) 99991 format (24h1check of error returns ///24h the following summarize, * 33hs calls with indicated parameters//23h ncalc different from n, * 30hb indicates some form of error//26h see documentation for rib, * 16hesl for details //7x, 3harg, 12x, 5halpha, 6x, 2hnb, 3x, 2hiz, * 7x, 3hres, 6x, 5hncalc//) 99990 format (2e15.7, 2i5, e15.7, i5) end subroutine machar(ibeta, it, irnd, ngrd, machep, negep, iexp, * minexp, maxexp, eps, epsneg, xmin, xmax) c*********************************************************************72 c cc MACHAR determines characteristics of the floating point arithmetic. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 January 2016 c C C THIS SUBROUTINE IS INTENDED TO DETERMINE THE CHARACTERISTICS C OF THE FLOATING-POINT ARITHMETIC SYSTEM THAT ARE SPECIFIED C BELOW. THE FIRST THREE ARE DETERMINED ACCORDING TO AN C ALGORITHM DUE TO M. MALCOLM, CACM 15 (1972), PP. 949-951, C INCORPORATING SOME, BUT NOT ALL, OF THE IMPROVEMENTS C SUGGESTED BY M. GENTLEMAN AND S. MAROVICH, CACM 17 (1974), C PP. 276-277. THE VERSION GIVEN HERE IS FOR SINGLE PRECISION. C CARDS CONTAINING CD IN COLUMNS 1 AND 2 CAN BE USED TO C CONVERT THE SUBROUTINE TO DOUBLE PRECISION BY REPLACING C EXISTING CARDS IN THE OBVIOUS MANNER. C C C IBETA - THE RADIX OF 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 C NGRD - THE NUMBER OF GUARD DIGITS FOR MULTIPLICATION. IT IS C 0 IF IRND=1, OR IF IRND=0 AND ONLY IT BASE IBETA C DIGITS PARTICIPATE IN THE POST NORMALIZATION SHIFT C OF THE FLOATING-POINT SIGNIFICAND IN MULTIPLICATION C 1 IF IRND=0 AND MORE THAN IT BASE IBETA DIGITS C PARTICIPATE IN THE POST NORMALIZATION SHIFT OF THE C FLOATING-POINT 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 A POSITIVE FLOATING-POINT C NUMBER C MAXEXP - THE LARGEST POSITIVE INTEGER EXPONENT FOR A FINITE C FLOATING-POINT NUMBER C EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH C THAT 1.0+EPS .NE. 1.0. IN PARTICULAR, IF EITHER C IBETA = 2 OR IRND = 0, EPS = FLOAT(IBETA)**MACHEP. C OTHERWISE, EPS = (FLOAT(IBETA)**MACHEP)/2 C EPSNEG - A SMALL POSITIVE FLOATING-POINT NUMBER SUCH THAT C 1.0-EPSNEG .NE. 1.0. IN PARTICULAR, IF IBETA = 2 C OR IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS. C OTHERWISE, EPSNEG = (IBETA**NEGEPS)/2. BECAUSE C NEGEPS IS BOUNDED BELOW BY -(IT+3), EPSNEG MAY NOT C BE THE SMALLEST NUMBER WHICH CAN ALTER 1.0 BY C SUBTRACTION. C XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT POWER OF THE C RADIX. IN PARTICULAR, 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 REVISION - OCTOBER 22, 1979 C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C----------------------------------------------------------------- implicit none integer i, ibeta, iexp, irnd, it, iz, j, k, machep, maxexp, * minexp, mx, negep, ngrd double precision a, b, beta, betain, betam1, eps, epsneg, one, * xmax, xmin, y, z, zero one = dble(float(1)) zero = 0.0d0 c----------------------------------------------------------------- c determine ibeta,beta ala malcolm c----------------------------------------------------------------- a = one 10 a = a + a if (((a+one)-a)-one.eq.zero) go to 10 b = one 20 b = b + b if ((a+b)-a.eq.zero) go to 20 ibeta = int(sngl((a+b)-a)) beta = dble(float(ibeta)) c----------------------------------------------------------------- c determine it, irnd c----------------------------------------------------------------- it = 0 b = one 30 it = it + 1 b = b*beta if (((b+one)-b)-one.eq.zero) go to 30 irnd = 0 betam1 = beta - one if ((a+betam1)-a.ne.zero) irnd = 1 c----------------------------------------------------------------- c determine negep, epsneg c----------------------------------------------------------------- negep = it + 3 betain = one/beta a = one c do 40 i=1,negep a = a*betain 40 continue c b = a 50 if ((one-a)-one.ne.zero) go to 60 a = a*beta negep = negep - 1 go to 50 60 negep = -negep epsneg = a if ((ibeta.eq.2) .or. (irnd.eq.0)) go to 70 a = (a*(one+a))/(one+one) if ((one-a)-one.ne.zero) epsneg = a c----------------------------------------------------------------- c determine machep, eps c----------------------------------------------------------------- 70 machep = -it - 3 a = b 80 if ((one+a)-one.ne.zero) go to 90 a = a*beta machep = machep + 1 go to 80 90 eps = a if ((ibeta.eq.2) .or. (irnd.eq.0)) go to 100 a = (a*(one+a))/(one+one) if ((one+a)-one.ne.zero) eps = a c----------------------------------------------------------------- c determine ngrd c----------------------------------------------------------------- 100 ngrd = 0 if ((irnd.eq.0) .and. ((one+eps)*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 110 y = z z = y*y c----------------------------------------------------------------- c check for underflow here c----------------------------------------------------------------- a = z*one if ((a+a.eq.zero) .or. (dabs(z).ge.y)) go to 120 i = i + 1 k = k + k go to 110 120 if (ibeta.eq.10) go to 130 iexp = i + 1 mx = k + k go to 160 c----------------------------------------------------------------- c for decimal machines only c----------------------------------------------------------------- 130 iexp = 2 iz = ibeta 140 if (k.lt.iz) go to 150 iz = iz*ibeta iexp = iexp + 1 go to 140 150 mx = iz + iz - 1 c----------------------------------------------------------------- c loop to determine minexp, xmin c exit from loop is signaled by an underflow. c----------------------------------------------------------------- 160 xmin = y y = y*betain c----------------------------------------------------------------- c check for underflow here c----------------------------------------------------------------- a = y*one if (((a+a).eq.zero) .or. (dabs(y).ge.xmin)) go to 170 k = k + 1 go to 160 170 minexp = -k c----------------------------------------------------------------- c determine maxexp, xmax c----------------------------------------------------------------- if ((mx.gt.k+k-3) .or. (ibeta.eq.10)) go to 180 mx = mx + mx iexp = iexp + 1 180 maxexp = mx + minexp c----------------------------------------------------------------- c adjust for machines with implicit leading c bit in binary significand and machines with c radix point at extreme 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 200 c do 190 j=1,i if (ibeta.eq.2) xmax = xmax + xmax if (ibeta.ne.2) xmax = xmax*beta 190 continue c 200 return end function ren ( k ) c*********************************************************************72 c cc ren() is a random number generator. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 January 2016 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 implicit none integer IY integer J integer K double precision ren save iy data iy /100001/ j = k iy = iy*125 iy = iy - (iy/2796203)*2796203 ren = dble(float(iy))/2796203.0d0*(1.0d0+1.0d-6+1.0d-12) return end subroutine timestamp ( ) c*********************************************************************72 c cc timestamp() prints out the current YMDHMS date as a timestamp. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 June 2014 c c Author: c c John Burkardt c implicit none character * ( 8 ) ampm integer d character * ( 8 ) date integer h integer m integer mm character * ( 9 ) month(12) integer n integer s character * ( 10 ) time integer y save month data month / & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' / call date_and_time ( date, time ) read ( date, '(i4,i2,i2)' ) y, m, d read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm if ( h .lt. 12 ) then ampm = 'AM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h .lt. 12 ) then ampm = 'PM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, & '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, & trim ( ampm ) return end