Re: [Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic

[Tobias Burnus <tobias dot burnus at physik dot fu-berlin dot de>]

Hi all, hi Steve,

Tobias Burnus wrote:
> Here is an updated version. [...] JN now works as it should.
> YN also works, though for some of the examples, the result
> differs quite a bit from the result of mpfr_yn.

I had another look using a C program and seemingly I had
chosen some values which are difficult to get correct. The
following result is using libm of GLIBC 2.5.

(Alternatively, apply the patch and use the Fortran program;
 as currently called in bessel_5.f90, the mpfr functions
 are called and never the library version.)

I think one can use the transformational form, but one might
want to add a big warning that the result of the
transformational intrinsic is obtained using a recurrence
algorithm which may lead to large errors - or at least to
different errors than the elemental function. (Not that the
implemention of non-recurrence algorithms as used in libm
are necessarily correct.)

Some results for YN for N=0 to 10 and different x values
using the C program below. The columns are: N, the libm's
ynf result, the result for the recurrence algorithm, the
difference, the difference divided by machine-precision
epsilon.

Recurrence algorithm for x =     +1.00000000
 1: +8.82569551e-02 | +8.82569551e-02 | +0.00000000e+00 |     +0.00000000
 2: -7.81212807e-01 | -7.81212807e-01 | +0.00000000e+00 |     +0.00000000
 2: -1.65068257e+00 | -1.65068257e+00 | +0.00000000e+00 |     +0.00000000
 3: -5.82151747e+00 | -5.82151747e+00 | +0.00000000e+00 |     +0.00000000
 4: -3.32784195e+01 | -3.32784195e+01 | +0.00000000e+00 |     +0.00000000
 5: -2.60405853e+02 | -2.60405853e+02 | +0.00000000e+00 |     +0.00000000
 6: -2.57078027e+03 | -2.57078027e+03 | +0.00000000e+00 |     +0.00000000
 7: -3.05889570e+04 | -3.05889570e+04 | +0.00000000e+00 |     +0.00000000
 8: -4.25674625e+05 | -4.25674625e+05 | +0.00000000e+00 |     +0.00000000
 9: -6.78020500e+06 | -6.78020500e+06 | +0.00000000e+00 |     +0.00000000
10: -1.21618016e+08 | -1.21618016e+08 | +0.00000000e+00 |     +0.00000000

Recurrence algorithm for x =     +2.00000000
 1: +5.10375679e-01 | +5.10375679e-01 | +0.00000000e+00 |     +0.00000000
 2: -1.07032441e-01 | -1.07032441e-01 | +0.00000000e+00 |     +0.00000000
 2: -6.17408097e-01 | -6.17408097e-01 | +0.00000000e+00 |     +0.00000000
 3: -1.12778378e+00 | -1.12778378e+00 | +0.00000000e+00 |     +0.00000000
 4: -2.76594329e+00 | -2.76594329e+00 | +0.00000000e+00 |     +0.00000000
 5: -9.93598938e+00 | -9.93598938e+00 | +0.00000000e+00 |     +0.00000000
 6: -4.69140053e+01 | -4.69140053e+01 | +0.00000000e+00 |     +0.00000000
 7: -2.71548035e+02 | -2.71548035e+02 | +0.00000000e+00 |     +0.00000000
 8: -1.85392212e+03 | -1.85392212e+03 | +0.00000000e+00 |     +0.00000000
 9: -1.45598291e+04 | -1.45598291e+04 | +0.00000000e+00 |     +0.00000000
10: -1.29184539e+05 | -1.29184539e+05 | +0.00000000e+00 |     +0.00000000

Recurrence algorithm for x =     +3.00000000
 1: +3.76850069e-01 | +3.76850069e-01 | +0.00000000e+00 |     +0.00000000
 2: +3.24674398e-01 | +3.24674398e-01 | +0.00000000e+00 |     +0.00000000
 2: -1.60400465e-01 | -1.60400465e-01 | +0.00000000e+00 |     +0.00000000
 3: -5.38541675e-01 | -5.38541675e-01 | +0.00000000e+00 |     +0.00000000
 4: -9.16682899e-01 | -9.16682899e-01 | +0.00000000e+00 |     +0.00000000
 5: -1.90594614e+00 | -1.90594614e+00 | +0.00000000e+00 |     +0.00000000
 6: -5.43647099e+00 | -5.43647146e+00 | +4.76837158e-07 |     +8.00062370
 7: -1.98399372e+01 | -1.98399391e+01 | +1.90734863e-06 |    +32.00249481
 8: -8.71499023e+01 | -8.71499176e+01 | +1.52587891e-05 |   +256.01995850
 9: -4.44959564e+02 | -4.44959625e+02 | +6.10351562e-05 |  +1024.07983398
10: -2.58260742e+03 | -2.58260791e+03 | +4.88281250e-04 |  +8192.63867188

 1: -1.81523263e-01 | -1.81523263e-01 | +0.00000000e+00 |     +0.00000000
 2: +3.12031418e-01 | +3.12031418e-01 | +0.00000000e+00 |     +0.00000000
 2: +3.21541846e-01 | +3.21541846e-01 | +0.00000000e+00 |     +0.00000000
 3: -2.34589577e-02 | -2.34589577e-02 | +0.00000000e+00 |     +0.00000000
 4: -3.53122234e-01 | -3.53122234e-01 | +0.00000000e+00 |     +0.00000000
 5: -6.10370517e-01 | -6.10370517e-01 | +0.00000000e+00 |     +0.00000000
 6: -1.01634288e+00 | -1.01634264e+00 | -2.38418579e-07 |     -4.00031185
 7: -2.12602520e+00 | -2.12602425e+00 | -9.53674316e-07 |    -16.00124741
 8: -5.66177082e+00 | -5.66176844e+00 | -2.38418579e-06 |    -40.00311661
 9: -1.81989326e+01 | -1.81989250e+01 | -7.62939453e-06 |   -128.00997925
10: -6.78362808e+01 | -6.78362503e+01 | -3.05175781e-05 |   -512.03991699

Recurrence algorithm for x =     +5.00000000
 1: -3.08517635e-01 | -3.08517635e-01 | +0.00000000e+00 |     +0.00000000
 2: +1.47863135e-01 | +1.47863135e-01 | +0.00000000e+00 |     +0.00000000
 2: +3.67662877e-01 | +3.67662877e-01 | +0.00000000e+00 |     +0.00000000
 3: +1.46267161e-01 | +1.46267161e-01 | +0.00000000e+00 |     +0.00000000
 4: -1.92142278e-01 | -1.92142278e-01 | +0.00000000e+00 |     +0.00000000
 5: -4.53694820e-01 | -4.53694820e-01 | +0.00000000e+00 |     +0.00000000
 6: -7.15247393e-01 | -7.15247393e-01 | +0.00000000e+00 |     +0.00000000
 7: -1.26289904e+00 | -1.26289904e+00 | +0.00000000e+00 |     +0.00000000
 8: -2.82086992e+00 | -2.82086992e+00 | +0.00000000e+00 |     +0.00000000
 9: -7.76388502e+00 | -7.76388502e+00 | +0.00000000e+00 |     +0.00000000
10: -2.51291161e+01 | -2.51291180e+01 | +1.90734863e-06 |    +32.00249481

Tobias

The C program, link with -lm; without _GNU_SOURCE "math.h"
does not provide the Bessel function as function prototype.

#define _GNU_SOURCE
#include <math.h>
#include <stdio.h>

int main ()
{
  float last1, last2, x, x2rev, cur, eps = 5.96e-08f;
  x = 1.00f;

  printf("Recurrence algorithm for x = %+15.8f\n", x);

  last1 = ynf(0, x);
  printf ("%2i: %+15.8e | %+15.8e | %+15.8e | %+15.8f\n", 1,
          ynf(0, x), last1, (ynf(0, x) - last1),
          (ynf(0, x) - last1)/eps);

  last2 = ynf(1, x);
  printf ("%2i: %+15.8e | %+15.8e | %+15.8e | %+15.8f\n", 2,
          ynf(1, x), last2, (ynf(1, x) - last2),
          (ynf(1, x) - last2)/eps);

  x2rev = 2.0f/x;

  for (int i = 2; i < 11; i++)
    {
      cur = x2rev * (i-1) * last2 - last1;
      printf ("%2i: %+15.8e | %+15.8e | %+15.8e | %+15.8f\n", i,
              ynf(i, x), cur, (ynf(i, x) - cur),
              (ynf(i, x) - cur)/eps);
      last1 = last2;
      last2 = cur;
    }
  return 0;
}