contrib/mpfr/jyn_asympt.c

   1 /* mpfr_jn_asympt, mpfr_yn_asympt -- shared code for mpfr_jn and mpfr_yn
   2
   3 Copyright 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 2.1 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #ifdef MPFR_JN
  24 # define FUNCTION mpfr_jn_asympt
  25 #else
  26 # ifdef MPFR_YN
  27 #  define FUNCTION mpfr_yn_asympt
  28 # else
  29 #  error "neither MPFR_JN nor MPFR_YN is defined"
  30 # endif
  31 #endif
  32
  33 /* Implements asymptotic expansion for jn or yn (formulae 9.2.5 and 9.2.6
  34    from Abramowitz & Stegun).
  35    Assumes |z| > p log(2)/2, where p is the target precision
  36    (z can be negative only for jn).
  37    Return 0 if the expansion does not converge enough (the value 0 as inexact
  38    flag should not happen for normal input).
  39 */
  40 static int
  41 FUNCTION (mpfr_ptr res, long n, mpfr_srcptr z, mp_rnd_t r)
  42 {
  43   mpfr_t s, c, P, Q, t, iz, err_t, err_s, err_u;
  44   mp_prec_t w;
  45   long k;
  46   int inex, stop, diverge = 0;
  47   mp_exp_t err2, err;
  48   MPFR_ZIV_DECL (loop);
  49
  50   mpfr_init (c);
  51
  52   w = MPFR_PREC(res) + MPFR_INT_CEIL_LOG2(MPFR_PREC(res)) + 4;
  53
  54   MPFR_ZIV_INIT (loop, w);
  55   for (;;)
  56     {
  57       mpfr_set_prec (c, w);
  58       mpfr_init2 (s, w);
  59       mpfr_init2 (P, w);
  60       mpfr_init2 (Q, w);
  61       mpfr_init2 (t, w);
  62       mpfr_init2 (iz, w);
  63       mpfr_init2 (err_t, 31);
  64       mpfr_init2 (err_s, 31);
  65       mpfr_init2 (err_u, 31);
  66
  67       /* Approximate sin(z) and cos(z). In the following, err <= k means that
  68          the approximate value y and the true value x are related by
  69          y = x * (1 + u)^k with |u| <= 2^(-w), following Higham's method. */
  70       mpfr_sin_cos (s, c, z, GMP_RNDN);
  71       if (MPFR_IS_NEG(z))
  72         mpfr_neg (s, s, GMP_RNDN); /* compute jn/yn(|z|), fix sign later */
  73       /* The absolute error on s/c is bounded by 1/2 ulp(1/2) <= 2^(-w-1). */
  74       mpfr_add (t, s, c, GMP_RNDN);
  75       mpfr_sub (c, s, c, GMP_RNDN);
  76       mpfr_swap (s, t);
  77       /* now s approximates sin(z)+cos(z), and c approximates sin(z)-cos(z),
  78          with total absolute error bounded by 2^(1-w). */
  79
  80       /* precompute 1/(8|z|) */
  81       mpfr_si_div (iz, MPFR_IS_POS(z) ? 1 : -1, z, GMP_RNDN);   /* err <= 1 */
  82       mpfr_div_2ui (iz, iz, 3, GMP_RNDN);
  83
  84       /* compute P and Q */
  85       mpfr_set_ui (P, 1, GMP_RNDN);
  86       mpfr_set_ui (Q, 0, GMP_RNDN);
  87       mpfr_set_ui (t, 1, GMP_RNDN); /* current term */
  88       mpfr_set_ui (err_t, 0, GMP_RNDN); /* error on t */
  89       mpfr_set_ui (err_s, 0, GMP_RNDN); /* error on P and Q (sum of errors) */
  90       for (k = 1, stop = 0; stop < 4; k++)
  91         {
  92           /* compute next term: t(k)/t(k-1) = (2n+2k-1)(2n-2k+1)/(8kz) */
  93           mpfr_mul_si (t, t, 2 * (n + k) - 1, GMP_RNDN); /* err <= err_k + 1 */
  94           mpfr_mul_si (t, t, 2 * (n - k) + 1, GMP_RNDN); /* err <= err_k + 2 */
  95           mpfr_div_ui (t, t, k, GMP_RNDN);               /* err <= err_k + 3 */
  96           mpfr_mul (t, t, iz, GMP_RNDN);                 /* err <= err_k + 5 */
  97           /* the relative error on t is bounded by (1+u)^(5k)-1, which is
  98              bounded by 6ku for 6ku <= 0.02: first |5 log(1+u)| <= |5.5u|
  99              for |u| <= 0.15, then |exp(5.5u)-1| <= 6u for |u| <= 0.02. */
 100           mpfr_mul_ui (err_t, t, 6 * k, MPFR_IS_POS(t) ? GMP_RNDU : GMP_RNDD);
 101           mpfr_abs (err_t, err_t, GMP_RNDN); /* exact */
 102           /* the absolute error on t is bounded by err_t * 2^(-w) */
 103           mpfr_abs (err_u, t, GMP_RNDU);
 104           mpfr_mul_2ui (err_u, err_u, w, GMP_RNDU); /* t * 2^w */
 105           mpfr_add (err_u, err_u, err_t, GMP_RNDU); /* max|t| * 2^w */
 106           if (stop >= 2)
 107             {
 108               /* take into account the neglected terms: t * 2^w */
 109               mpfr_div_2ui (err_s, err_s, w, GMP_RNDU);
 110               if (MPFR_IS_POS(t))
 111                 mpfr_add (err_s, err_s, t, GMP_RNDU);
 112               else
 113                 mpfr_sub (err_s, err_s, t, GMP_RNDU);
 114               mpfr_mul_2ui (err_s, err_s, w, GMP_RNDU);
 115               stop ++;
 116             }
 117           /* if k is odd, add to Q, otherwise to P */
 118           else if (k & 1)
 119             {
 120               /* if k = 1 mod 4, add, otherwise subtract */
 121               if ((k & 2) == 0)
 122                 mpfr_add (Q, Q, t, GMP_RNDN);
 123               else
 124                 mpfr_sub (Q, Q, t, GMP_RNDN);
 125               /* check if the next term is smaller than ulp(Q): if EXP(err_u)
 126                  <= EXP(Q), since the current term is bounded by
 127                  err_u * 2^(-w), it is bounded by ulp(Q) */
 128               if (MPFR_EXP(err_u) <= MPFR_EXP(Q))
 129                 stop ++;
 130               else
 131                 stop = 0;
 132             }
 133           else
 134             {
 135               /* if k = 0 mod 4, add, otherwise subtract */
 136               if ((k & 2) == 0)
 137                 mpfr_add (P, P, t, GMP_RNDN);
 138               else
 139                 mpfr_sub (P, P, t, GMP_RNDN);
 140               /* check if the next term is smaller than ulp(P) */
 141               if (MPFR_EXP(err_u) <= MPFR_EXP(P))
 142                 stop ++;
 143               else
 144                 stop = 0;
 145             }
 146           mpfr_add (err_s, err_s, err_t, GMP_RNDU);
 147           /* the sum of the rounding errors on P and Q is bounded by
 148              err_s * 2^(-w) */
 149
 150           /* stop when start to diverge */
 151           if (stop < 2 &&
 152               ((MPFR_IS_POS(z) && mpfr_cmp_ui (z, (k + 1) / 2) < 0) ||
 153                (MPFR_IS_NEG(z) && mpfr_cmp_si (z, - ((k + 1) / 2)) > 0)))
 154             {
 155               /* if we have to stop the series because it diverges, then
 156                  increasing the precision will most probably fail, since
 157                  we will stop to the same point, and thus compute a very
 158                  similar approximation */
 159               diverge = 1;
 160               stop = 2; /* force stop */
 161             }
 162         }
 163       /* the sum of the total errors on P and Q is bounded by err_s * 2^(-w) */
 164
 165       /* Now combine: the sum of the rounding errors on P and Q is bounded by
 166          err_s * 2^(-w), and the absolute error on s/c is bounded by 2^(1-w) */
 167       if ((n & 1) == 0) /* n even: P * (sin + cos) + Q (cos - sin) for jn
 168                                    Q * (sin + cos) + P (sin - cos) for yn */
 169         {
 170 #ifdef MPFR_JN
 171           mpfr_mul (c, c, Q, GMP_RNDN); /* Q * (sin - cos) */
 172           mpfr_mul (s, s, P, GMP_RNDN); /* P * (sin + cos) */
 173 #else
 174           mpfr_mul (c, c, P, GMP_RNDN); /* P * (sin - cos) */
 175           mpfr_mul (s, s, Q, GMP_RNDN); /* Q * (sin + cos) */
 176 #endif
 177           err = MPFR_EXP(c);
 178           if (MPFR_EXP(s) > err)
 179             err = MPFR_EXP(s);
 180 #ifdef MPFR_JN
 181           mpfr_sub (s, s, c, GMP_RNDN);
 182 #else
 183           mpfr_add (s, s, c, GMP_RNDN);
 184 #endif
 185         }
 186       else /* n odd: P * (sin - cos) + Q (cos + sin) for jn,
 187                      Q * (sin - cos) - P (cos + sin) for yn */
 188         {
 189 #ifdef MPFR_JN
 190           mpfr_mul (c, c, P, GMP_RNDN); /* P * (sin - cos) */
 191           mpfr_mul (s, s, Q, GMP_RNDN); /* Q * (sin + cos) */
 192 #else
 193           mpfr_mul (c, c, Q, GMP_RNDN); /* Q * (sin - cos) */
 194           mpfr_mul (s, s, P, GMP_RNDN); /* P * (sin + cos) */
 195 #endif
 196           err = MPFR_EXP(c);
 197           if (MPFR_EXP(s) > err)
 198             err = MPFR_EXP(s);
 199 #ifdef MPFR_JN
 200           mpfr_add (s, s, c, GMP_RNDN);
 201 #else
 202           mpfr_sub (s, c, s, GMP_RNDN);
 203 #endif
 204         }
 205       if ((n & 2) != 0)
 206         mpfr_neg (s, s, GMP_RNDN);
 207       if (MPFR_EXP(s) > err)
 208         err = MPFR_EXP(s);
 209       /* the absolute error on s is bounded by P*err(s/c) + Q*err(s/c)
 210          + err(P)*(s/c) + err(Q)*(s/c) + 3 * 2^(err - w - 1)
 211          <= (|P|+|Q|) * 2^(1-w) + err_s * 2^(1-w) + 2^err * 2^(1-w),
 212          since |c|, |old_s| <= 2. */
 213       err2 = (MPFR_EXP(P) >= MPFR_EXP(Q)) ? MPFR_EXP(P) + 2 : MPFR_EXP(Q) + 2;
 214       /* (|P| + |Q|) * 2^(1 - w) <= 2^(err2 - w) */
 215       err = MPFR_EXP(err_s) >= err ? MPFR_EXP(err_s) + 2 : err + 2;
 216       /* err_s * 2^(1-w) + 2^old_err * 2^(1-w) <= 2^err * 2^(-w) */
 217       err2 = (err >= err2) ? err + 1 : err2 + 1;
 218       /* now the absolute error on s is bounded by 2^(err2 - w) */
 219
 220       /* multiply by sqrt(1/(Pi*z)) */
 221       mpfr_const_pi (c, GMP_RNDN);     /* Pi, err <= 1 */
 222       mpfr_mul (c, c, z, GMP_RNDN);    /* err <= 2 */
 223       mpfr_si_div (c, MPFR_IS_POS(z) ? 1 : -1, c, GMP_RNDN); /* err <= 3 */
 224       mpfr_sqrt (c, c, GMP_RNDN);      /* err<=5/2, thus the absolute error is
 225                                           bounded by 3*u*|c| for |u| <= 0.25 */
 226       mpfr_mul (err_t, c, s, MPFR_SIGN(c)==MPFR_SIGN(s) ? GMP_RNDU : GMP_RNDD);
 227       mpfr_abs (err_t, err_t, GMP_RNDU);
 228       mpfr_mul_ui (err_t, err_t, 3, GMP_RNDU);
 229       /* 3*2^(-w)*|old_c|*|s| [see below] is bounded by err_t * 2^(-w) */
 230       err2 += MPFR_EXP(c);
 231       /* |old_c| * 2^(err2 - w) [see below] is bounded by 2^(err2-w) */
 232       mpfr_mul (c, c, s, GMP_RNDN);    /* the absolute error on c is bounded by
 233                                           1/2 ulp(c) + 3*2^(-w)*|old_c|*|s|
 234                                           + |old_c| * 2^(err2 - w) */
 235       /* compute err_t * 2^(-w) + 1/2 ulp(c) = (err_t + 2^EXP(c)) * 2^(-w) */
 236       err = (MPFR_EXP(err_t) > MPFR_EXP(c)) ? MPFR_EXP(err_t) + 1 : MPFR_EXP(c) + 1;
 237       /* err_t * 2^(-w) + 1/2 ulp(c) <= 2^(err - w) */
 238       /* now err_t * 2^(-w) bounds 1/2 ulp(c) + 3*2^(-w)*|old_c|*|s| */
 239       err = (err >= err2) ? err + 1 : err2 + 1;
 240       /* the absolute error on c is bounded by 2^(err - w) */
 241
 242       mpfr_clear (s);
 243       mpfr_clear (P);
 244       mpfr_clear (Q);
 245       mpfr_clear (t);
 246       mpfr_clear (iz);
 247       mpfr_clear (err_t);
 248       mpfr_clear (err_s);
 249       mpfr_clear (err_u);
 250
 251       err -= MPFR_EXP(c);
 252       if (MPFR_LIKELY (MPFR_CAN_ROUND (c, w - err, MPFR_PREC(res), r)))
 253         break;
 254       if (diverge != 0)
 255         {
 256           mpfr_set (c, z, r); /* will force inex=0 below, which means the
 257                                asymptotic expansion failed */
 258           break;
 259         }
 260       MPFR_ZIV_NEXT (loop, w);
 261     }
 262   MPFR_ZIV_FREE (loop);
 263
 264   inex = (MPFR_IS_POS(z) || ((n & 1) == 0)) ? mpfr_set (res, c, r)
 265     : mpfr_neg (res, c, r);
 266   mpfr_clear (c);
 267
 268   return inex;
 269 }