1 /* mpfr_jn_asympt, mpfr_yn_asympt -- shared code for mpfr_jn and mpfr_yn
2
3 Copyright 2007-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 3 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.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, 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 Note: for MPFR_RNDF, it returns 0 if the expansion failed, and a non-zero
40 value otherwise (with no other meaning).
41 */
42 static int
43 FUNCTION (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r)
44 {
45 mpfr_t s, c, P, Q, t, iz, err_t, err_s, err_u;
46 mpfr_prec_t w;
47 long k;
48 int inex, stop, diverge = 0;
49 mpfr_exp_t err2, err;
50 MPFR_ZIV_DECL (loop);
51
52 mpfr_init2 (c, 64);
53
54 /* The terms of the asymptotic expansion grow like mu^(2k)/(8z)^(2k), where
55 mu = 4n^2, thus we need mu < 8|z| so that it converges,
56 i.e., n^2/2 < |z| */
57 MPFR_ASSERTD (n >= 0);
58 mpfr_set_ui (c, n, MPFR_RNDU);
59 mpfr_mul_ui (c, c, n, MPFR_RNDU);
60 mpfr_div_2ui (c, c, 1, MPFR_RNDU);
61 if (mpfr_cmpabs (c, z) >= 0)
62 {
63 mpfr_clear (c);
64 return 0; /* asymptotic expansion failed */
65 }
66
67 w = MPFR_PREC(res) + MPFR_INT_CEIL_LOG2(MPFR_PREC(res)) + 4;
68
69 MPFR_ZIV_INIT (loop, w);
70 for (;;)
71 {
72 int ok = 0;
73
74 mpfr_set_prec (c, w);
75 mpfr_init2 (s, w);
76 mpfr_init2 (P, w);
77 mpfr_init2 (Q, w);
78 mpfr_init2 (t, w);
79 mpfr_init2 (iz, w);
80 mpfr_init2 (err_t, 31);
81 mpfr_init2 (err_s, 31);
82 mpfr_init2 (err_u, 31);
83
84 /* Approximate sin(z) and cos(z). In the following, err <= k means that
85 the approximate value y and the true value x are related by
86 y = x * (1 + u)^k with |u| <= 2^(-w), following Higham's method. */
87 mpfr_sin_cos (s, c, z, MPFR_RNDN);
88 if (MPFR_IS_NEG(z))
89 mpfr_neg (s, s, MPFR_RNDN); /* compute jn/yn(|z|), fix sign later */
90 /* The absolute error on s/c is bounded by 1/2 ulp(1/2) <= 2^(-w-1). */
91 mpfr_add (t, s, c, MPFR_RNDN);
92 mpfr_sub (c, s, c, MPFR_RNDN);
93 mpfr_swap (s, t);
94 /* now s approximates sin(z)+cos(z), and c approximates sin(z)-cos(z),
95 with total absolute error bounded by 2^(1-w). */
96
97 /* if s or c is zero, MPFR_GET_EXP will fail below */
98 if (MPFR_IS_ZERO(s) || MPFR_IS_ZERO(c))
99 goto clear; /* with ok=0 */
100
101 /* precompute 1/(8|z|) */
102 mpfr_si_div (iz, MPFR_IS_POS(z) ? 1 : -1, z, MPFR_RNDN); /* err <= 1 */
103 mpfr_div_2ui (iz, iz, 3, MPFR_RNDN);
104
105 /* compute P and Q */
106 mpfr_set_ui (P, 1, MPFR_RNDN);
107 mpfr_set_ui (Q, 0, MPFR_RNDN);
108 mpfr_set_ui (t, 1, MPFR_RNDN); /* current term */
109 mpfr_set_ui (err_t, 0, MPFR_RNDN); /* error on t */
110 mpfr_set_ui (err_s, 0, MPFR_RNDN); /* error on P and Q (sum of errors) */
111 for (k = 1, stop = 0; stop < 4; k++)
112 {
113 /* compute next term: t(k)/t(k-1) = (2n+2k-1)(2n-2k+1)/(8kz) */
114 MPFR_LOG_MSG (("loop (k,stop) = (%ld,%d)\n", k, stop));
115 mpfr_mul_si (t, t, 2 * (n + k) - 1, MPFR_RNDN); /* err <= err_k + 1 */
116 mpfr_mul_si (t, t, 2 * (n - k) + 1, MPFR_RNDN); /* err <= err_k + 2 */
117 mpfr_div_ui (t, t, k, MPFR_RNDN); /* err <= err_k + 3 */
118 mpfr_mul (t, t, iz, MPFR_RNDN); /* err <= err_k + 5 */
119 /* the relative error on t is bounded by (1+u)^(5k)-1, which is
120 bounded by 6ku for 6ku <= 0.02: first |5 log(1+u)| <= |5.5u|
121 for |u| <= 0.15, then |exp(5.5u)-1| <= 6u for |u| <= 0.02. */
122 mpfr_mul_ui (err_t, t, 6 * k, MPFR_IS_POS(t) ? MPFR_RNDU : MPFR_RNDD);
123 mpfr_abs (err_t, err_t, MPFR_RNDN); /* exact */
124 /* the absolute error on t is bounded by err_t * 2^(-w) */
125 mpfr_abs (err_u, t, MPFR_RNDU);
126 mpfr_mul_2ui (err_u, err_u, w, MPFR_RNDU); /* t * 2^w */
127 mpfr_add (err_u, err_u, err_t, MPFR_RNDU); /* max|t| * 2^w */
128 if (stop >= 2)
129 {
130 /* take into account the neglected terms: t * 2^w */
131 mpfr_div_2ui (err_s, err_s, w, MPFR_RNDU);
132 if (MPFR_IS_POS(t))
133 mpfr_add (err_s, err_s, t, MPFR_RNDU);
134 else
135 mpfr_sub (err_s, err_s, t, MPFR_RNDU);
136 mpfr_mul_2ui (err_s, err_s, w, MPFR_RNDU);
137 stop ++;
138 }
139 /* if k is odd, add to Q, otherwise to P */
140 else if (k & 1)
141 {
142 /* if k = 1 mod 4, add, otherwise subtract */
143 if ((k & 2) == 0)
144 mpfr_add (Q, Q, t, MPFR_RNDN);
145 else
146 mpfr_sub (Q, Q, t, MPFR_RNDN);
147 /* check if the next term is smaller than ulp(Q): if EXP(err_u)
148 <= EXP(Q), since the current term is bounded by
149 err_u * 2^(-w), it is bounded by ulp(Q) */
150 if (MPFR_GET_EXP (err_u) <= MPFR_GET_EXP (Q))
151 stop ++;
152 else
153 stop = 0;
154 }
155 else
156 {
157 /* if k = 0 mod 4, add, otherwise subtract */
158 if ((k & 2) == 0)
159 mpfr_add (P, P, t, MPFR_RNDN);
160 else
161 mpfr_sub (P, P, t, MPFR_RNDN);
162 /* check if the next term is smaller than ulp(P) */
163 if (MPFR_GET_EXP (err_u) <= MPFR_GET_EXP (P))
164 stop ++;
165 else
166 stop = 0;
167 }
168 mpfr_add (err_s, err_s, err_t, MPFR_RNDU);
169 /* the sum of the rounding errors on P and Q is bounded by
170 err_s * 2^(-w) */
171
172 /* stop when start to diverge */
173 if (stop < 2 &&
174 ((MPFR_IS_POS(z) && mpfr_cmp_ui (z, (k + 1) / 2) < 0) ||
175 (MPFR_IS_NEG(z) && mpfr_cmp_si (z, - ((k + 1) / 2)) > 0)))
176 {
177 /* if we have to stop the series because it diverges, then
178 increasing the precision will most probably fail, since
179 we will stop to the same point, and thus compute a very
180 similar approximation */
181 diverge = 1;
182 stop = 2; /* force stop */
183 }
184 }
185 /* the sum of the total errors on P and Q is bounded by err_s * 2^(-w) */
186
187 /* Now combine: the sum of the rounding errors on P and Q is bounded by
188 err_s * 2^(-w), and the absolute error on s/c is bounded by 2^(1-w) */
189 if ((n & 1) == 0) /* n even: P * (sin + cos) + Q (cos - sin) for jn
190 Q * (sin + cos) + P (sin - cos) for yn */
191 {
192 #ifdef MPFR_JN
193 mpfr_mul (c, c, Q, MPFR_RNDN); /* Q * (sin - cos) */
194 mpfr_mul (s, s, P, MPFR_RNDN); /* P * (sin + cos) */
195 #else
196 mpfr_mul (c, c, P, MPFR_RNDN); /* P * (sin - cos) */
197 mpfr_mul (s, s, Q, MPFR_RNDN); /* Q * (sin + cos) */
198 #endif
199 err = MPFR_GET_EXP (c);
200 if (MPFR_GET_EXP (s) > err)
201 err = MPFR_EXP (s);
202 #ifdef MPFR_JN
203 mpfr_sub (s, s, c, MPFR_RNDN);
204 #else
205 mpfr_add (s, s, c, MPFR_RNDN);
206 #endif
207 }
208 else /* n odd: P * (sin - cos) + Q (cos + sin) for jn,
209 Q * (sin - cos) - P (cos + sin) for yn */
210 {
211 #ifdef MPFR_JN
212 mpfr_mul (c, c, P, MPFR_RNDN); /* P * (sin - cos) */
213 mpfr_mul (s, s, Q, MPFR_RNDN); /* Q * (sin + cos) */
214 #else
215 mpfr_mul (c, c, Q, MPFR_RNDN); /* Q * (sin - cos) */
216 mpfr_mul (s, s, P, MPFR_RNDN); /* P * (sin + cos) */
217 #endif
218 err = MPFR_GET_EXP (c);
219 if (MPFR_GET_EXP (s) > err)
220 err = MPFR_EXP (s);
221 #ifdef MPFR_JN
222 mpfr_add (s, s, c, MPFR_RNDN);
223 #else
224 mpfr_sub (s, c, s, MPFR_RNDN);
225 #endif
226 }
227 if (MPFR_IS_ZERO(s))
228 goto clear; /* with ok=0 */
229 ok = 1;
230 if ((n & 2) != 0)
231 mpfr_neg (s, s, MPFR_RNDN);
232 if (MPFR_GET_EXP (s) > err)
233 err = MPFR_EXP (s);
234 /* the absolute error on s is bounded by P*err(s/c) + Q*err(s/c)
235 + err(P)*(s/c) + err(Q)*(s/c) + 3 * 2^(err - w - 1)
236 <= (|P|+|Q|) * 2^(1-w) + err_s * 2^(1-w) + 2^err * 2^(1-w),
237 since |c|, |old_s| <= 2. */
238 err2 = (MPFR_GET_EXP (P) >= MPFR_GET_EXP (Q))
239 ? MPFR_EXP (P) + 2 : MPFR_EXP (Q) + 2;
240 /* (|P| + |Q|) * 2^(1 - w) <= 2^(err2 - w) */
241 err = MPFR_GET_EXP (err_s) >= err ? MPFR_EXP (err_s) + 2 : err + 2;
242 /* err_s * 2^(1-w) + 2^old_err * 2^(1-w) <= 2^err * 2^(-w) */
243 err2 = (err >= err2) ? err + 1 : err2 + 1;
244 /* now the absolute error on s is bounded by 2^(err2 - w) */
245
246 /* multiply by sqrt(1/(Pi*z)) */
247 mpfr_const_pi (c, MPFR_RNDN); /* Pi, err <= 1 */
248 mpfr_mul (c, c, z, MPFR_RNDN); /* err <= 2 */
249 mpfr_si_div (c, MPFR_IS_POS(z) ? 1 : -1, c, MPFR_RNDN); /* err <= 3 */
250 mpfr_sqrt (c, c, MPFR_RNDN); /* err<=5/2, thus the absolute error is
251 bounded by 3*u*|c| for |u| <= 0.25 */
252 mpfr_mul (err_t, c, s, MPFR_SIGN(c)==MPFR_SIGN(s) ? MPFR_RNDU : MPFR_RNDD);
253 mpfr_abs (err_t, err_t, MPFR_RNDU);
254 mpfr_mul_ui (err_t, err_t, 3, MPFR_RNDU);
255 /* 3*2^(-w)*|old_c|*|s| [see below] is bounded by err_t * 2^(-w) */
256 err2 += MPFR_GET_EXP (c);
257 /* |old_c| * 2^(err2 - w) [see below] is bounded by 2^(err2-w) */
258 mpfr_mul (c, c, s, MPFR_RNDN); /* the absolute error on c is bounded by
259 1/2 ulp(c) + 3*2^(-w)*|old_c|*|s|
260 + |old_c| * 2^(err2 - w) */
261 /* compute err_t * 2^(-w) + 1/2 ulp(c) = (err_t + 2^EXP(c)) * 2^(-w) */
262 err = (MPFR_GET_EXP (err_t) > MPFR_GET_EXP (c)) ?
263 MPFR_EXP (err_t) + 1 : MPFR_EXP (c) + 1;
264 /* err_t * 2^(-w) + 1/2 ulp(c) <= 2^(err - w) */
265 /* now err_t * 2^(-w) bounds 1/2 ulp(c) + 3*2^(-w)*|old_c|*|s| */
266 err = (err >= err2) ? err + 1 : err2 + 1;
267 /* the absolute error on c is bounded by 2^(err - w) */
268
269 err -= MPFR_GET_EXP (c);
270
271 clear:
272 mpfr_clear (s);
273 mpfr_clear (P);
274 mpfr_clear (Q);
275 mpfr_clear (t);
276 mpfr_clear (iz);
277 mpfr_clear (err_t);
278 mpfr_clear (err_s);
279 mpfr_clear (err_u);
280
281 if (ok && MPFR_LIKELY (MPFR_CAN_ROUND (c, w - err, MPFR_PREC(res), r)))
282 break;
283 if (diverge != 0)
284 {
285 MPFR_ZIV_FREE (loop);
286 mpfr_clear (c);
287 return 0; /* means that the asymptotic expansion failed */
288 }
289 MPFR_ZIV_NEXT (loop, w);
290 }
291 MPFR_ZIV_FREE (loop);
292
293 inex = mpfr_set4 (res, c, r, MPFR_IS_POS (z) || (n & 1) == 0 ?
294 MPFR_SIGN (c) : - MPFR_SIGN (c));
295 mpfr_clear (c);
296
297 /* for RNDF, mpfr_set or mpfr_neg may return 0, but if we return 0, it
298 would mean the asymptotic expansion failed, thus we return 1 instead */
299 return (r != MPFR_RNDF) ? inex : 1;
300 }