contrib/mpfr/src/yn.c

   1 /* mpfr_y0, mpfr_y1, mpfr_yn -- Bessel functions of 2nd kind, integer order.
   2    http://www.opengroup.org/onlinepubs/009695399/functions/y0.html
   3
   4 Copyright 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   5 Contributed by the Arenaire 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 static int mpfr_yn_asympt (mpfr_ptr, long, mpfr_srcptr, mpfr_rnd_t);
  28
  29 int
  30 mpfr_y0 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r)
  31 {
  32   return mpfr_yn (res, 0, z, r);
  33 }
  34
  35 int
  36 mpfr_y1 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r)
  37 {
  38   return mpfr_yn (res, 1, z, r);
  39 }
  40
  41 /* compute in s an approximation of S1 = sum((n-k)!/k!*y^k,k=0..n)
  42    return e >= 0 the exponent difference between the maximal value of |s|
  43    during the for loop and the final value of |s|.
  44 */
  45 static mpfr_exp_t
  46 mpfr_yn_s1 (mpfr_ptr s, mpfr_srcptr y, unsigned long n)
  47 {
  48   unsigned long k;
  49   mpz_t f;
  50   mpfr_exp_t e, emax;
  51
  52   mpz_init_set_ui (f, 1);
  53   /* we compute n!*S1 = sum(a[k]*y^k,k=0..n) where a[k] = n!*(n-k)!/k!,
  54      a[0] = (n!)^2, a[1] = n!*(n-1)!, ..., a[n-1] = n, a[n] = 1 */
  55   mpfr_set_ui (s, 1, MPFR_RNDN); /* a[n] */
  56   emax = MPFR_EXP(s);
  57   for (k = n; k-- > 0;)
  58     {
  59       /* a[k]/a[k+1] = (n-k)!/k!/(n-(k+1))!*(k+1)! = (k+1)*(n-k) */
  60       mpfr_mul (s, s, y, MPFR_RNDN);
  61       mpz_mul_ui (f, f, n - k);
  62       mpz_mul_ui (f, f, k + 1);
  63       /* invariant: f = a[k] */
  64       mpfr_add_z (s, s, f, MPFR_RNDN);
  65       e = MPFR_EXP(s);
  66       if (e > emax)
  67         emax = e;
  68     }
  69   /* now we have f = (n!)^2 */
  70   mpz_sqrt (f, f);
  71   mpfr_div_z (s, s, f, MPFR_RNDN);
  72   mpz_clear (f);
  73   return emax - MPFR_EXP(s);
  74 }
  75
  76 /* compute in s an approximation of
  77    S3 = c*sum((h(k)+h(n+k))*y^k/k!/(n+k)!,k=0..infinity)
  78    where h(k) = 1 + 1/2 + ... + 1/k
  79    k=0: h(n)
  80    k=1: 1+h(n+1)
  81    k=2: 3/2+h(n+2)
  82    Returns e such that the error is bounded by 2^e ulp(s).
  83 */
  84 static mpfr_exp_t
  85 mpfr_yn_s3 (mpfr_ptr s, mpfr_srcptr y, mpfr_srcptr c, unsigned long n)
  86 {
  87   unsigned long k, zz;
  88   mpfr_t t, u;
  89   mpz_t p, q; /* p/q will store h(k)+h(n+k) */
  90   mpfr_exp_t exps, expU;
  91
  92   zz = mpfr_get_ui (y, MPFR_RNDU); /* y = z^2/4 */
  93   MPFR_ASSERTN (zz < ULONG_MAX - 2);
  94   zz += 2; /* z^2 <= 2^zz */
  95   mpz_init_set_ui (p, 0);
  96   mpz_init_set_ui (q, 1);
  97   /* initialize p/q to h(n) */
  98   for (k = 1; k <= n; k++)
  99     {
 100       /* p/q + 1/k = (k*p+q)/(q*k) */
 101       mpz_mul_ui (p, p, k);
 102       mpz_add (p, p, q);
 103       mpz_mul_ui (q, q, k);
 104     }
 105   mpfr_init2 (t, MPFR_PREC(s));
 106   mpfr_init2 (u, MPFR_PREC(s));
 107   mpfr_fac_ui (t, n, MPFR_RNDN);
 108   mpfr_div (t, c, t, MPFR_RNDN);    /* c/n! */
 109   mpfr_mul_z (u, t, p, MPFR_RNDN);
 110   mpfr_div_z (s, u, q, MPFR_RNDN);
 111   exps = MPFR_EXP (s);
 112   expU = exps;
 113   for (k = 1; ;k ++)
 114     {
 115       /* update t */
 116       mpfr_mul (t, t, y, MPFR_RNDN);
 117       mpfr_div_ui (t, t, k, MPFR_RNDN);
 118       mpfr_div_ui (t, t, n + k, MPFR_RNDN);
 119       /* update p/q:
 120          p/q + 1/k + 1/(n+k) = [p*k*(n+k) + q*(n+k) + q*k]/(q*k*(n+k)) */
 121       mpz_mul_ui (p, p, k);
 122       mpz_mul_ui (p, p, n + k);
 123       mpz_addmul_ui (p, q, n + 2 * k);
 124       mpz_mul_ui (q, q, k);
 125       mpz_mul_ui (q, q, n + k);
 126       mpfr_mul_z (u, t, p, MPFR_RNDN);
 127       mpfr_div_z (u, u, q, MPFR_RNDN);
 128       exps = MPFR_EXP (u);
 129       if (exps > expU)
 130         expU = exps;
 131       mpfr_add (s, s, u, MPFR_RNDN);
 132       exps = MPFR_EXP (s);
 133       if (exps > expU)
 134         expU = exps;
 135       if (MPFR_EXP (u) + (mpfr_exp_t) MPFR_PREC (u) < MPFR_EXP (s) &&
 136           zz / (2 * k) < k + n)
 137         break;
 138     }
 139   mpfr_clear (t);
 140   mpfr_clear (u);
 141   mpz_clear (p);
 142   mpz_clear (q);
 143   exps = expU - MPFR_EXP (s);
 144   /* the error is bounded by (6k^2+33/2k+11) 2^exps ulps
 145      <= 8*(k+2)^2 2^exps ulps */
 146   return 3 + 2 * MPFR_INT_CEIL_LOG2(k + 2) + exps;
 147 }
 148
 149 int
 150 mpfr_yn (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r)
 151 {
 152   int inex;
 153   unsigned long absn;
 154   MPFR_SAVE_EXPO_DECL (expo);
 155
 156   MPFR_LOG_FUNC
 157     (("n=%ld x[%Pu]=%.*Rg rnd=%d", n, mpfr_get_prec (z), mpfr_log_prec, z, r),
 158      ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (res), mpfr_log_prec, res, inex));
 159
 160   absn = SAFE_ABS (unsigned long, n);
 161
 162   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (z)))
 163     {
 164       if (MPFR_IS_NAN (z))
 165         {
 166           MPFR_SET_NAN (res); /* y(n,NaN) = NaN */
 167           MPFR_RET_NAN;
 168         }
 169       /* y(n,z) tends to zero when z goes to +Inf, oscillating around
 170          0. We choose to return +0 in that case. */
 171       else if (MPFR_IS_INF (z))
 172         {
 173           if (MPFR_SIGN(z) > 0)
 174             return mpfr_set_ui (res, 0, r);
 175           else /* y(n,-Inf) = NaN */
 176             {
 177               MPFR_SET_NAN (res);
 178               MPFR_RET_NAN;
 179             }
 180         }
 181       else /* y(n,z) tends to -Inf for n >= 0 or n even, to +Inf otherwise,
 182               when z goes to zero */
 183         {
 184           MPFR_SET_INF(res);
 185           if (n >= 0 || ((unsigned long) n & 1) == 0)
 186             MPFR_SET_NEG(res);
 187           else
 188             MPFR_SET_POS(res);
 189           mpfr_set_divby0 ();
 190           MPFR_RET(0);
 191         }
 192     }
 193
 194   /* for z < 0, y(n,z) is imaginary except when j(n,|z|) = 0, which we
 195      assume does not happen for a rational z. */
 196   if (MPFR_SIGN(z) < 0)
 197     {
 198       MPFR_SET_NAN (res);
 199       MPFR_RET_NAN;
 200     }
 201
 202   /* now z is not singular, and z > 0 */
 203
 204   MPFR_SAVE_EXPO_MARK (expo);
 205
 206   /* Deal with tiny arguments. We have:
 207      y0(z) = 2 log(z)/Pi + 2 (euler - log(2))/Pi + O(log(z)*z^2), more
 208      precisely for 0 <= z <= 1/2, with g(z) = 2/Pi + 2(euler-log(2))/Pi/log(z),
 209                 g(z) - 0.41*z^2 < y0(z)/log(z) < g(z)
 210      thus since log(z) is negative:
 211              g(z)*log(z) < y0(z) < (g(z) - z^2/2)*log(z)
 212      and since |g(z)| >= 0.63 for 0 <= z <= 1/2, the relative error on
 213      y0(z)/log(z) is bounded by 0.41*z^2/0.63 <= 0.66*z^2.
 214      Note: we use both the main term in log(z) and the constant term, because
 215      otherwise the relative error would be only in 1/log(|log(z)|).
 216   */
 217   if (n == 0 && MPFR_EXP(z) < - (mpfr_exp_t) (MPFR_PREC(res) / 2))
 218     {
 219       mpfr_t l, h, t, logz;
 220       mpfr_prec_t prec;
 221       int ok, inex2;
 222
 223       prec = MPFR_PREC(res) + 10;
 224       mpfr_init2 (l, prec);
 225       mpfr_init2 (h, prec);
 226       mpfr_init2 (t, prec);
 227       mpfr_init2 (logz, prec);
 228       /* first enclose log(z) + euler - log(2) = log(z/2) + euler */
 229       mpfr_log (logz, z, MPFR_RNDD);    /* lower bound of log(z) */
 230       mpfr_set (h, logz, MPFR_RNDU);    /* exact */
 231       mpfr_nextabove (h);              /* upper bound of log(z) */
 232       mpfr_const_euler (t, MPFR_RNDD);  /* lower bound of euler */
 233       mpfr_add (l, logz, t, MPFR_RNDD); /* lower bound of log(z) + euler */
 234       mpfr_nextabove (t);              /* upper bound of euler */
 235       mpfr_add (h, h, t, MPFR_RNDU);    /* upper bound of log(z) + euler */
 236       mpfr_const_log2 (t, MPFR_RNDU);   /* upper bound of log(2) */
 237       mpfr_sub (l, l, t, MPFR_RNDD);    /* lower bound of log(z/2) + euler */
 238       mpfr_nextbelow (t);              /* lower bound of log(2) */
 239       mpfr_sub (h, h, t, MPFR_RNDU);    /* upper bound of log(z/2) + euler */
 240       mpfr_const_pi (t, MPFR_RNDU);     /* upper bound of Pi */
 241       mpfr_div (l, l, t, MPFR_RNDD);    /* lower bound of (log(z/2)+euler)/Pi */
 242       mpfr_nextbelow (t);              /* lower bound of Pi */
 243       mpfr_div (h, h, t, MPFR_RNDD);    /* upper bound of (log(z/2)+euler)/Pi */
 244       mpfr_mul_2ui (l, l, 1, MPFR_RNDD); /* lower bound on g(z)*log(z) */
 245       mpfr_mul_2ui (h, h, 1, MPFR_RNDU); /* upper bound on g(z)*log(z) */
 246       /* we now have l <= g(z)*log(z) <= h, and we need to add -z^2/2*log(z)
 247          to h */
 248       mpfr_mul (t, z, z, MPFR_RNDU);     /* upper bound on z^2 */
 249       /* since logz is negative, a lower bound corresponds to an upper bound
 250          for its absolute value */
 251       mpfr_neg (t, t, MPFR_RNDD);
 252       mpfr_div_2ui (t, t, 1, MPFR_RNDD);
 253       mpfr_mul (t, t, logz, MPFR_RNDU); /* upper bound on z^2/2*log(z) */
 254       mpfr_add (h, h, t, MPFR_RNDU);
 255       inex = mpfr_prec_round (l, MPFR_PREC(res), r);
 256       inex2 = mpfr_prec_round (h, MPFR_PREC(res), r);
 257       /* we need h=l and inex=inex2 */
 258       ok = (inex == inex2) && mpfr_equal_p (l, h);
 259       if (ok)
 260         mpfr_set (res, h, r); /* exact */
 261       mpfr_clear (l);
 262       mpfr_clear (h);
 263       mpfr_clear (t);
 264       mpfr_clear (logz);
 265       if (ok)
 266         goto end;
 267     }
 268
 269   /* small argument check for y1(z) = -2/Pi/z + O(log(z)):
 270      for 0 <= z <= 1, |y1(z) + 2/Pi/z| <= 0.25 */
 271   if (n == 1 && MPFR_EXP(z) + 1 < - (mpfr_exp_t) MPFR_PREC(res))
 272     {
 273       mpfr_t y;
 274       mpfr_prec_t prec;
 275       mpfr_exp_t err1;
 276       int ok;
 277       MPFR_BLOCK_DECL (flags);
 278
 279       /* since 2/Pi > 0.5, and |y1(z)| >= |2/Pi/z|, if z <= 2^(-emax-1),
 280          then |y1(z)| > 2^emax */
 281       prec = MPFR_PREC(res) + 10;
 282       mpfr_init2 (y, prec);
 283       mpfr_const_pi (y, MPFR_RNDU); /* Pi*(1+u)^2, where here and below u
 284                                       represents a quantity <= 1/2^prec */
 285       mpfr_mul (y, y, z, MPFR_RNDU); /* Pi*z * (1+u)^4, upper bound */
 286       MPFR_BLOCK (flags, mpfr_ui_div (y, 2, y, MPFR_RNDZ));
 287       /* 2/Pi/z * (1+u)^6, lower bound, with possible overflow */
 288       if (MPFR_OVERFLOW (flags))
 289         {
 290           mpfr_clear (y);
 291           MPFR_SAVE_EXPO_FREE (expo);
 292           return mpfr_overflow (res, r, -1);
 293         }
 294       mpfr_neg (y, y, MPFR_RNDN);
 295       /* (1+u)^6 can be written 1+7u [for another value of u], thus the
 296          error on 2/Pi/z is less than 7ulp(y). The truncation error is less
 297          than 1/4, thus if ulp(y)>=1/4, the total error is less than 8ulp(y),
 298          otherwise it is less than 1/4+7/8 <= 2. */
 299       if (MPFR_EXP(y) + 2 >= MPFR_PREC(y)) /* ulp(y) >= 1/4 */
 300         err1 = 3;
 301       else /* ulp(y) <= 1/8 */
 302         err1 = (mpfr_exp_t) MPFR_PREC(y) - MPFR_EXP(y) + 1;
 303       ok = MPFR_CAN_ROUND (y, prec - err1, MPFR_PREC(res), r);
 304       if (ok)
 305         inex = mpfr_set (res, y, r);
 306       mpfr_clear (y);
 307       if (ok)
 308         goto end;
 309     }
 310
 311   /* we can use the asymptotic expansion as soon as z > p log(2)/2,
 312      but to get some margin we use it for z > p/2 */
 313   if (mpfr_cmp_ui (z, MPFR_PREC(res) / 2 + 3) > 0)
 314     {
 315       inex = mpfr_yn_asympt (res, n, z, r);
 316       if (inex != 0)
 317         goto end;
 318     }
 319
 320   /* General case */
 321   {
 322     mpfr_prec_t prec;
 323     mpfr_exp_t err1, err2, err3;
 324     mpfr_t y, s1, s2, s3;
 325     MPFR_ZIV_DECL (loop);
 326
 327     mpfr_init (y);
 328     mpfr_init (s1);
 329     mpfr_init (s2);
 330     mpfr_init (s3);
 331
 332     prec = MPFR_PREC(res) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (res)) + 13;
 333     MPFR_ZIV_INIT (loop, prec);
 334     for (;;)
 335       {
 336         mpfr_set_prec (y, prec);
 337         mpfr_set_prec (s1, prec);
 338         mpfr_set_prec (s2, prec);
 339         mpfr_set_prec (s3, prec);
 340
 341         mpfr_mul (y, z, z, MPFR_RNDN);
 342         mpfr_div_2ui (y, y, 2, MPFR_RNDN); /* z^2/4 */
 343
 344         /* store (z/2)^n temporarily in s2 */
 345         mpfr_pow_ui (s2, z, absn, MPFR_RNDN);
 346         mpfr_div_2si (s2, s2, absn, MPFR_RNDN);
 347
 348         /* compute S1 * (z/2)^(-n) */
 349         if (n == 0)
 350           {
 351             mpfr_set_ui (s1, 0, MPFR_RNDN);
 352             err1 = 0;
 353           }
 354         else
 355           err1 = mpfr_yn_s1 (s1, y, absn - 1);
 356         mpfr_div (s1, s1, s2, MPFR_RNDN); /* (z/2)^(-n) * S1 */
 357         /* See algorithms.tex: the relative error on s1 is bounded by
 358            (3n+3)*2^(e+1-prec). */
 359         err1 = MPFR_INT_CEIL_LOG2 (3 * absn + 3) + err1 + 1;
 360         /* rel_err(s1) <= 2^(err1-prec), thus err(s1) <= 2^err1 ulps */
 361
 362         /* compute (z/2)^n * S3 */
 363         mpfr_neg (y, y, MPFR_RNDN); /* -z^2/4 */
 364         err3 = mpfr_yn_s3 (s3, y, s2, absn); /* (z/2)^n * S3 */
 365         /* the error on s3 is bounded by 2^err3 ulps */
 366
 367         /* add s1+s3 */
 368         err1 += MPFR_EXP(s1);
 369         mpfr_add (s1, s1, s3, MPFR_RNDN);
 370         /* the error is bounded by 1/2 + 2^err1*2^(- EXP(s1))
 371            + 2^err3*2^(EXP(s3) - EXP(s1)) */
 372         err3 += MPFR_EXP(s3);
 373         err1 = (err3 > err1) ? err3 + 1 : err1 + 1;
 374         err1 -= MPFR_EXP(s1);
 375         err1 = (err1 >= 0) ? err1 + 1 : 1;
 376         /* now the error on s1 is bounded by 2^err1*ulp(s1) */
 377
 378         /* compute S2 */
 379         mpfr_div_2ui (s2, z, 1, MPFR_RNDN); /* z/2 */
 380         mpfr_log (s2, s2, MPFR_RNDN); /* log(z/2) */
 381         mpfr_const_euler (s3, MPFR_RNDN);
 382         err2 = MPFR_EXP(s2) > MPFR_EXP(s3) ? MPFR_EXP(s2) : MPFR_EXP(s3);
 383         mpfr_add (s2, s2, s3, MPFR_RNDN); /* log(z/2) + gamma */
 384         err2 -= MPFR_EXP(s2);
 385         mpfr_mul_2ui (s2, s2, 1, MPFR_RNDN); /* 2*(log(z/2) + gamma) */
 386         mpfr_jn (s3, absn, z, MPFR_RNDN); /* Jn(z) */
 387         mpfr_mul (s2, s2, s3, MPFR_RNDN); /* 2*(log(z/2) + gamma)*Jn(z) */
 388         err2 += 4; /* the error on s2 is bounded by 2^err2 ulps, see
 389                       algorithms.tex */
 390
 391         /* add all three sums */
 392         err1 += MPFR_EXP(s1); /* the error on s1 is bounded by 2^err1 */
 393         err2 += MPFR_EXP(s2); /* the error on s2 is bounded by 2^err2 */
 394         mpfr_sub (s2, s2, s1, MPFR_RNDN); /* s2 - (s1+s3) */
 395         err2 = (err1 > err2) ? err1 + 1 : err2 + 1;
 396         err2 -= MPFR_EXP(s2);
 397         err2 = (err2 >= 0) ? err2 + 1 : 1;
 398         /* now the error on s2 is bounded by 2^err2*ulp(s2) */
 399         mpfr_const_pi (y, MPFR_RNDN); /* error bounded by 1 ulp */
 400         mpfr_div (s2, s2, y, MPFR_RNDN); /* error bounded by
 401                                            2^(err2+1)*ulp(s2) */
 402         err2 ++;
 403
 404         if (MPFR_LIKELY (MPFR_CAN_ROUND (s2, prec - err2, MPFR_PREC(res), r)))
 405           break;
 406         MPFR_ZIV_NEXT (loop, prec);
 407       }
 408     MPFR_ZIV_FREE (loop);
 409
 410     /* Assume two's complement for the test n & 1 */
 411     inex = mpfr_set4 (res, s2, r, n >= 0 || (n & 1) == 0 ?
 412                       MPFR_SIGN (s2) : - MPFR_SIGN (s2));
 413
 414     mpfr_clear (y);
 415     mpfr_clear (s1);
 416     mpfr_clear (s2);
 417     mpfr_clear (s3);
 418   }
 419
 420  end:
 421   MPFR_SAVE_EXPO_FREE (expo);
 422   return mpfr_check_range (res, inex, r);
 423 }
 424
 425 #define MPFR_YN
 426 #include "jyn_asympt.c"