From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 28345 invoked by alias); 17 Aug 2010 15:36:22 -0000 Received: (qmail 28326 invoked by uid 22791); 17 Aug 2010 15:36:20 -0000 X-SWARE-Spam-Status: No, hits=-2.0 required=5.0 tests=AWL,BAYES_00,T_RP_MATCHES_RCVD X-Spam-Check-By: sourceware.org Received: from troutmask.apl.washington.edu (HELO troutmask.apl.washington.edu) (128.208.78.105) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Tue, 17 Aug 2010 15:36:16 +0000 Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.4/8.14.4) with ESMTP id o7HFaBZf035789; Tue, 17 Aug 2010 08:36:11 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.4/8.14.4/Submit) id o7HFaBGl035788; Tue, 17 Aug 2010 08:36:11 -0700 (PDT) (envelope-from sgk) Date: Tue, 17 Aug 2010 15:39:00 -0000 From: Steve Kargl To: Tobias Burnus Cc: Daniel Kraft , gcc patches , gfortran Subject: Re: [Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic Message-ID: <20100817153611.GA35578@troutmask.apl.washington.edu> References: <4C6A82A2.6020909@net-b.de> <4C6A8D42.2000609@domob.eu> <20100817141438.GA35172@troutmask.apl.washington.edu> <4C6AA1A1.2040303@net-b.de> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <4C6AA1A1.2040303@net-b.de> User-Agent: Mutt/1.4.2.3i Mailing-List: contact gcc-patches-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Id: List-Archive: List-Post: List-Help: Sender: gcc-patches-owner@gcc.gnu.org X-SW-Source: 2010-08/txt/msg01277.txt.bz2 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