w32tex
About: TeX Live provides a comprehensive TeX system including all the major TeX-related programs, macro packages, and fonts that are free software. Windows sources.
  Fossies Dox: w32tex-src.tar.xz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

ai.c
Go to the documentation of this file.
1 /* mpfr_ai -- Airy function Ai
2 
3 Copyright 2010-2020 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
52 {
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 */
73  ("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
74  ("y[%Pu]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y) );
75 
76  /* Save current exponents range */
77  MPFR_SAVE_EXPO_MARK (expo);
78 
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);
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 */
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. */
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  {
151 
152  MPFR_BLOCK (flags, cond = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
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;
168 
169  MPFR_BLOCK (flags, t = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
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_init (ti);
192  mpfr_init (tip1);
193  mpfr_init (temp1);
194  mpfr_init (temp2);
195  mpfr_init (x3);
196  mpfr_init (s);
197 
198  /* ZIV loop */
199  for (;;)
200  {
201  MPFR_LOG_MSG (("Working precision: %Pu\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))
211  x3u = mpfr_get_ui (x3, MPFR_RNDU); /* x3u >= ceil(x^3) */
212  if (MPFR_IS_NEG (x))
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: %Pu\n", err));
261  MPFR_LOG_MSG (("Approxim error: %Pu\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 
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
323 {
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[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
346  ("y[%Pu]=%.*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 
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  {
376 
377  MPFR_BLOCK (flags, cond = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
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;
393 
394  MPFR_BLOCK (flags, t = mpfr_get_ui (tmp2_sp, MPFR_RNDU));
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_init (z[j]);
424 
425  mpfr_init (s);
426  mpfr_init (u0); mpfr_init (u1);
427  mpfr_init (result);
428  mpfr_init (temp1);
429  mpfr_init (temp2);
430 
431  /* ZIV loop */
432  for (;;)
433  {
434  MPFR_LOG_MSG (("working precision: %Pu\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 
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);
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);
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: %Pu\n", err));
552  MPFR_LOG_MSG (("Approxim error: %Pu\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_init (z[j]);
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",
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 
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
645 {
646  mpfr_t temp1, temp2;
647  int use_ai2;
648  MPFR_SAVE_EXPO_DECL (expo);
649 
650  /* Special cases */
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 
668 
669  mpfr_set (temp1, x, 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 }
int mpfr_ai(mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
Definition: ai.c:644
#define MPFR_AI_SCALE
Definition: ai.c:641
static int mpfr_ai1(mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
Definition: ai.c:51
static int mpfr_ai2(mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
Definition: ai.c:322
#define MPFR_AI_THRESHOLD1
Definition: mparam.h:230
#define MPFR_AI_THRESHOLD2
Definition: mparam.h:231
#define MPFR_AI_THRESHOLD3
Definition: mparam.h:232
#define x3
int mpfr_cbrt(mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
Definition: cbrt.c:49
int z
Definition: dviconv.c:26
void mpfr_gamma_one_and_two_third(mpfr_ptr y1, mpfr_ptr y2, mpfr_prec_t prec)
void mpfr_div_ui2(mpfr_ptr y, mpfr_srcptr x, unsigned long int v1, unsigned long int v2, mpfr_rnd_t mode)
Definition: gammaonethird.c:66
#define s
Definition: afcover.h:80
#define t
Definition: afcover.h:96
unsigned long __gmpfr_isqrt(unsigned long n)
Definition: isqrt.c:27
int test1()
#define NULL
Definition: ftobjs.h:61
small capitals from c petite p scientific f u
Definition: afcover.h:88
small capitals from c petite p scientific i
Definition: afcover.h:80
kerning y
Definition: ttdriver.c:212
#define ULONG_MAX
Definition: c-minmax.h:66
#define loop
Definition: tie.c:8
void mpfr_init(mpfr_ptr x)
Definition: init.c:26
void mpfr_free_func(void *ptr, size_t size)
Definition: mpfr-gmp.c:330
void * mpfr_allocate_func(size_t alloc_size)
Definition: mpfr-gmp.c:308
#define MPFR_IS_INF(x)
Definition: mpfr-impl.h:1082
#define MPFR_SET_NAN(x)
Definition: mpfr-impl.h:1081
#define MPFR_UNLIKELY(x)
Definition: mpfr-impl.h:1490
#define MPFR_LIKELY(x)
Definition: mpfr-impl.h:1489
#define MPFR_IS_NAN(x)
Definition: mpfr-impl.h:1080
#define MPFR_INT_CEIL_LOG2(x)
Definition: mpfr-impl.h:1558
#define MPFR_ZIV_INIT(_x, _p)
Definition: mpfr-impl.h:2053
#define MPFR_SMALL_PRECISION
Definition: mpfr-impl.h:591
#define MPFR_IS_NEG(x)
Definition: mpfr-impl.h:1140
#define MPFR_SAVE_EXPO_MARK(x)
Definition: mpfr-impl.h:1730
#define MPFR_PREC(x)
Definition: mpfr-impl.h:958
#define MPFR_GET_EXP(x)
Definition: mpfr-impl.h:1058
#define MPFR_IS_POS(x)
Definition: mpfr-impl.h:1141
#define MPFR_EXP_FSPEC
Definition: mpfr-impl.h:1005
#define MPFR_SAVE_EXPO_FREE(x)
Definition: mpfr-impl.h:1736
#define MPFR_IS_ZERO(x)
Definition: mpfr-impl.h:1084
#define MPFR_BLOCK_DECL(_flags)
Definition: mpfr-impl.h:427
#define MPFR_ADD_PREC(P, X)
Definition: mpfr-impl.h:2045
long mpfr_eexp_t
Definition: mpfr-impl.h:1001
#define MPFR_ZIV_NEXT(_x, _p)
Definition: mpfr-impl.h:2054
#define MPFR_RET_NAN
Definition: mpfr-impl.h:1184
#define mpfr_check_range(x, t, r)
Definition: mpfr-impl.h:1744
#define MPFR_CAN_ROUND(b, err, prec, rnd)
Definition: mpfr-impl.h:1934
#define MPFR_SAVE_EXPO_DECL(x)
Definition: mpfr-impl.h:1729
#define MPFR_ERANGEFLAG(_flags)
Definition: mpfr-impl.h:451
#define MPFR_ASSERTN(expr)
Definition: mpfr-impl.h:495
#define MPFR_ZIV_DECL(_x)
Definition: mpfr-impl.h:2052
#define MPFR_ZIV_FREE(x)
Definition: mpfr-impl.h:2055
#define MPFR_LOG_FUNC(x, y)
Definition: mpfr-impl.h:2227
#define MPFR_ASSERTD(expr)
Definition: mpfr-impl.h:516
#define MPFR_LOG_MSG(x)
Definition: mpfr-impl.h:2226
#define MPFR_CHANGE_SIGN(x)
Definition: mpfr-impl.h:1146
#define MPFR_BLOCK(_flags, _op)
Definition: mpfr-impl.h:431
#define MPFR_IS_SINGULAR(x)
Definition: mpfr-impl.h:1100
int mpfr_add(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
Definition: add.c:26
void mpfr_clear(mpfr_ptr m)
Definition: clear.c:26
int mpfr_div(mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
Definition: div.c:748
unsigned long mpfr_get_ui(mpfr_srcptr f, mpfr_rnd_t rnd)
Definition: get_ui.c:26
void mpfr_init2(mpfr_ptr x, mpfr_prec_t p)
Definition: init2.c:26
int mpfr_mul_si(mpfr_ptr, mpfr_srcptr, long, mpfr_rnd_t)
int mpfr_neg(mpfr_ptr, mpfr_srcptr, mpfr_rnd_t)
Definition: neg.c:26
int mpfr_ui_div(mpfr_ptr, unsigned long, mpfr_srcptr, mpfr_rnd_t)
int mpfr_sqr(mpfr_ptr, mpfr_srcptr, mpfr_rnd_t)
Definition: sqr.c:508
int mpfr_pow_ui(mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t)
#define mpfr_cmp_si(b, i)
Definition: mpfr.h:863
long mpfr_prec_t
Definition: mpfr.h:168
#define mpfr_set(a, b, r)
Definition: mpfr.h:864
__mpfr_struct mpfr_t[1]
Definition: mpfr.h:250
int mpfr_set_si(mpfr_ptr, long, mpfr_rnd_t)
Definition: set_si.c:27
#define mpfr_get_prec(_x)
Definition: mpfr.h:848
int mpfr_set_ui(mpfr_ptr, unsigned long, mpfr_rnd_t)
Definition: set_ui.c:27
#define mpfr_abs(a, b, r)
Definition: mpfr.h:865
mpfr_rnd_t
Definition: mpfr.h:102
@ MPFR_RNDN
Definition: mpfr.h:103
@ MPFR_RNDU
Definition: mpfr.h:105
@ MPFR_RNDD
Definition: mpfr.h:106
int mpfr_sqrt(mpfr_ptr, mpfr_srcptr, mpfr_rnd_t)
Definition: sqrt.c:498
int mpfr_set_str(mpfr_ptr, const char *, int, mpfr_rnd_t)
long mpfr_exp_t
Definition: mpfr.h:195
int mpfr_mul_ui(mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t)
int mpfr_mul(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t)
Definition: mul.c:732
void mpfr_set_prec(mpfr_ptr, mpfr_prec_t)
Definition: set_prec.c:26
@ err
Definition: mtxline.h:24
float x
Definition: cordic.py:15
int k
Definition: otp-parser.c:70
set set set set set set set set set set set set set set set set set set set set *set set set macro pixldst cond
int r
Definition: ppmqvga.c:68
#define y1
#define y2
#define flags
static int test2(RN *rn, const char *set)
Definition: liolib.c:433
Definition: dvips.h:235
int j
Definition: t4ht.c:1589
@ L
Definition: ubidiimp.h:45
#define rnd(x)
Definition: xim.h:106