1 /* mpfr_jn_asympt, mpfr_yn_asympt -- shared code for mpfr_jn and mpfr_yn
3 Copyright 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
24 # define FUNCTION mpfr_jn_asympt
27 # define FUNCTION mpfr_yn_asympt
29 # error "neither MPFR_JN nor MPFR_YN is defined"
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).
41 FUNCTION (mpfr_ptr res
, long n
, mpfr_srcptr z
, mp_rnd_t r
)
43 mpfr_t s
, c
, P
, Q
, t
, iz
, err_t
, err_s
, err_u
;
46 int inex
, stop
, diverge
= 0;
52 w
= MPFR_PREC(res
) + MPFR_INT_CEIL_LOG2(MPFR_PREC(res
)) + 4;
54 MPFR_ZIV_INIT (loop
, w
);
63 mpfr_init2 (err_t
, 31);
64 mpfr_init2 (err_s
, 31);
65 mpfr_init2 (err_u
, 31);
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
);
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
);
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). */
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
);
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
++)
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 */
108 /* take into account the neglected terms: t * 2^w */
109 mpfr_div_2ui (err_s
, err_s
, w
, GMP_RNDU
);
111 mpfr_add (err_s
, err_s
, t
, GMP_RNDU
);
113 mpfr_sub (err_s
, err_s
, t
, GMP_RNDU
);
114 mpfr_mul_2ui (err_s
, err_s
, w
, GMP_RNDU
);
117 /* if k is odd, add to Q, otherwise to P */
120 /* if k = 1 mod 4, add, otherwise subtract */
122 mpfr_add (Q
, Q
, t
, GMP_RNDN
);
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
))
135 /* if k = 0 mod 4, add, otherwise subtract */
137 mpfr_add (P
, P
, t
, GMP_RNDN
);
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
))
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
150 /* stop when start to diverge */
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)))
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 */
160 stop
= 2; /* force stop */
163 /* the sum of the total errors on P and Q is bounded by err_s * 2^(-w) */
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 */
171 mpfr_mul (c
, c
, Q
, GMP_RNDN
); /* Q * (sin - cos) */
172 mpfr_mul (s
, s
, P
, GMP_RNDN
); /* P * (sin + cos) */
174 mpfr_mul (c
, c
, P
, GMP_RNDN
); /* P * (sin - cos) */
175 mpfr_mul (s
, s
, Q
, GMP_RNDN
); /* Q * (sin + cos) */
178 if (MPFR_EXP(s
) > err
)
181 mpfr_sub (s
, s
, c
, GMP_RNDN
);
183 mpfr_add (s
, s
, c
, GMP_RNDN
);
186 else /* n odd: P * (sin - cos) + Q (cos + sin) for jn,
187 Q * (sin - cos) - P (cos + sin) for yn */
190 mpfr_mul (c
, c
, P
, GMP_RNDN
); /* P * (sin - cos) */
191 mpfr_mul (s
, s
, Q
, GMP_RNDN
); /* Q * (sin + cos) */
193 mpfr_mul (c
, c
, Q
, GMP_RNDN
); /* Q * (sin - cos) */
194 mpfr_mul (s
, s
, P
, GMP_RNDN
); /* P * (sin + cos) */
197 if (MPFR_EXP(s
) > err
)
200 mpfr_add (s
, s
, c
, GMP_RNDN
);
202 mpfr_sub (s
, c
, s
, GMP_RNDN
);
206 mpfr_neg (s
, s
, GMP_RNDN
);
207 if (MPFR_EXP(s
) > err
)
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) */
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) */
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) */
252 if (MPFR_LIKELY (MPFR_CAN_ROUND (c
, w
- err
, MPFR_PREC(res
), r
)))
256 mpfr_set (c
, z
, r
); /* will force inex=0 below, which means the
257 asymptotic expansion failed */
260 MPFR_ZIV_NEXT (loop
, w
);
262 MPFR_ZIV_FREE (loop
);
264 inex
= (MPFR_IS_POS(z
) || ((n
& 1) == 0)) ? mpfr_set (res
, c
, r
)
265 : mpfr_neg (res
, c
, r
);