[Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic
Steve Kargl
sgk@troutmask.apl.washington.edu
Fri Aug 20 01:05:00 GMT 2010
On Thu, Aug 19, 2010 at 09:30:10AM +0200, Tobias Burnus wrote:
> On 08/19/2010 05:34 AM, Steve Kargl wrote:
> >I read your patch and it appears to be correct,
> >except that 2.0 in mpfr_ui_div should probably
> >be an integer 2.
>
> Fixed. I also added a link to DLMF for the recurrence algorithm.
>
> >I haven't had time to investigate your program and
> >findings, but I have no reason to suspect that your
> >results are flawed. I am somewhat puzzled by them.
> >After all the work you did, I hesitate to suggest
> >that I would rather have gfortran produce correct
> >results slowly (ie, n2 - n1 + 1 calls to mpfr_{jn,yn})
> >than to get possibly wrong results fast via
> >recurrence algorithm. If you want to change back
> >to your original algorithm that's fine with me;
> >otherwise, the current patch can be committed and
> >I'll see if I can find time to poke at it.
>
> I have committed the version with the recurrence algorithm but added a
> note to intrinsics.texi that the transformational version uses a
> recurrence algorithm which might lead to different results than invoking
> the elemental function. From what I have heard from the J3/WG5 committee
> the transformational version has been added explicitly because it allows
> the implementation of the recurrence algorithm. Thus, it makes sense to
> use it.
>
> The next step is to add an implementation for run-time evaluation, which
> will also use the recurrence algorithm - thus the mpfr and the run-time
> version will give the same result, which is useful by itself.
>
> I think it were be useful if you could find time to poke at it. One
> could also consider using the algorithm for JN but not for YN - but
> first one should do more checks to make sure that I haven't just been
> lucky that JN worked so well.
>
I compared the gfortran results against my math library. I think
you're see that your algorithm is using a more accurate approximation
for the initial two functions in the recursion algorithm. Note
my math library routines use y0[f], y1[f], and jn[f] from libm
to initialize the recursion and thus issues in libm may potentially
affect my results. When we observe a large relative error, I think
the algorithm may be hitting a castastrophic cancelation issue (ie,
subtraction of nearly equal large number).
There may be a method to improve the values computed, but it
would substantially increase the computation. Instead of
recursing on the actually estimates of the bessel function,
we could recurse on the significands by careful bookkeeping
of the exponents (which may require an additional array), e.g.,
Let jnm1 = j(n-1), jn = j(n), and jnp1 = j(n+1):
jnm1 = (2n/x) * jn - jnp1
= (2n/x) * fn * 2**en - fnp1 * 2**enp1
Assume en > enp1, then
= 2**en * ((2n/x) * fn - fnp1 * 2**(enp1-en))
and due to the recursion |enp1 - en| is normally small
so
fnm1 = (2n/x) * fn - fnp1 * 2**(enp1-en)
I used this to compute fairly large sequence, ie., n1 = 0
and n2 = 15000 to 20000 for spherical bessel functions.
I haven't proved that it avoids catastrophic cancellation
issues, it does allow long sequence that otherwise would
be +-Inf or NaN.
Here are some numbers from my comparison. The first column is
the gfortran result for bessel_jn(0,20,X) where X = 1., 2., and
10. The 2nd column is my single precision result and the
third column is my double precision result. The last two columns
are the relative errors, measured in epsilons, in comparison to
my SP and DP results, respectively. The last line in each
group is the result from bessel_jn(20,X), ie., the MPFR value
from the constant folding and the relative error.
j(1.)
0 7.651976E-01 7.651977E-01 7.65197687E-01 0.65 0.57
1 4.400505E-01 4.400506E-01 4.40050586E-01 1.14 0.83
2 1.149035E-01 1.149035E-01 1.14903485E-01 1.09 0.92
3 1.956335E-02 1.956335E-02 1.95633540E-02 0.80 0.55
4 2.476639E-03 2.476639E-03 2.47663896E-03 0.00 0.28
5 2.497577E-04 2.497577E-04 2.49757730E-04 0.00 0.26
6 2.093834E-05 2.093834E-05 2.09383380E-05 0.00 0.04
7 1.502326E-06 1.502326E-06 1.50232582E-06 0.00 0.23
8 9.422345E-08 9.422345E-08 9.42234417E-08 0.00 0.30
9 5.249250E-09 5.249250E-09 5.24925018E-09 0.00 0.30
10 2.630615E-10 2.630615E-10 2.63061512E-10 0.00 0.15
11 1.198007E-11 1.198007E-11 1.19800675E-11 0.00 0.25
12 4.999718E-13 4.999718E-13 4.99971818E-13 0.00 0.15
13 1.925617E-14 1.925617E-14 1.92561676E-14 0.74 0.80
14 6.885407E-16 6.885408E-16 6.88540820E-16 0.64 0.88
15 2.297531E-17 2.297532E-17 2.29753153E-17 0.60 0.58
16 7.186396E-19 7.186397E-19 7.18639659E-19 0.60 0.37
17 2.115375E-20 2.115376E-20 2.11537557E-20 0.64 0.36
18 5.880344E-22 5.880344E-22 5.88034457E-22 0.00 0.28
19 1.548478E-23 1.548478E-23 1.54847844E-23 0.00 0.11
20 3.873503E-25 3.873503E-25 3.87350301E-25 0.00 0.04
3.873503E-25 0.00
j(2.)
0 2.238908E-01 2.238908E-01 2.23890779E-01 1.67 1.03
1 5.767247E-01 5.767248E-01 5.76724808E-01 1.73 1.45
2 3.528340E-01 3.528340E-01 3.52834029E-01 2.13 1.72
3 1.289432E-01 1.289433E-01 1.28943249E-01 2.91 1.93
4 3.399571E-02 3.399572E-02 3.39957198E-02 2.76 2.34
5 7.039628E-03 7.039630E-03 7.03962976E-03 2.22 1.55
6 1.202429E-03 1.202429E-03 1.20242897E-03 2.44 1.80
7 1.749440E-04 1.749441E-04 1.74944075E-04 2.09 1.59
8 2.217955E-05 2.217955E-05 2.21795523E-05 2.06 1.25
9 2.492343E-06 2.492344E-06 2.49234344E-06 1.53 1.03
10 2.515386E-07 2.515386E-07 2.51538628E-07 0.95 0.66
11 2.304285E-08 2.304285E-08 2.30428476E-08 0.65 0.02
12 1.932695E-09 1.932695E-09 1.93269515E-09 0.96 0.15
13 1.494942E-10 1.494942E-10 1.49494201E-10 0.78 0.35
14 1.072946E-11 1.072947E-11 1.07294645E-11 0.68 0.22
15 7.183016E-13 7.183017E-13 7.18301636E-13 0.63 0.28
16 4.506006E-14 4.506006E-14 4.50600590E-14 0.63 0.01
17 2.659308E-15 2.659308E-15 2.65930781E-15 0.67 0.03
18 1.481737E-16 1.481737E-16 1.48173725E-16 0.75 0.09
19 7.819243E-18 7.819244E-18 7.81924327E-18 0.89 0.07
20 3.918973E-19 3.918973E-19 3.91897281E-19 0.55 0.16
3.918973E-19 0.00
j(10.)
0 -2.459358E-01 -2.459358E-01 -2.45935764E-01 2.03 2.15
1 4.347282E-02 4.347274E-02 4.34727462E-02 14.38 14.06
2 2.546304E-01 2.546303E-01 2.54630314E-01 2.95 2.41
3 5.837934E-02 5.837938E-02 5.83793793E-02 6.42 6.05
4 -2.196028E-01 -2.196027E-01 -2.19602686E-01 3.98 4.10
5 -2.340616E-01 -2.340615E-01 -2.34061528E-01 1.07 1.46
6 -1.445878E-02 -1.445884E-02 -1.44588421E-02 38.90 38.62
7 2.167110E-01 2.167109E-01 2.16710918E-01 4.61 4.90
8 3.178542E-01 3.178541E-01 3.17854127E-01 2.36 2.61
9 2.918558E-01 2.918557E-01 2.91855685E-01 1.71 1.93
10 2.074862E-01 2.074861E-01 2.07486107E-01 1.81 1.86
11 1.231165E-01 1.231165E-01 1.23116528E-01 1.02 1.18
12 6.337026E-02 6.337025E-02 6.33702550E-02 0.99 0.35
13 2.897208E-02 2.897208E-02 2.89720839E-02 0.54 0.60
14 1.195716E-02 1.195716E-02 1.19571632E-02 0.00 0.17
15 4.507973E-03 4.507973E-03 4.50797314E-03 0.87 0.34
16 1.566756E-03 1.566756E-03 1.56675619E-03 0.62 0.49
17 5.056467E-04 5.056466E-04 5.05646670E-04 0.97 0.35
18 1.524425E-04 1.524425E-04 1.52442485E-04 0.80 0.66
19 4.314628E-05 4.314628E-05 4.31462775E-05 0.00 0.32
20 1.151337E-05 1.151337E-05 1.15133692E-05 0.00 0.24
1.151337E-05 0.00
y(1.)
0 8.825696E-02 8.825696E-02 8.82569642E-02 0.00 0.15
1 -7.812128E-01 -7.812128E-01 -7.81212821E-01 0.00 0.16
2 -1.650683E+00 -1.650683E+00 -1.65068261E+00 0.00 0.19
3 -5.821517E+00 -5.821517E+00 -5.82151761E+00 0.00 0.20
4 -3.327842E+01 -3.327842E+01 -3.32784230E+01 0.96 0.89
5 -2.604059E+02 -2.604059E+02 -2.60405867E+02 0.98 0.43
6 -2.570780E+03 -2.570781E+03 -2.57078024E+03 0.80 0.10
7 -3.058896E+04 -3.058896E+04 -3.05889571E+04 1.07 0.01
8 -4.256746E+05 -4.256747E+05 -4.25674618E+05 1.23 0.13
9 -6.780205E+06 -6.780206E+06 -6.78020494E+06 1.24 0.08
10 -1.216180E+08 -1.216180E+08 -1.21618014E+08 1.10 0.12
11 -2.425580E+09 -2.425581E+09 -2.42558008E+09 0.89 0.72
12 -5.324115E+10 -5.324116E+10 -5.32411438E+10 1.29 0.57
13 -1.275362E+12 -1.275362E+12 -1.27536187E+12 1.72 0.11
14 -3.310617E+13 -3.310617E+13 -3.31061675E+13 1.59 0.10
15 -9.256974E+14 -9.256975E+14 -9.25697328E+14 1.22 0.55
16 -2.773782E+16 -2.773782E+16 -2.77378137E+16 1.30 0.64
17 -8.866844E+17 -8.866845E+17 -8.86684340E+17 1.30 0.42
18 -3.011953E+19 -3.011954E+19 -3.01195297E+19 1.22 0.39
19 -1.083416E+21 -1.083417E+21 -1.08341639E+21 1.63 0.08
20 -4.113970E+22 -4.113971E+22 -4.11397031E+22 1.84 0.12
-4.113970E+22 0.00
y(2.)
0 5.103757E-01 5.103757E-01 5.10375673E-01 0.00 0.10
1 -1.070324E-01 -1.070324E-01 -1.07032432E-01 0.58 0.13
2 -6.174081E-01 -6.174081E-01 -6.17408104E-01 0.00 0.10
3 -1.127784E+00 -1.127784E+00 -1.12778378E+00 0.00 0.01
4 -2.765943E+00 -2.765943E+00 -2.76594323E+00 0.00 0.19
5 -9.935989E+00 -9.935989E+00 -9.93598913E+00 0.00 0.21
6 -4.691401E+01 -4.691401E+01 -4.69140024E+01 0.00 0.51
7 -2.715480E+02 -2.715480E+02 -2.71548025E+02 0.00 0.29
8 -1.853922E+03 -1.853922E+03 -1.85392218E+03 0.55 0.25
9 -1.455983E+04 -1.455983E+04 -1.45598294E+04 0.56 0.16
10 -1.291845E+05 -1.291845E+05 -1.29184542E+05 0.51 0.20
11 -1.277286E+06 -1.277286E+06 -1.27728559E+06 0.82 0.61
12 -1.392096E+07 -1.392096E+07 -1.39209570E+07 1.21 1.19
13 -1.657742E+08 -1.657742E+08 -1.65774198E+08 0.81 1.12
14 -2.141143E+09 -2.141144E+09 -2.14114362E+09 0.50 0.76
15 -2.981024E+10 -2.981024E+10 -2.98102365E+10 0.00 0.30
16 -4.450124E+11 -4.450124E+11 -4.45012403E+11 0.00 0.49
17 -7.090388E+12 -7.090388E+12 -7.09038822E+12 0.00 0.79
18 -1.200916E+14 -1.200916E+14 -1.20091587E+14 0.00 1.06
19 -2.154558E+15 -2.154558E+15 -2.15455818E+15 0.00 1.25
20 -4.081651E+16 -4.081651E+16 -4.08165139E+16 0.00 1.06
-4.081651E+16 0.88
y(10.)
0 5.567117E-02 5.567117E-02 5.56711673E-02 0.56 0.10
1 2.490154E-01 2.490154E-01 2.49015424E-01 0.50 0.12
2 -5.868081E-03 -5.868082E-03 -5.86808244E-03 2.00 2.06
3 -2.513627E-01 -2.513627E-01 -2.51362657E-01 0.99 0.19
4 -1.449495E-01 -1.449495E-01 -1.44949512E-01 0.00 0.79
5 1.354030E-01 1.354031E-01 1.35403048E-01 0.92 0.66
6 2.803526E-01 2.803526E-01 2.80352560E-01 0.89 0.09
7 2.010201E-01 2.010201E-01 2.01020024E-01 0.00 1.59
8 1.075506E-03 1.075499E-03 1.07547373E-03 59.93 253.31
9 -1.992992E-01 -1.992993E-01 -1.99299266E-01 0.63 0.83
10 -3.598142E-01 -3.598142E-01 -3.59814152E-01 0.00 0.35
11 -5.203291E-01 -5.203291E-01 -5.20329039E-01 0.96 1.28
12 -7.849100E-01 -7.849098E-01 -7.84909733E-01 1.91 2.47
13 -1.363455E+00 -1.363454E+00 -1.36345432E+00 2.20 3.07
14 -2.760073E+00 -2.760072E+00 -2.76007150E+00 2.90 3.68
15 -6.364748E+00 -6.364747E+00 -6.36474588E+00 2.51 3.43
16 -1.633417E+01 -1.633417E+01 -1.63341661E+01 2.94 3.63
17 -4.590461E+01 -4.590459E+01 -4.59045857E+01 2.79 3.68
18 -1.397415E+02 -1.397414E+02 -1.39741425E+02 3.66 4.53
19 -4.571648E+02 -4.571646E+02 -4.57164546E+02 3.92 5.13
20 -1.597485E+03 -1.597484E+03 -1.59748385E+03 3.85 5.33
-1.597484E+03 5.13
--
steve
More information about the Fortran
mailing list