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

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

On Tue, Aug 17, 2010 at 04:50:09PM +0200, Tobias Burnus wrote:
>  On 08/17/2010 04:14 PM, Steve Kargl wrote:
> >Not OK. :-)
> >In check.c,
> >+  if (scalar_check (n2, 0) == FAILURE)
> >+    return FAILURE;
> >The 0 should 1.
> Thanks for spotting this.
> 
> >Also note that if n1 and n2 are to be nonnegative.
> >So, one should add,
> >
> >    if (nonnegative_check ('n1', n1) == FAILURE)
> >      return FAILURE;
> >    if (nonnegative_check ('n2', n2) == FAILURE)
> >      return FAILURE;
> 
> 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.

+  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. 

> (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.

> >In the simplification route
> >
> >+  if (order->expr_type == EXPR_CONSTANT)
> >+    {
> >+      n = mpz_get_si (order->value.integer);
> >+      if (n<  0&&  gfc_notify_std (GFC_STD_GNU, "Extension: Negativ 
> >argument "
> >+                                   "N at %L",&order->where) == FAILURE)
> >+       return&gfc_bad_expr;
> >+    }
> >
> >    if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
> >      return NULL;
> >
> >Please move the 2nd if-statement to the top of the function.  There's
> >no reason to possibly issue an error for n<  0, if once that is fixed
> >simplification terminates due to a non-constant x.
> 
> Well, the error is useful, unless one already puts the error into 
> check.c. But I agree that one can move the check into check.c.

Yes, the check belong in check.c.

> >Finally, note that the mpfr_{jn,yn} routines are called n2 - n1 + 1
> >times, which may be very slow for large n2 - n1 + 1.  Bessel functions
> >are stable for downward recursion, and Neumann functions are stable
> >for upward recursion.  It may be better to use recursion, which
> >involves two evaluations of mfpr_{jn,yn} and then some simple
> >multiplications and additions
> 
> 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 you have Abromowitz and Stegun, look up Miller's
algorithm for computing Bessel functions.  The new NIST
rewrite has

http://dlmf.nist.gov/10.74#iv

The guts of my bessel function routine has

   if (x == 0.e0_knd) then  ! Avoid division by zero
      j = 0
      j(0) = 1
   else
      !  
      ! Do the downward recursion.
      !
      jmx1 = besjn(n + 2, x)
      jmx  = besjn(n + 1, x)
      ix = 2 / x
      do i = n + 1, 1, -1
         j(i - 1) = i * ix * jmx - jmx1
         jmx1 = jmx
         jmx = j(i - 1)
      end do

If the initial values of J(N,x) and J(N+1, x) are accurate
approximations then the above generates the required sequence.
If the values are inaccurate, the generated sequence is 
porportional to the desired sequence of Bessel functions.
One then to normalize the sequence, for my routine I do

      jmx = cj0(x)   ! J(0)
      jmx1 = cj1(x)  ! J(1)
      if (abs(jmx) >= abs(jmx1)) then
         jmx = jmx  / j(0)
      else
         jmx = jmx1 / j(1)
      end if
      j = jmx * j
   end if


-- 
Steve