contrib/mpfr/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 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 static int mpfr_yn_asympt (mpfr_ptr, long, mpfr_srcptr, mp_rnd_t);
  28
  29 int
  30 mpfr_y0 (mpfr_ptr res, mpfr_srcptr z, mp_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, mp_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 mp_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   mp_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, GMP_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, GMP_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, GMP_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, GMP_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 mp_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   mp_exp_t exps, expU;
  91
  92   zz = mpfr_get_ui (y, GMP_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, GMP_RNDN);
 108   mpfr_div (t, c, t, GMP_RNDN);    /* c/n! */
 109   mpfr_mul_z (u, t, p, GMP_RNDN);
 110   mpfr_div_z (s, u, q, GMP_RNDN);
 111   exps = MPFR_EXP (s);
 112   expU = exps;
 113   for (k = 1; ;k ++)
 114     {
 115       /* update t */
 116       mpfr_mul (t, t, y, GMP_RNDN);
 117       mpfr_div_ui (t, t, k, GMP_RNDN);
 118       mpfr_div_ui (t, t, n + k, GMP_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, GMP_RNDN);
 127       mpfr_div_z (u, u, q, GMP_RNDN);
 128       exps = MPFR_EXP (u);
 129       if (exps > expU)
 130         expU = exps;
 131       mpfr_add (s, s, u, GMP_RNDN);
 132       exps = MPFR_EXP (s);
 133       if (exps > expU)
 134         expU = exps;
 135       if (MPFR_EXP (u) + (mp_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, mp_rnd_t r)
 151 {
 152   int inex;
 153   unsigned long absn;
 154
 155   MPFR_LOG_FUNC (("x[%#R]=%R n=%d rnd=%d", z, z, n, r),
 156                  ("y[%#R]=%R", res, res));
 157
 158   absn = SAFE_ABS (unsigned long, n);
 159
 160   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (z)))
 161     {
 162       if (MPFR_IS_NAN (z))
 163         {
 164           MPFR_SET_NAN (res); /* y(n,NaN) = NaN */
 165           MPFR_RET_NAN;
 166         }
 167       /* y(n,z) tends to zero when z goes to +Inf, oscillating around
 168          0. We choose to return +0 in that case. */
 169       else if (MPFR_IS_INF (z))
 170         {
 171           if (MPFR_SIGN(z) > 0)
 172             return mpfr_set_ui (res, 0, r);
 173           else /* y(n,-Inf) = NaN */
 174             {
 175               MPFR_SET_NAN (res);
 176               MPFR_RET_NAN;
 177             }
 178         }
 179       else /* y(n,z) tends to -Inf for n >= 0 or n even, to +Inf otherwise,
 180               when z goes to zero */
 181         {
 182           MPFR_SET_INF(res);
 183           if (n >= 0 || (n & 1) == 0)
 184             MPFR_SET_NEG(res);
 185           else
 186             MPFR_SET_POS(res);
 187           MPFR_RET(0);
 188         }
 189     }
 190
 191   /* for z < 0, y(n,z) is imaginary except when j(n,|z|) = 0, which we
 192      assume does not happen for a rational z. */
 193   if (MPFR_SIGN(z) < 0)
 194     {
 195       MPFR_SET_NAN (res);
 196       MPFR_RET_NAN;
 197     }
 198
 199   /* now z is not singular, and z > 0 */
 200
 201   /* Deal with tiny arguments. We have:
 202      y0(z) = 2 log(z)/Pi + 2 (euler - log(2))/Pi + O(log(z)*z^2), more
 203      precisely for 0 <= z <= 1/2, with g(z) = 2/Pi + 2(euler-log(2))/Pi/log(z),
 204                 g(z) - 0.41*z^2 < y0(z)/log(z) < g(z)
 205      thus since log(z) is negative:
 206              g(z)*log(z) < y0(z) < (g(z) - z^2/2)*log(z)
 207      and since |g(z)| >= 0.63 for 0 <= z <= 1/2, the relative error on
 208      y0(z)/log(z) is bounded by 0.41*z^2/0.63 <= 0.66*z^2.
 209      Note: we use both the main term in log(z) and the constant term, because
 210      otherwise the relative error would be only in 1/log(|log(z)|).
 211   */
 212   if (n == 0 && MPFR_EXP(z) < - (mp_exp_t) (MPFR_PREC(res) / 2))
 213     {
 214       mpfr_t l, h, t, logz;
 215       mp_prec_t prec;
 216       int ok, inex2;
 217
 218       prec = MPFR_PREC(res) + 10;
 219       mpfr_init2 (l, prec);
 220       mpfr_init2 (h, prec);
 221       mpfr_init2 (t, prec);
 222       mpfr_init2 (logz, prec);
 223       /* first enclose log(z) + euler - log(2) = log(z/2) + euler */
 224       mpfr_log (logz, z, GMP_RNDD);    /* lower bound of log(z) */
 225       mpfr_set (h, logz, GMP_RNDU);    /* exact */
 226       mpfr_nextabove (h);              /* upper bound of log(z) */
 227       mpfr_const_euler (t, GMP_RNDD);  /* lower bound of euler */
 228       mpfr_add (l, logz, t, GMP_RNDD); /* lower bound of log(z) + euler */
 229       mpfr_nextabove (t);              /* upper bound of euler */
 230       mpfr_add (h, h, t, GMP_RNDU);    /* upper bound of log(z) + euler */
 231       mpfr_const_log2 (t, GMP_RNDU);   /* upper bound of log(2) */
 232       mpfr_sub (l, l, t, GMP_RNDD);    /* lower bound of log(z/2) + euler */
 233       mpfr_nextbelow (t);              /* lower bound of log(2) */
 234       mpfr_sub (h, h, t, GMP_RNDU);    /* upper bound of log(z/2) + euler */
 235       mpfr_const_pi (t, GMP_RNDU);     /* upper bound of Pi */
 236       mpfr_div (l, l, t, GMP_RNDD);    /* lower bound of (log(z/2)+euler)/Pi */
 237       mpfr_nextbelow (t);              /* lower bound of Pi */
 238       mpfr_div (h, h, t, GMP_RNDD);    /* upper bound of (log(z/2)+euler)/Pi */
 239       mpfr_mul_2ui (l, l, 1, GMP_RNDD); /* lower bound on g(z)*log(z) */
 240       mpfr_mul_2ui (h, h, 1, GMP_RNDU); /* upper bound on g(z)*log(z) */
 241       /* we now have l <= g(z)*log(z) <= h, and we need to add -z^2/2*log(z)
 242          to h */
 243       mpfr_mul (t, z, z, GMP_RNDU);     /* upper bound on z^2 */
 244       /* since logz is negative, a lower bound corresponds to an upper bound
 245          for its absolute value */
 246       mpfr_neg (t, t, GMP_RNDD);
 247       mpfr_div_2ui (t, t, 1, GMP_RNDD);
 248       mpfr_mul (t, t, logz, GMP_RNDU); /* upper bound on z^2/2*log(z) */
 249       /* an underflow may happen in the above instructions, clear flag */
 250       mpfr_clear_underflow ();
 251       mpfr_add (h, h, t, GMP_RNDU);
 252       inex = mpfr_prec_round (l, MPFR_PREC(res), r);
 253       inex2 = mpfr_prec_round (h, MPFR_PREC(res), r);
 254       /* we need h=l and inex=inex2 */
 255       ok = (inex == inex2) && (mpfr_cmp (l, h) == 0);
 256       if (ok)
 257         mpfr_set (res, h, r); /* exact */
 258       mpfr_clear (l);
 259       mpfr_clear (h);
 260       mpfr_clear (t);
 261       mpfr_clear (logz);
 262       if (ok)
 263         return inex;
 264     }
 265
 266   /* small argument check for y1(z) = -2/Pi/z + O(log(z)):
 267      for 0 <= z <= 1, |y1(z) + 2/Pi/z| <= 0.25 */
 268   if (n == 1 && MPFR_EXP(z) + 1 < - (mp_exp_t) MPFR_PREC(res))
 269     {
 270       mpfr_t y;
 271       mp_prec_t prec;
 272       mp_exp_t err1;
 273       int ok;
 274       MPFR_BLOCK_DECL (flags);
 275
 276       /* since 2/Pi > 0.5, and |y1(z)| >= |2/Pi/z|, if z <= 2^(-emax-1),
 277          then |y1(z)| > 2^emax */
 278       prec = MPFR_PREC(res) + 10;
 279       mpfr_init2 (y, prec);
 280       mpfr_const_pi (y, GMP_RNDU); /* Pi*(1+u)^2, where here and below u
 281                                       represents a quantity <= 1/2^prec */
 282       mpfr_mul (y, y, z, GMP_RNDU); /* Pi*z * (1+u)^4, upper bound */
 283       MPFR_BLOCK (flags, mpfr_ui_div (y, 2, y, GMP_RNDZ));
 284       /* 2/Pi/z * (1+u)^6, lower bound, with possible overflow */
 285       if (MPFR_OVERFLOW (flags))
 286         {
 287           mpfr_clear (y);
 288           return mpfr_overflow (res, r, -1);
 289         }
 290       mpfr_neg (y, y, GMP_RNDN);
 291       /* (1+u)^6 can be written 1+7u [for another value of u], thus the
 292          error on 2/Pi/z is less than 7ulp(y). The truncation error is less
 293          than 1/4, thus if ulp(y)>=1/4, the total error is less than 8ulp(y),
 294          otherwise it is less than 1/4+7/8 <= 2. */
 295       if (MPFR_EXP(y) + 2 >= MPFR_PREC(y)) /* ulp(y) >= 1/4 */
 296         err1 = 3;
 297       else /* ulp(y) <= 1/8 */
 298         err1 = (mp_exp_t) MPFR_PREC(y) - MPFR_EXP(y) + 1;
 299       ok = MPFR_CAN_ROUND (y, prec - err1, MPFR_PREC(res), r);
 300       if (ok)
 301         inex = mpfr_set (res, y, r);
 302       mpfr_clear (y);
 303       if (ok)
 304         return inex;
 305     }
 306
 307   /* we can use the asymptotic expansion as soon as z > p log(2)/2,
 308      but to get some margin we use it for z > p/2 */
 309   if (mpfr_cmp_ui (z, MPFR_PREC(res) / 2 + 3) > 0)
 310     {
 311       inex = mpfr_yn_asympt (res, n, z, r);
 312       if (inex != 0)
 313         return inex;
 314     }
 315
 316   /* General case */
 317   {
 318     mp_prec_t prec;
 319     mp_exp_t err1, err2, err3;
 320     mpfr_t y, s1, s2, s3;
 321     MPFR_ZIV_DECL (loop);
 322
 323     mpfr_init (y);
 324     mpfr_init (s1);
 325     mpfr_init (s2);
 326     mpfr_init (s3);
 327
 328     prec = MPFR_PREC(res) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (res)) + 13;
 329     MPFR_ZIV_INIT (loop, prec);
 330     for (;;)
 331       {
 332         mpfr_set_prec (y, prec);
 333         mpfr_set_prec (s1, prec);
 334         mpfr_set_prec (s2, prec);
 335         mpfr_set_prec (s3, prec);
 336
 337         mpfr_mul (y, z, z, GMP_RNDN);
 338         mpfr_div_2ui (y, y, 2, GMP_RNDN); /* z^2/4 */
 339
 340         /* store (z/2)^n temporarily in s2 */
 341         mpfr_pow_ui (s2, z, absn, GMP_RNDN);
 342         mpfr_div_2si (s2, s2, absn, GMP_RNDN);
 343
 344         /* compute S1 * (z/2)^(-n) */
 345         if (n == 0)
 346           {
 347             mpfr_set_ui (s1, 0, GMP_RNDN);
 348             err1 = 0;
 349           }
 350         else
 351           err1 = mpfr_yn_s1 (s1, y, absn - 1);
 352         mpfr_div (s1, s1, s2, GMP_RNDN); /* (z/2)^(-n) * S1 */
 353         /* See algorithms.tex: the relative error on s1 is bounded by
 354            (3n+3)*2^(e+1-prec). */
 355         err1 = MPFR_INT_CEIL_LOG2 (3 * absn + 3) + err1 + 1;
 356         /* rel_err(s1) <= 2^(err1-prec), thus err(s1) <= 2^err1 ulps */
 357
 358         /* compute (z/2)^n * S3 */
 359         mpfr_neg (y, y, GMP_RNDN); /* -z^2/4 */
 360         err3 = mpfr_yn_s3 (s3, y, s2, absn); /* (z/2)^n * S3 */
 361         /* the error on s3 is bounded by 2^err3 ulps */
 362
 363         /* add s1+s3 */
 364         err1 += MPFR_EXP(s1);
 365         mpfr_add (s1, s1, s3, GMP_RNDN);
 366         /* the error is bounded by 1/2 + 2^err1*2^(- EXP(s1))
 367            + 2^err3*2^(EXP(s3) - EXP(s1)) */
 368         err3 += MPFR_EXP(s3);
 369         err1 = (err3 > err1) ? err3 + 1 : err1 + 1;
 370         err1 -= MPFR_EXP(s1);
 371         err1 = (err1 >= 0) ? err1 + 1 : 1;
 372         /* now the error on s1 is bounded by 2^err1*ulp(s1) */
 373
 374         /* compute S2 */
 375         mpfr_div_2ui (s2, z, 1, GMP_RNDN); /* z/2 */
 376         mpfr_log (s2, s2, GMP_RNDN); /* log(z/2) */
 377         mpfr_const_euler (s3, GMP_RNDN);
 378         err2 = MPFR_EXP(s2) > MPFR_EXP(s3) ? MPFR_EXP(s2) : MPFR_EXP(s3);
 379         mpfr_add (s2, s2, s3, GMP_RNDN); /* log(z/2) + gamma */
 380         err2 -= MPFR_EXP(s2);
 381         mpfr_mul_2ui (s2, s2, 1, GMP_RNDN); /* 2*(log(z/2) + gamma) */
 382         mpfr_jn (s3, absn, z, GMP_RNDN); /* Jn(z) */
 383         mpfr_mul (s2, s2, s3, GMP_RNDN); /* 2*(log(z/2) + gamma)*Jn(z) */
 384         err2 += 4; /* the error on s2 is bounded by 2^err2 ulps, see
 385                       algorithms.tex */
 386
 387         /* add all three sums */
 388         err1 += MPFR_EXP(s1); /* the error on s1 is bounded by 2^err1 */
 389         err2 += MPFR_EXP(s2); /* the error on s2 is bounded by 2^err2 */
 390         mpfr_sub (s2, s2, s1, GMP_RNDN); /* s2 - (s1+s3) */
 391         err2 = (err1 > err2) ? err1 + 1 : err2 + 1;
 392         err2 -= MPFR_EXP(s2);
 393         err2 = (err2 >= 0) ? err2 + 1 : 1;
 394         /* now the error on s2 is bounded by 2^err2*ulp(s2) */
 395         mpfr_const_pi (y, GMP_RNDN); /* error bounded by 1 ulp */
 396         mpfr_div (s2, s2, y, GMP_RNDN); /* error bounded by
 397                                            2^(err2+1)*ulp(s2) */
 398         err2 ++;
 399
 400         if (MPFR_LIKELY (MPFR_CAN_ROUND (s2, prec - err2, MPFR_PREC(res), r)))
 401           break;
 402         MPFR_ZIV_NEXT (loop, prec);
 403       }
 404     MPFR_ZIV_FREE (loop);
 405
 406     inex = (n >= 0 || (n & 1) == 0)
 407       ? mpfr_set (res, s2, r)
 408       : mpfr_neg (res, s2, r);
 409
 410     mpfr_clear (y);
 411     mpfr_clear (s1);
 412     mpfr_clear (s2);
 413     mpfr_clear (s3);
 414   }
 415
 416   return inex;
 417 }
 418
 419 #define MPFR_YN
 420 #include "jyn_asympt.c"