contrib/mpfr/src/digamma.c

   1 /* mpfr_digamma -- digamma function of a floating-point number
   2
   3 Copyright 2009, 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 #include "mpfr-impl.h"
  24
  25 /* Put in s an approximation of digamma(x).
  26    Assumes x >= 2.
  27    Assumes s does not overlap with x.
  28    Returns an integer e such that the error is bounded by 2^e ulps
  29    of the result s.
  30 */
  31 static mpfr_exp_t
  32 mpfr_digamma_approx (mpfr_ptr s, mpfr_srcptr x)
  33 {
  34   mpfr_prec_t p = MPFR_PREC (s);
  35   mpfr_t t, u, invxx;
  36   mpfr_exp_t e, exps, f, expu;
  37   mpz_t *INITIALIZED(B);  /* variable B declared as initialized */
  38   unsigned long n0, n; /* number of allocated B[] */
  39
  40   MPFR_ASSERTN(MPFR_IS_POS(x) && (MPFR_EXP(x) >= 2));
  41
  42   mpfr_init2 (t, p);
  43   mpfr_init2 (u, p);
  44   mpfr_init2 (invxx, p);
  45
  46   mpfr_log (s, x, MPFR_RNDN);         /* error <= 1/2 ulp */
  47   mpfr_ui_div (t, 1, x, MPFR_RNDN);   /* error <= 1/2 ulp */
  48   mpfr_div_2exp (t, t, 1, MPFR_RNDN); /* exact */
  49   mpfr_sub (s, s, t, MPFR_RNDN);
  50   /* error <= 1/2 + 1/2*2^(EXP(olds)-EXP(s)) + 1/2*2^(EXP(t)-EXP(s)).
  51      For x >= 2, log(x) >= 2*(1/(2x)), thus olds >= 2t, and olds - t >= olds/2,
  52      thus 0 <= EXP(olds)-EXP(s) <= 1, and EXP(t)-EXP(s) <= 0, thus
  53      error <= 1/2 + 1/2*2 + 1/2 <= 2 ulps. */
  54   e = 2; /* initial error */
  55   mpfr_mul (invxx, x, x, MPFR_RNDZ);     /* invxx = x^2 * (1 + theta)
  56                                             for |theta| <= 2^(-p) */
  57   mpfr_ui_div (invxx, 1, invxx, MPFR_RNDU); /* invxx = 1/x^2 * (1 + theta)^2 */
  58
  59   /* in the following we note err=xxx when the ratio between the approximation
  60      and the exact result can be written (1 + theta)^xxx for |theta| <= 2^(-p),
  61      following Higham's method */
  62   B = mpfr_bernoulli_internal ((mpz_t *) 0, 0);
  63   mpfr_set_ui (t, 1, MPFR_RNDN); /* err = 0 */
  64   for (n = 1;; n++)
  65     {
  66       /* compute next Bernoulli number */
  67       B = mpfr_bernoulli_internal (B, n);
  68       /* The main term is Bernoulli[2n]/(2n)/x^(2n) = B[n]/(2n+1)!(2n)/x^(2n)
  69          = B[n]*t[n]/(2n) where t[n]/t[n-1] = 1/(2n)/(2n+1)/x^2. */
  70       mpfr_mul (t, t, invxx, MPFR_RNDU);        /* err = err + 3 */
  71       mpfr_div_ui (t, t, 2 * n, MPFR_RNDU);     /* err = err + 1 */
  72       mpfr_div_ui (t, t, 2 * n + 1, MPFR_RNDU); /* err = err + 1 */
  73       /* we thus have err = 5n here */
  74       mpfr_div_ui (u, t, 2 * n, MPFR_RNDU);     /* err = 5n+1 */
  75       mpfr_mul_z (u, u, B[n], MPFR_RNDU);       /* err = 5n+2, and the
  76                                                    absolute error is bounded
  77                                                    by 10n+4 ulp(u) [Rule 11] */
  78       /* if the terms 'u' are decreasing by a factor two at least,
  79          then the error coming from those is bounded by
  80          sum((10n+4)/2^n, n=1..infinity) = 24 */
  81       exps = mpfr_get_exp (s);
  82       expu = mpfr_get_exp (u);
  83       if (expu < exps - (mpfr_exp_t) p)
  84         break;
  85       mpfr_sub (s, s, u, MPFR_RNDN); /* error <= 24 + n/2 */
  86       if (mpfr_get_exp (s) < exps)
  87         e <<= exps - mpfr_get_exp (s);
  88       e ++; /* error in mpfr_sub */
  89       f = 10 * n + 4;
  90       while (expu < exps)
  91         {
  92           f = (1 + f) / 2;
  93           expu ++;
  94         }
  95       e += f; /* total rouding error coming from 'u' term */
  96     }
  97
  98   n0 = ++n;
  99   while (n--)
 100     mpz_clear (B[n]);
 101   (*__gmp_free_func) (B, n0 * sizeof (mpz_t));
 102
 103   mpfr_clear (t);
 104   mpfr_clear (u);
 105   mpfr_clear (invxx);
 106
 107   f = 0;
 108   while (e > 1)
 109     {
 110       f++;
 111       e = (e + 1) / 2;
 112       /* Invariant: 2^f * e does not decrease */
 113     }
 114   return f;
 115 }
 116
 117 /* Use the reflection formula Digamma(1-x) = Digamma(x) + Pi * cot(Pi*x),
 118    i.e., Digamma(x) = Digamma(1-x) - Pi * cot(Pi*x).
 119    Assume x < 1/2. */
 120 static int
 121 mpfr_digamma_reflection (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
 122 {
 123   mpfr_prec_t p = MPFR_PREC(y) + 10, q;
 124   mpfr_t t, u, v;
 125   mpfr_exp_t e1, expv;
 126   int inex;
 127   MPFR_ZIV_DECL (loop);
 128
 129   /* we want that 1-x is exact with precision q: if 0 < x < 1/2, then
 130      q = PREC(x)-EXP(x) is ok, otherwise if -1 <= x < 0, q = PREC(x)-EXP(x)
 131      is ok, otherwise for x < -1, PREC(x) is ok if EXP(x) <= PREC(x),
 132      otherwise we need EXP(x) */
 133   if (MPFR_EXP(x) < 0)
 134     q = MPFR_PREC(x) + 1 - MPFR_EXP(x);
 135   else if (MPFR_EXP(x) <= MPFR_PREC(x))
 136     q = MPFR_PREC(x) + 1;
 137   else
 138     q = MPFR_EXP(x);
 139   mpfr_init2 (u, q);
 140   MPFR_ASSERTN(mpfr_ui_sub (u, 1, x, MPFR_RNDN) == 0);
 141
 142   /* if x is half an integer, cot(Pi*x) = 0, thus Digamma(x) = Digamma(1-x) */
 143   mpfr_mul_2exp (u, u, 1, MPFR_RNDN);
 144   inex = mpfr_integer_p (u);
 145   mpfr_div_2exp (u, u, 1, MPFR_RNDN);
 146   if (inex)
 147     {
 148       inex = mpfr_digamma (y, u, rnd_mode);
 149       goto end;
 150     }
 151
 152   mpfr_init2 (t, p);
 153   mpfr_init2 (v, p);
 154
 155   MPFR_ZIV_INIT (loop, p);
 156   for (;;)
 157     {
 158       mpfr_const_pi (v, MPFR_RNDN);  /* v = Pi*(1+theta) for |theta|<=2^(-p) */
 159       mpfr_mul (t, v, x, MPFR_RNDN); /* (1+theta)^2 */
 160       e1 = MPFR_EXP(t) - (mpfr_exp_t) p + 1; /* bound for t: err(t) <= 2^e1 */
 161       mpfr_cot (t, t, MPFR_RNDN);
 162       /* cot(t * (1+h)) = cot(t) - theta * (1 + cot(t)^2) with |theta|<=t*h */
 163       if (MPFR_EXP(t) > 0)
 164         e1 = e1 + 2 * MPFR_EXP(t) + 1;
 165       else
 166         e1 = e1 + 1;
 167       /* now theta * (1 + cot(t)^2) <= 2^e1 */
 168       e1 += (mpfr_exp_t) p - MPFR_EXP(t); /* error is now 2^e1 ulps */
 169       mpfr_mul (t, t, v, MPFR_RNDN);
 170       e1 ++;
 171       mpfr_digamma (v, u, MPFR_RNDN);   /* error <= 1/2 ulp */
 172       expv = MPFR_EXP(v);
 173       mpfr_sub (v, v, t, MPFR_RNDN);
 174       if (MPFR_EXP(v) < MPFR_EXP(t))
 175         e1 += MPFR_EXP(t) - MPFR_EXP(v); /* scale error for t wrt new v */
 176       /* now take into account the 1/2 ulp error for v */
 177       if (expv - MPFR_EXP(v) - 1 > e1)
 178         e1 = expv - MPFR_EXP(v) - 1;
 179       else
 180         e1 ++;
 181       e1 ++; /* rounding error for mpfr_sub */
 182       if (MPFR_CAN_ROUND (v, p - e1, MPFR_PREC(y), rnd_mode))
 183         break;
 184       MPFR_ZIV_NEXT (loop, p);
 185       mpfr_set_prec (t, p);
 186       mpfr_set_prec (v, p);
 187     }
 188   MPFR_ZIV_FREE (loop);
 189
 190   inex = mpfr_set (y, v, rnd_mode);
 191
 192   mpfr_clear (t);
 193   mpfr_clear (v);
 194  end:
 195   mpfr_clear (u);
 196
 197   return inex;
 198 }
 199
 200 /* we have x >= 1/2 here */
 201 static int
 202 mpfr_digamma_positive (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
 203 {
 204   mpfr_prec_t p = MPFR_PREC(y) + 10, q;
 205   mpfr_t t, u, x_plus_j;
 206   int inex;
 207   mpfr_exp_t errt, erru, expt;
 208   unsigned long j = 0, min;
 209   MPFR_ZIV_DECL (loop);
 210
 211   /* compute a precision q such that x+1 is exact */
 212   if (MPFR_PREC(x) < MPFR_EXP(x))
 213     q = MPFR_EXP(x);
 214   else
 215     q = MPFR_PREC(x) + 1;
 216   mpfr_init2 (x_plus_j, q);
 217
 218   mpfr_init2 (t, p);
 219   mpfr_init2 (u, p);
 220   MPFR_ZIV_INIT (loop, p);
 221   for(;;)
 222     {
 223       /* Lower bound for x+j in mpfr_digamma_approx call: since the smallest
 224          term of the divergent series for Digamma(x) is about exp(-2*Pi*x), and
 225          we want it to be less than 2^(-p), this gives x > p*log(2)/(2*Pi)
 226          i.e., x >= 0.1103 p.
 227          To be safe, we ensure x >= 0.25 * p.
 228       */
 229       min = (p + 3) / 4;
 230       if (min < 2)
 231         min = 2;
 232
 233       mpfr_set (x_plus_j, x, MPFR_RNDN);
 234       mpfr_set_ui (u, 0, MPFR_RNDN);
 235       j = 0;
 236       while (mpfr_cmp_ui (x_plus_j, min) < 0)
 237         {
 238           j ++;
 239           mpfr_ui_div (t, 1, x_plus_j, MPFR_RNDN); /* err <= 1/2 ulp */
 240           mpfr_add (u, u, t, MPFR_RNDN);
 241           inex = mpfr_add_ui (x_plus_j, x_plus_j, 1, MPFR_RNDZ);
 242           if (inex != 0) /* we lost one bit */
 243             {
 244               q ++;
 245               mpfr_prec_round (x_plus_j, q, MPFR_RNDZ);
 246               mpfr_nextabove (x_plus_j);
 247             }
 248           /* since all terms are positive, the error is bounded by j ulps */
 249         }
 250       for (erru = 0; j > 1; erru++, j = (j + 1) / 2);
 251       errt = mpfr_digamma_approx (t, x_plus_j);
 252       expt = MPFR_EXP(t);
 253       mpfr_sub (t, t, u, MPFR_RNDN);
 254       if (MPFR_EXP(t) < expt)
 255         errt += expt - MPFR_EXP(t);
 256       if (MPFR_EXP(t) < MPFR_EXP(u))
 257         erru += MPFR_EXP(u) - MPFR_EXP(t);
 258       if (errt > erru)
 259         errt = errt + 1;
 260       else if (errt == erru)
 261         errt = errt + 2;
 262       else
 263         errt = erru + 1;
 264       if (MPFR_CAN_ROUND (t, p - errt, MPFR_PREC(y), rnd_mode))
 265         break;
 266       MPFR_ZIV_NEXT (loop, p);
 267       mpfr_set_prec (t, p);
 268       mpfr_set_prec (u, p);
 269     }
 270   MPFR_ZIV_FREE (loop);
 271   inex = mpfr_set (y, t, rnd_mode);
 272   mpfr_clear (t);
 273   mpfr_clear (u);
 274   mpfr_clear (x_plus_j);
 275   return inex;
 276 }
 277
 278 int
 279 mpfr_digamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
 280 {
 281   int inex;
 282   MPFR_SAVE_EXPO_DECL (expo);
 283
 284   MPFR_LOG_FUNC
 285     (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode),
 286      ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y, inex));
 287
 288
 289   if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)))
 290     {
 291       if (MPFR_IS_NAN(x))
 292         {
 293           MPFR_SET_NAN(y);
 294           MPFR_RET_NAN;
 295         }
 296       else if (MPFR_IS_INF(x))
 297         {
 298           if (MPFR_IS_POS(x)) /* Digamma(+Inf) = +Inf */
 299             {
 300               MPFR_SET_SAME_SIGN(y, x);
 301               MPFR_SET_INF(y);
 302               MPFR_RET(0);
 303             }
 304           else                /* Digamma(-Inf) = NaN */
 305             {
 306               MPFR_SET_NAN(y);
 307               MPFR_RET_NAN;
 308             }
 309         }
 310       else /* Zero case */
 311         {
 312           /* the following works also in case of overlap */
 313           MPFR_SET_INF(y);
 314           MPFR_SET_OPPOSITE_SIGN(y, x);
 315           mpfr_set_divby0 ();
 316           MPFR_RET(0);
 317         }
 318     }
 319
 320   /* Digamma is undefined for negative integers */
 321   if (MPFR_IS_NEG(x) && mpfr_integer_p (x))
 322     {
 323       MPFR_SET_NAN(y);
 324       MPFR_RET_NAN;
 325     }
 326
 327   /* now x is a normal number */
 328
 329   MPFR_SAVE_EXPO_MARK (expo);
 330   /* for x very small, we have Digamma(x) = -1/x - gamma + O(x), more precisely
 331      -1 < Digamma(x) + 1/x < 0 for -0.2 < x < 0.2, thus:
 332      (i) either x is a power of two, then 1/x is exactly representable, and
 333          as long as 1/2*ulp(1/x) > 1, we can conclude;
 334      (ii) otherwise assume x has <= n bits, and y has <= n+1 bits, then
 335    |y + 1/x| >= 2^(-2n) ufp(y), where ufp means unit in first place.
 336    Since |Digamma(x) + 1/x| <= 1, if 2^(-2n) ufp(y) >= 2, then
 337    |y - Digamma(x)| >= 2^(-2n-1)ufp(y), and rounding -1/x gives the correct result.
 338    If x < 2^E, then y > 2^(-E), thus ufp(y) > 2^(-E-1).
 339    A sufficient condition is thus EXP(x) <= -2 MAX(PREC(x),PREC(Y)). */
 340   if (MPFR_EXP(x) < -2)
 341     {
 342       if (MPFR_EXP(x) <= -2 * (mpfr_exp_t) MAX(MPFR_PREC(x), MPFR_PREC(y)))
 343         {
 344           int signx = MPFR_SIGN(x);
 345           inex = mpfr_si_div (y, -1, x, rnd_mode);
 346           if (inex == 0) /* x is a power of two */
 347             { /* result always -1/x, except when rounding down */
 348               if (rnd_mode == MPFR_RNDA)
 349                 rnd_mode = (signx > 0) ? MPFR_RNDD : MPFR_RNDU;
 350               if (rnd_mode == MPFR_RNDZ)
 351                 rnd_mode = (signx > 0) ? MPFR_RNDU : MPFR_RNDD;
 352               if (rnd_mode == MPFR_RNDU)
 353                 inex = 1;
 354               else if (rnd_mode == MPFR_RNDD)
 355                 {
 356                   mpfr_nextbelow (y);
 357                   inex = -1;
 358                 }
 359               else /* nearest */
 360                 inex = 1;
 361             }
 362           MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 363           goto end;
 364         }
 365     }
 366
 367   if (MPFR_IS_NEG(x))
 368     inex = mpfr_digamma_reflection (y, x, rnd_mode);
 369   /* if x < 1/2 we use the reflection formula */
 370   else if (MPFR_EXP(x) < 0)
 371     inex = mpfr_digamma_reflection (y, x, rnd_mode);
 372   else
 373     inex = mpfr_digamma_positive (y, x, rnd_mode);
 374
 375  end:
 376   MPFR_SAVE_EXPO_FREE (expo);
 377   return mpfr_check_range (y, inex, rnd_mode);
 378 }