tape from Fullerton at IMSL sent 21 Feb 84, read onto bowell 24 Feb by ehg
tr A-Z a-z <tapefile1 | fsplit
21 Mar '86
duplicate declarations removed from dbetai, dgamma, dbinom, dshi by wmc
From dessert.cam.nist.gov!boisvert Wed Mar 27 15:59:16 EST 1991
Date: Wed, 27 Mar 91 15:59:16 EST
From: boisvert@dessert.cam.nist.gov (Ronald_Boisvert_x3812)
To: ehg@research.att.com, na.brophy@na-net.ornl.gov, slatec@lanl.gov
Subject: FNLIB Bug Alert
---------
Bug Alert
---------
A serious bug is present in the routine DBESJ1 of FNLIB. This routine
computes the J1 Bessel function in double precision. The function values
produced by this software have the wrong sign for all arguments less
than -4.
The error is present in the version of FNLIB currently included in :
* netlib
* SLATEC
* CMLIB
* IMSL SFUN/LIBRARY (routine name : DBSJ1)
The following simple program illustrates the error. In it we tabulate
J1(x) for x in the interval (-10,0) as computed by BESJ1 and DBESJ1, the
single and double precision FNLIB routines, respectively.
real x, besj1
double precision xd, dbesj1
print *,' x besj1(x) dbesj1(x)'
print *,' ------------------------------------------------'
do 10 j=1,50
xd = -j*0.20d0
x = xd
print '(1p,2x,e9.2,2x,e12.5,2x,d20.13)',
+ x,besj1(x),dbesj1(xd)
10 continue
end
This program produces the following output on a Sun SPARCstation 1+ :
x besj1(x) dbesj1(x)
------------------------------------------------
-2.00E-01 -9.95008E-02 -9.9500832639236D-02
-4.00E-01 -1.96027E-01 -1.9602657795532D-01
-6.00E-01 -2.86701E-01 -2.8670098806392D-01
-8.00E-01 -3.68842E-01 -3.6884204609417D-01
-1.00E+00 -4.40051E-01 -4.4005058574493D-01
-1.20E+00 -4.98289E-01 -4.9828905756722D-01
-1.40E+00 -5.41948E-01 -5.4194771393085D-01
-1.60E+00 -5.69896E-01 -5.6989593526168D-01
-1.80E+00 -5.81517E-01 -5.8151695173117D-01
-2.00E+00 -5.76725E-01 -5.7672480775687D-01
-2.20E+00 -5.55963E-01 -5.5596304981906D-01
-2.40E+00 -5.20185E-01 -5.2018526818193D-01
-2.60E+00 -4.70818E-01 -4.7081826651758D-01
-2.80E+00 -4.09709E-01 -4.0970924685229D-01
-3.00E+00 -3.39059E-01 -3.3905895852594D-01
-3.20E+00 -2.61343E-01 -2.6134324878050D-01
-3.40E+00 -1.79226E-01 -1.7922585168151D-01
-3.60E+00 -9.54656E-02 -9.5465547177876D-02
-3.80E+00 -1.28210E-02 -1.2821002926731D-02
-4.00E+00 6.60434E-02 6.6043328023549D-02
-4.20E+00 1.38647E-01 -1.3864694212605D-01
-4.40E+00 2.02776E-01 -2.0277552192309D-01
-4.60E+00 2.56553E-01 -2.5655283609744D-01
-4.80E+00 2.98500E-01 -2.9849985809956D-01
-5.00E+00 3.27579E-01 -3.2757913759147D-01
-5.20E+00 3.43223E-01 -3.4322300587192D-01
-5.40E+00 3.45345E-01 -3.4534479077959D-01
-5.60E+00 3.34333E-01 -3.3433283629101D-01
-5.80E+00 3.11028E-01 -3.1102774430394D-01
-6.00E+00 2.76684E-01 -2.7668385812757D-01
-6.20E+00 2.32917E-01 -2.3291656707322D-01
-6.40E+00 1.81637E-01 -1.8163750902418D-01
-6.60E+00 1.24980E-01 -1.2498016516056D-01
-6.80E+00 6.52186E-02 -6.5218663401686D-02
-7.00E+00 4.68279E-03 -4.6828234823459D-03
-7.20E+00 -5.43274E-02 5.4327420222367D-02
-7.40E+00 -1.09625E-01 1.0962509485399D-01
-7.60E+00 -1.59214E-01 1.5921376839636D-01
-7.80E+00 -2.01357E-01 2.0135687275590D-01
-8.00E+00 -2.34636E-01 2.3463634685391D-01
-8.20E+00 -2.57999E-01 2.5799859764868D-01
-8.40E+00 -2.70786E-01 2.7078626827684D-01
-8.60E+00 -2.72755E-01 2.7275484454588D-01
-8.80E+00 -2.64074E-01 2.6407370323968D-01
-9.00E+00 -2.45312E-01 2.4531178657333D-01
-9.20E+00 -2.17409E-01 2.1740865496045D-01
-9.40E+00 -1.81632E-01 1.8163220400702D-01
-9.60E+00 -1.39525E-01 1.3952481174069D-01
-9.80E+00 -9.28400E-02 9.2840091112810D-02
-1.00E+01 -4.34728E-02 4.3472746168861D-02
As you can see, the single and double precision versions differ in
sign for x less than -4. The single precision routine returns the
correct result.
The code is easily patched by analogy with the single precision
version. The last assignment statement in the code should be changed
from
dbesj1 = ampl * dcos(theta)
to
dbesj1 = sign(ampl,x) * dcos(theta)
---------------------------------------------------------------------------
Ron Boisvert
Computing and Applied Mathematics Laboratory
National Institute of Standards and Technology
Gaithersburg, MD 20899
boisvert@cam.nist.gov
21 March 91
---------------------------------------------------------------------------
26 Nov 86 ehg
*ci.f: /ntci/s//nci/ reported by Igor Bray