This is the mail archive of the
gcc-patches@gcc.gnu.org
mailing list for the GCC project.
Re: [Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic
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.
(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.)
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.
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)
Tobias