C ALGORITHM 713, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 19, NO. 1, MARCH, 1993, P. 131. # for Sun FC = f77 FFLAGS = -O2 LDFLAGS = # for Convex #FC = fc #FFLAGS = -O2 -fi #LDFLAGS = -fi # for Cray #FC = cf77 #FFLAGS = #LDFLAGS = OBJ = testall.o vfnlib.o fnlib.o util.o DOBJ = dtestall.o dvfnlib.o dfnlib.o util.o all : testall dtestall testall : $(OBJ) $(FC) -o $@ $(LDFLAGS) $(OBJ) dtestall : $(DOBJ) $(FC) -o $@ $(LDFLAGS) $(DOBJ) .f : $(FC) -c $(FFLAGS) $< clean : rm -f testall dtestall core *.o Portable Vectorized Software for Bessel Function Evaluation Ronald F. Boisvert and Bonita V. Saunders Computing and Applied Mathematics Laboratory National Institute of Standards and Technology Gaithersburg, MD 20899 boisvert@cam.nist.gov saunders@cam.nist.gov June 19, 1991 -------- Contents -------- Makefile ...... A Unix makefile which generates the test programs. dfnlib.f ...... FNLIB (double precision) dtestall.f .... Test program (double precision) dvfnlib.f ..... VFNLIB (double precision) fnlib.f ....... FNLIB (single precision) util.f ........ Machine-dependent utilities. testall.f ..... Test program (single precision) vfnlib.f ...... VFNLIB (single precision) ----- Notes ----- The VFNLIB package is contained in the files vfnlib.f, dvfnlib.f, and util.f. The test programs testall.f and dtestall.f compare the values computed by VFNLIB with the original FNLIB package. The required FNLIB codes are contained in the files fnlib.f and dfnlib.f. Parts of the package which are machine-dependent are stored in the file util.f. The routines in this file are as follows: I1MACH These return machine-dependent constants according R1MACH to the PORT framework. They are set up for machines D1MACH using IEEE arithmetic, like the Sun. To modify them, simply comment out the Sun constants and un-comment the lines corresponding to your local machine. SECOND Returns the elapsed CP time in seconds since the beginning of the run in seconds. This is set up to make standard Unix system calls. On the Cray, this routine is supplied by the system, so it should be deleted from this file. ------------ Dependencies ------------ testall.f requires vfnlib.f fnlib.f util.f dtestall.f requires dvfnlib.f dfnlib.f util.f vfnlib.f requires util.f dvfnlib.f requires util.f fnlib.f requires util.f dfnlib.f requires util.f function besi0 (x) c april 1977 version. w. fullerton, c3, los alamos scientific lab. dimension bi0cs(12) c external alog, besi0e, csevl, exp, inits, r1mach, sqrt c c series for bi0 on the interval 0. to 9.00000d+00 c with weighted error 2.46e-18 c log weighted error 17.61 c significant figures required 17.90 c decimal places required 18.15 c data bi0 cs( 1) / -.0766054725 2839144951e0 / data bi0 cs( 2) / 1.9273379539 93808270e0 / data bi0 cs( 3) / .2282644586 920301339e0 / data bi0 cs( 4) / .0130489146 6707290428e0 / data bi0 cs( 5) / .0004344270 9008164874e0 / data bi0 cs( 6) / .0000094226 5768600193e0 / data bi0 cs( 7) / .0000001434 0062895106e0 / data bi0 cs( 8) / .0000000016 1384906966e0 / data bi0 cs( 9) / .0000000000 1396650044e0 / data bi0 cs(10) / .0000000000 0009579451e0 / data bi0 cs(11) / .0000000000 0000053339e0 / data bi0 cs(12) / .0000000000 0000000245e0 / c data nti0, xsml, xmax / 0, 0., 0. / c if (nti0.ne.0) go to 10 nti0 = inits (bi0cs, 12, 0.1*r1mach(3)) xsml = sqrt (4.0*r1mach(3)) xmax = alog (r1mach(2)) c 10 y = abs(x) if (y.gt.3.0) go to 20 c besi0 = 1.0 if (y.gt.xsml) besi0 = 2.75 + csevl (y*y/4.5-1.0, bi0cs, nti0) return c 20 if (y.gt.xmax) call seteru (34hbesi0 abs(x) so big i0 overflows, 1 34, 1, 2) c besi0 = exp(y) * besi0e(x) c return end function besi1 (x) c oct 1983 version. w. fullerton, c3, los alamos scientific lab. dimension bi1cs(11) c external alog, besi1e, csevl, exp, inits, r1mach, sqrt c c series for bi1 on the interval 0. to 9.00000d+00 c with weighted error 2.40e-17 c log weighted error 16.62 c significant figures required 16.23 c decimal places required 17.14 c data bi1 cs( 1) / -.0019717132 61099859e0 / data bi1 cs( 2) / .4073488766 7546481e0 / data bi1 cs( 3) / .0348389942 99959456e0 / data bi1 cs( 4) / .0015453945 56300123e0 / data bi1 cs( 5) / .0000418885 21098377e0 / data bi1 cs( 6) / .0000007649 02676483e0 / data bi1 cs( 7) / .0000000100 42493924e0 / data bi1 cs( 8) / .0000000000 99322077e0 / data bi1 cs( 9) / .0000000000 00766380e0 / data bi1 cs(10) / .0000000000 00004741e0 / data bi1 cs(11) / .0000000000 00000024e0 / c c data nti1, xmin, xsml, xmax / 0, 3*0. / c if (nti1.ne.0) go to 10 nti1 = inits (bi1cs, 11, 0.1*r1mach(3)) xmin = 2.0*r1mach(1) xsml = sqrt (8.0*r1mach(3)) xmax = alog (r1mach(2)) c 10 y = abs(x) if (y.gt.3.0) go to 20 c besi1 = 0.0 if (y.eq.0.0) return c if (y.le.xmin) call seteru ( 1 37hbesi1 abs(x) so small i1 underflows, 37, 1, 0) if (y.gt.xmin) besi1 = 0.5*x if (y.gt.xsml) besi1 = x * (.875 + csevl(y*y/4.5-1., bi1cs, nti1)) return c 20 if (y.gt.xmax) call seteru ( 1 34hbesi1 abs(x) so big i1 overflows, 34, 2, 2) c besi1 = exp(y) * besi1e(x) c return end function besj0 (x) c april 1977 version. w. fullerton, c3, los alamos scientific lab. dimension bj0cs(13), bm0cs(21), bth0cs(24) c external cos, csevl, inits, r1mach, sqrt c c series for bj0 on the interval 0. to 1.60000d+01 c with weighted error 7.47e-18 c log weighted error 17.13 c significant figures required 16.98 c decimal places required 17.68 c data bj0 cs( 1) / .1002541619 68939137e0 / data bj0 cs( 2) / -.6652230077 64405132e0 / data bj0 cs( 3) / .2489837034 98281314e0 / data bj0 cs( 4) / -.0332527231 700357697e0 / data bj0 cs( 5) / .0023114179 304694015e0 / data bj0 cs( 6) / -.0000991127 741995080e0 / data bj0 cs( 7) / .0000028916 708643998e0 / data bj0 cs( 8) / -.0000000612 108586630e0 / data bj0 cs( 9) / .0000000009 838650793e0 / data bj0 cs(10) / -.0000000000 124235515e0 / data bj0 cs(11) / .0000000000 001265433e0 / data bj0 cs(12) / -.0000000000 000010619e0 / data bj0 cs(13) / .0000000000 000000074e0 / c c series for bm0 on the interval 0. to 6.25000d-02 c with weighted error 4.98e-17 c log weighted error 16.30 c significant figures required 14.97 c decimal places required 16.96 c data bm0 cs( 1) / .0928496163 7381644e0 / data bm0 cs( 2) / -.0014298770 7403484e0 / data bm0 cs( 3) / .0000283057 9271257e0 / data bm0 cs( 4) / -.0000014330 0611424e0 / data bm0 cs( 5) / .0000001202 8628046e0 / data bm0 cs( 6) / -.0000000139 7113013e0 / data bm0 cs( 7) / .0000000020 4076188e0 / data bm0 cs( 8) / -.0000000003 5399669e0 / data bm0 cs( 9) / .0000000000 7024759e0 / data bm0 cs(10) / -.0000000000 1554107e0 / data bm0 cs(11) / .0000000000 0376226e0 / data bm0 cs(12) / -.0000000000 0098282e0 / data bm0 cs(13) / .0000000000 0027408e0 / data bm0 cs(14) / -.0000000000 0008091e0 / data bm0 cs(15) / .0000000000 0002511e0 / data bm0 cs(16) / -.0000000000 0000814e0 / data bm0 cs(17) / .0000000000 0000275e0 / data bm0 cs(18) / -.0000000000 0000096e0 / data bm0 cs(19) / .0000000000 0000034e0 / data bm0 cs(20) / -.0000000000 0000012e0 / data bm0 cs(21) / .0000000000 0000004e0 / c c series for bth0 on the interval 0. to 6.25000d-02 c with weighted error 3.67e-17 c log weighted error 16.44 c significant figures required 15.53 c decimal places required 17.13 c data bth0cs( 1) / -.2463916377 4300119e0 / data bth0cs( 2) / .0017370983 07508963e0 / data bth0cs( 3) / -.0000621836 33402968e0 / data bth0cs( 4) / .0000043680 50165742e0 / data bth0cs( 5) / -.0000004560 93019869e0 / data bth0cs( 6) / .0000000621 97400101e0 / data bth0cs( 7) / -.0000000103 00442889e0 / data bth0cs( 8) / .0000000019 79526776e0 / data bth0cs( 9) / -.0000000004 28198396e0 / data bth0cs(10) / .0000000001 02035840e0 / data bth0cs(11) / -.0000000000 26363898e0 / data bth0cs(12) / .0000000000 07297935e0 / data bth0cs(13) / -.0000000000 02144188e0 / data bth0cs(14) / .0000000000 00663693e0 / data bth0cs(15) / -.0000000000 00215126e0 / data bth0cs(16) / .0000000000 00072659e0 / data bth0cs(17) / -.0000000000 00025465e0 / data bth0cs(18) / .0000000000 00009229e0 / data bth0cs(19) / -.0000000000 00003448e0 / data bth0cs(20) / .0000000000 00001325e0 / data bth0cs(21) / -.0000000000 00000522e0 / data bth0cs(22) / .0000000000 00000210e0 / data bth0cs(23) / -.0000000000 00000087e0 / data bth0cs(24) / .0000000000 00000036e0 / c data pi4 / 0.7853981633 9744831e0 / data ntj0, ntm0, ntth0, xsml, xmax / 3*0, 2*0./ c if (ntj0.ne.0) go to 10 ntj0 = inits (bj0cs, 13, 0.1*r1mach(3)) ntm0 = inits (bm0cs, 21, 0.1*r1mach(3)) ntth0 = inits (bth0cs, 24, 0.1*r1mach(3)) c xsml = sqrt (4.0*r1mach(3)) xmax = 1.0/r1mach(4) c 10 y = abs(x) if (y.gt.4.0) go to 20 c besj0 = 1.0 if (y.gt.xsml) besj0 = csevl (.125*y*y-1., bj0cs, ntj0) return c 20 if (y.gt.xmax) call seteru ( 1 42hbesj0 no precision because abs(x) is big, 42, 1, 2) c z = 32.0/y**2 - 1.0 ampl = (0.75 + csevl (z, bm0cs, ntm0)) / sqrt(y) theta = y - pi4 + csevl (z, bth0cs, ntth0) / y besj0 = ampl * cos (theta) c return end function besj1 (x) c sept 1983 edition. w. fullerton, c3, los alamos scientific lab. dimension bj1cs(12), bm1cs(21), bth1cs(24) c external cos, csevl, inits, r1mach, sqrt c c series for bj1 on the interval 0. to 1.60000d+01 c with weighted error 4.48e-17 c log weighted error 16.35 c significant figures required 15.77 c decimal places required 16.89 c data bj1 cs( 1) / -.1172614151 3332787e0 / data bj1 cs( 2) / -.2536152183 0790640e0 / data bj1 cs( 3) / .0501270809 84469569e0 / data bj1 cs( 4) / -.0046315148 09625081e0 / data bj1 cs( 5) / .0002479962 29415914e0 / data bj1 cs( 6) / -.0000086789 48686278e0 / data bj1 cs( 7) / .0000002142 93917143e0 / data bj1 cs( 8) / -.0000000039 36093079e0 / data bj1 cs( 9) / .0000000000 55911823e0 / data bj1 cs(10) / -.0000000000 00632761e0 / data bj1 cs(11) / .0000000000 00005840e0 / data bj1 cs(12) / -.0000000000 00000044e0 / c c series for bm1 on the interval 0. to 6.25000d-02 c with weighted error 5.61e-17 c log weighted error 16.25 c significant figures required 14.97 c decimal places required 16.91 c data bm1 cs( 1) / .1047362510 931285e0 / data bm1 cs( 2) / .0044244389 3702345e0 / data bm1 cs( 3) / -.0000566163 9504035e0 / data bm1 cs( 4) / .0000023134 9417339e0 / data bm1 cs( 5) / -.0000001737 7182007e0 / data bm1 cs( 6) / .0000000189 3209930e0 / data bm1 cs( 7) / -.0000000026 5416023e0 / data bm1 cs( 8) / .0000000004 4740209e0 / data bm1 cs( 9) / -.0000000000 8691795e0 / data bm1 cs(10) / .0000000000 1891492e0 / data bm1 cs(11) / -.0000000000 0451884e0 / data bm1 cs(12) / .0000000000 0116765e0 / data bm1 cs(13) / -.0000000000 0032265e0 / data bm1 cs(14) / .0000000000 0009450e0 / data bm1 cs(15) / -.0000000000 0002913e0 / data bm1 cs(16) / .0000000000 0000939e0 / data bm1 cs(17) / -.0000000000 0000315e0 / data bm1 cs(18) / .0000000000 0000109e0 / data bm1 cs(19) / -.0000000000 0000039e0 / data bm1 cs(20) / .0000000000 0000014e0 / data bm1 cs(21) / -.0000000000 0000005e0 / c c series for bth1 on the interval 0. to 6.25000d-02 c with weighted error 4.10e-17 c log weighted error 16.39 c significant figures required 15.96 c decimal places required 17.08 c data bth1cs( 1) / .7406014102 6313850e0 / data bth1cs( 2) / -.0045717556 59637690e0 / data bth1cs( 3) / .0001198185 10964326e0 / data bth1cs( 4) / -.0000069645 61891648e0 / data bth1cs( 5) / .0000006554 95621447e0 / data bth1cs( 6) / -.0000000840 66228945e0 / data bth1cs( 7) / .0000000133 76886564e0 / data bth1cs( 8) / -.0000000024 99565654e0 / data bth1cs( 9) / .0000000005 29495100e0 / data bth1cs(10) / -.0000000001 24135944e0 / data bth1cs(11) / .0000000000 31656485e0 / data bth1cs(12) / -.0000000000 08668640e0 / data bth1cs(13) / .0000000000 02523758e0 / data bth1cs(14) / -.0000000000 00775085e0 / data bth1cs(15) / .0000000000 00249527e0 / data bth1cs(16) / -.0000000000 00083773e0 / data bth1cs(17) / .0000000000 00029205e0 / data bth1cs(18) / -.0000000000 00010534e0 / data bth1cs(19) / .0000000000 00003919e0 / data bth1cs(20) / -.0000000000 00001500e0 / data bth1cs(21) / .0000000000 00000589e0 / data bth1cs(22) / -.0000000000 00000237e0 / data bth1cs(23) / .0000000000 00000097e0 / data bth1cs(24) / -.0000000000 00000040e0 / c c data pi4 / 0.7853981633 9744831e0 / data ntj1, ntm1, ntth1, xsml, xmin, xmax / 3*0, 3*0./ c if (ntj1.ne.0) go to 10 ntj1 = inits (bj1cs, 12, 0.1*r1mach(3)) ntm1 = inits (bm1cs, 21, 0.1*r1mach(3)) ntth1 = inits (bth1cs, 24, 0.1*r1mach(3)) c xsml = sqrt (8.0*r1mach(3)) xmin = 2.0*r1mach(1) xmax = 1.0/r1mach(4) c 10 y = abs(x) if (y.gt.4.0) go to 20 c besj1 = 0.0 if (y.eq.0.0) return if (y.le.xmin) call seteru ( 1 37hbesj1 abs(x) so small j1 underflows, 37, 1, 0) if (y.gt.xmin) besj1 = 0.5*x if (y.gt.xsml) besj1 = x * (.25 + csevl(.125*y*y-1., bj1cs, ntj1)) return c 20 if (y.gt.xmax) call seteru ( 1 42hbesj1 no precision because abs(x) is big, 42, 2, 2) z = 32.0/y**2 - 1.0 ampl = (0.75 + csevl (z, bm1cs, ntm1)) / sqrt(y) theta = y - 3.0*pi4 + csevl (z, bth1cs, ntth1) / y besj1 = sign (ampl, x) * cos (theta) c return end function besk0 (x) c april 1977 version. w. fullerton, c3, los alamos scientific lab. dimension bk0cs(11) c external alog, besi0, besk0e, csevl, exp, inits, r1mach, sqrt c c series for bk0 on the interval 0. to 4.00000d+00 c with weighted error 3.57e-19 c log weighted error 18.45 c significant figures required 17.99 c decimal places required 18.97 c data bk0 cs( 1) / -.0353273932 3390276872e0 / data bk0 cs( 2) / .3442898999 246284869e0 / data bk0 cs( 3) / .0359799365 1536150163e0 / data bk0 cs( 4) / .0012646154 1144692592e0 / data bk0 cs( 5) / .0000228621 2103119451e0 / data bk0 cs( 6) / .0000002534 7910790261e0 / data bk0 cs( 7) / .0000000019 0451637722e0 / data bk0 cs( 8) / .0000000000 1034969525e0 / data bk0 cs( 9) / .0000000000 0004259816e0 / data bk0 cs(10) / .0000000000 0000013744e0 / data bk0 cs(11) / .0000000000 0000000035e0 / c data ntk0, xsml, xmax / 0, 0., 0. / c if (ntk0.ne.0) go to 10 ntk0 = inits (bk0cs, 11, 0.1*r1mach(3)) xsml = sqrt (4.0*r1mach(3)) xmax = -alog(r1mach(1)) xmax = xmax - 0.5*xmax*alog(xmax)/(xmax+0.5) - 0.01 c 10 if (x.le.0.) call seteru (29hbesk0 x is zero or negative, 29, 1 2, 2) if (x.gt.2.) go to 20 c y = 0. if (x.gt.xsml) y = x*x besk0 = -alog(0.5*x)*besi0(x) - .25 + csevl (.5*y-1., bk0cs, ntk0) return c 20 besk0 = 0. if (x.gt.xmax) call seteru (30hbesk0 x so big k0 underflows, 30, 1 1, 0) if (x.gt.xmax) return c besk0 = exp(-x) * besk0e(x) c return end function besk1 (x) c july 1980 version. w. fullerton, c3, los alamos scientific lab. dimension bk1cs(11) c external alog, besi1, besk1e, csevl, exp, inits, r1mach, sqrt c c series for bk1 on the interval 0. to 4.00000d+00 c with weighted error 7.02e-18 c log weighted error 17.15 c significant figures required 16.73 c decimal places required 17.67 c data bk1 cs( 1) / .0253002273 389477705e0 / data bk1 cs( 2) / -.3531559607 76544876e0 / data bk1 cs( 3) / -.1226111808 22657148e0 / data bk1 cs( 4) / -.0069757238 596398643e0 / data bk1 cs( 5) / -.0001730288 957513052e0 / data bk1 cs( 6) / -.0000024334 061415659e0 / data bk1 cs( 7) / -.0000000221 338763073e0 / data bk1 cs( 8) / -.0000000001 411488392e0 / data bk1 cs( 9) / -.0000000000 006666901e0 / data bk1 cs(10) / -.0000000000 000024274e0 / data bk1 cs(11) / -.0000000000 000000070e0 / c data ntk1, xmin, xsml, xmax /0, 0., 0., 0. / c if (ntk1.ne.0) go to 10 ntk1 = inits (bk1cs, 11, 0.1*r1mach(3)) xmin = exp (amax1(alog(r1mach(1)), -alog(r1mach(2))) + .01) xsml = sqrt (4.0*r1mach(3)) xmax = -alog(r1mach(1)) xmax = xmax - 0.5*xmax*alog(xmax)/(xmax+0.5) - 0.01 c 10 if (x.le.0.) call seteru (29hbesk1 x is zero or negative, 29, 1 2, 2) if (x.gt.2.0) go to 20 c if (x.lt.xmin) call seteru (31hbesk1 x so small k1 overflows, 1 31, 3, 2) y = 0. if (x.gt.xsml) y = x*x besk1 = alog(0.5*x)*besi1(x) + 1 (0.75 + csevl (.5*y-1., bk1cs, ntk1))/x return c 20 besk1 = 0. if (x.gt.xmax) call seteru (30hbesk1 x so big k1 underflows, 1 30, 1, 0) if (x.gt.xmax) return c besk1 = exp(-x) * besk1e(x) c return end function besy0 (x) c august 1980 version. w. fullerton, c3, los alamos scientific lab. dimension by0cs(13), bm0cs(21), bth0cs(24) c external alog, besj0, csevl, inits, r1mach, sin, sqrt c c series for by0 on the interval 0. to 1.60000d+01 c with weighted error 1.20e-17 c log weighted error 16.92 c significant figures required 16.15 c decimal places required 17.48 c data by0 cs( 1) / -.0112778393 92865573e0 / data by0 cs( 2) / -.1283452375 6042035e0 / data by0 cs( 3) / -.1043788479 9794249e0 / data by0 cs( 4) / .0236627491 83969695e0 / data by0 cs( 5) / -.0020903916 47700486e0 / data by0 cs( 6) / .0001039754 53939057e0 / data by0 cs( 7) / -.0000033697 47162423e0 / data by0 cs( 8) / .0000000772 93842676e0 / data by0 cs( 9) / -.0000000013 24976772e0 / data by0 cs(10) / .0000000000 17648232e0 / data by0 cs(11) / -.0000000000 00188105e0 / data by0 cs(12) / .0000000000 00001641e0 / data by0 cs(13) / -.0000000000 00000011e0 / c c series for bm0 on the interval 0. to 6.25000d-02 c with weighted error 4.98e-17 c log weighted error 16.30 c significant figures required 14.97 c decimal places required 16.96 c data bm0 cs( 1) / .0928496163 7381644e0 / data bm0 cs( 2) / -.0014298770 7403484e0 / data bm0 cs( 3) / .0000283057 9271257e0 / data bm0 cs( 4) / -.0000014330 0611424e0 / data bm0 cs( 5) / .0000001202 8628046e0 / data bm0 cs( 6) / -.0000000139 7113013e0 / data bm0 cs( 7) / .0000000020 4076188e0 / data bm0 cs( 8) / -.0000000003 5399669e0 / data bm0 cs( 9) / .0000000000 7024759e0 / data bm0 cs(10) / -.0000000000 1554107e0 / data bm0 cs(11) / .0000000000 0376226e0 / data bm0 cs(12) / -.0000000000 0098282e0 / data bm0 cs(13) / .0000000000 0027408e0 / data bm0 cs(14) / -.0000000000 0008091e0 / data bm0 cs(15) / .0000000000 0002511e0 / data bm0 cs(16) / -.0000000000 0000814e0 / data bm0 cs(17) / .0000000000 0000275e0 / data bm0 cs(18) / -.0000000000 0000096e0 / data bm0 cs(19) / .0000000000 0000034e0 / data bm0 cs(20) / -.0000000000 0000012e0 / data bm0 cs(21) / .0000000000 0000004e0 / c c series for bth0 on the interval 0. to 6.25000d-02 c with weighted error 3.67e-17 c log weighted error 16.44 c significant figures required 15.53 c decimal places required 17.13 c data bth0cs( 1) / -.2463916377 4300119e0 / data bth0cs( 2) / .0017370983 07508963e0 / data bth0cs( 3) / -.0000621836 33402968e0 / data bth0cs( 4) / .0000043680 50165742e0 / data bth0cs( 5) / -.0000004560 93019869e0 / data bth0cs( 6) / .0000000621 97400101e0 / data bth0cs( 7) / -.0000000103 00442889e0 / data bth0cs( 8) / .0000000019 79526776e0 / data bth0cs( 9) / -.0000000004 28198396e0 / data bth0cs(10) / .0000000001 02035840e0 / data bth0cs(11) / -.0000000000 26363898e0 / data bth0cs(12) / .0000000000 07297935e0 / data bth0cs(13) / -.0000000000 02144188e0 / data bth0cs(14) / .0000000000 00663693e0 / data bth0cs(15) / -.0000000000 00215126e0 / data bth0cs(16) / .0000000000 00072659e0 / data bth0cs(17) / -.0000000000 00025465e0 / data bth0cs(18) / .0000000000 00009229e0 / data bth0cs(19) / -.0000000000 00003448e0 / data bth0cs(20) / .0000000000 00001325e0 / data bth0cs(21) / -.0000000000 00000522e0 / data bth0cs(22) / .0000000000 00000210e0 / data bth0cs(23) / -.0000000000 00000087e0 / data bth0cs(24) / .0000000000 00000036e0 / c c twodpi = 2.0/pi data twodpi / 0.6366197723 6758134e0 / data pi4 / 0.7853981633 9744831e0 / data alnhaf / -0.6931471805 59945309e0 / data nty0, ntm0, ntth0, xsml, xmax / 3*0, 2*0./ c if (nty0.ne.0) go to 10 nty0 = inits (by0cs, 13, 0.1*r1mach(3)) ntm0 = inits (bm0cs, 21, 0.1*r1mach(3)) ntth0 = inits (bth0cs, 24, 0.1*r1mach(3)) c xsml = sqrt (4.0*r1mach(3)) xmax = 1.0/r1mach(4) c 10 if (x.le.0.) call seteru (29hbesy0 x is zero or negative, 29, 1 1, 2) if (x.gt.4.0) go to 20 c y = 0. if (x.gt.xsml) y = x*x besy0 = twodpi*(alnhaf+alog(x))*besj0(x) + .375 + 1 csevl (.125*y-1., by0cs, nty0) return c 20 if (x.gt.xmax) call seteru ( 1 37hbesy0 no precision because x is big, 37, 2, 2) c z = 32.0/x**2 - 1.0 ampl = (0.75 + csevl (z, bm0cs, ntm0)) / sqrt(x) theta = x - pi4 + csevl (z, bth0cs, ntth0) / x besy0 = ampl * sin (theta) c return end function besy1 (x) c april 1977 version. w. fullerton, c3, los alamos scientific lab. dimension by1cs(14), bm1cs(21), bth1cs(24) c external alog, besj1, csevl, exp, inits, r1mach, sin, sqrt c c series for by1 on the interval 0. to 1.60000d+01 c with weighted error 1.87e-18 c log weighted error 17.73 c significant figures required 17.83 c decimal places required 18.30 c data by1 cs( 1) / .0320804710 0611908629e0 / data by1 cs( 2) / 1.2627078974 33500450e0 / data by1 cs( 3) / .0064999618 9992317500e0 / data by1 cs( 4) / -.0893616452 8860504117e0 / data by1 cs( 5) / .0132508812 2175709545e0 / data by1 cs( 6) / -.0008979059 1196483523e0 / data by1 cs( 7) / .0000364736 1487958306e0 / data by1 cs( 8) / -.0000010013 7438166600e0 / data by1 cs( 9) / .0000000199 4539657390e0 / data by1 cs(10) / -.0000000003 0230656018e0 / data by1 cs(11) / .0000000000 0360987815e0 / data by1 cs(12) / -.0000000000 0003487488e0 / data by1 cs(13) / .0000000000 0000027838e0 / data by1 cs(14) / -.0000000000 0000000186e0 / c c series for bm1 on the interval 0. to 6.25000d-02 c with weighted error 5.61e-17 c log weighted error 16.25 c significant figures required 14.97 c decimal places required 16.91 c data bm1 cs( 1) / .1047362510 931285e0 / data bm1 cs( 2) / .0044244389 3702345e0 / data bm1 cs( 3) / -.0000566163 9504035e0 / data bm1 cs( 4) / .0000023134 9417339e0 / data bm1 cs( 5) / -.0000001737 7182007e0 / data bm1 cs( 6) / .0000000189 3209930e0 / data bm1 cs( 7) / -.0000000026 5416023e0 / data bm1 cs( 8) / .0000000004 4740209e0 / data bm1 cs( 9) / -.0000000000 8691795e0 / data bm1 cs(10) / .0000000000 1891492e0 / data bm1 cs(11) / -.0000000000 0451884e0 / data bm1 cs(12) / .0000000000 0116765e0 / data bm1 cs(13) / -.0000000000 0032265e0 / data bm1 cs(14) / .0000000000 0009450e0 / data bm1 cs(15) / -.0000000000 0002913e0 / data bm1 cs(16) / .0000000000 0000939e0 / data bm1 cs(17) / -.0000000000 0000315e0 / data bm1 cs(18) / .0000000000 0000109e0 / data bm1 cs(19) / -.0000000000 0000039e0 / data bm1 cs(20) / .0000000000 0000014e0 / data bm1 cs(21) / -.0000000000 0000005e0 / c c series for bth1 on the interval 0. to 6.25000d-02 c with weighted error 4.10e-17 c log weighted error 16.39 c significant figures required 15.96 c decimal places required 17.08 c data bth1cs( 1) / .7406014102 6313850e0 / data bth1cs( 2) / -.0045717556 59637690e0 / data bth1cs( 3) / .0001198185 10964326e0 / data bth1cs( 4) / -.0000069645 61891648e0 / data bth1cs( 5) / .0000006554 95621447e0 / data bth1cs( 6) / -.0000000840 66228945e0 / data bth1cs( 7) / .0000000133 76886564e0 / data bth1cs( 8) / -.0000000024 99565654e0 / data bth1cs( 9) / .0000000005 29495100e0 / data bth1cs(10) / -.0000000001 24135944e0 / data bth1cs(11) / .0000000000 31656485e0 / data bth1cs(12) / -.0000000000 08668640e0 / data bth1cs(13) / .0000000000 02523758e0 / data bth1cs(14) / -.0000000000 00775085e0 / data bth1cs(15) / .0000000000 00249527e0 / data bth1cs(16) / -.0000000000 00083773e0 / data bth1cs(17) / .0000000000 00029205e0 / data bth1cs(18) / -.0000000000 00010534e0 / data bth1cs(19) / .0000000000 00003919e0 / data bth1cs(20) / -.0000000000 00001500e0 / data bth1cs(21) / .0000000000 00000589e0 / data bth1cs(22) / -.0000000000 00000237e0 / data bth1cs(23) / .0000000000 00000097e0 / data bth1cs(24) / -.0000000000 00000040e0 / c c twodpi = 2.0/pi data twodpi / 0.6366197723 6758134e0 / data pi4 / 0.7853981633 9744831e0 / data nty1, ntm1, ntth1, xmin, xsml, xmax / 3*0, 3*0./ c if (nty1.ne.0) go to 10 nty1 = inits (by1cs, 14, 0.1*r1mach(3)) ntm1 = inits (bm1cs, 21, 0.1*r1mach(3)) ntth1 = inits (bth1cs, 24, 0.1*r1mach(3)) c xmin = 1.571*exp ( amax1(alog(r1mach(1)), -alog(r1mach(2)))+.01) xsml = sqrt (4.0*r1mach(3)) xmax = 1.0/r1mach(4) c 10 if (x.le.0.) call seteru (29hbesy1 x is zero or negative, 29, 1 1, 2) if (x.gt.4.0) go to 20 c if (x.lt.xmin) call seteru (31hbesy1 x so small y1 overflows, 1 31, 3, 2) y = 0. if (x.gt.xsml) y = x*x besy1 = twodpi*alog(0.5*x)*besj1(x) + 1 (0.5 + csevl (.125*y-1., by1cs, nty1))/x return c 20 if (x.gt.xmax) call seteru ( 1 37hbesy1 no precision because x is big, 37, 2, 2) c z = 32.0/x**2 - 1.0 ampl = (0.75 + csevl (z, bm1cs, ntm1)) / sqrt(x) theta = x - 3.0*pi4 + csevl (z, bth1cs, ntth1) / x besy1 = ampl * sin (theta) c return end function inits (os, nos, eta) c april 1977 version. w. fullerton, c3, los alamos scientific lab. c c initialize the orthogonal series so that inits is the number of terms c needed to insure the error is no larger than eta. ordinarily, eta c will be chosen to be one-tenth machine precision. c c input arguments -- c os array of nos coefficients in an orthogonal series. c nos number of coefficients in os. c eta requested accuracy of series. c dimension os(nos) c if (nos.lt.1) call seteru ( 1 35hinits number of coefficients lt 1, 35, 2, 2) c err = 0. do 10 ii=1,nos i = nos + 1 - ii err = err + abs(os(i)) if (err.gt.eta) go to 20 10 continue c 20 if (i.eq.nos) call seteru (28hinits eta may be too small, 28, 1 1, 2) inits = i c return end function csevl (x, cs, n) c april 1977 version. w. fullerton, c3, los alamos scientific lab. c c evaluate the n-term chebyshev series cs at x. adapted from c r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973). also see fox c and parker, chebyshev polys in numerical analysis, oxford press, p.56. c c input arguments -- c x value at which the series is to be evaluated. c cs array of n terms of a chebyshev series. in eval- c uating cs, only half the first coef is summed. c n number of terms in array cs. c dimension cs(1) c if (n.lt.1) call seteru (28hcsevl number of terms le 0, 28, 2,2) if (n.gt.1000) call seteru (31hcsevl number of terms gt 1000, 1 31, 3, 2) if (x.lt.(-1.1) .or. x.gt.1.1) call seteru ( 1 25hcsevl x outside (-1,+1), 25, 1, 1) c b1 = 0. b0 = 0. twox = 2.*x do 10 i=1,n b2 = b1 b1 = b0 ni = n + 1 - i b0 = twox*b1 - b2 + cs(ni) 10 continue c csevl = 0.5 * (b0-b2) c return end subroutine seteru (messg, nmessg, nerr, iopt) common /cseter/ iunflo integer messg(1) data iunflo / 0 / c if (iopt.ne.0) call seterr (messg, nmessg, nerr, iopt) if (iopt.ne.0) return c if (iunflo.le.0) return call seterr (messg, nmessg, nerr, 1) c return end function besi0e (x) c july 1977 version. w. fullerton, c3, los alamos scientific lab. dimension bi0cs(12), ai0cs(21), ai02cs(22) c external csevl, exp, inits, r1mach, sqrt c c series for bi0 on the interval 0. to 9.00000d+00 c with weighted error 2.46e-18 c log weighted error 17.61 c significant figures required 17.90 c decimal places required 18.15 c data bi0 cs( 1) / -.0766054725 2839144951e0 / data bi0 cs( 2) / 1.9273379539 93808270e0 / data bi0 cs( 3) / .2282644586 920301339e0 / data bi0 cs( 4) / .0130489146 6707290428e0 / data bi0 cs( 5) / .0004344270 9008164874e0 / data bi0 cs( 6) / .0000094226 5768600193e0 / data bi0 cs( 7) / .0000001434 0062895106e0 / data bi0 cs( 8) / .0000000016 1384906966e0 / data bi0 cs( 9) / .0000000000 1396650044e0 / data bi0 cs(10) / .0000000000 0009579451e0 / data bi0 cs(11) / .0000000000 0000053339e0 / data bi0 cs(12) / .0000000000 0000000245e0 / c c series for ai0 on the interval 1.25000d-01 to 3.33333d-01 c with weighted error 7.87e-17 c log weighted error 16.10 c significant figures required 14.69 c decimal places required 16.76 c data ai0 cs( 1) / .0757599449 4023796e0 / data ai0 cs( 2) / .0075913808 1082334e0 / data ai0 cs( 3) / .0004153131 3389237e0 / data ai0 cs( 4) / .0000107007 6463439e0 / data ai0 cs( 5) / -.0000079011 7997921e0 / data ai0 cs( 6) / -.0000007826 1435014e0 / data ai0 cs( 7) / .0000002783 8499429e0 / data ai0 cs( 8) / .0000000082 5247260e0 / data ai0 cs( 9) / -.0000000120 4463945e0 / data ai0 cs(10) / .0000000015 5964859e0 / data ai0 cs(11) / .0000000002 2925563e0 / data ai0 cs(12) / -.0000000001 1916228e0 / data ai0 cs(13) / .0000000000 1757854e0 / data ai0 cs(14) / .0000000000 0112822e0 / data ai0 cs(15) / -.0000000000 0114684e0 / data ai0 cs(16) / .0000000000 0027155e0 / data ai0 cs(17) / -.0000000000 0002415e0 / data ai0 cs(18) / -.0000000000 0000608e0 / data ai0 cs(19) / .0000000000 0000314e0 / data ai0 cs(20) / -.0000000000 0000071e0 / data ai0 cs(21) / .0000000000 0000007e0 / c c series for ai02 on the interval 0. to 1.25000d-01 c with weighted error 3.79e-17 c log weighted error 16.42 c significant figures required 14.86 c decimal places required 17.09 c data ai02cs( 1) / .0544904110 1410882e0 / data ai02cs( 2) / .0033691164 7825569e0 / data ai02cs( 3) / .0000688975 8346918e0 / data ai02cs( 4) / .0000028913 7052082e0 / data ai02cs( 5) / .0000002048 9185893e0 / data ai02cs( 6) / .0000000226 6668991e0 / data ai02cs( 7) / .0000000033 9623203e0 / data ai02cs( 8) / .0000000004 9406022e0 / data ai02cs( 9) / .0000000000 1188914e0 / data ai02cs(10) / -.0000000000 3149915e0 / data ai02cs(11) / -.0000000000 1321580e0 / data ai02cs(12) / -.0000000000 0179419e0 / data ai02cs(13) / .0000000000 0071801e0 / data ai02cs(14) / .0000000000 0038529e0 / data ai02cs(15) / .0000000000 0001539e0 / data ai02cs(16) / -.0000000000 0004151e0 / data ai02cs(17) / -.0000000000 0000954e0 / data ai02cs(18) / .0000000000 0000382e0 / data ai02cs(19) / .0000000000 0000176e0 / data ai02cs(20) / -.0000000000 0000034e0 / data ai02cs(21) / -.0000000000 0000027e0 / data ai02cs(22) / .0000000000 0000003e0 / c data nti0, ntai0, ntai02, xsml / 3*0, 0. / c if (nti0.ne.0) go to 10 nti0 = inits (bi0cs, 12, 0.1*r1mach(3)) ntai0 = inits (ai0cs, 21, 0.1*r1mach(3)) ntai02 = inits (ai02cs, 22, 0.1*r1mach(3)) xsml = sqrt (4.0*r1mach(3)) c 10 y = abs(x) if (y.gt.3.0) go to 20 c besi0e = 1.0 if (y.gt.xsml) besi0e = exp(-y) * ( 2.75 + 1 csevl (y*y/4.5-1.0, bi0cs, nti0) ) return c 20 if (y.le.8.) besi0e = (.375 + csevl ((48./y-11.)/5., ai0cs, ntai0) 1 ) / sqrt(y) if (y.gt.8.) besi0e = (.375 + csevl (16./y-1., ai02cs, ntai02)) 1 / sqrt(y) c return end function besi1e (x) c oct 1983 version. w. fullerton, c3, los alamos scientific lab. dimension bi1cs(11), ai1cs(21), ai12cs(22) c external csevl, exp, inits, r1mach, sqrt c c series for bi1 on the interval 0. to 9.00000d+00 c with weighted error 2.40e-17 c log weighted error 16.62 c significant figures required 16.23 c decimal places required 17.14 c data bi1 cs( 1) / -.0019717132 61099859e0 / data bi1 cs( 2) / .4073488766 7546481e0 / data bi1 cs( 3) / .0348389942 99959456e0 / data bi1 cs( 4) / .0015453945 56300123e0 / data bi1 cs( 5) / .0000418885 21098377e0 / data bi1 cs( 6) / .0000007649 02676483e0 / data bi1 cs( 7) / .0000000100 42493924e0 / data bi1 cs( 8) / .0000000000 99322077e0 / data bi1 cs( 9) / .0000000000 00766380e0 / data bi1 cs(10) / .0000000000 00004741e0 / data bi1 cs(11) / .0000000000 00000024e0 / c c series for ai1 on the interval 1.25000d-01 to 3.33333d-01 c with weighted error 6.98e-17 c log weighted error 16.16 c significant figures required 14.53 c decimal places required 16.82 c data ai1 cs( 1) / -.0284674418 1881479e0 / data ai1 cs( 2) / -.0192295323 1443221e0 / data ai1 cs( 3) / -.0006115185 8579437e0 / data ai1 cs( 4) / -.0000206997 1253350e0 / data ai1 cs( 5) / .0000085856 1914581e0 / data ai1 cs( 6) / .0000010494 9824671e0 / data ai1 cs( 7) / -.0000002918 3389184e0 / data ai1 cs( 8) / -.0000000155 9378146e0 / data ai1 cs( 9) / .0000000131 8012367e0 / data ai1 cs(10) / -.0000000014 4842341e0 / data ai1 cs(11) / -.0000000002 9085122e0 / data ai1 cs(12) / .0000000001 2663889e0 / data ai1 cs(13) / -.0000000000 1664947e0 / data ai1 cs(14) / -.0000000000 0166665e0 / data ai1 cs(15) / .0000000000 0124260e0 / data ai1 cs(16) / -.0000000000 0027315e0 / data ai1 cs(17) / .0000000000 0002023e0 / data ai1 cs(18) / .0000000000 0000730e0 / data ai1 cs(19) / -.0000000000 0000333e0 / data ai1 cs(20) / .0000000000 0000071e0 / data ai1 cs(21) / -.0000000000 0000006e0 / c c series for ai12 on the interval 0. to 1.25000d-01 c with weighted error 3.55e-17 c log weighted error 16.45 c significant figures required 14.69 c decimal places required 17.12 c data ai12cs( 1) / .0285762350 1828014e0 / data ai12cs( 2) / -.0097610974 9136147e0 / data ai12cs( 3) / -.0001105889 3876263e0 / data ai12cs( 4) / -.0000038825 6480887e0 / data ai12cs( 5) / -.0000002512 2362377e0 / data ai12cs( 6) / -.0000000263 1468847e0 / data ai12cs( 7) / -.0000000038 3538039e0 / data ai12cs( 8) / -.0000000005 5897433e0 / data ai12cs( 9) / -.0000000000 1897495e0 / data ai12cs(10) / .0000000000 3252602e0 / data ai12cs(11) / .0000000000 1412580e0 / data ai12cs(12) / .0000000000 0203564e0 / data ai12cs(13) / -.0000000000 0071985e0 / data ai12cs(14) / -.0000000000 0040836e0 / data ai12cs(15) / -.0000000000 0002101e0 / data ai12cs(16) / .0000000000 0004273e0 / data ai12cs(17) / .0000000000 0001041e0 / data ai12cs(18) / -.0000000000 0000382e0 / data ai12cs(19) / -.0000000000 0000186e0 / data ai12cs(20) / .0000000000 0000033e0 / data ai12cs(21) / .0000000000 0000028e0 / data ai12cs(22) / -.0000000000 0000003e0 / c data nti1, ntai1, ntai12, xmin, xsml / 3*0, 2*0. / c if (nti1.ne.0) go to 10 nti1 = inits (bi1cs, 11, 0.1*r1mach(3)) ntai1 = inits (ai1cs, 21, 0.1*r1mach(3)) ntai12 = inits (ai12cs, 22, 0.1*r1mach(3)) c xmin = 2.0*r1mach(1) xsml = sqrt (8.0*r1mach(3)) c 10 y = abs(x) if (y.gt.3.0) go to 20 c besi1e = 0.0 if (y.eq.0.0) return c if (y.le.xmin) call seteru ( 1 37hbesi1e abs(x) so small i1 underflows, 37, 1, 0) besi1e = 0.0 if (y.gt.xmin) besi1e = 0.5*x if (y.gt.xsml) besi1e = x * (.875 + csevl(y*y/4.5-1., bi1cs,nti1)) besi1e = exp(-y) * besi1e return c 20 if (y.le.8.) besi1e = (.375 + csevl ((48./y-11.)/5., ai1cs, ntai1) 1 ) / sqrt(y) if (y.gt.8.) besi1e = (.375 + csevl (16./y-1.0, ai12cs, ntai12)) 1 / sqrt(y) besi1e = sign (besi1e, x) c return end function besk0e (x) c july 1980 version. w. fullerton, c3, los alamos scientific lab. dimension bk0cs(11), ak0cs(17), ak02cs(14) c external alog, besi0, csevl, exp, inits, r1mach, sqrt c c series for bk0 on the interval 0. to 4.00000d+00 c with weighted error 3.57e-19 c log weighted error 18.45 c significant figures required 17.99 c decimal places required 18.97 c data bk0 cs( 1) / -.0353273932 3390276872e0 / data bk0 cs( 2) / .3442898999 246284869e0 / data bk0 cs( 3) / .0359799365 1536150163e0 / data bk0 cs( 4) / .0012646154 1144692592e0 / data bk0 cs( 5) / .0000228621 2103119451e0 / data bk0 cs( 6) / .0000002534 7910790261e0 / data bk0 cs( 7) / .0000000019 0451637722e0 / data bk0 cs( 8) / .0000000000 1034969525e0 / data bk0 cs( 9) / .0000000000 0004259816e0 / data bk0 cs(10) / .0000000000 0000013744e0 / data bk0 cs(11) / .0000000000 0000000035e0 / c c series for ak0 on the interval 1.25000d-01 to 5.00000d-01 c with weighted error 5.34e-17 c log weighted error 16.27 c significant figures required 14.92 c decimal places required 16.89 c data ak0 cs( 1) / -.0764394790 3327941e0 / data ak0 cs( 2) / -.0223565260 5699819e0 / data ak0 cs( 3) / .0007734181 1546938e0 / data ak0 cs( 4) / -.0000428100 6688886e0 / data ak0 cs( 5) / .0000030817 0017386e0 / data ak0 cs( 6) / -.0000002639 3672220e0 / data ak0 cs( 7) / .0000000256 3713036e0 / data ak0 cs( 8) / -.0000000027 4270554e0 / data ak0 cs( 9) / .0000000003 1694296e0 / data ak0 cs(10) / -.0000000000 3902353e0 / data ak0 cs(11) / .0000000000 0506804e0 / data ak0 cs(12) / -.0000000000 0068895e0 / data ak0 cs(13) / .0000000000 0009744e0 / data ak0 cs(14) / -.0000000000 0001427e0 / data ak0 cs(15) / .0000000000 0000215e0 / data ak0 cs(16) / -.0000000000 0000033e0 / data ak0 cs(17) / .0000000000 0000005e0 / c c series for ak02 on the interval 0. to 1.25000d-01 c with weighted error 2.34e-17 c log weighted error 16.63 c significant figures required 14.67 c decimal places required 17.20 c data ak02cs( 1) / -.0120186982 6307592e0 / data ak02cs( 2) / -.0091748526 9102569e0 / data ak02cs( 3) / .0001444550 9317750e0 / data ak02cs( 4) / -.0000040136 1417543e0 / data ak02cs( 5) / .0000001567 8318108e0 / data ak02cs( 6) / -.0000000077 7011043e0 / data ak02cs( 7) / .0000000004 6111825e0 / data ak02cs( 8) / -.0000000000 3158592e0 / data ak02cs( 9) / .0000000000 0243501e0 / data ak02cs(10) / -.0000000000 0020743e0 / data ak02cs(11) / .0000000000 0001925e0 / data ak02cs(12) / -.0000000000 0000192e0 / data ak02cs(13) / .0000000000 0000020e0 / data ak02cs(14) / -.0000000000 0000002e0 / c data ntk0, ntak0, ntak02, xsml / 3*0, 0. / c if (ntk0.ne.0) go to 10 ntk0 = inits (bk0cs, 11, 0.1*r1mach(3)) ntak0 = inits (ak0cs, 17, 0.1*r1mach(3)) ntak02 = inits (ak02cs, 14, 0.1*r1mach(3)) xsml = sqrt (4.0*r1mach(3)) c 10 if (x.le.0.) call seteru (29hbesk0e x is zero or negative, 29, 1 1, 2) if (x.gt.2.) go to 20 c y = 0. if (x.gt.xsml) y = x*x besk0e = exp(x) * (-alog(0.5*x)*besi0(x) 1 - .25 + csevl (.5*y-1., bk0cs, ntk0) ) return c 20 if (x.le.8.) besk0e = (1.25 + csevl ((16./x-5.)/3., ak0cs, ntak0)) 1 / sqrt(x) if (x.gt.8.) besk0e = (1.25 + csevl (16./x-1., ak02cs, ntak02)) 1 / sqrt(x) c return end function besk1e (x) c july 1980 version. w. fullerton, c3, los alamos scientific lab. dimension bk1cs(11), ak1cs(17), ak12cs(14) c external alog, besi1, csevl, exp, inits, r1mach, sqrt c c series for bk1 on the interval 0. to 4.00000d+00 c with weighted error 7.02e-18 c log weighted error 17.15 c significant figures required 16.73 c decimal places required 17.67 c data bk1 cs( 1) / .0253002273 389477705e0 / data bk1 cs( 2) / -.3531559607 76544876e0 / data bk1 cs( 3) / -.1226111808 22657148e0 / data bk1 cs( 4) / -.0069757238 596398643e0 / data bk1 cs( 5) / -.0001730288 957513052e0 / data bk1 cs( 6) / -.0000024334 061415659e0 / data bk1 cs( 7) / -.0000000221 338763073e0 / data bk1 cs( 8) / -.0000000001 411488392e0 / data bk1 cs( 9) / -.0000000000 006666901e0 / data bk1 cs(10) / -.0000000000 000024274e0 / data bk1 cs(11) / -.0000000000 000000070e0 / c c series for ak1 on the interval 1.25000d-01 to 5.00000d-01 c with weighted error 6.06e-17 c log weighted error 16.22 c significant figures required 15.41 c decimal places required 16.83 c data ak1 cs( 1) / .2744313406 973883e0 / data ak1 cs( 2) / .0757198995 3199368e0 / data ak1 cs( 3) / -.0014410515 5647540e0 / data ak1 cs( 4) / .0000665011 6955125e0 / data ak1 cs( 5) / -.0000043699 8470952e0 / data ak1 cs( 6) / .0000003540 2774997e0 / data ak1 cs( 7) / -.0000000331 1163779e0 / data ak1 cs( 8) / .0000000034 4597758e0 / data ak1 cs( 9) / -.0000000003 8989323e0 / data ak1 cs(10) / .0000000000 4720819e0 / data ak1 cs(11) / -.0000000000 0604783e0 / data ak1 cs(12) / .0000000000 0081284e0 / data ak1 cs(13) / -.0000000000 0011386e0 / data ak1 cs(14) / .0000000000 0001654e0 / data ak1 cs(15) / -.0000000000 0000248e0 / data ak1 cs(16) / .0000000000 0000038e0 / data ak1 cs(17) / -.0000000000 0000006e0 / c c series for ak12 on the interval 0. to 1.25000d-01 c with weighted error 2.58e-17 c log weighted error 16.59 c significant figures required 15.22 c decimal places required 17.16 c data ak12cs( 1) / .0637930834 3739001e0 / data ak12cs( 2) / .0283288781 3049721e0 / data ak12cs( 3) / -.0002475370 6739052e0 / data ak12cs( 4) / .0000057719 7245160e0 / data ak12cs( 5) / -.0000002068 9392195e0 / data ak12cs( 6) / .0000000097 3998344e0 / data ak12cs( 7) / -.0000000005 5853361e0 / data ak12cs( 8) / .0000000000 3732996e0 / data ak12cs( 9) / -.0000000000 0282505e0 / data ak12cs(10) / .0000000000 0023720e0 / data ak12cs(11) / -.0000000000 0002176e0 / data ak12cs(12) / .0000000000 0000215e0 / data ak12cs(13) / -.0000000000 0000022e0 / data ak12cs(14) / .0000000000 0000002e0 / c data ntk1, ntak1, ntak12, xmin, xsml / 3*0, 2*0. / c if (ntk1.ne.0) go to 10 ntk1 = inits (bk1cs, 11, 0.1*r1mach(3)) ntak1 = inits (ak1cs, 17, 0.1*r1mach(3)) ntak12 = inits (ak12cs, 14, 0.1*r1mach(3)) c xmin = exp (amax1(alog(r1mach(1)), -alog(r1mach(2))) + .01) xsml = sqrt (4.0*r1mach(3)) c 10 if (x.le.0.) call seteru (29hbesk1e x is zero or negative, 29, 1 1, 2) if (x.gt.2.0) go to 20 c if (x.lt.xmin) call seteru (31hbesk1e x so small k1 overflows, 1 31, 2, 2) y = 0. if (x.gt.xsml) y = x*x besk1e = exp(x) * (alog(0.5*x)*besi1(x) + 1 (0.75 + csevl (.5*y-1., bk1cs, ntk1))/x ) return c 20 if (x.le.8.) besk1e = (1.25 + csevl ((16./x-5.)/3., ak1cs, ntak1)) 1 / sqrt(x) if (x.gt.8.) besk1e = (1.25 + csevl (16./x-1., ak12cs, ntak12)) 1 / sqrt(x) c return end integer function i8save(isw,ivalue,set) c c if (isw = 1) i8save returns the current error number and c sets it to ivalue if set = .true. . c c if (isw = 2) i8save returns the current recovery switch and c sets it to ivalue if set = .true. . c logical set c integer iparam(2) c iparam(1) is the error number and iparam(2) is the recovery switch. c c start execution error free and with recovery turned off. c data iparam(1) /0/, iparam(2) /2/ c i8save=iparam(isw) if (set) iparam(isw)=ivalue c return c end subroutine eprint c c this subroutine prints the last error message, if any. c integer messg(1) c call e9rint(messg,1,1,.false.) return c end subroutine s88fmt( n, w, ifmt ) c c s88fmt replaces ifmt(1), ... , ifmt(n) with c the characters corresponding to the n least significant c digits of w. c integer n,w,ifmt(n) c integer nt,wt,digits(10) c data digits( 1) / 1h0 / data digits( 2) / 1h1 / data digits( 3) / 1h2 / data digits( 4) / 1h3 / data digits( 5) / 1h4 / data digits( 6) / 1h5 / data digits( 7) / 1h6 / data digits( 8) / 1h7 / data digits( 9) / 1h8 / data digits(10) / 1h9 / c nt = n wt = w c 10 if (nt .le. 0) return idigit = mod( wt, 10 ) ifmt(nt) = digits(idigit+1) wt = wt/10 nt = nt - 1 go to 10 c end subroutine seterr (messg, nmessg, nerr, iopt) c c this version modified by w. fullerton to dump if iopt = 1 and c not recovering. c seterr sets lerror = nerr, optionally prints the message and dumps c according to the following rules... c c if iopt = 1 and recovering - just remember the error. c if iopt = 1 and not recovering - print, dump and stop. c if iopt = 2 - print, dump and stop. c c input c c messg - the error message. c nmessg - the length of the message, in characters. c nerr - the error number. must have nerr non-zero. c iopt - the option. must have iopt=1 or 2. c c error states - c c 1 - message length not positive. c 2 - cannot have nerr=0. c 3 - an unrecovered error followed by another error. c 4 - bad value for iopt. c c only the first 72 characters of the message are printed. c c the error handler calls a subroutine named fdump to produce a c symbolic dump. to complete the package, a dummy version of fdump c is supplied, but it should be replaced by a locally written version c which at least gives a trace-back. c integer messg(1) external i1mach, i8save c c the unit for error messages. c iwunit=i1mach(4) c if (nmessg.ge.1) go to 10 c c a message of non-positive length is fatal. c write(iwunit,9000) 9000 format(52h1error 1 in seterr - message length not positive.) go to 60 c c nw is the number of words the message occupies. c 10 nw=(min0(nmessg,72)-1)/i1mach(6)+1 c if (nerr.ne.0) go to 20 c c cannot turn the error state off using seterr. c write(iwunit,9001) 9001 format(42h1error 2 in seterr - cannot have nerr=0// 1 34h the current error message follows///) call e9rint(messg,nw,nerr,.true.) itemp=i8save(1,1,.true.) go to 50 c c set lerror and test for a previous unrecovered error. c 20 if (i8save(1,nerr,.true.).eq.0) go to 30 c write(iwunit,9002) 9002 format(23h1error 3 in seterr -, 1 48h an unrecovered error followed by another error.// 2 48h the previous and current error messages follow.///) call eprint call e9rint(messg,nw,nerr,.true.) go to 50 c c save this message in case it is not recovered from properly. c 30 call e9rint(messg,nw,nerr,.true.) c if (iopt.eq.1 .or. iopt.eq.2) go to 40 c c must have iopt = 1 or 2. c write(iwunit,9003) 9003 format(42h1error 4 in seterr - bad value for iopt// 1 34h the current error message follows///) go to 50 c c test for recovery. c 40 if (iopt.eq.2) go to 50 c if (i8save(2,0,.false.).eq.1) return c c call eprint c stop c 50 call eprint 60 call fdump stop c end subroutine e9rint(messg,nw,nerr,save) c c this routine stores the current error message or prints the old one, c if any, depending on whether or not save = .true. . c integer messg(nw) logical save external i1mach, i8save c c messgp stores at least the first 72 characters of the previous c message. its length is machine dependent and must be at least c c 1 + 71/(the number of characters stored per integer word). c integer messgp(36),fmt(14),ccplus c c start with no previous message. c data messgp(1)/1h1/, nwp/0/, nerrp/0/ c c set up the format for printing the error message. c the format is simply (a1,14x,72axx) where xx=i1mach(6) is the c number of characters stored per integer word. c data ccplus / 1h+ / c data fmt( 1) / 1h( / data fmt( 2) / 1ha / data fmt( 3) / 1h1 / data fmt( 4) / 1h, / data fmt( 5) / 1h1 / data fmt( 6) / 1h4 / data fmt( 7) / 1hx / data fmt( 8) / 1h, / data fmt( 9) / 1h7 / data fmt(10) / 1h2 / data fmt(11) / 1ha / data fmt(12) / 1hx / data fmt(13) / 1hx / data fmt(14) / 1h) / c if (.not.save) go to 20 c c save the message. c nwp=nw nerrp=nerr do 10 i=1,nw 10 messgp(i)=messg(i) c go to 30 c 20 if (i8save(1,0,.false.).eq.0) go to 30 c c print the message. c iwunit=i1mach(4) write(iwunit,9000) nerrp 9000 format(7h error ,i4,4h in ) c call s88fmt(2,i1mach(6),fmt(12)) write(iwunit,fmt) ccplus,(messgp(i),i=1,nwp) c 30 return c end subroutine fdump return end C C --------------------------------------- C Compare scalar fnlib routines to vfnlib C --------------------------------------- C C Ron Boisvert and Bonita Saunders C Computing and Applied Mathematics Laboratory C National Institute of Standards and Technology C Gaithersburg, MD 20899 C C boisvert@cam.nist.gov C saunders@cam.nist.gov C C 21 Feb 91 C C---------------------------------------------------------------------- C INTEGER M PARAMETER ( M=2000 ) C REAL XLO, XHI, X, F, FX, WORK INTEGER I, IWORK, J C DIMENSION XLO(3,8), XHI(3,8), X(M), F(M), FX(M), WORK(7*M), + IWORK(M) C SAVE XLO, XHI C EXTERNAL VI0, BESI0, VI1, BESI1, VJ0, BESJ0, VJ1, BESJ1, + VK0, BESK0, VK1, BESK1, VY0, BESY0, VY1, BESY1 C DATA ((XLO(I,J),XHI(I,J),I=1,3),J=1,8) / + 0.0E0, 3.0E0 , -8.0E0, -3.0E0, 10.0E0, 50.0E0, + 0.0E0, 3.0E0 , -8.0E0, -3.0E0, 10.0E0, 50.0E0, + 0.0E0, 4.0E0 , -8.0E0, -4.0E0, 8.0E0, 50.0E0, + 0.0E0, 4.0E0 , -8.0E0, -4.0E0, 8.0E0, 50.0E0, + 0.0E0, 2.0E0 , 2.0E0, 8.0E0, 10.0E0, 50.0E0, + 0.0E0, 2.0E0 , 2.0E0, 8.0E0, 10.0E0, 50.0E0, + 0.0E0, 1.0E-3, 1.0E0, 4.0E0, 4.0E0, 50.0E0, + 0.0E0, 1.0E-3, 1.0E0, 4.0E0, 4.0E0, 50.0E0 / C C---------------------------------------------------------------------- C C ... TEST EACH USER-CALLABLE ROUTINE C CALL TEST('I0',BESI0,VI0,XLO(1,1),XHI(1,1),M,X,F,FX,WORK,IWORK) CALL TEST('I1',BESI1,VI1,XLO(1,2),XHI(1,2),M,X,F,FX,WORK,IWORK) CALL TEST('J0',BESJ0,VJ0,XLO(1,3),XHI(1,3),M,X,F,FX,WORK,IWORK) CALL TEST('J1',BESJ1,VJ1,XLO(1,4),XHI(1,4),M,X,F,FX,WORK,IWORK) CALL TEST('K0',BESK0,VK0,XLO(1,5),XHI(1,5),M,X,F,FX,WORK,IWORK) CALL TEST('K1',BESK1,VK1,XLO(1,6),XHI(1,6),M,X,F,FX,WORK,IWORK) CALL TEST('Y0',BESY0,VY0,XLO(1,7),XHI(1,7),M,X,F,FX,WORK,IWORK) CALL TEST('Y1',BESY1,VY1,XLO(1,8),XHI(1,8),M,X,F,FX,WORK,IWORK) C END SUBROUTINE TEST (NAME, SFUN, VSUB, XLOW, XHIGH, M, X, F, FX, + WORK, IWORK) C C---------------------------------------------------------------------- C C compare scalar function sfun to vectorized subroutine vsub for m C arguments distributed in various ways in the intervals C (xlo(i),xhi(i)), i=1,2,3. C C---------------------------------------------------------------------- C C ---------- C PARAMETERS C ---------- C REAL SFUN, XLOW, XHIGH, X, F, FX, WORK INTEGER IWORK, M CHARACTER*(*) NAME C DIMENSION XLOW(3), XHIGH(3), X(M), F(M), FX(M), WORK(6*M), + IWORK(M) C EXTERNAL SFUN, VSUB C C --------------- C LOCAL VARIABLES C --------------- C LOGICAL DEBUG INTEGER NREP, NTRY PARAMETER ( DEBUG=.FALSE., NREP=15, NTRY=19 ) C INTEGER I, INFO, IREP, IOFF, IPC(NTRY), JPC(NTRY), K, KPC, + MA, MB, MC, NFAIL, NREPV REAL DIFF, DMAX, DX, EPS, R1MACH REAL SECOND, SU, SUMAX, SUMIN, SUAVE, T0, T1, TS, TSAVE, TV, + TVAVE C SAVE IPC, JPC C DATA IPC / 0, 0, 0, 0, 0, 25, 25, 25, 25, 33, 50, 50, 50, + 75, 75,100, 1, 1, 98/ DATA JPC / 0, 25, 50, 75,100, 0, 25, 50, 75, 33, 0, 25, 50, + 0, 25, 0, 1, 98, 1/ C C---------------------------------------------------------------------- C C ntry = number of trials (each with different argument distribution) C nrep = number of times test repeated to get reliable timings C C the distribution of arguments in trial k is C C ipc(k)% in (xlow(1),xhigh(1)) C jpc(k)% in (xlow(2),xhigh(2)) C 100-ipc(k)-jpc(k)% in (xlow(3),xhigh(3)) C C scalar results (fx) and vector results (f) are compared and accepted C if abs((f - fx)/fx) < eps for fx > 1 (relative difference) C or abs(f - fx) < eps for fx < 1 (absolute difference) C C---------------------------------------------------------------------- C EPS = 5.0E0*R1MACH(4) C ... TRIGGER INITIALIZATIONS C X(1) = 1.0E0 CALL VSUB(1,X,F,WORK,IWORK,INFO) FX(1) = SFUN(X(1)) C PRINT * PRINT * PRINT *,' ----------' PRINT *,' TEST OF ',NAME PRINT *,' ----------' PRINT * PRINT *,' Vector Length = ',M PRINT *,' Repetitions = ',NREP PRINT *,' Eps = ',EPS PRINT * PRINT *,' Range A = (',XLOW(1),',',XHIGH(1),')' PRINT *,' Range B = (',XLOW(2),',',XHIGH(2),')' PRINT *,' Range C = (',XLOW(3),',',XHIGH(3),')' PRINT * PRINT *,' %A = percent arguments in range A' PRINT *,' %B = percent arguments in range B' PRINT *,' %C = percent arguments in range C' PRINT *,' Scalar = time for scalar code (FNLIB)' PRINT *,' Vector = time for vector code (VFNLIB)' PRINT *,' SU = speedup (Scalar/Vector)' PRINT *,' Nfail = number of scalar-vector differences .gt. eps' PRINT *,' Dmax = maximum scalar-vector difference' PRINT * PRINT *,' All times in seconds/argument.' PRINT * PRINT *,' Scalar-vector differences measured by' PRINT *,' abs(s-v)/max(abs(s),1)' PRINT *,' where s,v are scalar,vector results.' PRINT * PRINT *,' %A %B %C Scalar Vector SU Nfail Dmax' PRINT *,' -------------------------------------------------------' C SUMAX = 0.0 SUMIN = 1000.0 SUAVE = 0.0 TSAVE = 0.0 TVAVE = 0.0 C C ... PERFORM EACH TEST C DO 100 K=1,NTRY C MA = IPC(K)*M/100 MB = JPC(K)*M/100 KPC = 100 - IPC(K) - JPC(K) MC = M - MA - MB C C ... SELECT X VECTOR C IOFF = 0 C IF (MA .GT. 0) THEN DX = (XHIGH(1) - XLOW(1))/MA DO 10 I=1,MA X(I) = XLOW(1) + I*DX 10 CONTINUE IOFF = IOFF + MA ENDIF C IF (MB .GT. 0) THEN DX = (XHIGH(2) - XLOW(2))/MB DO 20 I=1,MB X(I+IOFF) = XLOW(2) + I*DX 20 CONTINUE IOFF = IOFF + MB ENDIF C IF (MC .GT. 0) THEN DX = (XHIGH(3) - XLOW(3))/MC DO 30 I=1,MC X(I+IOFF) = XLOW(3) + I*DX 30 CONTINUE ENDIF C C ... RUN VECTOR CODE C NREPV = 5*NREP T0 = SECOND() DO 40 IREP=1,NREPV CALL VSUB(M,X,F,WORK,IWORK,INFO) 40 CONTINUE T1 = SECOND() TV = (T1 - T0)/NREPV IF (INFO .LT. 0) PRINT *,'Warning in ',NAME,' : Info = ',INFO IF (INFO .GT. 0) PRINT *,'Failure in ',NAME,' : Info = ',INFO, + ' Position = ',IWORK(1) C C ... RUN SCALAR CODE C T0 = SECOND() DO 60 IREP=1,NREP DO 50 I=1,M FX(I) = SFUN(X(I)) 50 CONTINUE 60 CONTINUE T1 = SECOND() TS = (T1 - T0)/NREP C C ... COMPUTE TIMINGS C SU = TS/TV SUMIN = MIN(SU,SUMIN) SUMAX = MAX(SU,SUMAX) SUAVE = SUAVE + SU TSAVE = TSAVE + TS TVAVE = TVAVE + TV C C ... EVALUATE RESULTS C DMAX = 0.0E0 NFAIL = 0 DO 70 I=1,M DIFF = ABS(FX(I)-F(I))/MAX(ABS(FX(I)),1.0E0) DMAX = MAX(DIFF,DMAX) IF (DIFF .GT. EPS) THEN NFAIL = NFAIL + 1 IF (DEBUG) PRINT *,'Failure ',I,X(I),F(I),FX(I),DIFF ENDIF 70 CONTINUE C C ... PRINT RESULTS C PRINT 1000, IPC(K),JPC(K),KPC,TS/M,TV/M,SU,NFAIL,DMAX C 100 CONTINUE C SUAVE = SUAVE/NTRY TSAVE = TSAVE/M/NTRY TVAVE = TVAVE/M/NTRY C PRINT * PRINT 1002,'Ave Scalar = ',TSAVE,' sec/argument' PRINT 1002,'Ave Vector = ',TVAVE,' sec/argument' PRINT * PRINT 1001,'Max speed up = ',SUMAX PRINT 1001,'Min speed up = ',SUMIN PRINT * PRINT 1001,'Ave speed up = ',SUAVE C RETURN 1000 FORMAT(3(1X,I3),1P,2(3X,E7.1),0P,1X,F5.1,1X,I4,3X,1P,E7.1) 1001 FORMAT(1X,A,F5.1,A) 1002 FORMAT(1X,A,1P,E7.1,A) END INTEGER FUNCTION I1MACH(I) C***BEGIN PROLOGUE I1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 910131 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE Returns integer machine dependent constants C***DESCRIPTION C C This is the CMLIB version of I1MACH, the integer machine C constants subroutine originally developed for the PORT library. C C I1MACH can be used to obtain machine-dependent parameters C for the local machine environment. It is a function C subroutine with one (input) argument, and can be called C as follows, for example C C K = I1MACH(I) C C where I=1,...,16. The (output) value of K above is C determined by the (input) value of I. The results for C various values of I are discussed below. C C I/O unit numbers. C I1MACH( 1) = the standard input unit. C I1MACH( 2) = the standard output unit. C I1MACH( 3) = the standard punch unit. C I1MACH( 4) = the standard error message unit. C C Words. C I1MACH( 5) = the number of bits per integer storage unit. C I1MACH( 6) = the number of characters per integer storage unit. C C Integers. C assume integers are represented in the S-digit, base-A form C C sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) C C where 0 .LE. X(I) .LT. A for I=0,...,S-1. C I1MACH( 7) = A, the base. C I1MACH( 8) = S, the number of base-A digits. C I1MACH( 9) = A**S - 1, the largest magnitude. C C Floating-Point Numbers. C Assume floating-point numbers are represented in the T-digit, C base-B form C sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C where 0 .LE. X(I) .LT. B for I=1,...,T, C 0 .LT. X(1), and EMIN .LE. E .LE. EMAX. C I1MACH(10) = B, the base. C C Single-Precision C I1MACH(11) = T, the number of base-B digits. C I1MACH(12) = EMIN, the smallest exponent E. C I1MACH(13) = EMAX, the largest exponent E. C C Double-Precision C I1MACH(14) = T, the number of base-B digits. C I1MACH(15) = EMIN, the smallest exponent E. C I1MACH(16) = EMAX, the largest exponent E. C C To alter this function for a particular environment, C the desired set of DATA statements should be activated by C removing the C from column 1. Also, the values of C I1MACH(1) - I1MACH(4) should be checked for consistency C with the local operating system. C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***ROUTINES CALLED (NONE) C***END PROLOGUE I1MACH C INTEGER IMACH(16),OUTPUT EQUIVALENCE (IMACH(4),OUTPUT) C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300). C C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST C === MACHINE = SUN C === MACHINE = 68000 C === MACHINE = 8087 C === MACHINE = IBM.PC C === MACHINE = ATT.3B C === MACHINE = ATT.7300 C === MACHINE = ATT.6300 DATA IMACH( 1) / 5 / DATA IMACH( 2) / 6 / DATA IMACH( 3) / 7 / DATA IMACH( 4) / 6 / DATA IMACH( 5) / 32 / DATA IMACH( 6) / 4 / DATA IMACH( 7) / 2 / DATA IMACH( 8) / 31 / DATA IMACH( 9) / 2147483647 / DATA IMACH(10) / 2 / DATA IMACH(11) / 24 / DATA IMACH(12) / -125 / DATA IMACH(13) / 128 / DATA IMACH(14) / 53 / DATA IMACH(15) / -1021 / DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C === MACHINE = AMDAHL C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C === MACHINE = BURROUGHS.1700 C DATA IMACH( 1) / 7 / C DATA IMACH( 2) / 2 / C DATA IMACH( 3) / 2 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 33 / C DATA IMACH( 9) / Z1FFFFFFFF / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -256 / C DATA IMACH(13) / 255 / C DATA IMACH(14) / 60 / C DATA IMACH(15) / -256 / C DATA IMACH(16) / 255 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C === MACHINE = BURROUGHS.5700 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -50 / C DATA IMACH(16) / 76 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C === MACHINE = BURROUGHS.6700 C === MACHINE = BURROUGHS.7700 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -32754 / C DATA IMACH(16) / 32780 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (NATIVE MODE) C C === MACHINE = CONVEX.C1 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1023 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE CONVEX (NATIVE MODE) C WITH -R8 OPTION C C === MACHINE = CONVEX.C1.R8 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1023 / C DATA IMACH(13) / 1023 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1023 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (IEEE MODE) C C === MACHINE = CONVEX.C1.IEEE C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE CONVEX (IEEE MODE) C WITH -R8 OPTION C C === MACHINE = CONVEX.C1.IEEE.R8 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1021 / C DATA IMACH(13) / 1024 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5). C C === MACHINE = CYBER.170.NOS C === MACHINE = CYBER.180.NOS C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / O"00007777777777777777" / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 48 / C DATA IMACH(12) / -974 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 96 / C DATA IMACH(15) / -927 / C DATA IMACH(16) / 1070 / C C MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE C C === MACHINE = CYBER.180.NOS/VE C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 9223372036854775807 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -4095 / C DATA IMACH(13) / 4094 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -4095 / C DATA IMACH(16) / 4094 / C C MACHINE CONSTANTS FOR THE CYBER 205 C C === MACHINE = CYBER.205 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 47 / C DATA IMACH( 9) / X'00007FFFFFFFFFFF' / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -28625 / C DATA IMACH(13) / 28718 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -28625 / C DATA IMACH(16) / 28718 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C === MACHINE = CDC.6000 C === MACHINE = CDC.7000 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / 00007777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 48 / C DATA IMACH(12) / -974 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 96 / C DATA IMACH(15) / -927 / C DATA IMACH(16) / 1070 / C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C USING THE 46 BIT INTEGER COMPILER OPTION C C === MACHINE = CRAY.46-BIT-INTEGER C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 46 / C DATA IMACH( 9) / 777777777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -8189 / C DATA IMACH(13) / 8190 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -8099 / C DATA IMACH(16) / 8190 / C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C USING THE 64 BIT INTEGER COMPILER OPTION C C === MACHINE = CRAY.64-BIT-INTEGER C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 777777777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -8189 / C DATA IMACH(13) / 8190 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -8099 / C DATA IMACH(16) / 8190 /C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C === MACHINE = DATA_GENERAL.ECLIPSE.S/200 C DATA IMACH( 1) / 11 / C DATA IMACH( 2) / 12 / C DATA IMACH( 3) / 8 / C DATA IMACH( 4) / 10 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) /32767 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C ELXSI 6400 C C === MACHINE = ELSXI.6400 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 32 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1022 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE HARRIS 220 C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 C C === MACHINE = HARRIS.220 C === MACHINE = HARRIS.SLASH6 C === MACHINE = HARRIS.SLASH7 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / 3 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 23 / C DATA IMACH( 9) / 8388607 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 38 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C C === MACHINE = HONEYWELL.600/6000 C === MACHINE = HONEYWELL.DPS.8/70 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 43 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 63 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C === MACHINE = HP.2100.3_WORD_DP C DATA IMACH(1) / 5/ C DATA IMACH(2) / 6 / C DATA IMACH(3) / 4 / C DATA IMACH(4) / 1 / C DATA IMACH(5) / 16 / C DATA IMACH(6) / 2 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 15 / C DATA IMACH(9) / 32767 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 23 / C DATA IMACH(12)/ -128 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 39 / C DATA IMACH(15)/ -128 / C DATA IMACH(16)/ 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C === MACHINE = HP.2100.4_WORD_DP C DATA IMACH(1) / 5 / C DATA IMACH(2) / 6 / C DATA IMACH(3) / 4 / C DATA IMACH(4) / 1 / C DATA IMACH(5) / 16 / C DATA IMACH(6) / 2 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 15 / C DATA IMACH(9) / 32767 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 23 / C DATA IMACH(12)/ -128 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 55 / C DATA IMACH(15)/ -128 / C DATA IMACH(16)/ 127 / C C HP 9000 C C === MACHINE = HP.9000 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 7 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 32 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1015 / C DATA IMACH(16) / 1017 / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86 AND C THE INTERDATA 3230 AND INTERDATA 7/32. C C === MACHINE = IBM.360 C === MACHINE = IBM.370 C === MACHINE = XEROX.SIGMA.5 C === MACHINE = XEROX.SIGMA.7 C === MACHINE = XEROX.SIGMA.9 C === MACHINE = SEL.85 C === MACHINE = SEL.86 C === MACHINE = INTERDATA.3230 C === MACHINE = INTERDATA.7/32 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z7FFFFFFF / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C === MACHINE = INTERDATA.8/32.UNIX C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z'7FFFFFFF' / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 62 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 62 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C === MACHINE = PDP-10.KA C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 54 / C DATA IMACH(15) / -101 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C === MACHINE = PDP-10.KI C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 62 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGER ARITHMETIC. C C === MACHINE = PDP-11.32-BIT C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 16-BIT INTEGER ARITHMETIC. C C === MACHINE = PDP-11.16-BIT C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. C C === MACHINE = SEQUENT.BALANCE.8000 C DATA IMACH( 1) / 0 / C DATA IMACH( 2) / 0 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 1 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C === MACHINE = UNIVAC.1100 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 1 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 60 / C DATA IMACH(15) /-1024 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE VAX 11/780 C C === MACHINE = VAX.11/780 C DATA IMACH(1) / 5 / C DATA IMACH(2) / 6 / C DATA IMACH(3) / 5 / C DATA IMACH(4) / 6 / C DATA IMACH(5) / 32 / C DATA IMACH(6) / 4 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 31 / C DATA IMACH(9) /2147483647 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 24 / C DATA IMACH(12)/ -127 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 56 / C DATA IMACH(15)/ -127 / C DATA IMACH(16)/ 127 / C C C***FIRST EXECUTABLE STATEMENT I1MACH C I1MACH=IMACH(I) RETURN C END REAL FUNCTION R1MACH(I) C***BEGIN PROLOGUE R1MACH C***DATE WRITTEN 790101 (YYMMDD) C***REVISION DATE 910131 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE Returns single precision machine dependent constants C***DESCRIPTION C C This is the CMLIB version of R1MACH, the real machine C constants subroutine originally developed for the PORT library. C C R1MACH can be used to obtain machine-dependent parameters C for the local machine environment. It is a function C subroutine with one (input) argument, and can be called C as follows, for example C C A = R1MACH(I) C C where I=1,...,5. The (output) value of A above is C determined by the (input) value of I. The results for C various values of I are discussed below. C C Single-Precision Machine Constants C R1MACH(1) = B**(EMIN-1), the smallest positive magnitude. C R1MACH(2) = B**EMAX*(1 - B**(-T)), the largest magnitude. C R1MACH(3) = B**(-T), the smallest relative spacing. C R1MACH(4) = B**(1-T), the largest relative spacing. C R1MACH(5) = LOG10(B) C***REFERENCES FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR C A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE- C MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978, C PP. 177-188. C***ROUTINES CALLED (NONE) C***END PROLOGUE R1MACH C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) C REAL RMACH(5) C EQUIVALENCE (RMACH(1),SMALL(1)) EQUIVALENCE (RMACH(2),LARGE(1)) EQUIVALENCE (RMACH(3),RIGHT(1)) EQUIVALENCE (RMACH(4),DIVER(1)) EQUIVALENCE (RMACH(5),LOG10(1)) C C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300). C C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST C === MACHINE = SUN C === MACHINE = 68000 C === MACHINE = 8087 C === MACHINE = IBM.PC C === MACHINE = ATT.3B C === MACHINE = ATT.6300 C === MACHINE = ATT.7300 DATA SMALL(1) / 8388608 / DATA LARGE(1) / 2139095039 / DATA RIGHT(1) / 864026624 / DATA DIVER(1) / 872415232 / DATA LOG10(1) / 1050288283 / C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C === MACHINE = AMDAHL C DATA SMALL(1) / 1048576 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 990904320 / C DATA DIVER(1) / 1007681536 / C DATA LOG10(1) / 1091781651 / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C === MACHINE = BURROUGHS.1700 C DATA RMACH(1) / Z400800000 / C DATA RMACH(2) / Z5FFFFFFFF / C DATA RMACH(3) / Z4E9800000 / C DATA RMACH(4) / Z4EA800000 / C DATA RMACH(5) / Z500E730E8 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS. C C === MACHINE = BURROUGHS.5700 C === MACHINE = BURROUGHS.6700 C === MACHINE = BURROUGHS.7700 C DATA RMACH(1) / O1771000000000000 / C DATA RMACH(2) / O0777777777777777 / C DATA RMACH(3) / O1311000000000000 / C DATA RMACH(4) / O1301000000000000 / C DATA RMACH(5) / O1157163034761675 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (NATIVE MODE) C C === MACHINE = CONVEX.C1 C DATA RMACH(1) / 2.9387360E-39 / C DATA RMACH(2) / 1.7014117E+38 / C DATA RMACH(3) / 5.9604645E-08 / C DATA RMACH(4) / 1.1920929E-07 / C DATA RMACH(5) / 3.0102999E-01 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (NATIVE MODE) C WITH -R8 OPTION C C === MACHINE = CONVEX.C1.R8 C DATA RMACH(1) / 5.562684646268007D-309 / C DATA RMACH(2) / 8.988465674311577D+307 / C DATA RMACH(3) / 1.110223024625157D-016 / C DATA RMACH(4) / 2.220446049250313D-016 / C DATA RMACH(5) / 3.010299956639812D-001 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (IEEE MODE) C C === MACHINE = CONVEX.C1.IEEE C DATA RMACH(1) / 1.1754945E-38 / C DATA RMACH(2) / 3.4028234E+38 / C DATA RMACH(3) / 5.9604645E-08 / C DATA RMACH(4) / 1.1920929E-07 / C DATA RMACH(5) / 3.0102999E-01 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (IEEE MODE) C WITH -R8 OPTION C C === MACHINE = CONVEX.C1.IEEE.R8 C DATA RMACH(1) / 2.225073858507202D-308 / C DATA RMACH(2) / 1.797693134862315D+308 / C DATA RMACH(3) / 1.110223024625157D-016 / C DATA RMACH(4) / 2.220446049250313D-016 / C DATA RMACH(5) / 3.010299956639812D-001 / C C MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5). C C === MACHINE = CYBER.170.NOS C === MACHINE = CYBER.180.NOS C DATA RMACH(1) / O"00014000000000000000" / C DATA RMACH(2) / O"37767777777777777777" / C DATA RMACH(3) / O"16404000000000000000" / C DATA RMACH(4) / O"16414000000000000000" / C DATA RMACH(5) / O"17164642023241175720" / C C MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE C C === MACHINE = CYBER.180.NOS/VE C DATA RMACH(1) / Z"3001800000000000" / C DATA RMACH(2) / Z"4FFEFFFFFFFFFFFE" / C DATA RMACH(3) / Z"3FD2800000000000" / C DATA RMACH(4) / Z"3FD3800000000000" / C DATA RMACH(5) / Z"3FFF9A209A84FBCF" / C C MACHINE CONSTANTS FOR THE CYBER 205 C C === MACHINE = CYBER.205 C DATA RMACH(1) / X'9000400000000000' / C DATA RMACH(2) / X'6FFF7FFFFFFFFFFF' / C DATA RMACH(3) / X'FFA3400000000000' / C DATA RMACH(4) / X'FFA4400000000000' / C DATA RMACH(5) / X'FFD04D104D427DE8' / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C === MACHINE = CDC.6000 C === MACHINE = CDC.7000 C DATA RMACH(1) / 00014000000000000000B / C DATA RMACH(2) / 37767777777777777777B / C DATA RMACH(3) / 16404000000000000000B / C DATA RMACH(4) / 16414000000000000000B / C DATA RMACH(5) / 17164642023241175720B / C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C === MACHINE = CRAY.46-BIT-INTEGER C === MACHINE = CRAY.64-BIT-INTEGER C DATA RMACH(1) / 200034000000000000000B / C DATA RMACH(2) / 577767777777777777776B / C DATA RMACH(3) / 377224000000000000000B / C DATA RMACH(4) / 377234000000000000000B / C DATA RMACH(5) / 377774642023241175720B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - C STATIC RMACH(5) C C === MACHINE = DATA_GENERAL.ECLIPSE.S/200 C DATA SMALL/20K,0/,LARGE/77777K,177777K/ C DATA RIGHT/35420K,0/,DIVER/36020K,0/ C DATA LOG10/40423K,42023K/ C C ELXSI 6400 C C === MACHINE = ELSXI.6400 C DATA SMALL(1) / '00800000'X / C DATA LARGE(1) / '7F7FFFFF'X / C DATA RIGHT(1) / '33800000'X / C DATA DIVER(1) / '34000000'X / C DATA LOG10(1) / '3E9A209B'X / C C MACHINE CONSTANTS FOR THE HARRIS 220 C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 C C === MACHINE = HARRIS.220 C === MACHINE = HARRIS.SLASH6 C === MACHINE = HARRIS.SLASH7 C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '00000177 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000352 / C DATA DIVER(1),DIVER(2) / '20000000, '00000353 / C DATA LOG10(1),LOG10(2) / '23210115, '00000377 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C C === MACHINE = HONEYWELL.600/6000 C === MACHINE = HONEYWELL.DPS.8/70 C DATA RMACH(1) / O402400000000 / C DATA RMACH(2) / O376777777777 / C DATA RMACH(3) / O714400000000 / C DATA RMACH(4) / O716400000000 / C DATA RMACH(5) / O776464202324 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION WITH FTN4 C C === MACHINE = HP.2100.3_WORD_DP C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION WITH FTN4 C C === MACHINE = HP.2100.4_WORD_DP C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE91), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C HP 9000 C C R1MACH(1) = 1.17549435E-38 C R1MACH(2) = 1.70141163E+38 C R1MACH(3) = 5.960464478E-8 C R1MACH(4) = 1.119209290E-7 C R1MACH(5) = 3.01030010E-1 C C === MACHINE = HP.9000 C DATA SMALL(1) / 00040000000B / C DATA LARGE(1) / 17677777777B / C DATA RIGHT(1) / 06340000000B / C DATA DIVER(1) / 06400000000B / C DATA LOG10(1) / 07646420233B / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86 AND C THE INTERDATA 3230 AND INTERDATA 7/32. C C === MACHINE = IBM.360 C === MACHINE = IBM.370 C === MACHINE = XEROX.SIGMA.5 C === MACHINE = XEROX.SIGMA.7 C === MACHINE = XEROX.SIGMA.9 C === MACHINE = SEL.85 C === MACHINE = SEL.86 C === MACHINE = INTERDATA.3230 C === MACHINE = INTERDATA.7/32 C DATA RMACH(1) / Z00100000 / C DATA RMACH(2) / Z7FFFFFFF / C DATA RMACH(3) / Z3B100000 / C DATA RMACH(4) / Z3C100000 / C DATA RMACH(5) / Z41134413 / C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C === MACHINE = INTERDATA.8/32.UNIX C DATA RMACH(1) / Z'00100000' / C DATA RMACH(2) / Z'7EFFFFFF' / C DATA RMACH(3) / Z'3B100000' / C DATA RMACH(4) / Z'3C100000' / C DATA RMACH(5) / Z'41134413' / C C MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR). C C === MACHINE = PDP-10.KA C === MACHINE = PDP-10.KI C DATA RMACH(1) / "000400000000 / C DATA RMACH(2) / "377777777777 / C DATA RMACH(3) / "146400000000 / C DATA RMACH(4) / "147400000000 / C DATA RMACH(5) / "177464202324 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C === MACHINE = PDP-11.32-BIT C DATA SMALL(1) / 8388608 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 880803840 / C DATA DIVER(1) / 889192448 / C DATA LOG10(1) / 1067065499 / C C DATA RMACH(1) / O00040000000 / C DATA RMACH(2) / O17777777777 / C DATA RMACH(3) / O06440000000 / C DATA RMACH(4) / O06500000000 / C DATA RMACH(5) / O07746420233 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C === MACHINE = PDP-11.16-BIT C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA RIGHT(1),RIGHT(2) / 13440, 0 / C DATA DIVER(1),DIVER(2) / 13568, 0 / C DATA LOG10(1),LOG10(2) / 16282, 8347 / C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA RIGHT(1),RIGHT(2) / O032200, O000000 / C DATA DIVER(1),DIVER(2) / O032400, O000000 / C DATA LOG10(1),LOG10(2) / O037632, O020233 / C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. C C === MACHINE = SEQUENT.BALANCE.8000 C DATA SMALL(1) / $00800000 / C DATA LARGE(1) / $7F7FFFFF / C DATA RIGHT(1) / $33800000 / C DATA DIVER(1) / $34000000 / C DATA LOG10(1) / $3E9A209B / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C === MACHINE = UNIVAC.1100 C DATA RMACH(1) / O000400000000 / C DATA RMACH(2) / O377777777777 / C DATA RMACH(3) / O146400000000 / C DATA RMACH(4) / O147400000000 / C DATA RMACH(5) / O177464202324 / C C MACHINE CONSTANTS FOR THE VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C === MACHINE = VAX.11/780 C DATA SMALL(1) / 128 / C DATA LARGE(1) / -32769 / C DATA RIGHT(1) / 13440 / C DATA DIVER(1) / 13568 / C DATA LOG10(1) / 547045274 / C C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS*** C C DATA SMALL(1) / Z00000080 / C DATA LARGE(1) / ZFFFF7FFF / C DATA RIGHT(1) / Z00003480 / C DATA DIVER(1) / Z00003500 / C DATA LOG10(1) / Z209B3F9A / C C C***FIRST EXECUTABLE STATEMENT R1MACH C R1MACH = RMACH(I) RETURN C END DOUBLE PRECISION FUNCTION D1MACH(I) C***BEGIN PROLOGUE D1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 910131 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE Returns double precision machine dependent constants C***DESCRIPTION C C This is the CMLIB version of D1MACH, the double precision machine C constants subroutine originally developed for the PORT library. C C D1MACH can be used to obtain machine-dependent parameters C for the local machine environment. It is a function C subprogram with one (input) argument, and can be called C as follows, for example C C D = D1MACH(I) C C where I=1,...,5. The (output) value of D above is C determined by the (input) value of I. The results for C various values of I are discussed below. C C Double-precision machine constants C D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude. C D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude. C D1MACH( 3) = B**(-T), the smallest relative spacing. C D1MACH( 4) = B**(1-T), the largest relative spacing. C D1MACH( 5) = LOG10(B) C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***ROUTINES CALLED (NONE) C***END PROLOGUE D1MACH C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) C DOUBLE PRECISION DMACH(5) C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES AND MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST. C C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST C === MACHINE = SUN C === MACHINE = 68000 C === MACHINE = ATT.3B C === MACHINE = ATT.7300 DATA SMALL(1),SMALL(2) / 1048576, 0 / DATA LARGE(1),LARGE(2) / 2146435071, -1 / DATA RIGHT(1),RIGHT(2) / 1017118720, 0 / DATA DIVER(1),DIVER(2) / 1018167296, 0 / DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 / C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES AND 8087-BASED C MICROS, SUCH AS THE IBM PC AND AT&T 6300, IN WHICH THE LEAST C SIGNIFICANT BYTE IS STORED FIRST. C C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST C === MACHINE = 8087 C === MACHINE = IBM.PC C === MACHINE = ATT.6300 C DATA SMALL(1),SMALL(2) / 0, 1048576 / C DATA LARGE(1),LARGE(2) / -1, 2146435071 / C DATA RIGHT(1),RIGHT(2) / 0, 1017118720 / C DATA DIVER(1),DIVER(2) / 0, 1018167296 / C DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 / C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C === MACHINE = AMDAHL C DATA SMALL(1),SMALL(2) / 1048576, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 856686592, 0 / C DATA DIVER(1),DIVER(2) / 873463808, 0 / C DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C === MACHINE = BURROUGHS.1700 C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C === MACHINE = BURROUGHS.5700 C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C === MACHINE = BURROUGHS.6700 C === MACHINE = BURROUGHS.7700 C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (NATIVE MODE) C WITH OR WITHOUT -R8 OPTION C C === MACHINE = CONVEX.C1 C === MACHINE = CONVEX.C1.R8 C DATA DMACH(1) / 5.562684646268007D-309 / C DATA DMACH(2) / 8.988465674311577D+307 / C DATA DMACH(3) / 1.110223024625157D-016 / C DATA DMACH(4) / 2.220446049250313D-016 / C DATA DMACH(5) / 3.010299956639812D-001 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (IEEE MODE) C WITH OR WITHOUT -R8 OPTION C C === MACHINE = CONVEX.C1.IEEE C === MACHINE = CONVEX.C1.IEEE.R8 C DATA DMACH(1) / 2.225073858507202D-308 / C DATA DMACH(2) / 1.797693134862315D+308 / C DATA DMACH(3) / 1.110223024625157D-016 / C DATA DMACH(4) / 2.220446049250313D-016 / C DATA DMACH(5) / 3.010299956639812D-001 / C C MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5). C C === MACHINE = CYBER.170.NOS C === MACHINE = CYBER.180.NOS C DATA SMALL(1) / O"00604000000000000000" / C DATA SMALL(2) / O"00000000000000000000" / C DATA LARGE(1) / O"37767777777777777777" / C DATA LARGE(2) / O"37167777777777777777" / C DATA RIGHT(1) / O"15604000000000000000" / C DATA RIGHT(2) / O"15000000000000000000" / C DATA DIVER(1) / O"15614000000000000000" / C DATA DIVER(2) / O"15010000000000000000" / C DATA LOG10(1) / O"17164642023241175717" / C DATA LOG10(2) / O"16367571421742254654" / C C MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE C C === MACHINE = CYBER.180.NOS/VE C DATA SMALL(1) / Z"3001800000000000" / C DATA SMALL(2) / Z"3001000000000000" / C DATA LARGE(1) / Z"4FFEFFFFFFFFFFFE" / C DATA LARGE(2) / Z"4FFE000000000000" / C DATA RIGHT(1) / Z"3FD2800000000000" / C DATA RIGHT(2) / Z"3FD2000000000000" / C DATA DIVER(1) / Z"3FD3800000000000" / C DATA DIVER(2) / Z"3FD3000000000000" / C DATA LOG10(1) / Z"3FFF9A209A84FBCF" / C DATA LOG10(2) / Z"3FFFF7988F8959AC" / C C MACHINE CONSTANTS FOR THE CYBER 205 C C === MACHINE = CYBER.205 C DATA SMALL(1) / X'9000400000000000' / C DATA SMALL(2) / X'8FD1000000000000' / C DATA LARGE(1) / X'6FFF7FFFFFFFFFFF' / C DATA LARGE(2) / X'6FD07FFFFFFFFFFF' / C DATA RIGHT(1) / X'FF74400000000000' / C DATA RIGHT(2) / X'FF45000000000000' / C DATA DIVER(1) / X'FF75400000000000' / C DATA DIVER(2) / X'FF46000000000000' / C DATA LOG10(1) / X'FFD04D104D427DE7' / C DATA LOG10(2) / X'FFA17DE623E2566A' / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C === MACHINE = CDC.6000 C === MACHINE = CDC.7000 C DATA SMALL(1) / 00604000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C DATA LARGE(1) / 37767777777777777777B / C DATA LARGE(2) / 37167777777777777777B / C DATA RIGHT(1) / 15604000000000000000B / C DATA RIGHT(2) / 15000000000000000000B / C DATA DIVER(1) / 15614000000000000000B / C DATA DIVER(2) / 15010000000000000000B / C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B / C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C === MACHINE = CRAY.46-BIT-INTEGER C === MACHINE = CRAY.64-BIT-INTEGER C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777776B / C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - C STATIC DMACH(5) C C === MACHINE = DATA_GENERAL.ECLIPSE.S/200 C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ C DATA LOG10/40423K,42023K,50237K,74776K/ C C ELXSI 6400 C C === MACHINE = ELSXI.6400 C DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X / C DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X / C DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X / C DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X / C DATA LOG10(1), DIVER(2) / '3FD34413'X,'509F79FF'X / C C MACHINE CONSTANTS FOR THE HARRIS 220 C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 C C === MACHINE = HARRIS.220 C === MACHINE = HARRIS.SLASH6 C === MACHINE = HARRIS.SLASH7 C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1),DIVER(2) / '20000000, '00000334 / C DATA LOG10(1),LOG10(2) / '23210115, '10237777 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C C === MACHINE = HONEYWELL.600/6000 C === MACHINE = HONEYWELL.DPS.8/70 C DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C === MACHINE = HP.2100.3_WORD_DP C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B / C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B / C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C === MACHINE = HP.2100.4_WORD_DP C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177777B / C DATA LARGE(3), LARGE(4) / 177777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 225B / C DATA DIVER(1), DIVER(2) / 40000B, 0 / C DATA DIVER(3), DIVER(4) / 0, 227B / C DATA LOG10(1), LOG10(2) / 46420B, 46502B / C DATA LOG10(3), LOG10(4) / 76747B, 176377B / C C HP 9000 C C D1MACH(1) = 2.8480954D-306 C D1MACH(2) = 1.40444776D+306 C D1MACH(3) = 2.22044605D-16 C D1MACH(4) = 4.44089210D-16 C D1MACH(5) = 3.01029996D-1 C C === MACHINE = HP.9000 C DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B / C DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B / C DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B / C DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B / C DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE INTERDATA 3230 AND INTERDATA 7/32. C C === MACHINE = IBM.360 C === MACHINE = IBM.370 C === MACHINE = XEROX.SIGMA.5 C === MACHINE = XEROX.SIGMA.7 C === MACHINE = XEROX.SIGMA.9 C === MACHINE = SEL.85 C === MACHINE = SEL.86 C === MACHINE = INTERDATA.3230 C === MACHINE = INTERDATA.7/32 C DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C === MACHINE = INTERDATA.8/32.UNIX C DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000' / C DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000' / C DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C === MACHINE = PDP-10.KA C DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C === MACHINE = PDP-10.KI C DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "047674776746 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C === MACHINE = PDP-11.32-BIT C DATA SMALL(1),SMALL(2) / 8388608, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 612368384, 0 / C DATA DIVER(1),DIVER(2) / 620756992, 0 / C DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 / C C DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C === MACHINE = PDP-11.16-BIT C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA SMALL(3),SMALL(4) / 0, 0 / C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA LARGE(3),LARGE(4) / -1, -1 / C DATA RIGHT(1),RIGHT(2) / 9344, 0 / C DATA RIGHT(3),RIGHT(4) / 0, 0 / C DATA DIVER(1),DIVER(2) / 9472, 0 / C DATA DIVER(3),DIVER(4) / 0, 0 / C DATA LOG10(1),LOG10(2) / 16282, 8346 / C DATA LOG10(3),LOG10(4) / -31493, -12296 / C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA SMALL(3),SMALL(4) / O000000, O000000 / C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA LARGE(3),LARGE(4) / O177777, O177777 / C DATA RIGHT(1),RIGHT(2) / O022200, O000000 / C DATA RIGHT(3),RIGHT(4) / O000000, O000000 / C DATA DIVER(1),DIVER(2) / O022400, O000000 / C DATA DIVER(3),DIVER(4) / O000000, O000000 / C DATA LOG10(1),LOG10(2) / O037632, O020232 / C DATA LOG10(3),LOG10(4) / O102373, O147770 / C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000 C C === MACHINE = SEQUENT.BALANCE.8000 C DATA SMALL(1),SMALL(2) / $00000000, $00100000 / C DATA LARGE(1),LARGE(2) / $FFFFFFFF, $7FEFFFFF / C DATA RIGHT(1),RIGHT(2) / $00000000, $3CA00000 / C DATA DIVER(1),DIVER(2) / $00000000, $3CB00000 / C DATA LOG10(1),LOG10(2) / $509F79FF, $3FD34413 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C === MACHINE = UNIVAC.1100 C DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 / C C MACHINE CONSTANTS FOR VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C === MACHINE = VAX.11/780 C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA LOG10(1), LOG10(2) / 546979738, -805796613 / C C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB / C C MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING) C (EXPRESSED IN INTEGER AND HEXADECIMAL) C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1), SMALL(2) / 16, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 15552, 0 / C DATA DIVER(1), DIVER(2) / 15568, 0 / C DATA LOG10(1), LOG10(2) / 1142112243, 2046775455 / C C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 / C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / C C C***FIRST EXECUTABLE STATEMENT D1MACH C D1MACH = DMACH(I) RETURN C END SUBROUTINE WFERR (SNAME, MESSAG, LEVEL) C***BEGIN PROLOGUE WFERR C***PURPOSE Processes an error message for the vectorized special C function package C***LIBRARY VFNLIB C***TYPE CHARACTER C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C WFERR processes an error message for the vectorized special C function package. C C P A R A M E T E R S C C SNAME (Input) Character*(*) C The name of the routine issuing the message. C C MESSAG (Input) Character*(*) C The error message. Must be less than 80 characters. C C LEVEL (Input) Integer C The message level. 1 for recoverable, 2 for fatal. C C This routine issues a STOP for fatal errors. C C C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WFERR C C ... PARAMETERS C CHARACTER*(*) MESSAG, SNAME INTEGER LEVEL C C ... LOCAL VARIABLES C INTEGER IUNIT, I1MACH, MLEN C C***FIRST EXECUTABLE STATEMENT WFERR C IUNIT = I1MACH(4) MLEN = MIN(LEN(MESSAG),80) C WRITE(IUNIT,*) IF (LEVEL .EQ. 2) THEN WRITE(IUNIT,*) 'FATAL ERROR IN ',SNAME,' ...' WRITE(IUNIT,*) MESSAG(1:MLEN) STOP ELSE WRITE(IUNIT,*) 'RECOVERABLE ERROR IN ',SNAME,' ...' WRITE(IUNIT,*) MESSAG(1:MLEN) RETURN ENDIF C END REAL FUNCTION SECOND () C C OBTAIN ELAPSED CPU TIME SINCE START OF RUN. C C ON CRAY -- THIS ROUTINE NOT NEEDED. C C ON CONVEX -- MORE ACCURATE TIMINGS ARE OBTAINED BY CALLING C THE VECLIB ROUTINE CPUTIME AS FOLLOWS C C SECOND = CPUTIME(0.0E0) C C ON UNIX -- (SUN, FOR EXAMPLE) THE FOLLOWING MAY BE USED C REAL TARRAY(2), ETIME, TOTAL TOTAL = ETIME(TARRAY) SECOND = TARRAY(1) + TARRAY(2) C RETURN END SUBROUTINE VI0 (M, X, F, WORK, IWORK, INFO) C***BEGIN PROLOGUE VI0 C***PURPOSE Computes the hyperbolic Bessel function of the first kind C of order zero (I0) for a vector of real arguments C***LIBRARY VFNLIB C***CATEGORY C10B1 C***TYPE SINGLE PRECISION (VI0-S, DVI0-D) C***KEYWORDS BESSEL FUNCTION,FIRST KIND,HYPERBOLIC BESSEL FUNCTION, C MODIFIED BESSEL FUNCTION, ORDER ZERO, VECTORIZED C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C VI0 computes the modified (hyperbolic) Bessel function of the C first kind of order zero (I0) for real arguments using uniform C approximation by Chebyshev polynomials. C C C P A R A M E T E R S C C M (Input) Integer (M .GT. 0) C The number of arguments at which the function is to be C evaluated. C C X (Input) Real array of length M C The arguments at which the function is to be evaluated are C stored in X(1) to X(M) in any order. C C F (Output) Real array of length M C F(i) contains the value of the function at X(i), i=1,..,M. C C WORK (Work) Real vector of length 7*M C C IWORK (Work) Integer vector of length M C C INFO (Output) Integer C Indicates status of computed result. The following table C lists possible values and their meanings. If OK=Yes then C all F(i) have been set, otherwise none have been set. C C INFO OK Description C ------------------------------------------------------------ C 0 Yes Successfull execution. C 1 No Error: M .LE. 0 C 2 No Error: Some abs(X(i)) so big I0 overflows. C The index of the first offending argument is C returned in IWORK(1). C C C ********************************************************************* C This routine is a modification of the function BESI0 developed by C W. Fullerton of LANL. C ********************************************************************* C C***REFERENCES (NONE) C***ROUTINES CALLED WI0 C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE VI0 C C ---------- C PARAMETERS C ---------- C INTEGER INFO, IWORK, M REAL F, X, WORK C DIMENSION X(M), F(M), WORK(7*M), IWORK(M) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER IWB0, IWB1, IWB2, IWY, IWTC, IWYC, IWZC, JIN C C***FIRST EXECUTABLE STATEMENT VI0 C C ... PARTITION WORK ARRAYS C IWY = 1 IWTC = IWY + M IWYC = IWTC + M IWZC = IWYC + M IWB0 = IWZC + M IWB1 = IWB0 + M IWB2 = IWB1 + M C Total = IWB2 + M C JIN = 1 C Total = JIN + M C C ... WI0 DOES ALL THE WORK C CALL WI0(M,X,F,WORK(IWY),WORK(IWTC),WORK(IWYC),WORK(IWZC), + IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO) C RETURN END SUBROUTINE WI0 (M, X, F, Y, TCMP, YCMP, ZCMP, INDX, B0, B1, B2, + INFO) C***BEGIN PROLOGUE WI0 C***SUBSIDIARY C***PURPOSE Computes the hyperbolic Bessel function of the first kind C of order zero (I0) for a vector of arguments C***LIBRARY VFNLIB C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***ROUTINES CALLED IWCS, R1MACH, WNGT, WGTHR, WGTLE, WGT, WLE, C WSCTR, WCS C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WI0 C C ---------- C PARAMETERS C ---------- C INTEGER INFO, INDX, M REAL B0, B1, B2, F, X, Y, TCMP, YCMP, ZCMP C DIMENSION B0(M), B1(M), B2(M), F(M), INDX(M), X(M), Y(M), + TCMP(M), YCMP(M), ZCMP(M) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER LAI0, LAI02, LBI0 PARAMETER ( LAI0=21, LAI02=22, LBI0=12 ) C INTEGER I, IWCS, J, KEY, N, NTI0, NTAI0, NTAI02 REAL AI0CS, AI02CS, BI0CS, EPMACH, EPS, R1MACH, XSML, + XMAX C DIMENSION AI0CS(LAI0), AI02CS(LAI02), BI0CS(LBI0) C SAVE AI0CS, AI02CS, BI0CS, N, NTAI0, NTAI02, NTI0, XSML, XMAX C C---------------------------------------------------------------------- C C Series for BI0 on the interval 0. to 9.00000D+00 C with weighted error 2.46E-18 C log weighted error 17.61 C significant figures required 17.90 C decimal places required 18.15 C DATA BI0 CS( 1) / -.0766054725 2839144951E0 / DATA BI0 CS( 2) / 1.9273379539 93808270E0 / DATA BI0 CS( 3) / .2282644586 920301339E0 / DATA BI0 CS( 4) / .0130489146 6707290428E0 / DATA BI0 CS( 5) / .0004344270 9008164874E0 / DATA BI0 CS( 6) / .0000094226 5768600193E0 / DATA BI0 CS( 7) / .0000001434 0062895106E0 / DATA BI0 CS( 8) / .0000000016 1384906966E0 / DATA BI0 CS( 9) / .0000000000 1396650044E0 / DATA BI0 CS(10) / .0000000000 0009579451E0 / DATA BI0 CS(11) / .0000000000 0000053339E0 / DATA BI0 CS(12) / .0000000000 0000000245E0 / C C---------------------------------------------------------------------- C C Series for AI0 on the interval 1.25000D-01 to 3.33333D-01 C with weighted error 7.87E-17 C log weighted error 16.10 C significant figures required 14.69 C decimal places required 16.76 C DATA AI0 CS( 1) / .0757599449 4023796E0 / DATA AI0 CS( 2) / .0075913808 1082334E0 / DATA AI0 CS( 3) / .0004153131 3389237E0 / DATA AI0 CS( 4) / .0000107007 6463439E0 / DATA AI0 CS( 5) / -.0000079011 7997921E0 / DATA AI0 CS( 6) / -.0000007826 1435014E0 / DATA AI0 CS( 7) / .0000002783 8499429E0 / DATA AI0 CS( 8) / .0000000082 5247260E0 / DATA AI0 CS( 9) / -.0000000120 4463945E0 / DATA AI0 CS(10) / .0000000015 5964859E0 / DATA AI0 CS(11) / .0000000002 2925563E0 / DATA AI0 CS(12) / -.0000000001 1916228E0 / DATA AI0 CS(13) / .0000000000 1757854E0 / DATA AI0 CS(14) / .0000000000 0112822E0 / DATA AI0 CS(15) / -.0000000000 0114684E0 / DATA AI0 CS(16) / .0000000000 0027155E0 / DATA AI0 CS(17) / -.0000000000 0002415E0 / DATA AI0 CS(18) / -.0000000000 0000608E0 / DATA AI0 CS(19) / .0000000000 0000314E0 / DATA AI0 CS(20) / -.0000000000 0000071E0 / DATA AI0 CS(21) / .0000000000 0000007E0 / C C---------------------------------------------------------------------- C C Series for AI02 on the interval 0. to 1.25000D-01 C with weighted error 3.79E-17 C log weighted error 16.42 C significant figures required 14.86 C decimal places required 17.09 C DATA AI02CS( 1) / .0544904110 1410882E0 / DATA AI02CS( 2) / .0033691164 7825569E0 / DATA AI02CS( 3) / .0000688975 8346918E0 / DATA AI02CS( 4) / .0000028913 7052082E0 / DATA AI02CS( 5) / .0000002048 9185893E0 / DATA AI02CS( 6) / .0000000226 6668991E0 / DATA AI02CS( 7) / .0000000033 9623203E0 / DATA AI02CS( 8) / .0000000004 9406022E0 / DATA AI02CS( 9) / .0000000000 1188914E0 / DATA AI02CS(10) / -.0000000000 3149915E0 / DATA AI02CS(11) / -.0000000000 1321580E0 / DATA AI02CS(12) / -.0000000000 0179419E0 / DATA AI02CS(13) / .0000000000 0071801E0 / DATA AI02CS(14) / .0000000000 0038529E0 / DATA AI02CS(15) / .0000000000 0001539E0 / DATA AI02CS(16) / -.0000000000 0004151E0 / DATA AI02CS(17) / -.0000000000 0000954E0 / DATA AI02CS(18) / .0000000000 0000382E0 / DATA AI02CS(19) / .0000000000 0000176E0 / DATA AI02CS(20) / -.0000000000 0000034E0 / DATA AI02CS(21) / -.0000000000 0000027E0 / DATA AI02CS(22) / .0000000000 0000003E0 / C C---------------------------------------------------------------------- C DATA NTI0 / 0 / C C***FIRST EXECUTABLE STATEMENT WI0 C IF (M .LE. 0) GO TO 910 C IF (NTI0.EQ.0) THEN EPMACH = R1MACH(3) EPS = 0.10E0*EPMACH NTI0 = IWCS(BI0CS, LBI0, EPS) NTAI0 = IWCS(AI0CS, LAI0, EPS) NTAI02 = IWCS(AI02CS, LAI02, EPS) XSML = SQRT(4.0E0*EPMACH) XMAX = LOG(R1MACH(2)) ENDIF C DO 10 I=1,M Y(I) = ABS(X(I)) 10 CONTINUE C CALL WNGT(M,Y,XMAX,KEY) IF (KEY .NE. 0) GO TO 920 C C ---------------- C CASE Y