1 /* Functions for evaluating Gamma(1/3) and Gamma(2/3). Used by mpfr_ai. 2 3 Copyright 2010, 2011, 2012, 2013 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 #define MPFR_ACC_OR_MUL(v) \ 27 do \ 28 { \ 29 if (v <= ULONG_MAX / acc) \ 30 acc *= v; \ 31 else \ 32 { \ 33 mpfr_mul_ui (y, y, acc, mode); acc = v; \ 34 } \ 35 } \ 36 while (0) 37 38 #define MPFR_ACC_OR_DIV(v) \ 39 do \ 40 { \ 41 if (v <= ULONG_MAX / acc) \ 42 acc *= v; \ 43 else \ 44 { \ 45 mpfr_div_ui (y, y, acc, mode); acc = v; \ 46 } \ 47 } \ 48 while (0) 49 50 static void 51 mpfr_mul_ui5 (mpfr_ptr y, mpfr_srcptr x, 52 unsigned long int v1, unsigned long int v2, 53 unsigned long int v3, unsigned long int v4, 54 unsigned long int v5, mpfr_rnd_t mode) 55 { 56 unsigned long int acc = v1; 57 mpfr_set (y, x, mode); 58 MPFR_ACC_OR_MUL (v2); 59 MPFR_ACC_OR_MUL (v3); 60 MPFR_ACC_OR_MUL (v4); 61 MPFR_ACC_OR_MUL (v5); 62 mpfr_mul_ui (y, y, acc, mode); 63 } 64 65 void 66 mpfr_div_ui2 (mpfr_ptr y, mpfr_srcptr x, 67 unsigned long int v1, unsigned long int v2, mpfr_rnd_t mode) 68 { 69 unsigned long int acc = v1; 70 mpfr_set (y, x, mode); 71 MPFR_ACC_OR_DIV (v2); 72 mpfr_div_ui (y, y, acc, mode); 73 } 74 75 static void 76 mpfr_div_ui8 (mpfr_ptr y, mpfr_srcptr x, 77 unsigned long int v1, unsigned long int v2, 78 unsigned long int v3, unsigned long int v4, 79 unsigned long int v5, unsigned long int v6, 80 unsigned long int v7, unsigned long int v8, mpfr_rnd_t mode) 81 { 82 unsigned long int acc = v1; 83 mpfr_set (y, x, mode); 84 MPFR_ACC_OR_DIV (v2); 85 MPFR_ACC_OR_DIV (v3); 86 MPFR_ACC_OR_DIV (v4); 87 MPFR_ACC_OR_DIV (v5); 88 MPFR_ACC_OR_DIV (v6); 89 MPFR_ACC_OR_DIV (v7); 90 MPFR_ACC_OR_DIV (v8); 91 mpfr_div_ui (y, y, acc, mode); 92 } 93 94 95 /* Gives an approximation of omega = Gamma(1/3)^6 * sqrt(10) / (12pi^4) */ 96 /* using C. H. Brown's formula. */ 97 /* The computed value s satisfies |s-omega| <= 2^{1-prec}*omega */ 98 /* As usual, the variable s is supposed to be initialized. */ 99 static void 100 mpfr_Browns_const (mpfr_ptr s, mpfr_prec_t prec) 101 { 102 mpfr_t uk; 103 unsigned long int k; 104 105 mpfr_prec_t working_prec = prec + 10 + MPFR_INT_CEIL_LOG2 (2 + prec / 10); 106 107 mpfr_init2 (uk, working_prec); 108 mpfr_set_prec (s, working_prec); 109 110 mpfr_set_ui (uk, 1, MPFR_RNDN); 111 mpfr_set (s, uk, MPFR_RNDN); 112 k = 1; 113 114 /* Invariants: uk ~ u(k-1) and s ~ sum(i=0..k-1, u(i)) */ 115 for (;;) 116 { 117 mpfr_mul_ui5 (uk, uk, 6 * k - 5, 6 * k - 4, 6 * k - 3, 6 * k - 2, 118 6 * k - 1, MPFR_RNDN); 119 mpfr_div_ui8 (uk, uk, k, k, 3 * k - 2, 3 * k - 1, 3 * k, 80, 160, 160, 120 MPFR_RNDN); 121 MPFR_CHANGE_SIGN (uk); 122 123 mpfr_add (s, s, uk, MPFR_RNDN); 124 k++; 125 if (MPFR_GET_EXP (uk) + prec <= MPFR_GET_EXP (s) + 7) 126 break; 127 } 128 129 mpfr_clear (uk); 130 return; 131 } 132 133 /* Returns y such that |Gamma(1/3)-y| <= 2^{1-prec}*Gamma(1/3) */ 134 static void 135 mpfr_gamma_one_third (mpfr_ptr y, mpfr_prec_t prec) 136 { 137 mpfr_t tmp, tmp2, tmp3; 138 139 mpfr_init2 (tmp, prec + 9); 140 mpfr_init2 (tmp2, prec + 9); 141 mpfr_init2 (tmp3, prec + 4); 142 mpfr_set_prec (y, prec + 2); 143 144 mpfr_const_pi (tmp, MPFR_RNDN); 145 mpfr_sqr (tmp, tmp, MPFR_RNDN); 146 mpfr_sqr (tmp, tmp, MPFR_RNDN); 147 mpfr_mul_ui (tmp, tmp, 12, MPFR_RNDN); 148 149 mpfr_Browns_const (tmp2, prec + 9); 150 mpfr_mul (tmp, tmp, tmp2, MPFR_RNDN); 151 152 mpfr_set_ui (tmp2, 10, MPFR_RNDN); 153 mpfr_sqrt (tmp2, tmp2, MPFR_RNDN); 154 mpfr_div (tmp, tmp, tmp2, MPFR_RNDN); 155 156 mpfr_sqrt (tmp3, tmp, MPFR_RNDN); 157 mpfr_cbrt (y, tmp3, MPFR_RNDN); 158 159 mpfr_clear (tmp); 160 mpfr_clear (tmp2); 161 mpfr_clear (tmp3); 162 return; 163 } 164 165 /* Computes y1 and y2 such that: */ 166 /* |y1-Gamma(1/3)| <= 2^{1-prec}Gamma(1/3) */ 167 /* and |y2-Gamma(2/3)| <= 2^{1-prec}Gamma(2/3) */ 168 /* */ 169 /* Uses the formula Gamma(z)Gamma(1-z) = pi / sin(pi*z) */ 170 /* to compute Gamma(2/3) from Gamma(1/3). */ 171 void 172 mpfr_gamma_one_and_two_third (mpfr_ptr y1, mpfr_ptr y2, mpfr_prec_t prec) 173 { 174 mpfr_t temp; 175 176 mpfr_init2 (temp, prec + 4); 177 mpfr_set_prec (y2, prec + 4); 178 179 mpfr_gamma_one_third (y1, prec + 4); 180 181 mpfr_set_ui (temp, 3, MPFR_RNDN); 182 mpfr_sqrt (temp, temp, MPFR_RNDN); 183 mpfr_mul (temp, y1, temp, MPFR_RNDN); 184 185 mpfr_const_pi (y2, MPFR_RNDN); 186 mpfr_mul_2ui (y2, y2, 1, MPFR_RNDN); 187 188 mpfr_div (y2, y2, temp, MPFR_RNDN); 189 190 mpfr_clear (temp); 191 } 192