118334Speter/* mpfr_li2 -- Dilogarithm. 218334Speter 372562SobrienCopyright 2007-2023 Free Software Foundation, Inc. 490075SobrienContributed by the AriC and Caramba projects, INRIA. 518334Speter 690075SobrienThis file is part of the GNU MPFR Library. 718334Speter 890075SobrienThe GNU MPFR Library is free software; you can redistribute it and/or modify 990075Sobrienit under the terms of the GNU Lesser General Public License as published by 1090075Sobrienthe Free Software Foundation; either version 3 of the License, or (at your 1190075Sobrienoption) any later version. 1218334Speter 1390075SobrienThe GNU MPFR Library is distributed in the hope that it will be useful, but 1490075SobrienWITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 1590075Sobrienor FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 1690075SobrienLicense for more details. 1718334Speter 1818334SpeterYou should have received a copy of the GNU Lesser General Public License 1990075Sobrienalong with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 2090075Sobrienhttps://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2190075Sobrien51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 2218334Speter 2318334Speter#define MPFR_NEED_LONGLONG_H 2418334Speter#include "mpfr-impl.h" 2550397Sobrien 2652284Sobrien/* Compute the alternating series 2718334Speter s = S(z) = \sum_{k=0}^infty B_{2k} (z))^{2k+1} / (2k+1)! 2818334Speter with 0 < z <= log(2) to the precision of s rounded in the direction 2990075Sobrien rnd_mode. 3018334Speter Return the maximum index of the truncature which is useful 3118334Speter for determinating the relative error. 3218334Speter*/ 3390075Sobrienstatic int 3418334Speterli2_series (mpfr_ptr sum, mpfr_srcptr z, mpfr_rnd_t rnd_mode) 3518334Speter{ 3618334Speter int i; 3790075Sobrien mpfr_t s, u, v, w; 3890075Sobrien mpfr_prec_t sump, p; 3990075Sobrien mpfr_exp_t se, err; 4090075Sobrien MPFR_ZIV_DECL (loop); 4190075Sobrien 4290075Sobrien /* The series converges for |z| < 2 pi, but in mpfr_li2 the argument is 4390075Sobrien reduced so that 0 < z <= log(2). Here is an additional check that z 4490075Sobrien is (nearly) correct. */ 4590075Sobrien MPFR_ASSERTD (MPFR_IS_STRICTPOS (z)); 4690075Sobrien MPFR_ASSERTD (mpfr_cmp_ui_2exp (z, 89, -7) <= 0); /* z <= 0.6953125 */ 4790075Sobrien 4890075Sobrien sump = MPFR_PREC (sum); /* target precision */ 4990075Sobrien p = sump + MPFR_INT_CEIL_LOG2 (sump) + 4; /* the working precision */ 5090075Sobrien mpfr_init2 (s, p); 5190075Sobrien mpfr_init2 (u, p); 5290075Sobrien mpfr_init2 (v, p); 5390075Sobrien mpfr_init2 (w, p); 5490075Sobrien 5518334Speter MPFR_ZIV_INIT (loop, p); 5618334Speter for (;;) 5796263Sobrien { 5818334Speter mpfr_sqr (u, z, MPFR_RNDU); 5918334Speter mpfr_set (v, z, MPFR_RNDU); 6018334Speter mpfr_set (s, z, MPFR_RNDU); 6118334Speter se = MPFR_GET_EXP (s); 6218334Speter err = 0; 6318334Speter 6490075Sobrien for (i = 1;; i++) 6518334Speter { 6618334Speter mpfr_mul (v, u, v, MPFR_RNDU); 6718334Speter mpfr_div_ui (v, v, 2 * i, MPFR_RNDU); 6818334Speter mpfr_div_ui (v, v, 2 * i, MPFR_RNDU); 6918334Speter mpfr_div_ui (v, v, 2 * i + 1, MPFR_RNDU); 7018334Speter mpfr_div_ui (v, v, 2 * i + 1, MPFR_RNDU); 7118334Speter /* here, v_2i = v_{2i-2} / (2i * (2i+1))^2 */ 7218334Speter 7318334Speter mpfr_mul_z (w, v, mpfr_bernoulli_cache(i), MPFR_RNDN); 7418334Speter /* here, w_2i = v_2i * B_2i * (2i+1)! with 7590075Sobrien error(w_2i) < 2^(5 * i + 8) ulp(w_2i) (see algorithms.tex) */ 7690075Sobrien 7790075Sobrien mpfr_add (s, s, w, MPFR_RNDN); 7890075Sobrien 7990075Sobrien err = MAX (err + se, 5 * i + 8 + MPFR_GET_EXP (w)) 8090075Sobrien - MPFR_GET_EXP (s); 8190075Sobrien err = 2 + MAX (-1, err); 8290075Sobrien se = MPFR_GET_EXP (s); 8390075Sobrien if (MPFR_GET_EXP (w) <= se - (mpfr_exp_t) p) 8490075Sobrien break; 8590075Sobrien } 8690075Sobrien 8790075Sobrien /* the previous value of err is the rounding error, 8890075Sobrien the truncation error is less than EXP(z) - 6 * i - 5 8990075Sobrien (see algorithms.tex) */ 9090075Sobrien err = MAX (err, MPFR_GET_EXP (z) - 6 * i - 5) + 1; 9190075Sobrien if (MPFR_CAN_ROUND (s, (mpfr_exp_t) p - err, sump, rnd_mode)) 9290075Sobrien break; 9390075Sobrien 9418334Speter MPFR_ZIV_NEXT (loop, p); 9518334Speter mpfr_set_prec (s, p); 9618334Speter mpfr_set_prec (u, p); 9718334Speter mpfr_set_prec (v, p); 9818334Speter mpfr_set_prec (w, p); 9918334Speter } 10018334Speter MPFR_ZIV_FREE (loop); 10118334Speter mpfr_set (sum, s, rnd_mode); 10218334Speter 10318334Speter mpfr_clears (s, u, v, w, (mpfr_ptr) 0); 10418334Speter 10518334Speter /* Let K be the returned value. 10618334Speter 1. As we compute an alternating series, the truncation error has the same 10718334Speter sign as the next term w_{K+2} which is positive iff K%4 == 0. 10850397Sobrien 2. Assume that error(z) <= (1+t) z', where z' is the actual value, then 10918334Speter error(s) <= 2 * (K+1) * t (see algorithms.tex). 11018334Speter */ 11118334Speter return 2 * i; 11218334Speter} 11318334Speter 11418334Speter/* try asymptotic expansion when x is large and positive: 11518334Speter Li2(x) = -log(x)^2/2 + Pi^2/3 - 1/x + O(1/x^2). 11690075Sobrien More precisely for x >= 2 we have for g(x) = -log(x)^2/2 + Pi^2/3: 11718334Speter -2 <= x * (Li2(x) - g(x)) <= -1 11818334Speter thus |Li2(x) - g(x)| <= 2/x. 11950397Sobrien Assumes x >= 38, which ensures log(x)^2/2 >= 2*Pi^2/3, and g(x) <= -3.3. 12018334Speter Return 0 if asymptotic expansion failed (unable to round), otherwise 12150397Sobrien returns 1 for RNDF, and correct ternary value otherwise. 12250397Sobrien*/ 12350397Sobrienstatic int 12418334Spetermpfr_li2_asympt_pos (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 12550397Sobrien{ 12650397Sobrien mpfr_t g, h; 12750397Sobrien mpfr_prec_t w = MPFR_PREC (y) + 20; 12850397Sobrien int inex = 0; 12950397Sobrien 13050397Sobrien MPFR_ASSERTN (mpfr_cmp_ui (x, 38) >= 0); 13118334Speter 13250397Sobrien mpfr_init2 (g, w); 13350397Sobrien mpfr_init2 (h, w); 13450397Sobrien mpfr_log (g, x, MPFR_RNDN); /* rel. error <= |(1 + theta) - 1| */ 13550397Sobrien mpfr_sqr (g, g, MPFR_RNDN); /* rel. error <= |(1 + theta)^3 - 1| <= 2^(2-w) */ 13650397Sobrien mpfr_div_2ui (g, g, 1, MPFR_RNDN); /* rel. error <= 2^(2-w) */ 13750397Sobrien mpfr_const_pi (h, MPFR_RNDN); /* error <= 2^(1-w) */ 13818334Speter mpfr_sqr (h, h, MPFR_RNDN); /* rel. error <= 2^(2-w) */ 13918334Speter mpfr_div_ui (h, h, 3, MPFR_RNDN); /* rel. error <= |(1 + theta)^4 - 1| 14018334Speter <= 5 * 2^(-w) */ 14118334Speter /* since x is chosen such that log(x)^2/2 >= 2 * (Pi^2/3), we should have 14218334Speter g >= 2*h, thus |g-h| >= |h|, and the relative error on g is at most 14318334Speter multiplied by 2 in the difference, and that by h is unchanged. */ 14450397Sobrien MPFR_ASSERTN (MPFR_EXP (g) > MPFR_EXP (h)); 14518334Speter mpfr_sub (g, h, g, MPFR_RNDN); /* err <= ulp(g)/2 + g*2^(3-w) + g*5*2^(-w) 14618334Speter <= ulp(g) * (1/2 + 8 + 5) < 14 ulp(g). 14718334Speter 14818334Speter If in addition 2/x <= 2 ulp(g), i.e., 14918334Speter 1/x <= ulp(g), then the total error is 15018334Speter bounded by 16 ulp(g). */ 15118334Speter if ((MPFR_EXP (x) >= (mpfr_exp_t) w - MPFR_EXP (g)) && 15218334Speter MPFR_CAN_ROUND (g, w - 4, MPFR_PREC (y), rnd_mode)) 15318334Speter { 15418334Speter inex = mpfr_set (y, g, rnd_mode); 15518334Speter if (rnd_mode == MPFR_RNDF) 15618334Speter inex = 1; 15718334Speter } 15818334Speter 15918334Speter mpfr_clear (g); 16018334Speter mpfr_clear (h); 16118334Speter 16218334Speter return inex; 16350397Sobrien} 16418334Speter 16518334Speter/* try asymptotic expansion when x is large and negative: 16650397Sobrien Li2(x) = -log(-x)^2/2 - Pi^2/6 - 1/x + O(1/x^2). 16718334Speter More precisely for x <= -2 we have for g(x) = -log(-x)^2/2 - Pi^2/6: 16818334Speter |Li2(x) - g(x)| <= 1/|x|. 16950397Sobrien Assumes x <= -7, which ensures |log(-x)^2/2| >= Pi^2/6, and g(x) <= -3.5. 17018334Speter Return 0 if asymptotic expansion failed (unable to round), otherwise 17118334Speter returns 1 for RNDF, and correct ternary value otherwise. 17218334Speter*/ 17318334Speterstatic int 17418334Spetermpfr_li2_asympt_neg (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 17518334Speter{ 17650397Sobrien mpfr_t g, h; 17750397Sobrien mpfr_prec_t w = MPFR_PREC (y) + 20; 17850397Sobrien int inex = 0; 17918334Speter 18018334Speter MPFR_ASSERTN (mpfr_cmp_si (x, -7) <= 0); 18118334Speter 18218334Speter mpfr_init2 (g, w); 18350397Sobrien mpfr_init2 (h, w); 18450397Sobrien mpfr_neg (g, x, MPFR_RNDN); 18550397Sobrien mpfr_log (g, g, MPFR_RNDN); /* rel. error <= |(1 + theta) - 1| */ 18618334Speter mpfr_sqr (g, g, MPFR_RNDN); /* rel. error <= |(1 + theta)^3 - 1| <= 2^(2-w) */ 18718334Speter mpfr_div_2ui (g, g, 1, MPFR_RNDN); /* rel. error <= 2^(2-w) */ 18850397Sobrien mpfr_const_pi (h, MPFR_RNDN); /* error <= 2^(1-w) */ 18950397Sobrien mpfr_sqr (h, h, MPFR_RNDN); /* rel. error <= 2^(2-w) */ 19090075Sobrien mpfr_div_ui (h, h, 6, MPFR_RNDN); /* rel. error <= |(1 + theta)^4 - 1| 19190075Sobrien <= 5 * 2^(-w) */ 19290075Sobrien MPFR_ASSERTN (MPFR_EXP (g) >= MPFR_EXP (h)); 19390075Sobrien mpfr_add (g, g, h, MPFR_RNDN); /* err <= ulp(g)/2 + g*2^(2-w) + g*5*2^(-w) 19490075Sobrien <= ulp(g) * (1/2 + 4 + 5) < 10 ulp(g). 19590075Sobrien 19690075Sobrien If in addition |1/x| <= 4 ulp(g), then the 19718334Speter total error is bounded by 16 ulp(g). */ 19818334Speter if ((MPFR_EXP (x) >= (mpfr_exp_t) (w - 2) - MPFR_EXP (g)) && 19918334Speter MPFR_CAN_ROUND (g, w - 4, MPFR_PREC (y), rnd_mode)) 20018334Speter { 20118334Speter inex = mpfr_neg (y, g, rnd_mode); 20218334Speter if (rnd_mode == MPFR_RNDF) 20318334Speter inex = 1; 20418334Speter } 20518334Speter 20618334Speter mpfr_clear (g); 20718334Speter mpfr_clear (h); 20818334Speter 20918334Speter return inex; 21018334Speter} 21118334Speter 21218334Speter/* Compute the real part of the dilogarithm defined by 21350397Sobrien Li2(x) = -\Int_{t=0}^x log(1-t)/t dt */ 21450397Sobrienint 21550397Sobrienmpfr_li2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 21650397Sobrien{ 21750397Sobrien int inexact; 21850397Sobrien mpfr_exp_t err; 21918334Speter mpfr_prec_t yp, m; 22090075Sobrien MPFR_ZIV_DECL (loop); 22190075Sobrien MPFR_SAVE_EXPO_DECL (expo); 22290075Sobrien 22390075Sobrien MPFR_LOG_FUNC 22490075Sobrien (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 22590075Sobrien ("y[%Pu]=%.*Rg inexact=%d", 22690075Sobrien mpfr_get_prec (y), mpfr_log_prec, y, inexact)); 22790075Sobrien 22890075Sobrien if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 22990075Sobrien { 23090075Sobrien if (MPFR_IS_NAN (x)) 23190075Sobrien { 23290075Sobrien MPFR_SET_NAN (y); 23390075Sobrien MPFR_RET_NAN; 23490075Sobrien } 23590075Sobrien else if (MPFR_IS_INF (x)) 23690075Sobrien { 23790075Sobrien MPFR_SET_NEG (y); 23890075Sobrien MPFR_SET_INF (y); 23990075Sobrien MPFR_RET (0); 24090075Sobrien } 24190075Sobrien else /* x is zero */ 24290075Sobrien { 24390075Sobrien MPFR_ASSERTD (MPFR_IS_ZERO (x)); 24490075Sobrien MPFR_SET_SAME_SIGN (y, x); 24590075Sobrien MPFR_SET_ZERO (y); 24690075Sobrien MPFR_RET (0); 24790075Sobrien } 24890075Sobrien } 24990075Sobrien 25090075Sobrien /* Li2(x) = x + x^2/4 + x^3/9 + ..., more precisely for 0 < x <= 1/2 25190075Sobrien we have |Li2(x) - x| < x^2/2 <= 2^(2EXP(x)-1) and for -1/2 <= x < 0 25290075Sobrien we have |Li2(x) - x| < x^2/4 <= 2^(2EXP(x)-2) */ 25390075Sobrien if (MPFR_IS_POS (x)) 25490075Sobrien MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -MPFR_GET_EXP (x), 1, 1, rnd_mode, 25590075Sobrien {}); 25690075Sobrien else 25790075Sobrien MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -MPFR_GET_EXP (x), 2, 0, rnd_mode, 25890075Sobrien {}); 25990075Sobrien 26090075Sobrien MPFR_SAVE_EXPO_MARK (expo); 26190075Sobrien yp = MPFR_PREC (y); 26290075Sobrien m = yp + MPFR_INT_CEIL_LOG2 (yp) + 13; 26390075Sobrien 26490075Sobrien if (MPFR_LIKELY ((mpfr_cmp_ui (x, 0) > 0) && (mpfr_cmp_ui_2exp (x, 1, -1) <= 0))) 26590075Sobrien /* 0 < x <= 1/2: Li2(x) = S(-log(1-x))-log^2(1-x)/4 */ 26690075Sobrien { 26790075Sobrien mpfr_t s, u; 26890075Sobrien mpfr_exp_t expo_l; 26990075Sobrien int k; 27090075Sobrien 27190075Sobrien mpfr_init2 (u, m); 27218334Speter mpfr_init2 (s, m); 27318334Speter 27418334Speter MPFR_ZIV_INIT (loop, m); 27518334Speter for (;;) 27618334Speter { 27790075Sobrien mpfr_ui_sub (u, 1, x, MPFR_RNDN); 27818334Speter mpfr_log (u, u, MPFR_RNDU); 27918334Speter if (MPFR_IS_ZERO(u)) 28018334Speter goto next_m; 28118334Speter mpfr_neg (u, u, MPFR_RNDN); /* u = -log(1-x) */ 28218334Speter expo_l = MPFR_GET_EXP (u); 28352284Sobrien k = li2_series (s, u, MPFR_RNDU); 28418334Speter err = 1 + MPFR_INT_CEIL_LOG2 (k + 1); 28552284Sobrien 28652284Sobrien mpfr_sqr (u, u, MPFR_RNDU); 28752284Sobrien mpfr_div_2ui (u, u, 2, MPFR_RNDU); /* u = log^2(1-x) / 4 */ 28818334Speter mpfr_sub (s, s, u, MPFR_RNDN); 28990075Sobrien 29018334Speter /* error(s) <= (0.5 + 2^(d-EXP(s)) 29190075Sobrien + 2^(3 + MAX(1, - expo_l) - EXP(s))) ulp(s) */ 29290075Sobrien err = MAX (err, MAX (1, - expo_l) - 1) - MPFR_GET_EXP (s); 29318334Speter err = 2 + MAX (-1, err); 29418334Speter if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 29590075Sobrien break; 29618334Speter 29790075Sobrien next_m: 29890075Sobrien MPFR_ZIV_NEXT (loop, m); 29990075Sobrien mpfr_set_prec (u, m); 30090075Sobrien mpfr_set_prec (s, m); 30190075Sobrien } 30296263Sobrien MPFR_ZIV_FREE (loop); 30318334Speter inexact = mpfr_set (y, s, rnd_mode); 30490075Sobrien 30552284Sobrien mpfr_clear (u); 30618334Speter mpfr_clear (s); 30718334Speter MPFR_SAVE_EXPO_FREE (expo); 30818334Speter return mpfr_check_range (y, inexact, rnd_mode); 30918334Speter } 31018334Speter else if (!mpfr_cmp_ui (x, 1)) 31118334Speter /* Li2(1)= pi^2 / 6 */ 31218334Speter { 31318334Speter mpfr_t u; 31418334Speter mpfr_init2 (u, m); 31518334Speter 31618334Speter MPFR_ZIV_INIT (loop, m); 31718334Speter for (;;) 31818334Speter { 31990075Sobrien mpfr_const_pi (u, MPFR_RNDU); 32018334Speter mpfr_sqr (u, u, MPFR_RNDN); 32118334Speter mpfr_div_ui (u, u, 6, MPFR_RNDN); 32218334Speter 32318334Speter err = m - 4; /* error(u) <= 19/2 ulp(u) */ 32418334Speter if (MPFR_CAN_ROUND (u, err, yp, rnd_mode)) 32518334Speter break; 32618334Speter 32718334Speter MPFR_ZIV_NEXT (loop, m); 32896263Sobrien mpfr_set_prec (u, m); 32996263Sobrien } 33096263Sobrien MPFR_ZIV_FREE (loop); 33196263Sobrien inexact = mpfr_set (y, u, rnd_mode); 33296263Sobrien 33396263Sobrien mpfr_clear (u); 33418334Speter MPFR_SAVE_EXPO_FREE (expo); 33590075Sobrien return mpfr_check_range (y, inexact, rnd_mode); 33690075Sobrien } 33790075Sobrien else if (mpfr_cmp_ui (x, 2) >= 0) 33890075Sobrien /* x >= 2: Li2(x) = -S(-log(1-1/x))-log^2(x)/2+log^2(1-1/x)/4+pi^2/3 */ 33990075Sobrien { 34090075Sobrien int k; 34196263Sobrien mpfr_exp_t expo_l; 34296263Sobrien mpfr_t s, u, xx; 34396263Sobrien 34496263Sobrien if (mpfr_cmp_ui (x, 38) >= 0) 34590075Sobrien { 34690075Sobrien inexact = mpfr_li2_asympt_pos (y, x, rnd_mode); 34718334Speter if (inexact != 0) 34896263Sobrien goto end_of_case_gt2; 34990075Sobrien } 35096263Sobrien 35190075Sobrien mpfr_init2 (u, m); 35290075Sobrien mpfr_init2 (s, m); 35390075Sobrien mpfr_init2 (xx, m); 35418334Speter 35518334Speter MPFR_ZIV_INIT (loop, m); 35618334Speter for (;;) 35752284Sobrien { 35852284Sobrien mpfr_ui_div (xx, 1, x, MPFR_RNDN); 35952284Sobrien mpfr_neg (xx, xx, MPFR_RNDN); 36052284Sobrien mpfr_log1p (u, xx, MPFR_RNDD); 36152284Sobrien mpfr_neg (u, u, MPFR_RNDU); /* u = -log(1-1/x) */ 36252284Sobrien expo_l = MPFR_GET_EXP (u); 36352284Sobrien k = li2_series (s, u, MPFR_RNDN); 36452284Sobrien mpfr_neg (s, s, MPFR_RNDN); 36552284Sobrien err = MPFR_INT_CEIL_LOG2 (k + 1) + 1; /* error(s) <= 2^err ulp(s) */ 36652284Sobrien 36752284Sobrien mpfr_sqr (u, u, MPFR_RNDN); 36852284Sobrien mpfr_div_2ui (u, u, 2, MPFR_RNDN); /* u= log^2(1-1/x)/4 */ 36918334Speter mpfr_add (s, s, u, MPFR_RNDN); 37096263Sobrien err = 37118334Speter MAX (err, 37290075Sobrien 3 + MAX (1, -expo_l) + MPFR_GET_EXP (u)) - MPFR_GET_EXP (s); 37318334Speter err = 2 + MAX (-1, err); /* error(s) <= 2^err ulp(s) */ 37418334Speter err += MPFR_GET_EXP (s); 37518334Speter 37618334Speter mpfr_log (u, x, MPFR_RNDU); 37718334Speter mpfr_sqr (u, u, MPFR_RNDN); 37890075Sobrien mpfr_div_2ui (u, u, 1, MPFR_RNDN); /* u = log^2(x)/2 */ 37990075Sobrien mpfr_sub (s, s, u, MPFR_RNDN); 38090075Sobrien err = MAX (err, 3 + MPFR_GET_EXP (u)) - MPFR_GET_EXP (s); 38190075Sobrien err = 2 + MAX (-1, err); /* error(s) <= 2^err ulp(s) */ 38290075Sobrien err += MPFR_GET_EXP (s); 38390075Sobrien 38490075Sobrien mpfr_const_pi (u, MPFR_RNDU); 38590075Sobrien mpfr_sqr (u, u, MPFR_RNDN); 38690075Sobrien mpfr_div_ui (u, u, 3, MPFR_RNDN); /* u = pi^2/3 */ 38790075Sobrien mpfr_add (s, s, u, MPFR_RNDN); 38890075Sobrien err = MAX (err, 2) - MPFR_GET_EXP (s); 38990075Sobrien err = 2 + MAX (-1, err); /* error(s) <= 2^err ulp(s) */ 39090075Sobrien if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 39190075Sobrien break; 39290075Sobrien 39390075Sobrien MPFR_ZIV_NEXT (loop, m); 39490075Sobrien mpfr_set_prec (u, m); 39596263Sobrien mpfr_set_prec (s, m); 39696263Sobrien mpfr_set_prec (xx, m); 39796263Sobrien } 39896263Sobrien MPFR_ZIV_FREE (loop); 39996263Sobrien inexact = mpfr_set (y, s, rnd_mode); 40096263Sobrien mpfr_clears (s, u, xx, (mpfr_ptr) 0); 40196263Sobrien 40296263Sobrien end_of_case_gt2: 40396263Sobrien MPFR_SAVE_EXPO_FREE (expo); 40490075Sobrien return mpfr_check_range (y, inexact, rnd_mode); 40590075Sobrien } 40690075Sobrien else if (mpfr_cmp_ui (x, 1) > 0) 40790075Sobrien /* 2 > x > 1: Li2(x) = S(log(x))+log^2(x)/4-log(x)log(x-1)+pi^2/6 */ 40890075Sobrien { 40990075Sobrien int k; 41090075Sobrien mpfr_exp_t e1, e2; 41190075Sobrien mpfr_t s, u, v, xx; 41218334Speter mpfr_init2 (s, m); 41318334Speter mpfr_init2 (u, m); 41418334Speter mpfr_init2 (v, m); 41518334Speter mpfr_init2 (xx, m); 41618334Speter 41718334Speter MPFR_ZIV_INIT (loop, m); 41890075Sobrien for (;;) 41990075Sobrien { 42018334Speter mpfr_log (v, x, MPFR_RNDU); 42190075Sobrien k = li2_series (s, v, MPFR_RNDN); 42290075Sobrien e1 = MPFR_GET_EXP (s); 42318334Speter 42418334Speter mpfr_sqr (u, v, MPFR_RNDN); 42518334Speter mpfr_div_2ui (u, u, 2, MPFR_RNDN); /* u = log^2(x)/4 */ 42618334Speter mpfr_add (s, s, u, MPFR_RNDN); 42718334Speter 42818334Speter mpfr_sub_ui (xx, x, 1, MPFR_RNDN); 42918334Speter mpfr_log (u, xx, MPFR_RNDU); 43018334Speter e2 = MPFR_GET_EXP (u); 43190075Sobrien mpfr_mul (u, v, u, MPFR_RNDN); /* u = log(x) * log(x-1) */ 43290075Sobrien mpfr_sub (s, s, u, MPFR_RNDN); 43390075Sobrien 43490075Sobrien mpfr_const_pi (u, MPFR_RNDU); 43518334Speter mpfr_sqr (u, u, MPFR_RNDN); 43690075Sobrien mpfr_div_ui (u, u, 6, MPFR_RNDN); /* u = pi^2/6 */ 43790075Sobrien mpfr_add (s, s, u, MPFR_RNDN); 43890075Sobrien /* error(s) <= (31 + (k+1) * 2^(1-e1) + 2^(1-e2)) ulp(s) 43990075Sobrien see algorithms.tex */ 44090075Sobrien err = MAX (MPFR_INT_CEIL_LOG2 (k + 1) + 1 - e1, 1 - e2); 44190075Sobrien err = 2 + MAX (5, err); 44290075Sobrien if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 44390075Sobrien break; 44490075Sobrien 44590075Sobrien MPFR_ZIV_NEXT (loop, m); 44652284Sobrien mpfr_set_prec (s, m); 44790075Sobrien mpfr_set_prec (u, m); 44890075Sobrien mpfr_set_prec (v, m); 44990075Sobrien mpfr_set_prec (xx, m); 45090075Sobrien } 45190075Sobrien MPFR_ZIV_FREE (loop); 45252284Sobrien inexact = mpfr_set (y, s, rnd_mode); 45318334Speter 45418334Speter mpfr_clears (s, u, v, xx, (mpfr_ptr) 0); 45518334Speter MPFR_SAVE_EXPO_FREE (expo); 45618334Speter return mpfr_check_range (y, inexact, rnd_mode); 45718334Speter } 45818334Speter else if (mpfr_cmp_ui_2exp (x, 1, -1) > 0) /* 1/2 < x < 1 */ 45918334Speter /* 1 > x > 1/2: Li2(x) = -S(-log(x))+log^2(x)/4-log(x)log(1-x)+pi^2/6 */ 46018334Speter { 46118334Speter int k; 46218334Speter mpfr_t s, u, v, xx; 46318334Speter mpfr_init2 (s, m); 46418334Speter mpfr_init2 (u, m); 46518334Speter mpfr_init2 (v, m); 46690075Sobrien mpfr_init2 (xx, m); 46790075Sobrien 46890075Sobrien 46918334Speter MPFR_ZIV_INIT (loop, m); 47018334Speter for (;;) 47118334Speter { 47218334Speter mpfr_log (u, x, MPFR_RNDD); 47318334Speter mpfr_neg (u, u, MPFR_RNDN); 47418334Speter k = li2_series (s, u, MPFR_RNDN); 47590075Sobrien mpfr_neg (s, s, MPFR_RNDN); 47618334Speter err = 1 + MPFR_INT_CEIL_LOG2 (k + 1) - MPFR_GET_EXP (s); 47718334Speter 47818334Speter mpfr_ui_sub (xx, 1, x, MPFR_RNDN); 47918334Speter mpfr_log (v, xx, MPFR_RNDU); 48018334Speter mpfr_mul (v, v, u, MPFR_RNDN); /* v = - log(x) * log(1-x) */ 48190075Sobrien mpfr_add (s, s, v, MPFR_RNDN); 48290075Sobrien err = MAX (err, 1 - MPFR_GET_EXP (v)); 48390075Sobrien err = 2 + MAX (3, err) - MPFR_GET_EXP (s); 48490075Sobrien 48590075Sobrien mpfr_sqr (u, u, MPFR_RNDN); 48690075Sobrien mpfr_div_2ui (u, u, 2, MPFR_RNDN); /* u = log^2(x)/4 */ 48790075Sobrien mpfr_add (s, s, u, MPFR_RNDN); 48818334Speter err = MAX (err, 2 + MPFR_GET_EXP (u)) - MPFR_GET_EXP (s); 48918334Speter err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); 49018334Speter 49118334Speter mpfr_const_pi (u, MPFR_RNDU); 49218334Speter mpfr_sqr (u, u, MPFR_RNDN); 49318334Speter mpfr_div_ui (u, u, 6, MPFR_RNDN); /* u = pi^2/6 */ 49418334Speter mpfr_add (s, s, u, MPFR_RNDN); 49590075Sobrien err = MAX (err, 3) - MPFR_GET_EXP (s); 49618334Speter err = 2 + MAX (-1, err); 49718334Speter 49818334Speter if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 49918334Speter break; 50018334Speter 50118334Speter MPFR_ZIV_NEXT (loop, m); 50218334Speter mpfr_set_prec (s, m); 50318334Speter mpfr_set_prec (u, m); 50418334Speter mpfr_set_prec (v, m); 50518334Speter mpfr_set_prec (xx, m); 50652284Sobrien } 50718334Speter MPFR_ZIV_FREE (loop); 50818334Speter inexact = mpfr_set (y, s, rnd_mode); 50918334Speter 51052284Sobrien mpfr_clears (s, u, v, xx, (mpfr_ptr) 0); 51152284Sobrien MPFR_SAVE_EXPO_FREE (expo); 51252284Sobrien return mpfr_check_range (y, inexact, rnd_mode); 51352284Sobrien } 51452284Sobrien else if (mpfr_cmp_si (x, -1) >= 0) 51552284Sobrien /* 0 > x >= -1: Li2(x) = -S(log(1-x))-log^2(1-x)/4 */ 51652284Sobrien { 51752284Sobrien int k; 51852284Sobrien mpfr_exp_t expo_l; 51952284Sobrien mpfr_t s, u, xx; 52052284Sobrien mpfr_init2 (s, m); 52152284Sobrien mpfr_init2 (u, m); 52252284Sobrien mpfr_init2 (xx, m); 52352284Sobrien 52452284Sobrien MPFR_ZIV_INIT (loop, m); 52590075Sobrien for (;;) 52652284Sobrien { 52718334Speter mpfr_neg (xx, x, MPFR_RNDN); 52818334Speter mpfr_log1p (u, xx, MPFR_RNDN); 52918334Speter k = li2_series (s, u, MPFR_RNDN); 53090075Sobrien mpfr_neg (s, s, MPFR_RNDN); 53118334Speter expo_l = MPFR_GET_EXP (u); 53218334Speter err = 1 + MPFR_INT_CEIL_LOG2 (k + 1) - MPFR_GET_EXP (s); 53318334Speter 53418334Speter mpfr_sqr (u, u, MPFR_RNDN); 53518334Speter mpfr_div_2ui (u, u, 2, MPFR_RNDN); /* u = log^2(1-x)/4 */ 53696263Sobrien mpfr_sub (s, s, u, MPFR_RNDN); 53796263Sobrien err = MAX (err, - expo_l); 53896263Sobrien err = 2 + MAX (err, 3); 53918334Speter if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 54018334Speter break; 54118334Speter 54218334Speter MPFR_ZIV_NEXT (loop, m); 54318334Speter mpfr_set_prec (s, m); 54450397Sobrien mpfr_set_prec (u, m); 54518334Speter mpfr_set_prec (xx, m); 54618334Speter } 54790075Sobrien MPFR_ZIV_FREE (loop); 54818334Speter inexact = mpfr_set (y, s, rnd_mode); 54990075Sobrien 55018334Speter mpfr_clears (s, u, xx, (mpfr_ptr) 0); 55118334Speter MPFR_SAVE_EXPO_FREE (expo); 55218334Speter return mpfr_check_range (y, inexact, rnd_mode); 55318334Speter } 55418334Speter else 55518334Speter /* x < -1: Li2(x) 55690075Sobrien = S(log(1-1/x))-log^2(-x)/4-log(1-x)log(-x)/2+log^2(1-x)/4-pi^2/6 */ 55752284Sobrien { 55818334Speter int k; 55918334Speter mpfr_t s, u, v, w, xx; 56018334Speter 56150397Sobrien if (mpfr_cmp_si (x, -7) <= 0) 56250397Sobrien { 56350397Sobrien inexact = mpfr_li2_asympt_neg (y, x, rnd_mode); 56450397Sobrien if (inexact != 0) 56518334Speter goto end_of_case_ltm1; 56690075Sobrien } 56750397Sobrien 56818334Speter mpfr_init2 (s, m); 56918334Speter mpfr_init2 (u, m); 57018334Speter mpfr_init2 (v, m); 57118334Speter mpfr_init2 (w, m); 57218334Speter mpfr_init2 (xx, m); 57318334Speter 57418334Speter MPFR_ZIV_INIT (loop, m); 57518334Speter for (;;) 57618334Speter { 57718334Speter mpfr_ui_div (xx, 1, x, MPFR_RNDN); 57818334Speter mpfr_neg (xx, xx, MPFR_RNDN); 57918334Speter mpfr_log1p (u, xx, MPFR_RNDN); 58090075Sobrien k = li2_series (s, u, MPFR_RNDN); 58118334Speter 58218334Speter mpfr_ui_sub (xx, 1, x, MPFR_RNDN); 58318334Speter mpfr_log (u, xx, MPFR_RNDU); 58418334Speter mpfr_neg (xx, x, MPFR_RNDN); 58518334Speter mpfr_log (v, xx, MPFR_RNDU); 58690075Sobrien mpfr_mul (w, v, u, MPFR_RNDN); 58790075Sobrien mpfr_div_2ui (w, w, 1, MPFR_RNDN); /* w = log(-x) * log(1-x) / 2 */ 58818334Speter mpfr_sub (s, s, w, MPFR_RNDN); 58918334Speter err = 1 + MAX (3, MPFR_INT_CEIL_LOG2 (k+1) + 1 - MPFR_GET_EXP (s)) 59096263Sobrien + MPFR_GET_EXP (s); 59190075Sobrien 59218334Speter mpfr_sqr (w, v, MPFR_RNDN); 59318334Speter mpfr_div_2ui (w, w, 2, MPFR_RNDN); /* w = log^2(-x) / 4 */ 59418334Speter mpfr_sub (s, s, w, MPFR_RNDN); 59590075Sobrien err = MAX (err, 3 + MPFR_GET_EXP(w)) - MPFR_GET_EXP (s); 59618334Speter err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); 59790075Sobrien 59890075Sobrien mpfr_sqr (w, u, MPFR_RNDN); 59918334Speter mpfr_div_2ui (w, w, 2, MPFR_RNDN); /* w = log^2(1-x) / 4 */ 60018334Speter mpfr_add (s, s, w, MPFR_RNDN); 60190075Sobrien err = MAX (err, 3 + MPFR_GET_EXP (w)) - MPFR_GET_EXP (s); 60218334Speter err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); 60318334Speter 60418334Speter mpfr_const_pi (w, MPFR_RNDU); 60518334Speter mpfr_sqr (w, w, MPFR_RNDN); 60618334Speter mpfr_div_ui (w, w, 6, MPFR_RNDN); /* w = pi^2 / 6 */ 60718334Speter mpfr_sub (s, s, w, MPFR_RNDN); 60818334Speter err = MAX (err, 3) - MPFR_GET_EXP (s); 60990075Sobrien err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); 61018334Speter 61118334Speter if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 61218334Speter break; 61318334Speter 61418334Speter MPFR_ZIV_NEXT (loop, m); 61518334Speter mpfr_set_prec (s, m); 61690075Sobrien mpfr_set_prec (u, m); 61718334Speter mpfr_set_prec (v, m); 61850397Sobrien mpfr_set_prec (w, m); 61918334Speter mpfr_set_prec (xx, m); 62018334Speter } 62118334Speter MPFR_ZIV_FREE (loop); 62218334Speter inexact = mpfr_set (y, s, rnd_mode); 62318334Speter mpfr_clears (s, u, v, w, xx, (mpfr_ptr) 0); 62418334Speter 62518334Speter end_of_case_ltm1: 62618334Speter MPFR_SAVE_EXPO_FREE (expo); 62718334Speter return mpfr_check_range (y, inexact, rnd_mode); 62818334Speter } 62918334Speter 63018334Speter MPFR_RET_NEVER_GO_HERE (); 63118334Speter} 63218334Speter