1 /* mpfr_ai -- Airy function Ai
3 Copyright 2010 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 #define MPFR_SMALL_PRECISION 32
29 /* Reminder and notations:
30 -----------------------
32 Ai is the solution of:
34 { Ai(0) = 1/ ( 9^(1/3)*Gamma(2/3) )
35 \ Ai'(0) = -1/ ( 3^(1/3)*Gamma(1/3) )
42 a_(i+3) = a_i / ((i+2)*(i+3))
43 t_(i+3) = t_i * x^3 / ((i+2)*(i+3))
51 /* Airy function Ai evaluated by the most naive algorithm */
53 mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
56 MPFR_SAVE_EXPO_DECL (expo);
57 mpfr_prec_t wprec; /* working precision */
58 mpfr_prec_t prec; /* target precision */
59 mpfr_prec_t err; /* used to estimate the evaluation error */
60 mpfr_prec_t correct_bits; /* estimates the number of correct bits*/
62 unsigned long int cond; /* condition number of the series */
63 unsigned long int assumed_exponent; /* used as a lowerbound of |EXP(Ai(x))| */
65 mpfr_t s; /* used to store the partial sum */
66 mpfr_t ti, tip1; /* used to store successive values of t_i */
67 mpfr_t x3; /* used to store x^3 */
68 mpfr_t tmp_sp, tmp2_sp; /* small precision variables */
69 unsigned long int x3u; /* used to store ceil(x^3) */
74 MPFR_LOG_FUNC ( ("x[%#R]=%R rnd=%d", x, x, rnd), ("y[%#R]=%R", y, y) );
77 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
84 else if (MPFR_IS_INF (x))
85 return mpfr_set_ui (y, 0, rnd);
88 /* FIXME: handle the case x == 0 (and in a consistent way for +0 and -0) */
90 /* Save current exponents range */
91 MPFR_SAVE_EXPO_MARK (expo);
93 /* FIXME: underflow for large values of |x| ? */
96 /* Set initial precision */
97 /* If we compute sum(i=0, N-1, t_i), the relative error is bounded by */
98 /* 2*(4N)*2^(1-wprec)*C(|x|)/Ai(x) */
99 /* where C(|x|) = 1 if 0<=x<=1 */
100 /* and C(|x|) = (1/2)*x^(-1/4)*exp(2/3 x^(3/2)) if x >= 1 */
102 /* A priori, we do not know N, so we estimate it to ~ prec */
103 /* If 0<=x<=1, we estimate Ai(x) ~ 1/8 */
104 /* if 1<=x, we estimate Ai(x) ~ (1/4)*x^(-1/4)*exp(-2/3 * x^(3/2)) */
107 /* We begin with 11 guard bits */
108 prec = MPFR_PREC (y)+11;
109 MPFR_ZIV_INIT (loop, prec);
111 /* The working precision is heuristically chosen in order to obtain */
112 /* approximately prec correct bits in the sum. To sum up: the sum */
113 /* is stopped when the *exact* sum gives ~ prec correct bit. And */
114 /* when it is stopped, the accuracy of the computed sum, with respect*/
115 /* to the exact one should be ~prec bits. */
116 mpfr_init2 (tmp_sp, MPFR_SMALL_PRECISION);
117 mpfr_init2 (tmp2_sp, MPFR_SMALL_PRECISION);
118 mpfr_abs (tmp_sp, x, MPFR_RNDU);
119 mpfr_pow_ui (tmp_sp, tmp_sp, 3, MPFR_RNDU);
120 mpfr_sqrt (tmp_sp, tmp_sp, MPFR_RNDU); /* tmp_sp ~ x^3/2 */
122 /* 0.96179669392597567 >~ 2/3 * log2(e). See algorithms.tex */
123 mpfr_set_str(tmp2_sp, "0.96179669392597567", 10, MPFR_RNDU);
124 mpfr_mul (tmp2_sp, tmp_sp, tmp2_sp, MPFR_RNDU);
126 /* cond represents the number of lost bits in the evaluation of the sum */
127 if ( (MPFR_IS_ZERO(x)) || (MPFR_GET_EXP (x) <= 0) )
130 cond = mpfr_get_ui (tmp2_sp, MPFR_RNDU) - (MPFR_GET_EXP (x)-1)/4 - 1;
132 /* The variable assumed_exponent is used to store the maximal assumed */
133 /* exponent of Ai(x). More precisely, we assume that |Ai(x)| will be */
134 /* greater than 2^{-assumed_exponent}. */
136 assumed_exponent = 2;
141 if (MPFR_GET_EXP (x) <= 0)
142 assumed_exponent = 3;
144 assumed_exponent = (2 + (MPFR_GET_EXP (x)/4 + 1)
145 + mpfr_get_ui (tmp2_sp, MPFR_RNDU));
147 /* We do not know Ai (x) yet */
148 /* We cover the case when EXP (Ai (x))>=-10 */
150 assumed_exponent = 10;
153 wprec = prec + MPFR_INT_CEIL_LOG2 (prec) + 5 + cond + assumed_exponent;
165 MPFR_LOG_MSG (("Working precision: %Pu\n", wprec));
166 mpfr_set_prec (ti, wprec);
167 mpfr_set_prec (tip1, wprec);
168 mpfr_set_prec (x3, wprec);
169 mpfr_set_prec (s, wprec);
171 mpfr_sqr (x3, x, MPFR_RNDU);
172 mpfr_mul (x3, x3, x, (MPFR_IS_POS (x)?MPFR_RNDU:MPFR_RNDD)); /* x3=x^3 */
174 MPFR_CHANGE_SIGN (x3);
175 x3u = mpfr_get_ui (x3, MPFR_RNDU); /* x3u >= ceil(x^3) */
177 MPFR_CHANGE_SIGN (x3);
179 mpfr_gamma_one_and_two_third (temp1, temp2, wprec);
180 mpfr_set_ui (ti, 9, MPFR_RNDN);
181 mpfr_cbrt (ti, ti, MPFR_RNDN);
182 mpfr_mul (ti, ti, temp2, MPFR_RNDN);
183 mpfr_ui_div (ti, 1, ti , MPFR_RNDN); /* ti = 1/( Gamma (2/3)*9^(1/3) ) */
185 mpfr_set_ui (tip1, 3, MPFR_RNDN);
186 mpfr_cbrt (tip1, tip1, MPFR_RNDN);
187 mpfr_mul (tip1, tip1, temp1, MPFR_RNDN);
188 mpfr_neg (tip1, tip1, MPFR_RNDN);
189 mpfr_div (tip1, x, tip1, MPFR_RNDN); /* tip1 = -x/(Gamma (1/3)*3^(1/3)) */
191 mpfr_add (s, ti, tip1, MPFR_RNDN);
194 /* Evaluation of the series */
198 mpfr_mul (ti, ti, x3, MPFR_RNDN);
199 mpfr_mul (tip1, tip1, x3, MPFR_RNDN);
201 mpfr_div_ui2 (ti, ti, k, (k+1), MPFR_RNDN);
202 mpfr_div_ui2 (tip1, tip1, (k+1), (k+2), MPFR_RNDN);
205 mpfr_add (s, s, ti, MPFR_RNDN);
206 mpfr_add (s, s, tip1, MPFR_RNDN);
209 test1 = MPFR_IS_ZERO (ti)
210 || (MPFR_GET_EXP (ti) + (mpfr_exp_t)prec + 3 <= MPFR_GET_EXP (s));
211 test2 = MPFR_IS_ZERO (tip1)
212 || (MPFR_GET_EXP (tip1) + (mpfr_exp_t)prec + 3 <= MPFR_GET_EXP (s));
214 if ( test1 && test2 && (x3u <= k*(k+1)/2) )
215 break; /* FIXME: if k*(k+1) overflows */
218 MPFR_LOG_MSG (("Truncation rank: %lu\n", k));
220 err = 4 + MPFR_INT_CEIL_LOG2 (k) + cond - MPFR_GET_EXP (s);
222 /* err is the number of bits lost due to the evaluation error */
223 /* wprec-(prec+1): number of bits lost due to the approximation error */
224 MPFR_LOG_MSG (("Roundoff error: %Pu\n", err));
225 MPFR_LOG_MSG (("Approxim error: %Pu\n", wprec-prec-1));
231 if (wprec < err+prec+1)
232 correct_bits = wprec - err - 1;
237 if (MPFR_LIKELY (MPFR_CAN_ROUND (s, correct_bits, MPFR_PREC (y), rnd)))
240 if (correct_bits == 0)
242 assumed_exponent *= 2;
243 MPFR_LOG_MSG (("Not a single bit correct (assumed_exponent=%lu)\n",
245 wprec = prec + 5 + MPFR_INT_CEIL_LOG2 (k) + cond + assumed_exponent;
249 if (correct_bits < prec)
250 { /* The precision was badly chosen */
251 MPFR_LOG_MSG (("Bad assumption on the exponent of Ai(x)", 0));
252 MPFR_LOG_MSG ((" (E=%ld)\n", (long) MPFR_GET_EXP (s)));
253 wprec = prec + err + 1;
256 { /* We are really in a bad case of the TMD */
257 MPFR_ZIV_NEXT (loop, prec);
259 /* We update wprec */
260 /* We assume that K will not be multiplied by more than 4 */
261 wprec = prec + (MPFR_INT_CEIL_LOG2 (k)+2) + 5 + cond
266 } /* End of ZIV loop */
268 MPFR_ZIV_FREE (loop);
270 r = mpfr_set (y, s, rnd);
279 mpfr_clear (tmp2_sp);
281 MPFR_SAVE_EXPO_FREE (expo);
282 return mpfr_check_range (y, r, rnd);