1 /* mpfr_ai -- Airy function Ai
2
3 Copyright 2010-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 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 /* Reminder and notations:
27 -----------------------
28
29 Ai is the solution of:
30 / y'' - x*y = 0
31 { Ai(0) = 1/ ( 9^(1/3)*Gamma(2/3) )
32 \ Ai'(0) = -1/ ( 3^(1/3)*Gamma(1/3) )
33
34 Series development:
35 Ai(x) = sum (a_i*x^i)
36 = sum (t_i)
37
38 Recurrences:
39 a_(i+3) = a_i / ((i+2)*(i+3))
40 t_(i+3) = t_i * x^3 / ((i+2)*(i+3))
41
42 Values:
43 a_0 = Ai(0) ~ 0.355
44 a_1 = Ai'(0) ~ -0.259
45 */
46
47
48 /* Airy function Ai evaluated by the most naive algorithm.
49 Assume that x is a finite number. */
50 static int
51 mpfr_ai1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
52 {
53 MPFR_ZIV_DECL (loop);
54 MPFR_SAVE_EXPO_DECL (expo);
55 mpfr_prec_t wprec; /* working precision */
56 mpfr_prec_t prec; /* target precision */
57 mpfr_prec_t err; /* used to estimate the evaluation error */
58 mpfr_prec_t correct_bits; /* estimates the number of correct bits*/
59 unsigned long int k;
60 unsigned long int cond; /* condition number of the series */
61 unsigned long int assumed_exponent; /* used as a lowerbound of |EXP(Ai(x))| */
62 int r;
63 mpfr_t s; /* used to store the partial sum */
64 mpfr_t ti, tip1; /* used to store successive values of t_i */
65 mpfr_t x3; /* used to store x^3 */
66 mpfr_t tmp_sp, tmp2_sp; /* small precision variables */
67 unsigned long int x3u; /* used to store ceil(x^3) */
68 mpfr_t temp1, temp2;
69 int test1, test2;
70
71 /* Logging */
72 MPFR_LOG_FUNC (
73 ("x[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
74 ("y[%Pd]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y) );
75
76 /* Save current exponents range */
77 MPFR_SAVE_EXPO_MARK (expo);
78
79 if (MPFR_UNLIKELY (MPFR_IS_ZERO (x)))
80 {
81 mpfr_t y1, y2;
82 prec = MPFR_ADD_PREC (MPFR_PREC (y), 3);
83 mpfr_init2 (y1, prec);
84 mpfr_init2 (y2, prec);
85 MPFR_ZIV_INIT (loop, prec);
86
87 /* ZIV loop */
88 for (;;)
89 {
90 mpfr_gamma_one_and_two_third (y1, y2, prec); /* y2 = Gamma(2/3)(1 + delta1), |delta1| <= 2^{1-prec}. */
91
92 r = mpfr_set_ui (y1, 9, MPFR_RNDN);
93 MPFR_ASSERTD (r == 0);
94 mpfr_cbrt (y1, y1, MPFR_RNDN); /* y1 = cbrt(9)(1 + delta2), |delta2| <= 2^{-prec}. */
95 mpfr_mul (y1, y1, y2, MPFR_RNDN);
96 mpfr_ui_div (y1, 1, y1, MPFR_RNDN);
97 if (MPFR_LIKELY (MPFR_CAN_ROUND (y1, prec - 3, MPFR_PREC (y), rnd)))
98 break;
99 MPFR_ZIV_NEXT (loop, prec);
100 }
101 r = mpfr_set (y, y1, rnd);
102 MPFR_ZIV_FREE (loop);
103 MPFR_SAVE_EXPO_FREE (expo);
104 mpfr_clear (y1);
105 mpfr_clear (y2);
106 return mpfr_check_range (y, r, rnd);
107 }
108
109 /* now x is not zero */
110 MPFR_ASSERTD(!MPFR_IS_ZERO(x));
111
112 /* FIXME: underflow for large values of |x| ? */
113
114
115 /* Set initial precision */
116 /* If we compute sum(i=0, N-1, t_i), the relative error is bounded by */
117 /* 2*(4N)*2^(1-wprec)*C(|x|)/Ai(x) */
118 /* where C(|x|) = 1 if 0<=x<=1 */
119 /* and C(|x|) = (1/2)*x^(-1/4)*exp(2/3 x^(3/2)) if x >= 1 */
120
121 /* A priori, we do not know N, so we estimate it to ~ prec */
122 /* If 0<=x<=1, we estimate Ai(x) ~ 1/8 */
123 /* if 1<=x, we estimate Ai(x) ~ (1/4)*x^(-1/4)*exp(-2/3 * x^(3/2)) */
124 /* if x<=0, ????? */
125
126 /* We begin with 11 guard bits */
127 prec = MPFR_ADD_PREC (MPFR_PREC (y), 11);
128 MPFR_ZIV_INIT (loop, prec);
129
130 /* The working precision is heuristically chosen in order to obtain */
131 /* approximately prec correct bits in the sum. To sum up: the sum */
132 /* is stopped when the *exact* sum gives ~ prec correct bit. And */
133 /* when it is stopped, the accuracy of the computed sum, with respect*/
134 /* to the exact one should be ~prec bits. */
135 mpfr_init2 (tmp_sp, MPFR_SMALL_PRECISION);
136 mpfr_init2 (tmp2_sp, MPFR_SMALL_PRECISION);
137 mpfr_abs (tmp_sp, x, MPFR_RNDU);
138 mpfr_pow_ui (tmp_sp, tmp_sp, 3, MPFR_RNDU);
139 mpfr_sqrt (tmp_sp, tmp_sp, MPFR_RNDU); /* tmp_sp ~ |x|^(3/2) */
140
141 /* 0.96179669392597567 >~ 2/3 * log2(e). See algorithms.tex */
142 mpfr_set_str (tmp2_sp, "0.96179669392597567", 10, MPFR_RNDU);
143 mpfr_mul (tmp2_sp, tmp_sp, tmp2_sp, MPFR_RNDU);
144
145 /* cond represents the number of lost bits in the evaluation of the sum */
146 if (MPFR_GET_EXP (x) <= 0)
147 cond = 0;
148 else
149 {
150 MPFR_BLOCK_DECL (flags);
151
152 MPFR_BLOCK (flags, cond = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
153 MPFR_ASSERTN (! MPFR_ERANGEFLAG (flags));
154 cond -= (MPFR_GET_EXP (x) - 1) / 4 + 1;
155 }
156
157 /* The variable assumed_exponent is used to store the maximal assumed */
158 /* exponent of Ai(x). More precisely, we assume that |Ai(x)| will be */
159 /* greater than 2^{-assumed_exponent}. */
160 if (MPFR_IS_POS (x))
161 {
162 if (MPFR_GET_EXP (x) <= 0)
163 assumed_exponent = 3;
164 else
165 {
166 unsigned long int t;
167 MPFR_BLOCK_DECL (flags);
168
169 MPFR_BLOCK (flags, t = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
170 MPFR_ASSERTN (! MPFR_ERANGEFLAG (flags));
171 assumed_exponent = t + 2 + (MPFR_GET_EXP (x) / 4 + 1);
172 MPFR_ASSERTN (assumed_exponent > t);
173 }
174 }
175 /* We do not know Ai (x) yet */
176 /* We cover the case when EXP (Ai (x))>=-10 */
177 else
178 assumed_exponent = 10;
179
180 {
181 unsigned long int t, u;
182
183 t = assumed_exponent + cond;
184 MPFR_ASSERTN (t >= cond);
185 u = MPFR_INT_CEIL_LOG2 (prec) + 5;
186 t += u;
187 MPFR_ASSERTN (t >= u);
188 wprec = MPFR_ADD_PREC (prec, t);
189 }
190
191 mpfr_init2 (ti, wprec);
192 mpfr_init2 (tip1, wprec);
193 mpfr_init2 (temp1, wprec);
194 mpfr_init2 (temp2, wprec);
195 mpfr_init2 (x3, wprec);
196 mpfr_init2 (s, wprec);
197
198 /* ZIV loop */
199 for (;;)
200 {
201 MPFR_LOG_MSG (("Working precision: %Pd\n", wprec));
202 mpfr_set_prec (ti, wprec);
203 mpfr_set_prec (tip1, wprec);
204 mpfr_set_prec (x3, wprec);
205 mpfr_set_prec (s, wprec);
206
207 mpfr_sqr (x3, x, MPFR_RNDU);
208 mpfr_mul (x3, x3, x, (MPFR_IS_POS (x)?MPFR_RNDU:MPFR_RNDD)); /* x3=x^3 */
209 if (MPFR_IS_NEG (x))
210 MPFR_CHANGE_SIGN (x3);
211 x3u = mpfr_get_ui (x3, MPFR_RNDU); /* x3u >= ceil(x^3) */
212 if (MPFR_IS_NEG (x))
213 MPFR_CHANGE_SIGN (x3);
214
215 mpfr_gamma_one_and_two_third (temp1, temp2, wprec);
216 mpfr_set_ui (ti, 9, MPFR_RNDN);
217 mpfr_cbrt (ti, ti, MPFR_RNDN);
218 mpfr_mul (ti, ti, temp2, MPFR_RNDN);
219 mpfr_ui_div (ti, 1, ti, MPFR_RNDN); /* ti = 1/( Gamma (2/3)*9^(1/3) ) */
220
221 mpfr_set_ui (tip1, 3, MPFR_RNDN);
222 mpfr_cbrt (tip1, tip1, MPFR_RNDN);
223 mpfr_mul (tip1, tip1, temp1, MPFR_RNDN);
224 mpfr_neg (tip1, tip1, MPFR_RNDN);
225 mpfr_div (tip1, x, tip1, MPFR_RNDN); /* tip1 = -x/(Gamma (1/3)*3^(1/3)) */
226
227 mpfr_add (s, ti, tip1, MPFR_RNDN);
228
229
230 /* Evaluation of the series */
231 k = 2;
232 for (;;)
233 {
234 mpfr_mul (ti, ti, x3, MPFR_RNDN);
235 mpfr_mul (tip1, tip1, x3, MPFR_RNDN);
236
237 mpfr_div_ui2 (ti, ti, k, (k+1), MPFR_RNDN);
238 mpfr_div_ui2 (tip1, tip1, (k+1), (k+2), MPFR_RNDN);
239
240 k += 3;
241 mpfr_add (s, s, ti, MPFR_RNDN);
242 mpfr_add (s, s, tip1, MPFR_RNDN);
243
244 /* FIXME: if s==0 */
245 test1 = MPFR_IS_ZERO (ti)
246 || (MPFR_GET_EXP (ti) + (mpfr_exp_t)prec + 3 <= MPFR_GET_EXP (s));
247 test2 = MPFR_IS_ZERO (tip1)
248 || (MPFR_GET_EXP (tip1) + (mpfr_exp_t)prec + 3 <= MPFR_GET_EXP (s));
249
250 if ( test1 && test2 && (x3u <= k*(k+1)/2) )
251 break; /* FIXME: if k*(k+1) overflows */
252 }
253
254 MPFR_LOG_MSG (("Truncation rank: %lu\n", k));
255
256 err = 4 + MPFR_INT_CEIL_LOG2 (k) + cond - MPFR_GET_EXP (s);
257
258 /* err is the number of bits lost due to the evaluation error */
259 /* wprec-(prec+1): number of bits lost due to the approximation error */
260 MPFR_LOG_MSG (("Roundoff error: %Pd\n", err));
261 MPFR_LOG_MSG (("Approxim error: %Pd\n", wprec-prec-1));
262
263 if (wprec < err + 1)
264 correct_bits = 0;
265 else if (wprec < err + prec +1)
266 correct_bits = wprec - err - 1; /* since wprec > err + 1,
267 correct_bits > 0 */
268 else
269 correct_bits = prec;
270
271 if (MPFR_LIKELY (MPFR_CAN_ROUND (s, correct_bits, MPFR_PREC (y), rnd)))
272 break;
273
274 if (correct_bits == 0)
275 {
276 assumed_exponent *= 2;
277 MPFR_LOG_MSG (("Not a single bit correct (assumed_exponent=%lu)\n",
278 assumed_exponent));
279 wprec = prec + 5 + MPFR_INT_CEIL_LOG2 (prec) + cond +
280 assumed_exponent;
281 }
282 else if (correct_bits < prec)
283 { /* The precision was badly chosen */
284 MPFR_LOG_MSG (("Bad assumption on the exponent of Ai(x)"
285 " (E=%" MPFR_EXP_FSPEC "d)\n",
286 (mpfr_eexp_t) MPFR_GET_EXP (s)));
287 wprec = prec + err + 1;
288 }
289 else
290 { /* We are really in a bad case of the TMD */
291 MPFR_ZIV_NEXT (loop, prec);
292
293 /* We update wprec */
294 /* We assume that K will not be multiplied by more than 4 */
295 wprec = prec + (MPFR_INT_CEIL_LOG2 (k) + 2) + 5 + cond
296 - MPFR_GET_EXP (s);
297 }
298
299 } /* End of ZIV loop */
300
301 MPFR_ZIV_FREE (loop);
302
303 r = mpfr_set (y, s, rnd);
304
305 mpfr_clear (ti);
306 mpfr_clear (tip1);
307 mpfr_clear (temp1);
308 mpfr_clear (temp2);
309 mpfr_clear (x3);
310 mpfr_clear (s);
311 mpfr_clear (tmp_sp);
312 mpfr_clear (tmp2_sp);
313
314 MPFR_SAVE_EXPO_FREE (expo);
315 return mpfr_check_range (y, r, rnd);
316 }
317
318
319 /* Airy function Ai evaluated by Smith algorithm.
320 Assume that x is a finite non-zero number. */
321 static int
322 mpfr_ai2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
323 {
324 MPFR_ZIV_DECL (loop);
325 MPFR_SAVE_EXPO_DECL (expo);
326 mpfr_prec_t wprec; /* working precision */
327 mpfr_prec_t prec; /* target precision */
328 mpfr_prec_t err; /* used to estimate the evaluation error */
329 mpfr_prec_t correctBits; /* estimates the number of correct bits*/
330 unsigned long int i, j, L, t;
331 unsigned long int cond; /* condition number of the series */
332 unsigned long int assumed_exponent; /* used as a lowerbound of |EXP(Ai(x))| */
333 int r; /* returned ternary value */
334 mpfr_t s; /* used to store the partial sum */
335 mpfr_t u0, u1;
336 mpfr_t *z; /* used to store the (x^3j) */
337 mpfr_t result;
338 mpfr_t tmp_sp, tmp2_sp; /* small precision variables */
339 unsigned long int x3u; /* used to store ceil (x^3) */
340 mpfr_t temp1, temp2;
341 int test0, test1;
342
343 /* Logging */
344 MPFR_LOG_FUNC (
345 ("x[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
346 ("y[%Pd]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y));
347
348 /* Save current exponents range */
349 MPFR_SAVE_EXPO_MARK (expo);
350
351 /* FIXME: underflow for large values of |x| */
352
353 /* Set initial precision */
354 /* See the analysis for the naive evaluation */
355
356 /* We begin with 11 guard bits */
357 prec = MPFR_PREC (y) + 11;
358 MPFR_ZIV_INIT (loop, prec);
359
360 mpfr_init2 (tmp_sp, MPFR_SMALL_PRECISION);
361 mpfr_init2 (tmp2_sp, MPFR_SMALL_PRECISION);
362 mpfr_abs (tmp_sp, x, MPFR_RNDU);
363 mpfr_pow_ui (tmp_sp, tmp_sp, 3, MPFR_RNDU);
364 mpfr_sqrt (tmp_sp, tmp_sp, MPFR_RNDU); /* tmp_sp ~ |x|^(3/2) */
365
366 /* 0.96179669392597567 >~ 2/3 * log2(e). See algorithms.tex */
367 mpfr_set_str (tmp2_sp, "0.96179669392597567", 10, MPFR_RNDU);
368 mpfr_mul (tmp2_sp, tmp_sp, tmp2_sp, MPFR_RNDU);
369
370 /* cond represents the number of lost bits in the evaluation of the sum */
371 if (MPFR_GET_EXP (x) <= 0)
372 cond = 0;
373 else
374 {
375 MPFR_BLOCK_DECL (flags);
376
377 MPFR_BLOCK (flags, cond = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
378 MPFR_ASSERTN (! MPFR_ERANGEFLAG (flags));
379 cond -= (MPFR_GET_EXP (x) - 1) / 4 + 1;
380 }
381
382 /* This variable is used to store the maximal assumed exponent of */
383 /* Ai(x). More precisely, we assume that |Ai(x)| will be greater than */
384 /* 2^{-assumed_exponent}. */
385 if (MPFR_IS_POS (x))
386 {
387 if (MPFR_GET_EXP (x) <= 0)
388 assumed_exponent = 3;
389 else
390 {
391 unsigned long int t;
392 MPFR_BLOCK_DECL (flags);
393
394 MPFR_BLOCK (flags, t = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
395 MPFR_ASSERTN (! MPFR_ERANGEFLAG (flags));
396 assumed_exponent = t + 2 + (MPFR_GET_EXP (x) / 4 + 1);
397 MPFR_ASSERTN (assumed_exponent > t);
398 }
399 }
400 /* We do not know Ai(x) yet */
401 /* We cover the case when EXP(Ai(x))>=-10 */
402 else
403 assumed_exponent = 10;
404
405 {
406 unsigned long int t, u;
407
408 t = assumed_exponent + cond;
409 MPFR_ASSERTN (t >= cond);
410 u = MPFR_INT_CEIL_LOG2 (prec) + 6;
411 t += u;
412 MPFR_ASSERTN (t >= u);
413 wprec = MPFR_ADD_PREC (prec, t);
414 }
415
416 /* We assume that the truncation rank will be ~ prec */
417 L = __gmpfr_isqrt (prec);
418 MPFR_LOG_MSG (("size of blocks L = %lu\n", L));
419
420 z = (mpfr_t *) mpfr_allocate_func ( (L + 1) * sizeof (mpfr_t) );
421 MPFR_ASSERTN (z != NULL);
422 for (j=0; j<=L; j++)
423 mpfr_init2 (z[j], wprec);
424
425 mpfr_init2 (s, wprec);
426 mpfr_init2 (u0, wprec); mpfr_init2 (u1, wprec);
427 mpfr_init2 (result, wprec);
428 mpfr_init2 (temp1, wprec);
429 mpfr_init2 (temp2, wprec);
430
431 /* ZIV loop */
432 for (;;)
433 {
434 MPFR_LOG_MSG (("working precision: %Pd\n", wprec));
435
436 for (j=0; j<=L; j++)
437 mpfr_set_prec (z[j], wprec);
438 mpfr_set_prec (s, wprec);
439 mpfr_set_prec (u0, wprec); mpfr_set_prec (u1, wprec);
440 mpfr_set_prec (result, wprec);
441
442 mpfr_set_ui (u0, 1, MPFR_RNDN);
443 mpfr_set (u1, x, MPFR_RNDN);
444
445 mpfr_set_ui (z[0], 1, MPFR_RNDU);
446 mpfr_sqr (z[1], u1, MPFR_RNDU);
447 mpfr_mul (z[1], z[1], x, (MPFR_IS_POS (x) ? MPFR_RNDU : MPFR_RNDD) );
448
449 if (MPFR_IS_NEG (x))
450 MPFR_CHANGE_SIGN (z[1]);
451 x3u = mpfr_get_ui (z[1], MPFR_RNDU); /* x3u >= ceil (x^3) */
452 if (MPFR_IS_NEG (x))
453 MPFR_CHANGE_SIGN (z[1]);
454
455 for (j=2; j<=L ;j++)
456 {
457 if (j%2 == 0)
458 mpfr_sqr (z[j], z[j/2], MPFR_RNDN);
459 else
460 mpfr_mul (z[j], z[j-1], z[1], MPFR_RNDN);
461 }
462
463 mpfr_gamma_one_and_two_third (temp1, temp2, wprec);
464 mpfr_set_ui (u0, 9, MPFR_RNDN);
465 mpfr_cbrt (u0, u0, MPFR_RNDN);
466 mpfr_mul (u0, u0, temp2, MPFR_RNDN);
467 mpfr_ui_div (u0, 1, u0, MPFR_RNDN); /* u0 = 1/( Gamma (2/3)*9^(1/3) ) */
468
469 mpfr_set_ui (u1, 3, MPFR_RNDN);
470 mpfr_cbrt (u1, u1, MPFR_RNDN);
471 mpfr_mul (u1, u1, temp1, MPFR_RNDN);
472 mpfr_neg (u1, u1, MPFR_RNDN);
473 mpfr_div (u1, x, u1, MPFR_RNDN); /* u1 = -x/(Gamma (1/3)*3^(1/3)) */
474
475 mpfr_set_ui (result, 0, MPFR_RNDN);
476 t = 0;
477
478 /* Evaluation of the series by Smith' method */
479 for (i=0; ; i++)
480 {
481 t += 3 * L;
482
483 /* k = 0 */
484 t -= 3;
485 mpfr_set (s, z[L-1], MPFR_RNDN);
486 for (j=L-2; ; j--)
487 {
488 t -= 3;
489 mpfr_div_ui2 (s, s, (t+2), (t+3), MPFR_RNDN);
490 mpfr_add (s, s, z[j], MPFR_RNDN);
491 if (j==0)
492 break;
493 }
494 mpfr_mul (s, s, u0, MPFR_RNDN);
495 mpfr_add (result, result, s, MPFR_RNDN);
496
497 mpfr_mul (u0, u0, z[L], MPFR_RNDN);
498 for (j=0; j<=L-1; j++)
499 {
500 mpfr_div_ui2 (u0, u0, (t + 2), (t + 3), MPFR_RNDN);
501 t += 3;
502 }
503
504 t++;
505
506 /* k = 1 */
507 t -= 3;
508 mpfr_set (s, z[L-1], MPFR_RNDN);
509 for (j=L-2; ; j--)
510 {
511 t -= 3;
512 mpfr_div_ui2 (s, s, (t + 2), (t + 3), MPFR_RNDN);
513 mpfr_add (s, s, z[j], MPFR_RNDN);
514 if (j==0)
515 break;
516 }
517 mpfr_mul (s, s, u1, MPFR_RNDN);
518 mpfr_add (result, result, s, MPFR_RNDN);
519
520 mpfr_mul (u1, u1, z[L], MPFR_RNDN);
521 for (j=0; j<=L-1; j++)
522 {
523 mpfr_div_ui2 (u1, u1, (t + 2), (t + 3), MPFR_RNDN);
524 t += 3;
525 }
526
527 t++;
528
529 /* k = 2 */
530 t++;
531
532 /* End of the loop over k */
533 t -= 3;
534
535 test0 = MPFR_IS_ZERO (u0) ||
536 MPFR_GET_EXP (u0) + (mpfr_exp_t)prec + 4 <= MPFR_GET_EXP (result);
537 test1 = MPFR_IS_ZERO (u1) ||
538 MPFR_GET_EXP (u1) + (mpfr_exp_t)prec + 4 <= MPFR_GET_EXP (result);
539
540 if ( test0 && test1 && (x3u <= (t + 2) * (t + 3) / 2) )
541 break;
542 }
543
544 MPFR_LOG_MSG (("Truncation rank: %lu\n", t));
545
546 err = (5 + MPFR_INT_CEIL_LOG2 (L+1) + MPFR_INT_CEIL_LOG2 (i+1)
547 + cond - MPFR_GET_EXP (result));
548
549 /* err is the number of bits lost due to the evaluation error */
550 /* wprec-(prec+1): number of bits lost due to the approximation error */
551 MPFR_LOG_MSG (("Roundoff error: %Pd\n", err));
552 MPFR_LOG_MSG (("Approxim error: %Pd\n", wprec - prec - 1));
553
554 if (wprec < err+1)
555 correctBits = 0;
556 else
557 {
558 if (wprec < err+prec+1)
559 correctBits = wprec - err - 1;
560 else
561 correctBits = prec;
562 }
563
564 if (MPFR_LIKELY (MPFR_CAN_ROUND (result, correctBits,
565 MPFR_PREC (y), rnd)))
566 break;
567
568 for (j=0; j<=L; j++)
569 mpfr_clear (z[j]);
570 mpfr_free_func (z, (L + 1) * sizeof (mpfr_t));
571 L = __gmpfr_isqrt (t);
572 MPFR_LOG_MSG (("size of blocks L = %lu\n", L));
573 z = (mpfr_t *) mpfr_allocate_func ( (L + 1) * sizeof (mpfr_t));
574 MPFR_ASSERTN (z != NULL);
575 for (j=0; j<=L; j++)
576 mpfr_init2 (z[j], wprec);
577
578 if (correctBits == 0)
579 {
580 assumed_exponent *= 2;
581 MPFR_LOG_MSG (("Not a single bit correct (assumed_exponent=%lu)\n",
582 assumed_exponent));
583 wprec = prec + 6 + MPFR_INT_CEIL_LOG2 (t) + cond + assumed_exponent;
584 }
585 else
586 {
587 if (correctBits < prec)
588 { /* The precision was badly chosen */
589 MPFR_LOG_MSG (("Bad assumption on the exponent of Ai(x)"
590 " (E=%" MPFR_EXP_FSPEC "d)\n",
591 (mpfr_eexp_t) MPFR_GET_EXP (result)));
592 wprec = prec + err + 1;
593 }
594 else
595 { /* We are really in a bad case of the TMD */
596 MPFR_ZIV_NEXT (loop, prec);
597
598 /* We update wprec */
599 /* We assume that t will not be multiplied by more than 4 */
600 wprec = (prec + (MPFR_INT_CEIL_LOG2 (t) + 2) + 6 + cond
601 - MPFR_GET_EXP (result));
602 }
603 }
604 } /* End of ZIV loop */
605
606 MPFR_ZIV_FREE (loop);
607
608 r = mpfr_set (y, result, rnd);
609
610 mpfr_clear (tmp_sp);
611 mpfr_clear (tmp2_sp);
612 for (j=0; j<=L; j++)
613 mpfr_clear (z[j]);
614 mpfr_free_func (z, (L + 1) * sizeof (mpfr_t));
615
616 mpfr_clear (s);
617 mpfr_clear (u0); mpfr_clear (u1);
618 mpfr_clear (result);
619 mpfr_clear (temp1);
620 mpfr_clear (temp2);
621
622 MPFR_SAVE_EXPO_FREE (expo);
623 return mpfr_check_range (y, r, rnd);
624 }
625
626 /* We consider that the boundary between the area where the naive method
627 should preferably be used and the area where Smith' method should preferably
628 be used has the following form:
629 it is a triangle defined by two lines (one for the negative values of x, and
630 one for the positive values of x) crossing at x=0.
631
632 More precisely,
633
634 * If x<0 and MPFR_AI_THRESHOLD1*x + MPFR_AI_THRESHOLD2*prec > MPFR_AI_SCALE,
635 use Smith' algorithm;
636 * If x>0 and MPFR_AI_THRESHOLD3*x + MPFR_AI_THRESHOLD2*prec > MPFR_AI_SCALE,
637 use Smith' algorithm;
638 * otherwise, use the naive method.
639 */
640
641 #define MPFR_AI_SCALE 1048576
642
643 int
644 mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
645 {
646 mpfr_t temp1, temp2;
647 int use_ai2;
648 MPFR_SAVE_EXPO_DECL (expo);
649
650 /* Special cases */
651 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
652 {
653 if (MPFR_IS_NAN (x))
654 {
655 MPFR_SET_NAN (y);
656 MPFR_RET_NAN;
657 }
658 else if (MPFR_IS_INF (x))
659 return mpfr_set_ui (y, 0, rnd);
660 /* the cases x = +0 or -0 will be treated below */
661 }
662
663 /* The exponent range must be large enough for the computation of temp1. */
664 MPFR_SAVE_EXPO_MARK (expo);
665
666 mpfr_init2 (temp1, MPFR_SMALL_PRECISION);
667 mpfr_init2 (temp2, MPFR_SMALL_PRECISION);
668
669 mpfr_set (temp1, x, MPFR_RNDN);
670 mpfr_set_si (temp2, MPFR_AI_THRESHOLD2, MPFR_RNDN);
671 mpfr_mul_ui (temp2, temp2, MPFR_PREC (y) > ULONG_MAX ?
672 ULONG_MAX : (unsigned long) MPFR_PREC (y), MPFR_RNDN);
673
674 if (MPFR_IS_NEG (x))
675 mpfr_mul_si (temp1, temp1, MPFR_AI_THRESHOLD1, MPFR_RNDN);
676 else
677 mpfr_mul_si (temp1, temp1, MPFR_AI_THRESHOLD3, MPFR_RNDN);
678
679 mpfr_add (temp1, temp1, temp2, MPFR_RNDN);
680 mpfr_clear (temp2);
681
682 use_ai2 = mpfr_cmp_si (temp1, MPFR_AI_SCALE) > 0;
683 mpfr_clear (temp1);
684
685 MPFR_SAVE_EXPO_FREE (expo); /* Ignore all previous exceptions. */
686
687 /* we use ai2 if |x|*AI_THRESHOLD1/3 + PREC(y)*AI_THRESHOLD2 > AI_SCALE,
688 which means x cannot be zero in mpfr_ai2 */
689 return use_ai2 ? mpfr_ai2 (y, x, rnd) : mpfr_ai1 (y, x, rnd);
690 }