1/* mpfr_ai -- Airy function Ai 2 3Copyright 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#define MPFR_NEED_LONGLONG_H 24#include "mpfr-impl.h" 25 26#define MPFR_SMALL_PRECISION 32 27 28 29/* Reminder and notations: 30 ----------------------- 31 32 Ai is the solution of: 33 / y'' - x*y = 0 34 { Ai(0) = 1/ ( 9^(1/3)*Gamma(2/3) ) 35 \ Ai'(0) = -1/ ( 3^(1/3)*Gamma(1/3) ) 36 37 Series development: 38 Ai(x) = sum (a_i*x^i) 39 = sum (t_i) 40 41 Recurrences: 42 a_(i+3) = a_i / ((i+2)*(i+3)) 43 t_(i+3) = t_i * x^3 / ((i+2)*(i+3)) 44 45 Values: 46 a_0 = Ai(0) ~ 0.355 47 a_1 = Ai'(0) ~ -0.259 48*/ 49 50 51/* Airy function Ai evaluated by the most naive algorithm */ 52int 53mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) 54{ 55 MPFR_ZIV_DECL (loop); 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*/ 61 unsigned long int k; 62 unsigned long int cond; /* condition number of the series */ 63 unsigned long int assumed_exponent; /* used as a lowerbound of |EXP(Ai(x))| */ 64 int r; 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) */ 70 mpfr_t temp1, temp2; 71 int test1, test2; 72 73 /* Logging */ 74 MPFR_LOG_FUNC ( ("x[%#R]=%R rnd=%d", x, x, rnd), ("y[%#R]=%R", y, y) ); 75 76 /* Special cases */ 77 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 78 { 79 if (MPFR_IS_NAN (x)) 80 { 81 MPFR_SET_NAN (y); 82 MPFR_RET_NAN; 83 } 84 else if (MPFR_IS_INF (x)) 85 return mpfr_set_ui (y, 0, rnd); 86 } 87 88 /* FIXME: handle the case x == 0 (and in a consistent way for +0 and -0) */ 89 90 /* Save current exponents range */ 91 MPFR_SAVE_EXPO_MARK (expo); 92 93 /* FIXME: underflow for large values of |x| ? */ 94 95 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 */ 101 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)) */ 105 /* if x<=0, ????? */ 106 107 /* We begin with 11 guard bits */ 108 prec = MPFR_PREC (y)+11; 109 MPFR_ZIV_INIT (loop, prec); 110 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 */ 121 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); 125 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) ) 128 cond = 0; 129 else 130 cond = mpfr_get_ui (tmp2_sp, MPFR_RNDU) - (MPFR_GET_EXP (x)-1)/4 - 1; 131 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}. */ 135 if (MPFR_IS_ZERO(x)) 136 assumed_exponent = 2; 137 else 138 { 139 if (MPFR_IS_POS (x)) 140 { 141 if (MPFR_GET_EXP (x) <= 0) 142 assumed_exponent = 3; 143 else 144 assumed_exponent = (2 + (MPFR_GET_EXP (x)/4 + 1) 145 + mpfr_get_ui (tmp2_sp, MPFR_RNDU)); 146 } 147 /* We do not know Ai (x) yet */ 148 /* We cover the case when EXP (Ai (x))>=-10 */ 149 else 150 assumed_exponent = 10; 151 } 152 153 wprec = prec + MPFR_INT_CEIL_LOG2 (prec) + 5 + cond + assumed_exponent; 154 155 mpfr_init (ti); 156 mpfr_init (tip1); 157 mpfr_init (temp1); 158 mpfr_init (temp2); 159 mpfr_init (x3); 160 mpfr_init (s); 161 162 /* ZIV loop */ 163 for (;;) 164 { 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); 170 171 mpfr_sqr (x3, x, MPFR_RNDU); 172 mpfr_mul (x3, x3, x, (MPFR_IS_POS (x)?MPFR_RNDU:MPFR_RNDD)); /* x3=x^3 */ 173 if (MPFR_IS_NEG (x)) 174 MPFR_CHANGE_SIGN (x3); 175 x3u = mpfr_get_ui (x3, MPFR_RNDU); /* x3u >= ceil(x^3) */ 176 if (MPFR_IS_NEG (x)) 177 MPFR_CHANGE_SIGN (x3); 178 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) ) */ 184 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)) */ 190 191 mpfr_add (s, ti, tip1, MPFR_RNDN); 192 193 194 /* Evaluation of the series */ 195 k = 2; 196 for (;;) 197 { 198 mpfr_mul (ti, ti, x3, MPFR_RNDN); 199 mpfr_mul (tip1, tip1, x3, MPFR_RNDN); 200 201 mpfr_div_ui2 (ti, ti, k, (k+1), MPFR_RNDN); 202 mpfr_div_ui2 (tip1, tip1, (k+1), (k+2), MPFR_RNDN); 203 204 k += 3; 205 mpfr_add (s, s, ti, MPFR_RNDN); 206 mpfr_add (s, s, tip1, MPFR_RNDN); 207 208 /* FIXME: if s==0 */ 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)); 213 214 if ( test1 && test2 && (x3u <= k*(k+1)/2) ) 215 break; /* FIXME: if k*(k+1) overflows */ 216 } 217 218 MPFR_LOG_MSG (("Truncation rank: %lu\n", k)); 219 220 err = 4 + MPFR_INT_CEIL_LOG2 (k) + cond - MPFR_GET_EXP (s); 221 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)); 226 227 if (wprec < err+1) 228 correct_bits=0; 229 else 230 { 231 if (wprec < err+prec+1) 232 correct_bits = wprec - err - 1; 233 else 234 correct_bits = prec; 235 } 236 237 if (MPFR_LIKELY (MPFR_CAN_ROUND (s, correct_bits, MPFR_PREC (y), rnd))) 238 break; 239 240 if (correct_bits == 0) 241 { 242 assumed_exponent *= 2; 243 MPFR_LOG_MSG (("Not a single bit correct (assumed_exponent=%lu)\n", 244 assumed_exponent)); 245 wprec = prec + 5 + MPFR_INT_CEIL_LOG2 (k) + cond + assumed_exponent; 246 } 247 else 248 { 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; 254 } 255 else 256 { /* We are really in a bad case of the TMD */ 257 MPFR_ZIV_NEXT (loop, prec); 258 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 262 - MPFR_GET_EXP (s); 263 } 264 } 265 266 } /* End of ZIV loop */ 267 268 MPFR_ZIV_FREE (loop); 269 270 r = mpfr_set (y, s, rnd); 271 272 mpfr_clear (ti); 273 mpfr_clear (tip1); 274 mpfr_clear (temp1); 275 mpfr_clear (temp2); 276 mpfr_clear (x3); 277 mpfr_clear (s); 278 mpfr_clear (tmp_sp); 279 mpfr_clear (tmp2_sp); 280 281 MPFR_SAVE_EXPO_FREE (expo); 282 return mpfr_check_range (y, r, rnd); 283} 284