contrib/mpfr/src/gammaonethird.c

   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 }