"Fossies" - the Fresh Open Source Software Archive

Member "mpfr-4.0.2/src/digamma.c" (7 Jan 2019, 12906 Bytes) of package /linux/misc/mpfr-4.0.2.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "digamma.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 4.0.1_vs_4.0.2.

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