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