program main !*********************************************************************72 ! !! toms597_test() tests toms597(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 June 2021 ! implicit none call timestamp ( ) write ( *, '(a)' ) '' write ( *, '(a)' ) 'toms597_test():' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' Test toms597():' call ren_test ( ) call ribesl_test ( ) ! ! Terminate. ! write ( *, '(a)' ) '' write ( *, '(a)' ) 'toms597_test():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) '' call timestamp ( ) stop 0 end subroutine ren_test ( ) !*****************************************************************************80 ! !! ren_test() tests ren(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 June 2021 ! ! Author: ! ! John Burkardt ! implicit none integer i integer k double precision ren double precision u write ( *, '(a)' ) '' write ( *, '(a)' ) 'ren_test():' write ( *, '(a)' ) ' FORTRAN90 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 ( ) !*****************************************************************************80 ! !! ribesl_test() tests ribesl(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 09 January 2016 ! ! PROGRAM TO TEST RIBESL ! ! DATA REQUIRED ! ! NONE ! ! SUBPROGRAMS REQUIRED FROM THIS PACKAGE ! ! MACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING ! INFORMATION ON THE FLOATING-POINT ARITHMETIC ! SYSTEM. NOTE THAT THE CALL TO MACHAR CAN ! BE DELETED PROVIDED THE FOLLOWING FIVE ! PARAMETERS ARE ASSIGNED THE VALUES INDICATED ! ! IBETA - THE RADIX OF THE FLOATING-POINT SYSTEM ! IT - THE NUMBER OF BASE-IBETA DIGITS IN THE ! SIGNIFICAND OF A FLOATING-POINT NUMBER ! MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE ! INTEGER SUCH THAT FLOAT(IBETA)**MINEXP ! IS A POSITIVE FLOATING-POINT NUMBER ! EPS - THE SMALLEST POSITIVE FLOATING-POINT ! NUMBER SUCH THAT 1.0+EPS .NE. 1.0 ! EPSNEG - THE SMALLEST POSITIVE FLOATING-POINT ! NUMBER SUCH THAT 1.0-EPSNEG .NE. 1.0 ! ! REN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL ! NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) ! ! ! STANDARD FORTRAN SUBPROGRAMS REQUIRED FOR SINGLE PRECISION ! ! ABS, ALOG, AMAX1, FLOAT, IFIX, SQRT ! ! STANDARD FORTRAN SUBPROGRAMS REQUIRED FOR DOUBLE PRECISION ! ! DABS, DLOG, DMAX1, DBLE, FLOAT, IFIX, SNGL, DSQRT ! ! ! LATEST REVISION - MAY 18, 1982 ! ! AUTHOR - W. J. CODY ! ARGONNE NATIONAL LABORATORY ! !------------------------------------------------------------------ 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/ !------------------------------------------------------------------ ! ! determine machine parameters and set constants ! !------------------------------------------------------------------ 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 !----------------------------------------------------------------- ! random argument accuracy tests !----------------------------------------------------------------- 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 ! do 50 i=1,n x = del*ren(jt) + xl if (j.eq.4) go to 20 !------------------------------------------------------------------ ! first three tests are based on maclaurin series !------------------------------------------------------------------ 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 !------------------------------------------------------------------ ! last test is based on asymptotic form for alpha = 0.5 !------------------------------------------------------------------ 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 !------------------------------------------------------------------ ! gather and print statistics for test !------------------------------------------------------------------ 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 !------------------------------------------------------------------ ! initialize for next test !------------------------------------------------------------------ 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 !----------------------------------------------------------------- ! test of error returns ! ! first, test with bad parameters !----------------------------------------------------------------- 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 !----------------------------------------------------------------- ! last tests are with extreme parameters !----------------------------------------------------------------- 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 !----------------------------------------------------------------- 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) !*********************************************************************72 ! !! MACHAR determines characteristics of the floating point arithmetic. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 09 January 2016 ! ! ! THIS SUBROUTINE IS INTENDED TO DETERMINE THE CHARACTERISTICS ! OF THE FLOATING-POINT ARITHMETIC SYSTEM THAT ARE SPECIFIED ! BELOW. THE FIRST THREE ARE DETERMINED ACCORDING TO AN ! ALGORITHM DUE TO M. MALCOLM, CACM 15 (1972), PP. 949-951, ! INCORPORATING SOME, BUT NOT ALL, OF THE IMPROVEMENTS ! SUGGESTED BY M. GENTLEMAN AND S. MAROVICH, CACM 17 (1974), ! PP. 276-277. THE VERSION GIVEN HERE IS FOR SINGLE PRECISION. ! CARDS CONTAINING CD IN COLUMNS 1 AND 2 CAN BE USED TO ! CONVERT THE SUBROUTINE TO DOUBLE PRECISION BY REPLACING ! EXISTING CARDS IN THE OBVIOUS MANNER. ! ! ! IBETA - THE RADIX OF THE FLOATING-POINT REPRESENTATION ! IT - THE NUMBER OF BASE IBETA DIGITS IN THE FLOATING-POINT ! SIGNIFICAND ! IRND - 0 IF FLOATING-POINT ADDITION CHOPS, ! 1 IF FLOATING-POINT ADDITION ROUNDS ! NGRD - THE NUMBER OF GUARD DIGITS FOR MULTIPLICATION. IT IS ! 0 IF IRND=1, OR IF IRND=0 AND ONLY IT BASE IBETA ! DIGITS PARTICIPATE IN THE POST NORMALIZATION SHIFT ! OF THE FLOATING-POINT SIGNIFICAND IN MULTIPLICATION ! 1 IF IRND=0 AND MORE THAN IT BASE IBETA DIGITS ! PARTICIPATE IN THE POST NORMALIZATION SHIFT OF THE ! FLOATING-POINT SIGNIFICAND IN MULTIPLICATION ! MACHEP - THE LARGEST NEGATIVE INTEGER SUCH THAT ! 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, EXCEPT THAT ! MACHEP IS BOUNDED BELOW BY -(IT+3) ! NEGEPS - THE LARGEST NEGATIVE INTEGER SUCH THAT ! 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, EXCEPT THAT ! NEGEPS IS BOUNDED BELOW BY -(IT+3) ! IEXP - THE NUMBER OF BITS (DECIMAL PLACES IF IBETA = 10) ! RESERVED FOR THE REPRESENTATION OF THE EXPONENT ! (INCLUDING THE BIAS OR SIGN) OF A FLOATING-POINT ! NUMBER ! MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE INTEGER SUCH THAT ! FLOAT(IBETA)**MINEXP IS A POSITIVE FLOATING-POINT ! NUMBER ! MAXEXP - THE LARGEST POSITIVE INTEGER EXPONENT FOR A FINITE ! FLOATING-POINT NUMBER ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH ! THAT 1.0+EPS .NE. 1.0. IN PARTICULAR, IF EITHER ! IBETA = 2 OR IRND = 0, EPS = FLOAT(IBETA)**MACHEP. ! OTHERWISE, EPS = (FLOAT(IBETA)**MACHEP)/2 ! EPSNEG - A SMALL POSITIVE FLOATING-POINT NUMBER SUCH THAT ! 1.0-EPSNEG .NE. 1.0. IN PARTICULAR, IF IBETA = 2 ! OR IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS. ! OTHERWISE, EPSNEG = (IBETA**NEGEPS)/2. BECAUSE ! NEGEPS IS BOUNDED BELOW BY -(IT+3), EPSNEG MAY NOT ! BE THE SMALLEST NUMBER WHICH CAN ALTER 1.0 BY ! SUBTRACTION. ! XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT POWER OF THE ! RADIX. IN PARTICULAR, XMIN = FLOAT(IBETA)**MINEXP ! XMAX - THE LARGEST FINITE FLOATING-POINT NUMBER. IN ! PARTICULAR XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP ! NOTE - ON SOME MACHINES XMAX WILL BE ONLY THE ! SECOND, OR PERHAPS THIRD, LARGEST NUMBER, BEING ! TOO SMALL BY 1 OR 2 UNITS IN THE LAST DIGIT OF ! THE SIGNIFICAND. ! ! LATEST REVISION - OCTOBER 22, 1979 ! ! AUTHOR - W. J. CODY ! ARGONNE NATIONAL LABORATORY ! !----------------------------------------------------------------- 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 !----------------------------------------------------------------- ! determine ibeta,beta ala malcolm !----------------------------------------------------------------- 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)) !----------------------------------------------------------------- ! determine it, irnd !----------------------------------------------------------------- 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 !----------------------------------------------------------------- ! determine negep, epsneg !----------------------------------------------------------------- negep = it + 3 betain = one/beta a = one ! do 40 i=1,negep a = a*betain 40 continue ! 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 !----------------------------------------------------------------- ! determine machep, eps !----------------------------------------------------------------- 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 !----------------------------------------------------------------- ! determine ngrd !----------------------------------------------------------------- 100 ngrd = 0 if ((irnd.eq.0) .and. ((one+eps)*one-one).ne.zero) ngrd = 1 !----------------------------------------------------------------- ! determine iexp, minexp, xmin ! ! loop to determine largest i and k = 2**i such that ! (1/beta) ** (2**(i)) ! does not underflow ! exit from loop is signaled by an underflow. !----------------------------------------------------------------- i = 0 k = 1 z = betain 110 y = z z = y*y !----------------------------------------------------------------- ! check for underflow here !----------------------------------------------------------------- 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 !----------------------------------------------------------------- ! for decimal machines only !----------------------------------------------------------------- 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 !----------------------------------------------------------------- ! loop to determine minexp, xmin ! exit from loop is signaled by an underflow. !----------------------------------------------------------------- 160 xmin = y y = y*betain !----------------------------------------------------------------- ! check for underflow here !----------------------------------------------------------------- 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 !----------------------------------------------------------------- ! determine maxexp, xmax !----------------------------------------------------------------- if ((mx.gt.k+k-3) .or. (ibeta.eq.10)) go to 180 mx = mx + mx iexp = iexp + 1 180 maxexp = mx + minexp !----------------------------------------------------------------- ! adjust for machines with implicit leading ! bit in binary significand and machines with ! radix point at extreme right of significand !----------------------------------------------------------------- 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 ! do 190 j=1,i if (ibeta.eq.2) xmax = xmax + xmax if (ibeta.ne.2) xmax = xmax*beta 190 continue ! 200 return end double precision function ren(k) !*********************************************************************72 ! !! REN is a random number generator. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 09 January 2016 ! ! RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 BY PIKE AND ! HILL (MODIFIED BY HANSSON), COMMUNICATIONS OF THE ACM, ! VOL. 8, NO. 10, OCTOBER 1965. ! ! THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH ! FIXED POINT WORDLENGTH OF AT LEAST 29 BITS. IT IS ! BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST ! 29 BITS. ! implicit none integer iy, j, k 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 ( ) !*********************************************************************72 ! !! timestamp() prints out the current YMDHMS date as a timestamp. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 June 2014 ! ! Author: ! ! John Burkardt ! 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