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

[Steve Kargl <sgk at troutmask dot apl dot washington dot edu>]

On Wed, Aug 18, 2010 at 12:37:29PM +0200, Tobias Burnus wrote:
> 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.)
> 

I read your patch and it appears to be correct,
except that 2.0 in mpfr_ui_div should probably
be an integer 2.

Thanks for moving the checks into check.c.

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.

-- 
Steve