diff options
| -rw-r--r-- | TODO | 2 | ||||
| -rw-r--r-- | erfc.c | 2 | ||||
| -rw-r--r-- | li2.c | 50 | ||||
| -rw-r--r-- | lngamma.c | 2 | ||||
| -rw-r--r-- | mpfr.texi | 4 | ||||
| -rw-r--r-- | pow.c | 6 | ||||
| -rw-r--r-- | rem1.c | 54 | ||||
| -rw-r--r-- | tests/tfmod.c | 42 | ||||
| -rw-r--r-- | tests/tgamma.c | 2 | ||||
| -rw-r--r-- | tests/tpow.c | 4 | ||||
| -rw-r--r-- | tests/tzeta.c | 2 |
11 files changed, 85 insertions, 85 deletions
@@ -112,7 +112,7 @@ Table of contents: - dilog() [the dilogarithm: dilog(x) = int(ln(t)/(1-t), t=1..x)] See Numerical evaluation of multiple polylogarithms Jens Vollinga and Stefan Weinzierl, Computer Physics Communications, - Volume 167, Issue 3, 1 May 2005, Pages 177-194 + Volume 167, Issue 3, 1 May 2005, Pages 177-194 Cf also http://calgo.acm.org/ (490.gz) - mpfr_printf [Sisyphus <kalinabears@iinet.net.au> Tue, 04 Jan 2005] for example mpfr_printf ("%.2Ff\n", x) @@ -160,7 +160,7 @@ mpfr_erfc (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd) /* For x < 0 going to -infinity, erfc(x) tends to 2 by below. More precisely, we have 2 + 1/sqrt(Pi)/x/exp(x^2) < erfc(x) < 2. Thus log2 |2 - erfc(x)| <= -log2|x| - x^2 / log(2). - If |2 - erfc(x)| < 2^(-PREC(y)) then the result is either 2 or + If |2 - erfc(x)| < 2^(-PREC(y)) then the result is either 2 or nextbelow(2). For x <= -27282, -log2|x| - x^2 / log(2) <= -2^30. */ @@ -81,10 +81,10 @@ bernoulli (mpz_t * b, unsigned long n) } /* Compute the alternating series - s = S(z) = \sum_{k=0}^infty B_{2k} (z))^{2k+1} / (2k+1)! + s = S(z) = \sum_{k=0}^infty B_{2k} (z))^{2k+1} / (2k+1)! with 0 < z <= log(2) to the precision of s rounded in the direction - rnd_mode. - Return the maximum index of the truncature which is useful + rnd_mode. + Return the maximum index of the truncature which is useful for determinating the relative error. */ static int @@ -98,7 +98,7 @@ li2_series (mpfr_t sum, mpfr_srcptr z, mpfr_rnd_t rnd_mode) mpz_t *B; MPFR_ZIV_DECL (loop); - /* The series converges for |z| < 2 pi, but in mpfr_li2 the argument is + /* The series converges for |z| < 2 pi, but in mpfr_li2 the argument is reduced so that 0 < z <= log(2). Here is additionnal check that z is (nearly) correct */ MPFR_ASSERTD (MPFR_IS_STRICTPOS (z)); @@ -142,15 +142,15 @@ li2_series (mpfr_t sum, mpfr_srcptr z, mpfr_rnd_t rnd_mode) mpfr_add (s, s, w, GMP_RNDN); err = MAX (err + se, 5 * i + 8 + MPFR_GET_EXP (w)) - - MPFR_GET_EXP (s); - err = 2 + MAX (-1, err); + - MPFR_GET_EXP (s); + err = 2 + MAX (-1, err); se = MPFR_GET_EXP (s); if (MPFR_GET_EXP (w) <= se - (mp_exp_t) p) break; } - /* the previous value of err is the rounding error, - the truncation error is less than EXP(z) - 4 * i - 4 + /* the previous value of err is the rounding error, + the truncation error is less than EXP(z) - 4 * i - 4 (see algorithms.tex) */ err = MAX (err, MPFR_GET_EXP (z) - 4 * i - 4) + 1; if (MPFR_CAN_ROUND (s, p - err, sump, rnd_mode)) @@ -313,7 +313,7 @@ mpfr_li2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) {}); else MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -MPFR_GET_EXP (x), 2, 0, rnd_mode, - {}); + {}); MPFR_SAVE_EXPO_MARK (expo); yp = MPFR_PREC (y); @@ -345,8 +345,8 @@ mpfr_li2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) /* error(s) <= (0.5 + 2^(d-EXP(s)) + 2^(3 + MAX(1, - expo_l) - EXP(s))) ulp(s) */ - err = MAX (err, MAX (1, - expo_l) - 1) - MPFR_GET_EXP (s); - err = 2 + MAX (-1, err); + err = MAX (err, MAX (1, - expo_l) - 1) - MPFR_GET_EXP (s); + err = 2 + MAX (-1, err); if (MPFR_CAN_ROUND (s, m - err, yp, rnd_mode)) break; @@ -490,8 +490,8 @@ mpfr_li2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) mpfr_sqr (u, u, GMP_RNDN); mpfr_div_ui (u, u, 6, GMP_RNDN); /* u = pi^2/6 */ mpfr_add (s, s, u, GMP_RNDN); - /* error(s) <= (31 + (k+1) * 2^(1-e1) + 2^(1-e2)) ulp(s) - see algorithms.tex */ + /* error(s) <= (31 + (k+1) * 2^(1-e1) + 2^(1-e2)) ulp(s) + see algorithms.tex */ err = MAX (MPFR_INT_CEIL_LOG2 (k + 1) + 1 - e1, 1 - e2); err = 2 + MAX (5, err); if (MPFR_CAN_ROUND (s, m - err, yp, rnd_mode)) @@ -535,20 +535,20 @@ mpfr_li2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) mpfr_mul (v, v, u, GMP_RNDN); /* v = - log(x) * log(1-x) */ mpfr_add (s, s, v, GMP_RNDN); err = MAX (err, 1 - MPFR_GET_EXP (v)); - err = 2 + MAX (3, err) - MPFR_GET_EXP (s); + err = 2 + MAX (3, err) - MPFR_GET_EXP (s); mpfr_sqr (u, u, GMP_RNDN); mpfr_div_2ui (u, u, 2, GMP_RNDN); /* u = log^2(x)/4 */ mpfr_add (s, s, u, GMP_RNDN); err = MAX (err, 2 + MPFR_GET_EXP (u)) - MPFR_GET_EXP (s); - err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); + err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); mpfr_const_pi (u, GMP_RNDU); mpfr_sqr (u, u, GMP_RNDN); mpfr_div_ui (u, u, 6, GMP_RNDN); /* u = pi^2/6 */ mpfr_add (s, s, u, GMP_RNDN); err = MAX (err, 3) - MPFR_GET_EXP (s); - err = 2 + MAX (-1, err); + err = 2 + MAX (-1, err); if (MPFR_CAN_ROUND (s, m - err, yp, rnd_mode)) break; @@ -590,7 +590,7 @@ mpfr_li2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) mpfr_div_2ui (u, u, 2, GMP_RNDN); /* u = log^2(1-x)/4 */ mpfr_sub (s, s, u, GMP_RNDN); err = MAX (err, - expo_l); - err = 2 + MAX (err, 3); + err = 2 + MAX (err, 3); if (MPFR_CAN_ROUND (s, m - err, yp, rnd_mode)) break; @@ -607,7 +607,7 @@ mpfr_li2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) return mpfr_check_range (y, inexact, rnd_mode); } else - /* x < -1: Li2(x) + /* x < -1: Li2(x) = S(log(1-1/x))-log^2(-x)/4-log(1-x)log(-x)/2+log^2(1-x)/4-pi^2/6 */ { int k; @@ -641,27 +641,27 @@ mpfr_li2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) mpfr_mul (w, v, u, GMP_RNDN); mpfr_div_2ui (w, w, 1, GMP_RNDN); /* w = log(-x) * log(1-x) / 2 */ mpfr_sub (s, s, w, GMP_RNDN); - err = 1 + MAX (3, MPFR_INT_CEIL_LOG2 (k+1) + 1 - MPFR_GET_EXP (s)) - + MPFR_GET_EXP (s); + err = 1 + MAX (3, MPFR_INT_CEIL_LOG2 (k+1) + 1 - MPFR_GET_EXP (s)) + + MPFR_GET_EXP (s); mpfr_sqr (w, v, GMP_RNDN); mpfr_div_2ui (w, w, 2, GMP_RNDN); /* w = log^2(-x) / 4 */ mpfr_sub (s, s, w, GMP_RNDN); - err = MAX (err, 3 + MPFR_GET_EXP(w)) - MPFR_GET_EXP (s); - err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); + err = MAX (err, 3 + MPFR_GET_EXP(w)) - MPFR_GET_EXP (s); + err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); mpfr_sqr (w, u, GMP_RNDN); mpfr_div_2ui (w, w, 2, GMP_RNDN); /* w = log^2(1-x) / 4 */ mpfr_add (s, s, w, GMP_RNDN); - err = MAX (err, 3 + MPFR_GET_EXP (w)) - MPFR_GET_EXP (s); - err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); + err = MAX (err, 3 + MPFR_GET_EXP (w)) - MPFR_GET_EXP (s); + err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); mpfr_const_pi (w, GMP_RNDU); mpfr_sqr (w, w, GMP_RNDN); mpfr_div_ui (w, w, 6, GMP_RNDN); /* w = pi^2 / 6 */ mpfr_sub (s, s, w, GMP_RNDN); err = MAX (err, 3) - MPFR_GET_EXP (s); - err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); + err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); if (MPFR_CAN_ROUND (s, m - err, yp, rnd_mode)) break; @@ -184,7 +184,7 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mp_rnd_t rnd) int ok, inex2; mp_prec_t prec = MPFR_PREC(y) + 14; MPFR_ZIV_DECL (loop); - + MPFR_ZIV_INIT (loop, prec); do { @@ -1556,7 +1556,7 @@ For negative @var{x}, the returned value is NaN. @deftypefun int mpfr_li2 (mpfr_t @var{rop}, mpfr_t @var{op}, mp_rnd_t @var{rnd_mode}) Set @var{rop} to real part of the dilogarithm of @var{op}, rounded in the -direction @var{rnd_mode}. The dilogarithm function is defined here as +direction @var{rnd_mode}. The dilogarithm function is defined here as @m{-\int_{t=0}^x log(1-t)/t dt,the integral of -log(1-t)/t from 0 to x}. @end deftypefun @@ -1818,7 +1818,7 @@ corresponding precision of @var{iop} and @var{fop} (equivalent to @deftypefunx int mpfr_remainder (mpfr_t @var{r}, mpfr_t @var{x}, mpfr_t @var{y}, mp_rnd_t @var{rnd}) @deftypefunx int mpfr_remquo (mpfr_t @var{r}, long* @var{q}, mpfr_t @var{x}, mpfr_t @var{y}, mp_rnd_t @var{rnd}) Set @var{r} to the value of @math{x - n y}, rounded according to the direction -@var{rnd}, where @math{n} is the integer quotient of @var{x} divided +@var{rnd}, where @math{n} is the integer quotient of @var{x} divided by @var{y}, defined as follows: @math{n} is rounded towards zero for @code{mpfr_fmod}, and to the nearest integer (ties rounded to even) for @code{mpfr_remainder} and @code{mpfr_remquo}. @@ -295,7 +295,7 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) /* now we have: (1) either x > 0 - (2) or x < 0 and y is an integer + (2) or x < 0 and y is an integer and in addition |x| <> 1. */ @@ -378,7 +378,7 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) return inexact; } - /* Special case (+/-2^b)^Y which could be exact. If x is negative, then + /* Special case (+/-2^b)^Y which could be exact. If x is negative, then necessarily y is a large integer. */ if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_GET_EXP (x) - 1) == 0) { @@ -386,7 +386,7 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) mp_exp_t b; int sgnx = MPFR_SIGN(x); - /* now x = +/-2^b, so x^y = (+/-)^y*2^(b*y) is exact whenever b*y is + /* now x = +/-2^b, so x^y = (+/-)^y*2^(b*y) is exact whenever b*y is an integer */ b = MPFR_GET_EXP (x) - 1; /* x = +/-2^b */ mpfr_init2 (tmp, MPFR_PREC (y) + BITS_PER_MP_LIMB); @@ -29,13 +29,13 @@ MA 02110-1301, USA. */ /* rem1 works as follows: - The first rounding mode rnd_q indicate if we are actually computing + The first rounding mode rnd_q indicate if we are actually computing a fmod (GMP_RNDZ) or a remainder/remquo (GMP_RNDN). - + Let q = x/y rounded to an integer in the direction rnd_q. Put x - q*y in rem, rounded according to rnd. If quo is not null, the value stored in *quo has the sign of q, - and agrees with q with the 2^n low order bits. + and agrees with q with the 2^n low order bits. In other words, *quo = q (mod 2^n) and *quo q >= 0. If rem is zero, then it has the sign of x. The returned 'int' is the inexact flag giving the place of rem wrt x - q*y. @@ -111,12 +111,12 @@ mpfr_rem1 (mpfr_ptr rem, long *quo, mp_rnd_t rnd_q, /* q = x/y = mx/(my*2^(ey-ex)) */ mpz_mul_2exp (my, my, ey - ex); /* divide mx by my*2^(ey-ex) */ if (rnd_q == GMP_RNDZ) - /* 0 <= |r| <= |my|, r has the same sign as mx */ - mpz_tdiv_qr (mx, r, mx, my); + /* 0 <= |r| <= |my|, r has the same sign as mx */ + mpz_tdiv_qr (mx, r, mx, my); else - /* 0 <= |r| <= |my|, r has the same sign as my */ - mpz_fdiv_qr (mx, r, mx, my); - + /* 0 <= |r| <= |my|, r has the same sign as my */ + mpz_fdiv_qr (mx, r, mx, my); + if (rnd_q == GMP_RNDN) q_is_odd = mpz_tstbit (mx, 0); if (quo) /* mx is the quotient */ @@ -137,7 +137,7 @@ mpfr_rem1 (mpfr_ptr rem, long *quo, mp_rnd_t rnd_q, /* Let X = mx*2^(ex-ey) and Y = my. Then both X and Y are integers. Assume X = R mod Y, then x = X*2^ey = R*2^ey mod (Y*2^ey=y). To be able to perform the rounding, we need the least significant - bit of the quotient, i.e., one more bit in the remainder, + bit of the quotient, i.e., one more bit in the remainder, which is obtained by dividing by 2Y. */ mpz_mul_2exp (my, my, 1); /* 2Y */ @@ -165,31 +165,31 @@ mpfr_rem1 (mpfr_ptr rem, long *quo, mp_rnd_t rnd_q, mpz_sub (r, r, my); } } - /* now 0 <= |r| < |my|, and if needed, + /* now 0 <= |r| < |my|, and if needed, q_is_odd is the least significant bit of q */ } if (mpz_cmp_ui (r, 0) == 0) inex = mpfr_set_ui (rem, 0, GMP_RNDN); - else + else { if (rnd_q == GMP_RNDN) - { - /* FIXME: the comparison 2*r < my could be done more efficiently - at the mpn level */ - mpz_mul_2exp (r, r, 1); - compare = mpz_cmpabs (r, my); - mpz_div_2exp (r, r, 1); - compare = ((compare > 0) || - ((rnd_q == GMP_RNDN) && (compare == 0) && q_is_odd)); - /* if compare != 0, we need to subtract my to r, and add 1 to quo */ - if (compare) - { - mpz_sub (r, r, my); - if (quo && (rnd_q == GMP_RNDN)) - *quo += 1; - } - } + { + /* FIXME: the comparison 2*r < my could be done more efficiently + at the mpn level */ + mpz_mul_2exp (r, r, 1); + compare = mpz_cmpabs (r, my); + mpz_div_2exp (r, r, 1); + compare = ((compare > 0) || + ((rnd_q == GMP_RNDN) && (compare == 0) && q_is_odd)); + /* if compare != 0, we need to subtract my to r, and add 1 to quo */ + if (compare) + { + mpz_sub (r, r, my); + if (quo && (rnd_q == GMP_RNDN)) + *quo += 1; + } + } inex = mpfr_set_z (rem, r, rnd); /* if ex > ey, rem should be multiplied by 2^ey, else by 2^ex */ MPFR_EXP (rem) += (ex > ey) ? ey : ex; diff --git a/tests/tfmod.c b/tests/tfmod.c index 8b115d5cf..31d973124 100644 --- a/tests/tfmod.c +++ b/tests/tfmod.c @@ -40,25 +40,25 @@ slow_fmod (mpfr_ptr r, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd) if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x) || MPFR_IS_SINGULAR (y))) { if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y) || MPFR_IS_INF (x) - || MPFR_IS_ZERO (y)) - { - MPFR_SET_NAN (r); - MPFR_RET_NAN; - } - else /* either y is Inf and x is 0 or non-special, - or x is 0 and y is non-special, - in both cases the quotient is zero. */ - return mpfr_set (r, x, rnd); + || MPFR_IS_ZERO (y)) + { + MPFR_SET_NAN (r); + MPFR_RET_NAN; + } + else /* either y is Inf and x is 0 or non-special, + or x is 0 and y is non-special, + in both cases the quotient is zero. */ + return mpfr_set (r, x, rnd); } /* regular cases */ /* if 2^(ex-1) <= |x| < 2^ex, and 2^(ey-1) <= |y| < 2^ey, then |x/y| < 2^(ex-ey+1) */ mpfr_init2 (q, - MAX (MPFR_PREC_MIN, mpfr_get_exp (x) - mpfr_get_exp (y) + 1)); + MAX (MPFR_PREC_MIN, mpfr_get_exp (x) - mpfr_get_exp (y) + 1)); mpfr_div (q, x, y, GMP_RNDZ); mpfr_trunc (q, q); /* may change inexact flag */ mpfr_prec_round (q, mpfr_get_prec (q) + mpfr_get_prec (y), GMP_RNDZ); - inexact = mpfr_mul (q, q, y, GMP_RNDZ); /* exact */ + inexact = mpfr_mul (q, q, y, GMP_RNDZ); /* exact */ inexact = mpfr_sub (r, x, q, rnd); mpfr_clear (q); return inexact; @@ -97,7 +97,7 @@ check (mpfr_t r0, mpfr_t x, mpfr_t y, mp_rnd_t rnd) mpfr_clear (r1); } -static void +static void special (void) { int inexact; @@ -114,7 +114,7 @@ special (void) { fprintf (stderr, "error : mpfr_fmod (NaN, NaN) should be exact\n"); goto error; - } + } /* NaN mod +0 is NaN */ mpfr_set_ui (y, 0, GMP_RNDN); @@ -125,8 +125,8 @@ special (void) { fprintf (stderr, "error : mpfr_fmod (NaN, +0) should be exact\n"); goto error; - } - + } + /* 3.1415 mod +0 is NaN */ mpfr_set_d (x, 3.1415, GMP_RNDN); inexact = mpfr_fmod (r, x, y, GMP_RNDN); @@ -136,7 +136,7 @@ special (void) { fprintf (stderr, "error : mpfr_fmod (3.1415, NaN) should be exact\n"); goto error; - } + } /* 3.1415 mod +Inf is 3.1415 */ mpfr_set_inf (y, 1); @@ -147,16 +147,16 @@ special (void) { fprintf (stderr, "error : mpfr_fmod (3.1415, +Inf) should be exact\n"); goto error; - } - + } + mpfr_clears (x, y, r, (void *)0); return; error: mpfr_clears (x, y, r, (void *)0); exit (1); - - + + } int @@ -191,7 +191,7 @@ main (int argc, char *argv[]) mpfr_set_prec (y, 8); mpfr_set_ui (y, 1, GMP_RNDD); mpfr_div_2ui (y, y, 3, GMP_RNDD); /* y = 1/8 */ - mpfr_set_prec (r, 123); + mpfr_set_prec (r, 123); check (r, x, y, GMP_RNDD); mpfr_clears (x, y, r, (void *)0); diff --git a/tests/tgamma.c b/tests/tgamma.c index b4b4a7ab4..a04c1205c 100644 --- a/tests/tgamma.c +++ b/tests/tgamma.c @@ -398,7 +398,7 @@ special_overflow (void) /* bug found by Kevin Rauch, 26 Oct 2007 */ mpfr_set_str (x, "1e19", 10, GMP_RNDN); inex = mpfr_gamma (x, x, GMP_RNDN); - MPFR_ASSERTN(mpfr_inf_p (x) && inex > 0); + MPFR_ASSERTN(mpfr_inf_p (x) && inex > 0); mpfr_clear (y); mpfr_clear (x); diff --git a/tests/tpow.c b/tests/tpow.c index 46098fce2..0d7bac438 100644 --- a/tests/tpow.c +++ b/tests/tpow.c @@ -665,7 +665,7 @@ special () mpfr_set_str (x, "-1.5", 10, GMP_RNDN); inex = mpfr_pow (z, x, y, GMP_RNDN); MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS (z) && inex > 0); - + /* (-n)^1.5 = NaN for n even */ mpfr_neg (y, y, GMP_RNDN); mpfr_set_str (x, "1.5", 10, GMP_RNDN); @@ -678,7 +678,7 @@ special () mpfr_set_str (y, "-1.5", 10, GMP_RNDN); inex = mpfr_pow (z, x, y, GMP_RNDN); MPFR_ASSERTN(mpfr_nan_p (z)); - + mpfr_set_emin (emin); mpfr_set_emax (emax); mpfr_clear (x); diff --git a/tests/tzeta.c b/tests/tzeta.c index f451977d9..c978bee59 100644 --- a/tests/tzeta.c +++ b/tests/tzeta.c @@ -381,7 +381,7 @@ main (int argc, char *argv[]) printf ("got "); mpfr_dump (z); exit (1); } - + /* bug reported by Kevin Rauch on 26 Oct 2007 */ mpfr_set_prec (s, 128); mpfr_set_prec (z, 128); |