xref: /netbsd-src/external/lgpl3/mpfr/dist/src/jyn_asympt.c (revision 8ecbf5f02b752fcb7debe1a8fab1dc82602bc760)
1 /* mpfr_jn_asympt, mpfr_yn_asympt -- shared code for mpfr_jn and mpfr_yn
2 
3 Copyright 2007-2018 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 http://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_init (c);
53 
54   w = MPFR_PREC(res) + MPFR_INT_CEIL_LOG2(MPFR_PREC(res)) + 4;
55 
56   MPFR_ZIV_INIT (loop, w);
57   for (;;)
58     {
59       mpfr_set_prec (c, w);
60       mpfr_init2 (s, w);
61       mpfr_init2 (P, w);
62       mpfr_init2 (Q, w);
63       mpfr_init2 (t, w);
64       mpfr_init2 (iz, w);
65       mpfr_init2 (err_t, 31);
66       mpfr_init2 (err_s, 31);
67       mpfr_init2 (err_u, 31);
68 
69       /* Approximate sin(z) and cos(z). In the following, err <= k means that
70          the approximate value y and the true value x are related by
71          y = x * (1 + u)^k with |u| <= 2^(-w), following Higham's method. */
72       mpfr_sin_cos (s, c, z, MPFR_RNDN);
73       if (MPFR_IS_NEG(z))
74         mpfr_neg (s, s, MPFR_RNDN); /* compute jn/yn(|z|), fix sign later */
75       /* The absolute error on s/c is bounded by 1/2 ulp(1/2) <= 2^(-w-1). */
76       mpfr_add (t, s, c, MPFR_RNDN);
77       mpfr_sub (c, s, c, MPFR_RNDN);
78       mpfr_swap (s, t);
79       /* now s approximates sin(z)+cos(z), and c approximates sin(z)-cos(z),
80          with total absolute error bounded by 2^(1-w). */
81 
82       /* precompute 1/(8|z|) */
83       mpfr_si_div (iz, MPFR_IS_POS(z) ? 1 : -1, z, MPFR_RNDN);   /* err <= 1 */
84       mpfr_div_2ui (iz, iz, 3, MPFR_RNDN);
85 
86       /* compute P and Q */
87       mpfr_set_ui (P, 1, MPFR_RNDN);
88       mpfr_set_ui (Q, 0, MPFR_RNDN);
89       mpfr_set_ui (t, 1, MPFR_RNDN); /* current term */
90       mpfr_set_ui (err_t, 0, MPFR_RNDN); /* error on t */
91       mpfr_set_ui (err_s, 0, MPFR_RNDN); /* error on P and Q (sum of errors) */
92       for (k = 1, stop = 0; stop < 4; k++)
93         {
94           /* compute next term: t(k)/t(k-1) = (2n+2k-1)(2n-2k+1)/(8kz) */
95           mpfr_mul_si (t, t, 2 * (n + k) - 1, MPFR_RNDN); /* err <= err_k + 1 */
96           mpfr_mul_si (t, t, 2 * (n - k) + 1, MPFR_RNDN); /* err <= err_k + 2 */
97           mpfr_div_ui (t, t, k, MPFR_RNDN);               /* err <= err_k + 3 */
98           mpfr_mul (t, t, iz, MPFR_RNDN);                 /* err <= err_k + 5 */
99           /* the relative error on t is bounded by (1+u)^(5k)-1, which is
100              bounded by 6ku for 6ku <= 0.02: first |5 log(1+u)| <= |5.5u|
101              for |u| <= 0.15, then |exp(5.5u)-1| <= 6u for |u| <= 0.02. */
102           mpfr_mul_ui (err_t, t, 6 * k, MPFR_IS_POS(t) ? MPFR_RNDU : MPFR_RNDD);
103           mpfr_abs (err_t, err_t, MPFR_RNDN); /* exact */
104           /* the absolute error on t is bounded by err_t * 2^(-w) */
105           mpfr_abs (err_u, t, MPFR_RNDU);
106           mpfr_mul_2ui (err_u, err_u, w, MPFR_RNDU); /* t * 2^w */
107           mpfr_add (err_u, err_u, err_t, MPFR_RNDU); /* max|t| * 2^w */
108           if (stop >= 2)
109             {
110               /* take into account the neglected terms: t * 2^w */
111               mpfr_div_2ui (err_s, err_s, w, MPFR_RNDU);
112               if (MPFR_IS_POS(t))
113                 mpfr_add (err_s, err_s, t, MPFR_RNDU);
114               else
115                 mpfr_sub (err_s, err_s, t, MPFR_RNDU);
116               mpfr_mul_2ui (err_s, err_s, w, MPFR_RNDU);
117               stop ++;
118             }
119           /* if k is odd, add to Q, otherwise to P */
120           else if (k & 1)
121             {
122               /* if k = 1 mod 4, add, otherwise subtract */
123               if ((k & 2) == 0)
124                 mpfr_add (Q, Q, t, MPFR_RNDN);
125               else
126                 mpfr_sub (Q, Q, t, MPFR_RNDN);
127               /* check if the next term is smaller than ulp(Q): if EXP(err_u)
128                  <= EXP(Q), since the current term is bounded by
129                  err_u * 2^(-w), it is bounded by ulp(Q) */
130               if (MPFR_EXP(err_u) <= MPFR_EXP(Q))
131                 stop ++;
132               else
133                 stop = 0;
134             }
135           else
136             {
137               /* if k = 0 mod 4, add, otherwise subtract */
138               if ((k & 2) == 0)
139                 mpfr_add (P, P, t, MPFR_RNDN);
140               else
141                 mpfr_sub (P, P, t, MPFR_RNDN);
142               /* check if the next term is smaller than ulp(P) */
143               if (MPFR_EXP(err_u) <= MPFR_EXP(P))
144                 stop ++;
145               else
146                 stop = 0;
147             }
148           mpfr_add (err_s, err_s, err_t, MPFR_RNDU);
149           /* the sum of the rounding errors on P and Q is bounded by
150              err_s * 2^(-w) */
151 
152           /* stop when start to diverge */
153           if (stop < 2 &&
154               ((MPFR_IS_POS(z) && mpfr_cmp_ui (z, (k + 1) / 2) < 0) ||
155                (MPFR_IS_NEG(z) && mpfr_cmp_si (z, - ((k + 1) / 2)) > 0)))
156             {
157               /* if we have to stop the series because it diverges, then
158                  increasing the precision will most probably fail, since
159                  we will stop to the same point, and thus compute a very
160                  similar approximation */
161               diverge = 1;
162               stop = 2; /* force stop */
163             }
164         }
165       /* the sum of the total errors on P and Q is bounded by err_s * 2^(-w) */
166 
167       /* Now combine: the sum of the rounding errors on P and Q is bounded by
168          err_s * 2^(-w), and the absolute error on s/c is bounded by 2^(1-w) */
169       if ((n & 1) == 0) /* n even: P * (sin + cos) + Q (cos - sin) for jn
170                                    Q * (sin + cos) + P (sin - cos) for yn */
171         {
172 #ifdef MPFR_JN
173           mpfr_mul (c, c, Q, MPFR_RNDN); /* Q * (sin - cos) */
174           mpfr_mul (s, s, P, MPFR_RNDN); /* P * (sin + cos) */
175 #else
176           mpfr_mul (c, c, P, MPFR_RNDN); /* P * (sin - cos) */
177           mpfr_mul (s, s, Q, MPFR_RNDN); /* Q * (sin + cos) */
178 #endif
179           err = MPFR_EXP(c);
180           if (MPFR_EXP(s) > err)
181             err = MPFR_EXP(s);
182 #ifdef MPFR_JN
183           mpfr_sub (s, s, c, MPFR_RNDN);
184 #else
185           mpfr_add (s, s, c, MPFR_RNDN);
186 #endif
187         }
188       else /* n odd: P * (sin - cos) + Q (cos + sin) for jn,
189                      Q * (sin - cos) - P (cos + sin) for yn */
190         {
191 #ifdef MPFR_JN
192           mpfr_mul (c, c, P, MPFR_RNDN); /* P * (sin - cos) */
193           mpfr_mul (s, s, Q, MPFR_RNDN); /* Q * (sin + cos) */
194 #else
195           mpfr_mul (c, c, Q, MPFR_RNDN); /* Q * (sin - cos) */
196           mpfr_mul (s, s, P, MPFR_RNDN); /* P * (sin + cos) */
197 #endif
198           err = MPFR_EXP(c);
199           if (MPFR_EXP(s) > err)
200             err = MPFR_EXP(s);
201 #ifdef MPFR_JN
202           mpfr_add (s, s, c, MPFR_RNDN);
203 #else
204           mpfr_sub (s, c, s, MPFR_RNDN);
205 #endif
206         }
207       if ((n & 2) != 0)
208         mpfr_neg (s, s, MPFR_RNDN);
209       if (MPFR_EXP(s) > err)
210         err = MPFR_EXP(s);
211       /* the absolute error on s is bounded by P*err(s/c) + Q*err(s/c)
212          + err(P)*(s/c) + err(Q)*(s/c) + 3 * 2^(err - w - 1)
213          <= (|P|+|Q|) * 2^(1-w) + err_s * 2^(1-w) + 2^err * 2^(1-w),
214          since |c|, |old_s| <= 2. */
215       err2 = (MPFR_EXP(P) >= MPFR_EXP(Q)) ? MPFR_EXP(P) + 2 : MPFR_EXP(Q) + 2;
216       /* (|P| + |Q|) * 2^(1 - w) <= 2^(err2 - w) */
217       err = MPFR_EXP(err_s) >= err ? MPFR_EXP(err_s) + 2 : err + 2;
218       /* err_s * 2^(1-w) + 2^old_err * 2^(1-w) <= 2^err * 2^(-w) */
219       err2 = (err >= err2) ? err + 1 : err2 + 1;
220       /* now the absolute error on s is bounded by 2^(err2 - w) */
221 
222       /* multiply by sqrt(1/(Pi*z)) */
223       mpfr_const_pi (c, MPFR_RNDN);     /* Pi, err <= 1 */
224       mpfr_mul (c, c, z, MPFR_RNDN);    /* err <= 2 */
225       mpfr_si_div (c, MPFR_IS_POS(z) ? 1 : -1, c, MPFR_RNDN); /* err <= 3 */
226       mpfr_sqrt (c, c, MPFR_RNDN);      /* err<=5/2, thus the absolute error is
227                                           bounded by 3*u*|c| for |u| <= 0.25 */
228       mpfr_mul (err_t, c, s, MPFR_SIGN(c)==MPFR_SIGN(s) ? MPFR_RNDU : MPFR_RNDD);
229       mpfr_abs (err_t, err_t, MPFR_RNDU);
230       mpfr_mul_ui (err_t, err_t, 3, MPFR_RNDU);
231       /* 3*2^(-w)*|old_c|*|s| [see below] is bounded by err_t * 2^(-w) */
232       err2 += MPFR_EXP(c);
233       /* |old_c| * 2^(err2 - w) [see below] is bounded by 2^(err2-w) */
234       mpfr_mul (c, c, s, MPFR_RNDN);    /* the absolute error on c is bounded by
235                                           1/2 ulp(c) + 3*2^(-w)*|old_c|*|s|
236                                           + |old_c| * 2^(err2 - w) */
237       /* compute err_t * 2^(-w) + 1/2 ulp(c) = (err_t + 2^EXP(c)) * 2^(-w) */
238       err = (MPFR_EXP(err_t) > MPFR_EXP(c)) ? MPFR_EXP(err_t) + 1 : MPFR_EXP(c) + 1;
239       /* err_t * 2^(-w) + 1/2 ulp(c) <= 2^(err - w) */
240       /* now err_t * 2^(-w) bounds 1/2 ulp(c) + 3*2^(-w)*|old_c|*|s| */
241       err = (err >= err2) ? err + 1 : err2 + 1;
242       /* the absolute error on c is bounded by 2^(err - w) */
243 
244       mpfr_clear (s);
245       mpfr_clear (P);
246       mpfr_clear (Q);
247       mpfr_clear (t);
248       mpfr_clear (iz);
249       mpfr_clear (err_t);
250       mpfr_clear (err_s);
251       mpfr_clear (err_u);
252 
253       err -= MPFR_EXP(c);
254       if (MPFR_LIKELY (MPFR_CAN_ROUND (c, w - err, MPFR_PREC(res), r)))
255         break;
256       if (diverge != 0)
257         {
258           MPFR_ZIV_FREE (loop);
259           mpfr_clear (c);
260           return 0; /* means that the asymptotic expansion failed */
261         }
262       MPFR_ZIV_NEXT (loop, w);
263     }
264   MPFR_ZIV_FREE (loop);
265 
266   inex = (MPFR_IS_POS(z) || ((n & 1) == 0)) ? mpfr_set (res, c, r)
267     : mpfr_neg (res, c, r);
268   mpfr_clear (c);
269 
270   /* for RNDF, mpfr_set or mpfr_neg may return 0, but if we return 0, it
271      would mean the asymptotic expansion failed, thus we return 1 instead */
272   return (r != MPFR_RNDF) ? inex : 1;
273 }
274