contrib/mpfr/src/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, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   5 Contributed by the AriC and Caramel 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 3 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.LESSER.  If not, see
  21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  22 51 Franklin St, Fifth Floor, Boston, 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, mpfr_rnd_t);
  32
  33 int
  34 mpfr_j0 (mpfr_ptr res, mpfr_srcptr z, mpfr_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, mpfr_rnd_t r)
  41 {
  42   return mpfr_jn (res, 1, z, r);
  43 }
  44
  45 /* Estimate k1 such that z^2/4 = k1 * (k1 + n)
  46    i.e., k1 = (sqrt(n^2+z^2)-n)/2 = n/2 * (sqrt(1+(z/n)^2) - 1) if n != 0.
  47    Return k0 = min(2*k1/log(2), ULONG_MAX).
  48 */
  49 static unsigned long
  50 mpfr_jn_k0 (unsigned 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   if (n == 0)
  58     {
  59       mpfr_abs (t, z, MPFR_RNDN);  /* t = 2*k1 */
  60     }
  61   else
  62     {
  63       mpfr_div_ui (t, z, n, MPFR_RNDN);
  64       mpfr_sqr (t, t, MPFR_RNDN);
  65       mpfr_add_ui (t, t, 1, MPFR_RNDN);
  66       mpfr_sqrt (t, t, MPFR_RNDN);
  67       mpfr_sub_ui (t, t, 1, MPFR_RNDN);
  68       mpfr_mul_ui (t, t, n, MPFR_RNDN);  /* t = 2*k1 */
  69     }
  70   /* the following is a 32-bit approximation to nearest to 1/log(2) */
  71   mpfr_set_str_binary (u, "1.0111000101010100011101100101001");
  72   mpfr_mul (t, t, u, MPFR_RNDN);
  73   if (mpfr_fits_ulong_p (t, MPFR_RNDN))
  74     k0 = mpfr_get_ui (t, MPFR_RNDN);
  75   else
  76     k0 = ULONG_MAX;
  77   mpfr_clear (t);
  78   mpfr_clear (u);
  79   return k0;
  80 }
  81
  82 int
  83 mpfr_jn (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r)
  84 {
  85   int inex;
  86   int exception = 0;
  87   unsigned long absn;
  88   mpfr_prec_t prec, pbound, err;
  89   mpfr_uprec_t uprec;
  90   mpfr_exp_t exps, expT, diffexp;
  91   mpfr_t y, s, t, absz;
  92   unsigned long k, zz, k0;
  93   MPFR_GROUP_DECL(g);
  94   MPFR_SAVE_EXPO_DECL (expo);
  95   MPFR_ZIV_DECL (loop);
  96
  97   MPFR_LOG_FUNC
  98     (("n=%d x[%Pu]=%.*Rg rnd=%d", n, mpfr_get_prec (z), mpfr_log_prec, z, r),
  99      ("res[%Pu]=%.*Rg inexact=%d",
 100       mpfr_get_prec (res), mpfr_log_prec, res, inex));
 101
 102   absn = SAFE_ABS (unsigned long, n);
 103
 104   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (z)))
 105     {
 106       if (MPFR_IS_NAN (z))
 107         {
 108           MPFR_SET_NAN (res);
 109           MPFR_RET_NAN;
 110         }
 111       /* j(n,z) tends to zero when z goes to +Inf or -Inf, oscillating around
 112          0. We choose to return +0 in that case. */
 113       else if (MPFR_IS_INF (z)) /* FIXME: according to j(-n,z) = (-1)^n j(n,z)
 114                                    we might want to give a sign depending on
 115                                    z and n */
 116         return mpfr_set_ui (res, 0, r);
 117       else /* z=0: j(0,0)=1, j(n odd,+/-0) = +/-0 if n > 0, -/+0 if n < 0,
 118               j(n even,+/-0) = +0 */
 119         {
 120           if (n == 0)
 121             return mpfr_set_ui (res, 1, r);
 122           else if (absn & 1) /* n odd */
 123             return (n > 0) ? mpfr_set (res, z, r) : mpfr_neg (res, z, r);
 124           else /* n even */
 125             return mpfr_set_ui (res, 0, r);
 126         }
 127     }
 128
 129   MPFR_SAVE_EXPO_MARK (expo);
 130
 131   /* check for tiny input for j0: j0(z) = 1 - z^2/4 + ..., more precisely
 132      |j0(z) - 1| <= z^2/4 for -1 <= z <= 1. */
 133   if (n == 0)
 134     MPFR_FAST_COMPUTE_IF_SMALL_INPUT (res, __gmpfr_one, -2 * MPFR_GET_EXP (z),
 135                                       2, 0, r, inex = _inexact; goto end);
 136
 137   /* idem for j1: j1(z) = z/2 - z^3/16 + ..., more precisely
 138      |j1(z) - z/2| <= |z^3|/16 for -1 <= z <= 1, with the sign of j1(z) - z/2
 139      being the opposite of that of z. */
 140   /* TODO: add a test to trigger an error when
 141        inex = _inexact; goto end
 142      is forgotten in MPFR_FAST_COMPUTE_IF_SMALL_INPUT below. */
 143   if (n == 1)
 144     {
 145       /* We first compute 2j1(z) = z - z^3/8 + ..., then divide by 2 using
 146          the "extra" argument of MPFR_FAST_COMPUTE_IF_SMALL_INPUT. But we
 147          must also handle the underflow case (an overflow is not possible
 148          for small inputs). If an underflow occurred in mpfr_round_near_x,
 149          the rounding was to zero or equivalent, and the result is 0, so
 150          that the division by 2 will give the wanted result. Otherwise...
 151          The rounded result in unbounded exponent range is res/2. If the
 152          division by 2 doesn't underflow, it is exact, and we can return
 153          this result. And an underflow in the division is a real underflow.
 154          In case of directed rounding mode, the result is correct. But in
 155          case of rounding to nearest, there is a double rounding problem,
 156          and the result is 0 iff the result before the division is the
 157          minimum positive number and _inexact has the same sign as z;
 158          but in rounding to nearest, res/2 will yield 0 iff |res| is the
 159          minimum positive number, so that we just need to test the result
 160          of the division and the sign of _inexact. */
 161       mpfr_clear_flags ();
 162       MPFR_FAST_COMPUTE_IF_SMALL_INPUT
 163         (res, z, -2 * MPFR_GET_EXP (z), 3, 0, r, {
 164           int inex2 = mpfr_div_2ui (res, res, 1, r);
 165           if (MPFR_UNLIKELY (r == MPFR_RNDN && MPFR_IS_ZERO (res)) &&
 166               (MPFR_ASSERTN (inex2 != 0), SIGN (_inexact) != MPFR_SIGN (z)))
 167             {
 168               mpfr_nexttoinf (res);
 169               inex = - inex2;
 170             }
 171           else
 172             inex = inex2 != 0 ? inex2 : _inexact;
 173           MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 174           goto end;
 175         });
 176     }
 177
 178   /* we can use the asymptotic expansion as soon as |z| > p log(2)/2,
 179      but to get some margin we use it for |z| > p/2 */
 180   pbound = MPFR_PREC (res) / 2 + 3;
 181   MPFR_ASSERTN (pbound <= ULONG_MAX);
 182   MPFR_ALIAS (absz, z, 1, MPFR_EXP (z));
 183   if (mpfr_cmp_ui (absz, pbound) > 0)
 184     {
 185       inex = mpfr_jn_asympt (res, n, z, r);
 186       if (inex != 0)
 187         goto end;
 188     }
 189
 190   MPFR_GROUP_INIT_3 (g, 32, y, s, t);
 191
 192   /* check underflow case: |j(n,z)| <= 1/sqrt(2 Pi n) (ze/2n)^n
 193      (see algorithms.tex) */
 194   /* FIXME: the code below doesn't detect all the underflow cases. Either
 195      this should be done, or the generic code should detect underflows. */
 196   if (absn > 0)
 197     {
 198       /* the following is an upper 32-bit approximation to exp(1)/2 */
 199       mpfr_set_str_binary (y, "1.0101101111110000101010001011001");
 200       if (MPFR_SIGN(z) > 0)
 201         mpfr_mul (y, y, z, MPFR_RNDU);
 202       else
 203         {
 204           mpfr_mul (y, y, z, MPFR_RNDD);
 205           mpfr_neg (y, y, MPFR_RNDU);
 206         }
 207       mpfr_div_ui (y, y, absn, MPFR_RNDU);
 208       /* now y is an upper approximation to |ze/2n|: y < 2^EXP(y),
 209          thus |j(n,z)| < 1/2*y^n < 2^(n*EXP(y)-1).
 210          If n*EXP(y) < emin then we have an underflow.
 211          Note that if emin = MPFR_EMIN_MIN and j = 1, this inequality
 212          will never be satisfied.
 213          Warning: absn is an unsigned long. */
 214       if ((MPFR_GET_EXP (y) < 0 && absn > - expo.saved_emin)
 215           || (absn <= - MPFR_EMIN_MIN &&
 216               MPFR_GET_EXP (y) < expo.saved_emin / (mpfr_exp_t) absn))
 217         {
 218           MPFR_GROUP_CLEAR (g);
 219           MPFR_SAVE_EXPO_FREE (expo);
 220           return mpfr_underflow (res, (r == MPFR_RNDN) ? MPFR_RNDZ : r,
 221                          (n % 2) ? ((n > 0) ? MPFR_SIGN(z) : -MPFR_SIGN(z))
 222                                  : MPFR_SIGN_POS);
 223         }
 224     }
 225
 226   /* the logarithm of the ratio between the largest term in the series
 227      and the first one is roughly bounded by k0, which we add to the
 228      working precision to take into account this cancellation */
 229   /* The following operations avoid integer overflow and ensure that
 230      prec <= MPFR_PREC_MAX (prec = MPFR_PREC_MAX won't prevent an abort,
 231      but the failure should be handled cleanly). */
 232   k0 = mpfr_jn_k0 (absn, z);
 233   MPFR_LOG_MSG (("k0 = %lu\n", k0));
 234   uprec = MPFR_PREC_MAX - 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC_MAX) - 3;
 235   if (k0 < uprec)
 236     uprec = k0;
 237   uprec += MPFR_PREC (res) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (res)) + 3;
 238   prec = uprec < MPFR_PREC_MAX ? (mpfr_prec_t) uprec : MPFR_PREC_MAX;
 239
 240   MPFR_ZIV_INIT (loop, prec);
 241   for (;;)
 242     {
 243       MPFR_BLOCK_DECL (flags);
 244
 245       MPFR_GROUP_REPREC_3 (g, prec, y, s, t);
 246       MPFR_BLOCK (flags, {
 247       mpfr_pow_ui (t, z, absn, MPFR_RNDN); /* z^|n| */
 248       mpfr_mul (y, z, z, MPFR_RNDN);       /* z^2 */
 249       mpfr_clear_erangeflag ();
 250       zz = mpfr_get_ui (y, MPFR_RNDU);
 251       /* FIXME: The error analysis is incorrect in case of range error. */
 252       MPFR_ASSERTN (! mpfr_erangeflag_p ()); /* since mpfr_clear_erangeflag */
 253       mpfr_div_2ui (y, y, 2, MPFR_RNDN);   /* z^2/4 */
 254       mpfr_fac_ui (s, absn, MPFR_RNDN);    /* |n|! */
 255       mpfr_div (t, t, s, MPFR_RNDN);
 256       if (absn > 0)
 257         mpfr_div_2ui (t, t, absn, MPFR_RNDN);
 258       mpfr_set (s, t, MPFR_RNDN);
 259       /* note: we assume here that the maximal error bound is proportional to
 260          2^exps, which is true also in the case where s=0 */
 261       exps = MPFR_IS_ZERO (s) ? MPFR_EMIN_MIN : MPFR_GET_EXP (s);
 262       expT = exps;
 263       for (k = 1; ; k++)
 264         {
 265           MPFR_LOG_MSG (("loop on k, k = %lu\n", k));
 266           mpfr_mul (t, t, y, MPFR_RNDN);
 267           mpfr_neg (t, t, MPFR_RNDN);
 268           /* Mathematically: absn <= LONG_MAX + 1 <= (ULONG_MAX + 1) / 2,
 269              and in practice, k is not very large, so that one should have
 270              k + absn <= ULONG_MAX. */
 271           MPFR_ASSERTN (absn <= ULONG_MAX - k);
 272           if (k + absn <= ULONG_MAX / k)
 273             mpfr_div_ui (t, t, k * (k + absn), MPFR_RNDN);
 274           else
 275             {
 276               mpfr_div_ui (t, t, k, MPFR_RNDN);
 277               mpfr_div_ui (t, t, k + absn, MPFR_RNDN);
 278             }
 279           /* see above note */
 280           exps = MPFR_IS_ZERO (s) ? MPFR_EMIN_MIN : MPFR_GET_EXP (t);
 281           if (exps > expT)
 282             expT = exps;
 283           mpfr_add (s, s, t, MPFR_RNDN);
 284           exps = MPFR_IS_ZERO (s) ? MPFR_EMIN_MIN : MPFR_GET_EXP (s);
 285           if (exps > expT)
 286             expT = exps;
 287           /* Above it has been checked that k + absn <= ULONG_MAX. */
 288           if (MPFR_GET_EXP (t) + (mpfr_exp_t) prec <= exps &&
 289               zz / (2 * k) < k + absn)
 290             break;
 291         }
 292       });
 293       /* the error is bounded by (4k^2+21/2k+7) ulp(s)*2^(expT-exps)
 294          <= (k+2)^2 ulp(s)*2^(2+expT-exps) */
 295       diffexp = expT - exps;
 296       err = 2 * MPFR_INT_CEIL_LOG2(k + 2) + 2;
 297       /* FIXME: Can an overflow occur in the following sum? */
 298       MPFR_ASSERTN (diffexp >= 0 && err >= 0 &&
 299                     diffexp <= MPFR_PREC_MAX - err);
 300       err += diffexp;
 301       if (MPFR_LIKELY (MPFR_CAN_ROUND (s, prec - err, MPFR_PREC(res), r)))
 302         {
 303           if (MPFR_LIKELY (! (MPFR_UNDERFLOW (flags) ||
 304                               MPFR_OVERFLOW (flags))))
 305             break;
 306           /* The error analysis is incorrect in case of exception.
 307              If an underflow or overflow occurred, try once more in
 308              a larger precision, and if this happens a second time,
 309              then abort to avoid a probable infinite loop. This is
 310              a problem that must be fixed! */
 311           MPFR_ASSERTN (! exception);
 312           exception = 1;
 313         }
 314       MPFR_ZIV_NEXT (loop, prec);
 315     }
 316   MPFR_ZIV_FREE (loop);
 317
 318   inex = ((n >= 0) || ((n & 1) == 0)) ? mpfr_set (res, s, r)
 319                                       : mpfr_neg (res, s, r);
 320
 321   MPFR_GROUP_CLEAR (g);
 322
 323  end:
 324   MPFR_SAVE_EXPO_FREE (expo);
 325   return mpfr_check_range (res, inex, r);
 326 }
 327
 328 #define MPFR_JN
 329 #include "jyn_asympt.c"