source/libs/mpfr/mpfr-3.1.3/src/ai.c

   1 /* mpfr_ai -- Airy function Ai
   2
   3 Copyright 2010-2015 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 }