diff options
| -rw-r--r-- | Makefile.am | 2 | ||||
| -rw-r--r-- | ai.c | 48 | ||||
| -rw-r--r-- | ai2.c | 323 | ||||
| -rw-r--r-- | gammaonethird.c | 174 | ||||
| -rw-r--r-- | mpfr.h | 1 |
5 files changed, 518 insertions, 30 deletions
diff --git a/Makefile.am b/Makefile.am index c85ae7cc6..dc8440bee 100644 --- a/Makefile.am +++ b/Makefile.am @@ -55,7 +55,7 @@ round_near_x.c constant.c abort_prec_max.c stack_interface.c lngamma.c \ zeta_ui.c set_d64.c get_d64.c jn.c yn.c rem1.c get_patches.c add_d.c \ sub_d.c d_sub.c mul_d.c div_d.c d_div.c li2.c rec_sqrt.c min_prec.c \ buildopt.c digamma.c bernoulli.c isregular.c set_flt.c get_flt.c \ -scale2.c set_z_exp.c ai.c +scale2.c set_z_exp.c ai.c ai2.c gammaonethird.c libmpfr_la_LIBADD = @LIBOBJS@ @@ -23,7 +23,12 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" -#define MPFR_SMALL_PRECISION MPFR_PREC_MIN +#define MPFR_SMALL_PRECISION 32 + +/* These functions are provided by file gammaonethird.c */ +int mpfr_gamma_one_and_two_third (mpfr_ptr, mpfr_ptr, mp_prec_t); +void mpfr_div_ui2 (mpfr_ptr, mpfr_srcptr, unsigned long int, unsigned long int, mpfr_rnd_t); + /* Reminder and notations: ----------------------- @@ -66,7 +71,7 @@ mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) mpfr_t x3; /* used to store x^3 */ mpfr_t tmp_sp, tmp2_sp; /* small precision variables */ unsigned long int x3u; /* used to store ceil(x^3) */ - mpfr_t temp; + mpfr_t temp1, temp2; int test1, test2; /* Logging */ @@ -142,7 +147,8 @@ mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) mpfr_init (ti); mpfr_init (tip1); - mpfr_init (temp); + mpfr_init (temp1); + mpfr_init (temp2); mpfr_init (x3); mpfr_init (s); @@ -152,7 +158,6 @@ mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) MPFR_LOG_MSG (("Working precision: %d\n", wprec, 0)); mpfr_set_prec (ti, wprec); mpfr_set_prec (tip1, wprec); - mpfr_set_prec (temp, wprec+4); mpfr_set_prec (x3, wprec); mpfr_set_prec (s, wprec); @@ -162,23 +167,17 @@ mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) x3u = mpfr_get_ui (x3, MPFR_RNDU); /* x3u >= ceil(x^3) */ if (MPFR_IS_NEG (x)) MPFR_CHANGE_SIGN (x3); - + mpfr_gamma_one_and_two_third (temp1, temp2, wprec); mpfr_set_ui (ti, 9, MPFR_RNDN); mpfr_cbrt (ti, ti, MPFR_RNDN); - mpfr_set_ui (temp, 2, MPFR_RNDN); - mpfr_div_ui (temp, temp, 3, MPFR_RNDN); - mpfr_gamma (temp, temp, MPFR_RNDN); - mpfr_mul (ti, ti, temp, MPFR_RNDN); - mpfr_ui_div (ti, 1, ti , MPFR_RNDN); /* ti = 1 / ( Gamma(2/3)*9^(1/3) ) */ + mpfr_mul (ti, ti, temp2, MPFR_RNDN); + mpfr_ui_div (ti, 1, ti , MPFR_RNDN); /* ti = 1 / ( Gamma (2/3)*9^(1/3) ) */ mpfr_set_ui (tip1, 3, MPFR_RNDN); mpfr_cbrt (tip1, tip1, MPFR_RNDN); - mpfr_set_ui (temp, 1, MPFR_RNDN); - mpfr_div_ui (temp, temp, 3, MPFR_RNDN); - mpfr_gamma (temp, temp, MPFR_RNDN); - mpfr_mul (tip1, tip1, temp, MPFR_RNDN); + mpfr_mul (tip1, tip1, temp1, MPFR_RNDN); mpfr_neg (tip1, tip1, MPFR_RNDN); - mpfr_div (tip1, x, tip1, MPFR_RNDN); /* tip1 = -x/(Gamma(1/3)*3^(1/3)) */ + mpfr_div (tip1, x, tip1, MPFR_RNDN); /* tip1 = -x/(Gamma (1/3)*3^(1/3)) */ mpfr_add (s, ti, tip1, MPFR_RNDN); @@ -190,18 +189,8 @@ mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) mpfr_mul (ti, ti, x3, MPFR_RNDN); mpfr_mul (tip1, tip1, x3, MPFR_RNDN); - if (k+2 <= ULONG_MAX/(k+1)) - { - mpfr_div_ui (ti, ti, k*(k+1), MPFR_RNDN); - mpfr_div_ui (tip1, tip1, (k+1)*(k+2), MPFR_RNDN); - } - else - { - mpfr_div_ui (ti, ti, k, MPFR_RNDN); - mpfr_div_ui (ti, ti, k+1, MPFR_RNDN); - mpfr_div_ui (tip1, tip1, k+1, MPFR_RNDN); - mpfr_div_ui (tip1, tip1, k+2, MPFR_RNDN); - } + mpfr_div_ui2 (ti, ti, k, (k+1), MPFR_RNDN); + mpfr_div_ui2 (tip1, tip1, (k+1), (k+2), MPFR_RNDN); k += 3; mpfr_add (s, s, ti, MPFR_RNDN); @@ -249,7 +238,7 @@ mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) { if (correct_bits < prec) { /* The precision was badly chosen */ - MPFR_LOG_MSG (("Bad assumption on the exponent of Ai(x) (E=%d)\n",MPFR_EXP (s), 0)); + MPFR_LOG_MSG (("Bad assumption on the exponent of Ai(x) (E=%d)\n", MPFR_EXP (s), 0)); wprec = prec + err + 1; } else @@ -272,7 +261,8 @@ mpfr_ai (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) mpfr_clear (ti); mpfr_clear (tip1); - mpfr_clear (temp); + mpfr_clear (temp1); + mpfr_clear (temp2); mpfr_clear (x3); mpfr_clear (s); mpfr_clear (tmp_sp); @@ -0,0 +1,323 @@ +/* mpfr_ai2 -- implementation of Airy Ai function by Smith algorithm. + +Copyright 2010 Free Software Foundation, Inc. +Contributed by the Arenaire and Cacao projects, INRIA. + +This file is part of the GNU MPFR Library. + +The GNU MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see +http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., +51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ + +#define MPFR_NEED_LONGLONG_H +#include <mpfr-impl.h> + +#define SMALL_PRECISION 32 + +/* These functions are provided by file gammaonethird.c */ +int mpfr_gamma_one_and_two_third (mpfr_ptr, mpfr_ptr, mp_prec_t); +void mpfr_div_ui2 (mpfr_ptr, mpfr_srcptr, unsigned long int, unsigned long int, mpfr_rnd_t); + + +unsigned long isqrt (unsigned long n) +{ + mpz_t tmp; + unsigned int res; + mpz_init (tmp); + mpz_set_ui (tmp, n); + mpz_sqrt (tmp, tmp); + res = mpz_get_ui (tmp); + mpz_clear (tmp); + return res; +} + +/* Reminder and notations: + ----------------------- + + Ai is the solution of: + / y'' - x*y = 0 + { Ai (0) = 1/ ( 9^(1/3) * Gamma (2/3) ) + \ Ai' (0) = -1/ ( 3^(1/3) * Gamma (1/3) ) + + Series development: + Ai (x) = sum (a_i*x^i) + = sum (t_i) + + Recurrences: + a_(i+3) = a_i / ((i+2)*(i+3)) + t_(i+3) = t_i * x^3 / ((i+2)*(i+3)) + + Values: + a_0 = Ai (0) ~ 0.355 + a_1 = Ai' (0) ~ -0.259 +*/ + + +/* Airy function Ai evaluated by Smith algorithm */ +int +mpfr_ai2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) +{ + MPFR_ZIV_DECL (loop); + MPFR_SAVE_EXPO_DECL (expo); + mpfr_prec_t wprec; /* working precision */ + mpfr_prec_t prec; /* target precision */ + mpfr_prec_t err; /* used to estimate the evaluation error */ + mpfr_prec_t correctBits; /* estimates the number of correct bits*/ + unsigned long int i, j, k, L, t; + unsigned long int cond; /* condition number of the series */ + unsigned assumedExp; /* used as a lowerbound of |EXP (Ai (x))| */ + int r; /* returned ternary value */ + mpfr_t s; /* used to store the partial sum */ + mpfr_t u0, u1; + mpfr_t *z; /* used to store the (x^3j) */ + mpfr_t result; + mpfr_t tmp_sp, tmp2_sp; /* small precision variables */ + unsigned long int x3u; /* used to store ceil (x^3) */ + mpfr_t temp1, temp2; + int test0, test1; + + /* Logging */ + MPFR_LOG_FUNC ( ("x[%#R]=%R rnd=%d", x, x, rnd), ("y[%#R]=%R", y, y) ); + + /* Special cases */ + if (MPFR_UNLIKELY (MPFR_IS_NAN (x))) { MPFR_SET_NAN (y); MPFR_RET_NAN; } + if (MPFR_UNLIKELY (MPFR_IS_INF (x))) { return mpfr_set_ui (y, 0, rnd); } + + /* Save current exponents range */ + MPFR_SAVE_EXPO_MARK (expo); + + /* FIXME: underflow for large values of |x| */ + + + /* Set initial precision */ + /* See the analysis for the naive evaluation */ + + /* We begin with 11 guard bits */ + prec = MPFR_PREC (y) + 11; + MPFR_ZIV_INIT (loop, prec); + + mpfr_init2 (tmp_sp, SMALL_PRECISION); + mpfr_init2 (tmp2_sp, SMALL_PRECISION); + mpfr_abs (tmp_sp, x, MPFR_RNDU); + mpfr_pow_ui (tmp_sp, tmp_sp, 3, MPFR_RNDU); + mpfr_sqrt (tmp_sp, tmp_sp, MPFR_RNDU); /* tmp_sp ~ x^3/2 */ + + /* 0.96179669392597567 >~ 2/3 * log2(e). See algorithms.tex */ + mpfr_mul_d (tmp2_sp, tmp_sp, 0.96179669392597567, MPFR_RNDU); + + /* cond represents the number of lost bits in the evaluation of the sum */ + if (MPFR_EXP (x) <= 0 ) cond = 0; + else cond = mpfr_get_ui (tmp2_sp, MPFR_RNDU) - (MPFR_EXP (x) - 1)/4 - 1; + + /* This variable is used to store the maximal assumed exponent of */ + /* Ai (x). More precisely, we assume that |Ai (x)| will be greater than */ + /* 2^{-assumedExp}. */ + if (MPFR_IS_POS (x)) + { + if (MPFR_EXP (x) <= 0) assumedExp = 3; + else assumedExp = 2 + (MPFR_EXP (x)/4 + 1) + mpfr_get_ui (tmp2_sp, MPFR_RNDU); + } + /* We do not know Ai (x) yet */ + /* We cover the case when EXP (Ai (x))>=-10 */ + else assumedExp = 10; + + wprec = prec + MPFR_INT_CEIL_LOG2 (prec) + 6 + cond + assumedExp; + + /* We assume that the truncation rank will be ~ prec */ + L = isqrt (prec); + MPFR_LOG_MSG (("size of blocks L = %d\n", L, 0)); + + z = (mpfr_t *) (*__gmp_allocate_func) ( (L + 1) * sizeof (mpfr_t) ); + MPFR_ASSERTN (z != NULL); + for (j=0; j<=L; j++) mpfr_init (z[j]); + + mpfr_init (s); + mpfr_init (u0); mpfr_init (u1); + mpfr_init (result); + mpfr_init (temp1); + mpfr_init (temp2); + + /* ZIV loop */ + for (;;) + { + MPFR_LOG_MSG (("working precision: %d\n", wprec, 0)); + + for (j=0; j<=L; j++) mpfr_set_prec (z[j], wprec); + mpfr_set_prec (s, wprec); + mpfr_set_prec (u0, wprec); mpfr_set_prec (u1, wprec); + mpfr_set_prec (result, wprec); + + mpfr_set_ui (u0, 1, MPFR_RNDN); + mpfr_set (u1, x, MPFR_RNDN); + + mpfr_set_ui (z[0], 1, MPFR_RNDU); + mpfr_sqr (z[1], u1, MPFR_RNDU); + mpfr_mul (z[1], z[1], x, (MPFR_IS_POS (x) ? MPFR_RNDU : MPFR_RNDD) ); + + if (MPFR_IS_NEG (x)) MPFR_CHANGE_SIGN (z[1]); + x3u = mpfr_get_ui (z[1], MPFR_RNDU); /* x3u >= ceil (x^3) */ + if (MPFR_IS_NEG (x)) MPFR_CHANGE_SIGN (z[1]); + + for (j=2; j<=L ;j++) + { + if (j%2 == 0) mpfr_sqr (z[j], z[j/2], MPFR_RNDN); + else mpfr_mul (z[j], z[j-1], z[1], MPFR_RNDN); + } + + mpfr_gamma_one_and_two_third (temp1, temp2, wprec); + mpfr_set_ui (u0, 9, MPFR_RNDN); + mpfr_cbrt (u0, u0, MPFR_RNDN); + mpfr_mul (u0, u0, temp2, MPFR_RNDN); + mpfr_ui_div (u0, 1, u0 , MPFR_RNDN); /* u0 = 1 / ( Gamma (2/3)*9^(1/3) ) */ + + mpfr_set_ui (u1, 3, MPFR_RNDN); + mpfr_cbrt (u1, u1, MPFR_RNDN); + mpfr_mul (u1, u1, temp1, MPFR_RNDN); + mpfr_neg (u1, u1, MPFR_RNDN); + mpfr_div (u1, x, u1, MPFR_RNDN); /* u1 = -x/(Gamma (1/3)*3^(1/3)) */ + + mpfr_set_ui (result, 0, MPFR_RNDN); + t = 0; + + /* Evaluation of the series by Smith' method */ + for (i=0; ; i++) + { + t += 3 * L; + + /* k = 0 */ + t -= 3; + mpfr_set (s, z[L-1], MPFR_RNDN); + for (j=L-2; ; j--) + { + t -= 3; + mpfr_div_ui2 (s, s, (t+2), (t+3), MPFR_RNDN); + mpfr_add (s, s, z[j], MPFR_RNDN); + if (j==0) break; + } + mpfr_mul (s, s, u0, MPFR_RNDN); + mpfr_add (result, result, s, MPFR_RNDN); + + mpfr_mul (u0, u0, z[L], MPFR_RNDN); + for (j=0; j<=L-1; j++) + { + mpfr_div_ui2 (u0, u0, (t + 2), (t + 3), MPFR_RNDN); + t += 3; + } + + t++; + + /* k = 1 */ + t -= 3; + mpfr_set (s, z[L-1], MPFR_RNDN); + for (j=L-2; ; j--) + { + t -= 3; + mpfr_div_ui2 (s, s, (t + 2), (t + 3), MPFR_RNDN); + mpfr_add (s, s, z[j], MPFR_RNDN); + if (j==0) break; + } + mpfr_mul (s, s, u1, MPFR_RNDN); + mpfr_add (result, result, s, MPFR_RNDN); + + mpfr_mul (u1, u1, z[L], MPFR_RNDN); + for (j=0; j<=L-1; j++) + { + mpfr_div_ui2 (u1, u1, (t + 2), (t + 3), MPFR_RNDN); + t += 3; + } + + t++; + + /* k = 2 */ + t++; + + /* End of the loop over k */ + t -= 3; + + test0 = MPFR_IS_ZERO (u0) || + (MPFR_EXP (u0) + (mp_exp_t)prec + 4 <= MPFR_EXP (result)); + test1 = MPFR_IS_ZERO (u1) || + (MPFR_EXP (u1) + (mp_exp_t)prec + 4 <= MPFR_EXP (result)); + + if ( test0 && test1 && (x3u <= (t + 2) * (t + 3) / 2) ) break; + } + + MPFR_LOG_MSG (("Truncation rank: %d\n", t, 0)); + + err = 5 + MPFR_INT_CEIL_LOG2 (L+1) + MPFR_INT_CEIL_LOG2 (i+1) + cond - MPFR_EXP (result); + + /* err is the number of bits lost due to the evaluation error */ + /* wprec-(prec+1): the number of bits lost due to the approximation error */ + MPFR_LOG_MSG (("Roundoff error: %d\n", err, 0)); + MPFR_LOG_MSG (("Approxim error: %d\n", wprec - prec - 1, 0)); + + if (wprec < err+1) correctBits = 0; + else + { + if (wprec < err+prec+1) correctBits = (unsigned) ( (signed)wprec - (signed)err - 1); + else correctBits = prec; + } + + if (MPFR_LIKELY (MPFR_CAN_ROUND (result, correctBits, MPFR_PREC (y), rnd))) break; + + for (j=0; j<=L; j++) mpfr_clear (z[j]); + (*__gmp_free_func) (z, (L + 1) * sizeof (mpfr_t)); + L = isqrt (t); + MPFR_LOG_MSG (("size of blocks L = %d\n", L, 0)); + z = (mpfr_t *) (*__gmp_allocate_func) ( (L + 1) * sizeof (mpfr_t)); + MPFR_ASSERTN (z != NULL); + for (j=0; j<=L; j++) mpfr_init (z[j]); + + if (correctBits == 0) + { + assumedExp *= 2; + MPFR_LOG_MSG (("Not a single bit correct (assumedExp=%d)\n", 0)); + wprec = prec + 6 + MPFR_INT_CEIL_LOG2 (t) + cond + assumedExp; + } + else + { + if (correctBits < prec) + { /* The precision was badly chosen */ + MPFR_LOG_MSG (("Bad assumption on the exponent of Ai (x) (E=%d)\n", MPFR_EXP (result), 0)); + wprec = prec + err + 1; + } + else + { /* We are really in a bad case of the TMD */ + MPFR_ZIV_NEXT (loop, prec); + + /* We update wprec */ + /* We assume that t will not be multiplied by more than 4 */ + wprec = prec + (MPFR_INT_CEIL_LOG2 (t) + 2) + 6 + cond - MPFR_EXP (result); + } + } + } /* End of ZIV loop */ + + MPFR_ZIV_FREE (loop); + MPFR_SAVE_EXPO_FREE (expo); + + r = mpfr_set (y, result, rnd); + + mpfr_clear (tmp_sp); + mpfr_clear (tmp2_sp); + for (j=0; j<=L; j++) mpfr_clear (z[j]); + (*__gmp_free_func) (z, (L + 1) * sizeof (mpfr_t)); + + mpfr_clear (s); + mpfr_clear (u0); mpfr_clear (u1); + mpfr_clear (result); + mpfr_clear (temp1); + mpfr_clear (temp2); + + return r; +} diff --git a/gammaonethird.c b/gammaonethird.c new file mode 100644 index 000000000..628471e71 --- /dev/null +++ b/gammaonethird.c @@ -0,0 +1,174 @@ +/* Functions for evaluating Gamma(1/3) and Gamma(2/3). Used by mpfr_ai. + +Copyright 2010 Free Software Foundation, Inc. +Contributed by the Arenaire and Cacao projects, INRIA. + +This file is part of the GNU MPFR Library. + +The GNU MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see +http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., +51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ + +#define MPFR_NEED_LONGLONG_H +#include <mpfr-impl.h> + +#define MPFR_ACC_OR_MUL(v) \ + if (v <= ULONG_MAX / acc) acc *= v; \ + else { mpfr_mul_ui (y, y, acc, mode); acc = v; } + +#define MPFR_ACC_OR_DIV(v) \ + if (v <= ULONG_MAX / acc) acc *= v; \ + else { mpfr_div_ui (y, y, acc, mode); acc = v; } + +void +mpfr_mul_ui5 (mpfr_ptr y, mpfr_srcptr x, + unsigned long int v1, unsigned long int v2, + unsigned long int v3, unsigned long int v4, + unsigned long int v5, mpfr_rnd_t mode) +{ + unsigned long int acc = v1; + mpfr_set (y, x, mode); + MPFR_ACC_OR_MUL (v2); + MPFR_ACC_OR_MUL (v3); + MPFR_ACC_OR_MUL (v4); + MPFR_ACC_OR_MUL (v5); + mpfr_mul_ui (y, y, acc, mode); +} + +void +mpfr_div_ui2 (mpfr_ptr y, mpfr_srcptr x, + unsigned long int v1, unsigned long int v2, mpfr_rnd_t mode) +{ + unsigned long int acc = v1; + mpfr_set (y, x, mode); + MPFR_ACC_OR_DIV (v2); + mpfr_div_ui (y, y, acc, mode); +} + +void +mpfr_div_ui8 (mpfr_ptr y, mpfr_srcptr x, + unsigned long int v1, unsigned long int v2, + unsigned long int v3, unsigned long int v4, + unsigned long int v5, unsigned long int v6, + unsigned long int v7, unsigned long int v8, mpfr_rnd_t mode) +{ + unsigned long int acc = v1; + mpfr_set (y, x, mode); + MPFR_ACC_OR_DIV (v2); + MPFR_ACC_OR_DIV (v3); + MPFR_ACC_OR_DIV (v4); + MPFR_ACC_OR_DIV (v5); + MPFR_ACC_OR_DIV (v6); + MPFR_ACC_OR_DIV (v7); + MPFR_ACC_OR_DIV (v8); + mpfr_div_ui (y, y, acc, mode); +} + + +/* Gives an approximation of omega = Gamma(1/3)^6 * sqrt(10) / (12pi^4) */ +/* using C. H. Brown's formula. */ +/* The computed value s satisfies |s-omega| <= 2^{1-prec}*omega */ +/* As usual, the variable s is supposed to be initialized. */ +mpfr_Browns_const (mpfr_ptr s, mp_prec_t prec) +{ + mpfr_t uk; + unsigned long int k; + + mp_prec_t working_prec = prec + 10 + MPFR_INT_CEIL_LOG2 (2 + prec / 10); + + mpfr_init2 (uk, working_prec); + mpfr_set_prec (s, working_prec); + + mpfr_set_ui (uk, 1, MPFR_RNDN); + mpfr_set (s, uk, MPFR_RNDN); + k = 1; + + /* Invariants: uk ~ u(k-1) and s ~ sum(i=0..k-1, u(i)) */ + for (;;) + { + mpfr_mul_ui5 (uk, uk, 6 * k - 5, 6 * k - 4, 6 * k - 3, 6 * k - 2, + 6 * k - 1, MPFR_RNDN); + mpfr_div_ui8 (uk, uk, k, k, 3 * k - 2, 3 * k - 1, 3 * k, 80, 160, 160, + MPFR_RNDN); + MPFR_CHANGE_SIGN (uk); + + mpfr_add (s, s, uk, MPFR_RNDN); + k++; + if (MPFR_EXP (uk) + (signed) prec <= MPFR_EXP (s) + 7) + break; + } + + mpfr_clear (uk); + return; +} + +/* Returns y such that |Gamma(1/3)-y| <= 2^{1-prec}*Gamma(1/3) */ +int +mpfr_gamma_one_third (mpfr_ptr y, mp_prec_t prec) +{ + mpfr_t tmp, tmp2, tmp3; + + mpfr_init2 (tmp, prec + 9); + mpfr_init2 (tmp2, prec + 9); + mpfr_init2 (tmp3, prec + 4); + mpfr_set_prec (y, prec + 2); + + mpfr_const_pi (tmp, MPFR_RNDN); + mpfr_sqr (tmp, tmp, MPFR_RNDN); + mpfr_sqr (tmp, tmp, MPFR_RNDN); + mpfr_mul_ui (tmp, tmp, 12, MPFR_RNDN); + + mpfr_Browns_const (tmp2, prec + 9); + mpfr_mul (tmp, tmp, tmp2, MPFR_RNDN); + + mpfr_set_ui (tmp2, 10, MPFR_RNDN); + mpfr_sqrt (tmp2, tmp2, MPFR_RNDN); + mpfr_div (tmp, tmp, tmp2, MPFR_RNDN); + + mpfr_sqrt (tmp3, tmp, MPFR_RNDN); + mpfr_cbrt (y, tmp3, MPFR_RNDN); + + mpfr_clear (tmp); + mpfr_clear (tmp2); + mpfr_clear (tmp3); + return; +} + +/* Computes y1 and y2 such that: */ +/* |y1-Gamma(1/3)| <= 2^{1-prec}Gamma(1/3) */ +/* and |y2-Gamma(2/3)| <= 2^{1-prec}Gamma(2/3) */ +/* */ +/* Uses the formula Gamma(z)Gamma(1-z) = pi / sin(pi*z) */ +/* to compute Gamma(2/3) from Gamma(1/3). */ +int +mpfr_gamma_one_and_two_third (mpfr_ptr y1, mpfr_ptr y2, mp_prec_t prec) +{ + mpfr_t temp; + + mpfr_init2 (temp, prec + 4); + mpfr_set_prec (y2, prec + 4); + + mpfr_gamma_one_third (y1, prec + 4); + + mpfr_set_ui (temp, 3, MPFR_RNDN); + mpfr_sqrt (temp, temp, MPFR_RNDN); + mpfr_mul (temp, y1, temp, MPFR_RNDN); + + mpfr_const_pi (y2, MPFR_RNDN); + mpfr_mul_2ui (y2, y2, 1, MPFR_RNDN); + + mpfr_div (y2, y2, temp, MPFR_RNDN); + + mpfr_clear (temp); +} @@ -659,6 +659,7 @@ __MPFR_DECLSPEC int mpfr_yn _MPFR_PROTO ((mpfr_ptr, long, mpfr_srcptr, mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_ai _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_rnd_t)); +__MPFR_DECLSPEC int mpfr_ai2 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_min _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t)); |