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

[Steve Kargl]

On Tue, Aug 17, 2010 at 06:59:12PM +0200, Tobias Burnus wrote:
>  On 08/17/2010 05:36 PM, Steve Kargl wrote:
> >>I disagree. Negative values are perfectly valid thus I would prefer
> >>adding such a restriction only for -std=f2003.
> >Then you may need to fix your simplicification function.
> 
> Let's do it differently: For the existing elemental function, we allow 
> non-negative values but for the transformational one, we don't.
> 
> >+  mpz_init_set_ui (result->shape[0], MAX (n2-n1+1, 0));
> >+
> >+  for (i = n1; i<= n2; ++i)
> >+    {
> >
> >What happens if n2 = -4 and n1 = 1?  Note, the standard
> >does not specify that n2>  n1.  It is implied by
> >
> >     Case (ii): The result of BESSEL JN (N1, N2, X) is a rank-one
> >       array with extent MAX (N2-N1+1, 0).
> >
> >but there is no requirement for this ordering.  It may be prudent to
> >add a check that n2>  n1.
> 
> I think that's wrong. First, n2 == n1 is perfectly valid; secondly, if 
> n2 < n1, one gets a zero-sized array, which might be not terribly 
> useful, but is perfectly valid.

Whoops.  I misread the condition in the for-loop,
and thought that you would enter an infinite loop
stomping on memory.


> >>(J(-m,x) = (-1)**m * J(m,x) -- and analogously for YN - at least for
> >>integral m; using the Gamma function for negative m won't work, as it
> >>becomes +/-infinite for negative integer values). Thus, there is no need
> >>for negative "m", but also no need to reject them. As gfortran allowed
> >>negative "m" so far, I think it should continue to do so for -std=gnu.)
> >I would prefer to remove this extension.
> Why?

Experience.  I've been computing/using J() and Y() for more
than 23 years.  Getting a correct and accurate value is
is difficult.  The functions in libm have horrible accuracy 
near their zeros.

> >>Do you mean the following algorithm?
> >>x2rev = 2.0/x
> >>J(N-1, x) = x2rev * N * J(N, x) - J(N+1, x)
> >>Y(N+1, x) = x2rev * N * Y(N, x) - Y(N-1, x
> >Yes.
> >    if (x == 0.e0_knd) then  ! Avoid division by zero
> >       j = 0
> The x == 0 is a good point!
> 
> How about the attached patch?

I'll look over the patch in more detail in a few hours.

> PS: Steve, are you interested in implementing a TREE or libgfortran 
> version of this algorithm? I think using a libgfortran version makes 
> more sense than implementing it as TREE; what do you think?

I think we want to use TREE only because mpfr_{jn,yn} will
give accurate values for the two initial values in the
recursion scheme.  As noted above, the libm routines give
horrible accuracy near their zeros.  For j0(), I found 
over a million ULP at its 1st zero.  If you're interested
in the details

http://www.freebsd.org/cgi/query-pr.cgi?pr=standards/142803

-- 
Steve