contrib/mpfr/jn.c

   1 /* mpfr_j0, mpfr_j1, mpfr_jn -- Bessel functions of 1st kind, integer order.
   2    http://www.opengroup.org/onlinepubs/009695399/functions/j0.html
   3
   4 Copyright 2007, 2008, 2009 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Cacao projects, INRIA.
   6
   7 This file is part of the GNU MPFR Library.
   8
   9 The GNU MPFR Library is free software; you can redistribute it and/or modify
  10 it under the terms of the GNU Lesser General Public License as published by
  11 the Free Software Foundation; either version 2.1 of the License, or (at your
  12 option) any later version.
  13
  14 The GNU MPFR Library is distributed in the hope that it will be useful, but
  15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17 License for more details.
  18
  19 You should have received a copy of the GNU Lesser General Public License
  20 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  21 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  22 MA 02110-1301, USA. */
  23
  24 #define MPFR_NEED_LONGLONG_H
  25 #include "mpfr-impl.h"
  26
  27 /* Relations: j(-n,z) = (-1)^n j(n,z)
  28               j(n,-z) = (-1)^n j(n,z)
  29 */
  30
  31 static int mpfr_jn_asympt (mpfr_ptr, long, mpfr_srcptr, mp_rnd_t);
  32
  33 int
  34 mpfr_j0 (mpfr_ptr res, mpfr_srcptr z, mp_rnd_t r)
  35 {
  36   return mpfr_jn (res, 0, z, r);
  37 }
  38
  39 int
  40 mpfr_j1 (mpfr_ptr res, mpfr_srcptr z, mp_rnd_t r)
  41 {
  42   return mpfr_jn (res, 1, z, r);
  43 }
  44
  45 /* Estimate k0 such that z^2/4 = k0 * (k0 + n)
  46    i.e., (sqrt(n^2+z^2)-n)/2 = n/2 * (sqrt(1+(z/n)^2) - 1).
  47    Return min(2*k0/log(2), ULONG_MAX).
  48 */
  49 static unsigned long
  50 mpfr_jn_k0 (long n, mpfr_srcptr z)
  51 {
  52   mpfr_t t, u;
  53   unsigned long k0;
  54
  55   mpfr_init2 (t, 32);
  56   mpfr_init2 (u, 32);
  57   mpfr_div_si (t, z, n, GMP_RNDN);
  58   mpfr_sqr (t, t, GMP_RNDN);
  59   mpfr_add_ui (t, t, 1, GMP_RNDN);
  60   mpfr_sqrt (t, t, GMP_RNDN);
  61   mpfr_sub_ui (t, t, 1, GMP_RNDN);
  62   mpfr_mul_si (t, t, n, GMP_RNDN);
  63   /* the following is a 32-bit approximation to nearest of log(2) */
  64   mpfr_set_str_binary (u, "0.10110001011100100001011111111");
  65   mpfr_div (t, t, u, GMP_RNDN);
  66   if (mpfr_fits_ulong_p (t, GMP_RNDN))
  67     k0 = mpfr_get_ui (t, GMP_RNDN);
  68   else
  69     k0 = ULONG_MAX;
  70   mpfr_clear (t);
  71   mpfr_clear (u);
  72   return k0;
  73 }
  74
  75 int
  76 mpfr_jn (mpfr_ptr res, long n, mpfr_srcptr z, mp_rnd_t r)
  77 {
  78   int inex;
  79   unsigned long absn;
  80   mp_prec_t prec, pbound, err;
  81   mp_exp_t exps, expT;
  82   mpfr_t y, s, t, absz;
  83   unsigned long k, zz, k0;
  84   MPFR_ZIV_DECL (loop);
  85
  86   MPFR_LOG_FUNC (("x[%#R]=%R n=%d rnd=%d", z, z, n, r),
  87                  ("y[%#R]=%R", res, res));
  88
  89   absn = SAFE_ABS (unsigned long, n);
  90
  91   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (z)))
  92     {
  93       if (MPFR_IS_NAN (z))
  94         {
  95           MPFR_SET_NAN (res);
  96           MPFR_RET_NAN;
  97         }
  98       /* j(n,z) tends to zero when z goes to +Inf or -Inf, oscillating around
  99          0. We choose to return +0 in that case. */
 100       else if (MPFR_IS_INF (z)) /* FIXME: according to j(-n,z) = (-1)^n j(n,z)
 101                                    we might want to give a sign depending on
 102                                    z and n */
 103         return mpfr_set_ui (res, 0, r);
 104       else /* z=0: j(0,0)=1, j(n odd,+/-0) = +/-0 if n > 0, -/+0 if n < 0,
 105               j(n even,+/-0) = +0 */
 106         {
 107           if (n == 0)
 108             return mpfr_set_ui (res, 1, r);
 109           else if (absn & 1) /* n odd */
 110             return (n > 0) ? mpfr_set (res, z, r) : mpfr_neg (res, z, r);
 111           else /* n even */
 112             return mpfr_set_ui (res, 0, r);
 113         }
 114     }
 115
 116   /* check for tiny input for j0: j0(z) = 1 - z^2/4 + ..., more precisely
 117      |j0(z) - 1| <= z^2/4 for -1 <= z <= 1. */
 118   if (n == 0)
 119     MPFR_FAST_COMPUTE_IF_SMALL_INPUT (res, __gmpfr_one, -2 * MPFR_GET_EXP (z),
 120                                       2, 0, r, return _inexact);
 121
 122   /* idem for j1: j1(z) = z/2 - z^3/16 + ..., more precisely
 123      |j1(z) - z/2| <= |z^3|/16 for -1 <= z <= 1, with the sign of j1(z) - z/2
 124      being the opposite of that of z. */
 125   if (n == 1)
 126     /* we first compute 2j1(z) = z - z^3/8 + ..., then divide by 2 using
 127        the "extra" argument of MPFR_FAST_COMPUTE_IF_SMALL_INPUT. */
 128     MPFR_FAST_COMPUTE_IF_SMALL_INPUT (res, z, -2 * MPFR_GET_EXP (z), 3,
 129                                       0, r, mpfr_div_2ui (res, res, 1, r));
 130
 131   /* we can use the asymptotic expansion as soon as |z| > p log(2)/2,
 132      but to get some margin we use it for |z| > p/2 */
 133   pbound = MPFR_PREC (res) / 2 + 3;
 134   MPFR_ASSERTN (pbound <= ULONG_MAX);
 135   MPFR_ALIAS (absz, z, 1, MPFR_EXP (z));
 136   if (mpfr_cmp_ui (absz, pbound) > 0)
 137     {
 138       inex = mpfr_jn_asympt (res, n, z, r);
 139       if (inex != 0)
 140         return inex;
 141     }
 142
 143   mpfr_init2 (y, 32);
 144
 145   /* check underflow case: |j(n,z)| <= 1/sqrt(2 Pi n) (ze/2n)^n
 146      (see algorithms.tex) */
 147   if (absn > 0)
 148     {
 149       /* the following is an upper 32-bit approximation of exp(1)/2 */
 150       mpfr_set_str_binary (y, "1.0101101111110000101010001011001");
 151       if (MPFR_SIGN(z) > 0)
 152         mpfr_mul (y, y, z, GMP_RNDU);
 153       else
 154         {
 155           mpfr_mul (y, y, z, GMP_RNDD);
 156           mpfr_neg (y, y, GMP_RNDU);
 157         }
 158       mpfr_div_ui (y, y, absn, GMP_RNDU);
 159       /* now y is an upper approximation of |ze/2n|: y < 2^EXP(y),
 160          thus |j(n,z)| < 1/2*y^n < 2^(n*EXP(y)-1).
 161          If n*EXP(y) < __gmpfr_emin then we have an underflow.
 162          Warning: absn is an unsigned long. */
 163       if ((MPFR_EXP(y) < 0 && absn > (unsigned long) (-__gmpfr_emin))
 164           || (absn <= (unsigned long) (-MPFR_EMIN_MIN) &&
 165               MPFR_EXP(y) < __gmpfr_emin / (mp_exp_t) absn))
 166         {
 167           mpfr_clear (y);
 168           return mpfr_underflow (res, (r == GMP_RNDN) ? GMP_RNDZ : r,
 169                          (n % 2) ? ((n > 0) ? MPFR_SIGN(z) : -MPFR_SIGN(z))
 170                                  : MPFR_SIGN_POS);
 171         }
 172     }
 173
 174   mpfr_init (s);
 175   mpfr_init (t);
 176
 177   /* the logarithm of the ratio between the largest term in the series
 178      and the first one is roughly bounded by k0, which we add to the
 179      working precision to take into account this cancellation */
 180   k0 = mpfr_jn_k0 (absn, z);
 181   prec = MPFR_PREC (res) + k0 + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (res)) + 3;
 182
 183   MPFR_ZIV_INIT (loop, prec);
 184   for (;;)
 185     {
 186       mpfr_set_prec (y, prec);
 187       mpfr_set_prec (s, prec);
 188       mpfr_set_prec (t, prec);
 189       mpfr_pow_ui (t, z, absn, GMP_RNDN); /* z^|n| */
 190       mpfr_mul (y, z, z, GMP_RNDN);       /* z^2 */
 191       zz = mpfr_get_ui (y, GMP_RNDU);
 192       MPFR_ASSERTN (zz < ULONG_MAX);
 193       mpfr_div_2ui (y, y, 2, GMP_RNDN);   /* z^2/4 */
 194       mpfr_fac_ui (s, absn, GMP_RNDN);    /* |n|! */
 195       mpfr_div (t, t, s, GMP_RNDN);
 196       if (absn > 0)
 197         mpfr_div_2ui (t, t, absn, GMP_RNDN);
 198       mpfr_set (s, t, GMP_RNDN);
 199       exps = MPFR_EXP (s);
 200       expT = exps;
 201       for (k = 1; ; k++)
 202         {
 203           mpfr_mul (t, t, y, GMP_RNDN);
 204           mpfr_neg (t, t, GMP_RNDN);
 205           if (k + absn <= ULONG_MAX / k)
 206             mpfr_div_ui (t, t, k * (k + absn), GMP_RNDN);
 207           else
 208             {
 209               mpfr_div_ui (t, t, k, GMP_RNDN);
 210               mpfr_div_ui (t, t, k + absn, GMP_RNDN);
 211             }
 212           exps = MPFR_EXP (t);
 213           if (exps > expT)
 214             expT = exps;
 215           mpfr_add (s, s, t, GMP_RNDN);
 216           exps = MPFR_EXP (s);
 217           if (exps > expT)
 218             expT = exps;
 219           if (MPFR_EXP (t) + (mp_exp_t) prec <= MPFR_EXP (s) &&
 220               zz / (2 * k) < k + n)
 221             break;
 222         }
 223       /* the error is bounded by (4k^2+21/2k+7) ulp(s)*2^(expT-exps)
 224          <= (k+2)^2 ulp(s)*2^(2+expT-exps) */
 225       err = 2 * MPFR_INT_CEIL_LOG2(k + 2) + 2 + expT - MPFR_EXP (s);
 226       if (MPFR_LIKELY (MPFR_CAN_ROUND (s, prec - err, MPFR_PREC(res), r)))
 227           break;
 228       MPFR_ZIV_NEXT (loop, prec);
 229     }
 230   MPFR_ZIV_FREE (loop);
 231
 232   inex = ((n >= 0) || ((n & 1) == 0)) ? mpfr_set (res, s, r)
 233                                       : mpfr_neg (res, s, r);
 234
 235   mpfr_clear (y);
 236   mpfr_clear (s);
 237   mpfr_clear (t);
 238
 239   return inex;
 240 }
 241
 242 #define MPFR_JN
 243 #include "jyn_asympt.c"