diff options
| -rw-r--r-- | gmp_op.c | 25 | ||||
| -rw-r--r-- | log10.c | 102 | ||||
| -rw-r--r-- | log1p.c | 78 | ||||
| -rw-r--r-- | log2.c | 96 | ||||
| -rw-r--r-- | pow.c | 78 | ||||
| -rw-r--r-- | pow_si.c | 63 | ||||
| -rw-r--r-- | pow_ui.c | 64 | ||||
| -rw-r--r-- | pow_z.c | 57 | ||||
| -rw-r--r-- | strtofr.c | 5 | ||||
| -rw-r--r-- | sum.c | 47 | ||||
| -rw-r--r-- | ui_pow_ui.c | 36 | ||||
| -rw-r--r-- | zeta.c | 221 |
12 files changed, 451 insertions, 421 deletions
@@ -143,9 +143,10 @@ int mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode) { mpfr_t t,q; - mp_prec_t p = MPFR_PREC(y)+10; + mp_prec_t p; mp_exp_t err; int res; + MPFR_ZIV_DECL (loop); if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) { @@ -171,8 +172,11 @@ mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode) } } + p = MPFR_PREC (y) + 10; mpfr_init2 (t, p); mpfr_init2 (q, p); + + MPFR_ZIV_INIT (loop, p); for (;;) { res = mpfr_set_q (q, z, GMP_RNDN); /* Error <= 1/2 ulp(q) */ @@ -192,18 +196,18 @@ mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode) if (MPFR_LIKELY (!MPFR_IS_ZERO (t))) { err = (mp_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0); - res = mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, - MPFR_PREC(y) + (rnd_mode == GMP_RNDN) ); - if (MPFR_LIKELY (res != 0)) /* We can round! */ + if (MPFR_LIKELY (mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + MPFR_PREC(y) + (rnd_mode == GMP_RNDN)))) { - res = mpfr_set(y, t, rnd_mode); + res = mpfr_set (y, t, rnd_mode); break; } } - p += BITS_PER_MP_LIMB; /* Next precision if we continue */ + MPFR_ZIV_NEXT (loop, p); mpfr_set_prec (t, p); mpfr_set_prec (q, p); } + MPFR_ZIV_FREE (loop); mpfr_clear (t); mpfr_clear (q); return res; @@ -213,9 +217,10 @@ int mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode) { mpfr_t t,q; - mp_prec_t p = MPFR_PREC (y)+10; + mp_prec_t p; int res; mp_exp_t err; + MPFR_ZIV_DECL (loop); if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) { @@ -246,8 +251,11 @@ mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode) } } + p = MPFR_PREC (y) + 10; mpfr_init2 (t, p); mpfr_init2 (q, p); + + MPFR_ZIV_INIT (loop, p); for(;;) { res = mpfr_set_q(q, z, GMP_RNDN); /* Error <= 1/2 ulp(q) */ @@ -275,10 +283,11 @@ mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode) break; } } - p += BITS_PER_MP_LIMB; /* Next precision if we continue */ + MPFR_ZIV_NEXT (loop, p); mpfr_set_prec (t, p); mpfr_set_prec (q, p); } + MPFR_ZIV_FREE (loop); mpfr_clear (t); mpfr_clear (q); return res; @@ -1,6 +1,6 @@ /* mpfr_log10 -- logarithm in base 10. -Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -30,89 +30,81 @@ MA 02111-1307, USA. */ int mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) { - int inexact = 0; + int inexact; + MPFR_SAVE_EXPO_DECL (expo); /* If a is NaN, the result is NaN */ - if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(a) )) + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a))) { - if (MPFR_IS_NAN(a)) + if (MPFR_IS_NAN (a)) { - MPFR_SET_NAN(r); + MPFR_SET_NAN (r); MPFR_RET_NAN; } /* check for infinity before zero */ - else if (MPFR_IS_INF(a)) + else if (MPFR_IS_INF (a)) { - if (MPFR_IS_NEG(a)) + if (MPFR_IS_NEG (a)) /* log10(-Inf) = NaN */ { - MPFR_SET_NAN(r); + MPFR_SET_NAN (r); MPFR_RET_NAN; } else /* log10(+Inf) = +Inf */ { - MPFR_SET_INF(r); - MPFR_SET_POS(r); - MPFR_RET(0); /* exact */ + MPFR_SET_INF (r); + MPFR_SET_POS (r); + MPFR_RET (0); /* exact */ } } else /* a = 0 */ { - MPFR_ASSERTD(MPFR_IS_ZERO(a)); - MPFR_SET_INF(r); - MPFR_SET_NEG(r); - MPFR_RET(0); /* log10(0) is an exact -infinity */ + MPFR_ASSERTD (MPFR_IS_ZERO (a)); + MPFR_SET_INF (r); + MPFR_SET_NEG (r); + MPFR_RET (0); /* log10(0) is an exact -infinity */ } } /* If a is negative, the result is NaN */ - if (MPFR_UNLIKELY( MPFR_IS_NEG(a) )) + if (MPFR_UNLIKELY (MPFR_IS_NEG (a))) { - MPFR_SET_NAN(r); + MPFR_SET_NAN (r); MPFR_RET_NAN; } /* If a is 1, the result is 0 */ if (mpfr_cmp_ui (a, 1) == 0) { - MPFR_SET_ZERO(r); - MPFR_SET_POS(r); - MPFR_RET(0); /* result is exact */ + MPFR_SET_ZERO (r); + MPFR_SET_POS (r); + MPFR_RET (0); /* result is exact */ } - - /* Useless due to mpfr_set - MPFR_CLEAR_FLAGS(r);*/ + MPFR_SAVE_EXPO_MARK (expo); /* General case */ { /* Declaration of the intermediary variable */ mpfr_t t, tt; - int ok; - + MPFR_ZIV_DECL (loop); /* Declaration of the size variable */ - mp_prec_t Nx = MPFR_PREC(a); /* Precision of input variable */ + mp_prec_t Nx = MPFR_PREC(a); /* Precision of input variable */ mp_prec_t Ny = MPFR_PREC(r); /* Precision of output variable */ - - mp_prec_t Nt; /* Precision of the intermediary variable */ - long int err; /* Precision of error */ + mp_prec_t Nt; /* Precision of the intermediary variable */ + mp_exp_t err; /* Precision of error */ /* compute the precision of intermediary variable */ - Nt = MAX(Nx, Ny); + Nt = MAX (Nx, Ny); /* the optimal number of bits : see algorithms.ps */ Nt = Nt + 4+ MPFR_INT_CEIL_LOG2 (Nt); /* initialise of intermediary variables */ - mpfr_init (t); - mpfr_init (tt); - + mpfr_init2 (t, Nt); + mpfr_init2 (tt, Nt); /* First computation of log10 */ - do { - - /* reactualisation of the precision */ - mpfr_set_prec (t, Nt); - mpfr_set_prec (tt, Nt); - + MPFR_ZIV_INIT (loop, Nt); + for(;;) { /* compute log10 */ mpfr_set_ui (t, 10, GMP_RNDN); /* 10 */ mpfr_log (t, t, GMP_RNDD); /* log(10) */ @@ -121,20 +113,25 @@ mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) /* estimation of the error */ err = Nt - 4; - - ok = mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, - Ny + (rnd_mode == GMP_RNDN)); - - /* log10(10^n) is exact */ - if ((MPFR_SIGN(t) > 0) && mpfr_integer_p (t)) - if (mpfr_ui_pow_ui (tt, 10, (unsigned long) (mpfr_get_d1 (t) + 0.5), - GMP_RNDN) == 0) - if (mpfr_cmp (a, tt) == 0) - ok = 1; + if (MPFR_LIKELY (mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN)))) + break; + + /* log10(10^n) is exact: + FIXME: Can we have 10^n exactly representable as a mpfr_t + but n can't fit an unsigned long? */ + if (MPFR_IS_POS (t) + && mpfr_integer_p (t) && mpfr_fits_ulong_p (t, GMP_RNDN) + && !mpfr_ui_pow_ui (tt, 10, mpfr_get_ui (t, GMP_RNDN), GMP_RNDN) + && mpfr_cmp (a, tt) == 0) + break; /* actualisation of the precision */ - Nt += 10; - } while ((err < 0) || !ok); + MPFR_ZIV_NEXT (loop, Nt); + mpfr_set_prec (t, Nt); + mpfr_set_prec (tt, Nt); + } + MPFR_ZIV_FREE (loop); inexact = mpfr_set (r, t, rnd_mode); @@ -142,5 +139,6 @@ mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) mpfr_clear (tt); } - return inexact; + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_check_range (r, inexact, rnd_mode); } @@ -1,6 +1,6 @@ /* mpfr_log1p -- Compute log(1+x) -Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -28,36 +28,37 @@ MA 02111-1307, USA. */ int mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) { - int comp, inexact = 0; + int comp, inexact; + MPFR_SAVE_EXPO_DECL (expo); - if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x))) + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) { - if (MPFR_IS_NAN(x)) + if (MPFR_IS_NAN (x)) { - MPFR_SET_NAN(y); + MPFR_SET_NAN (y); MPFR_RET_NAN; } /* check for inf or -inf (result is not defined) */ - else if (MPFR_IS_INF(x)) + else if (MPFR_IS_INF (x)) { - if (MPFR_IS_POS(x)) + if (MPFR_IS_POS (x)) { - MPFR_SET_INF(y); - MPFR_SET_POS(y); - MPFR_RET(0); + MPFR_SET_INF (y); + MPFR_SET_POS (y); + MPFR_RET (0); } else { - MPFR_SET_NAN(y); + MPFR_SET_NAN (y); MPFR_RET_NAN; } } else /* x is zero */ { - MPFR_ASSERTD(MPFR_IS_ZERO(x)); - MPFR_SET_ZERO(y); /* log1p(+/- 0) = +/- 0 */ - MPFR_SET_SAME_SIGN(y, x); - MPFR_RET(0); + MPFR_ASSERTD (MPFR_IS_ZERO (x)); + MPFR_SET_ZERO (y); /* log1p(+/- 0) = +/- 0 */ + MPFR_SET_SAME_SIGN (y, x); + MPFR_RET (0); } } @@ -68,60 +69,61 @@ mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) if (comp == 0) /* x=0: log1p(-1)=-inf (division by zero) */ { - MPFR_SET_INF(y); - MPFR_SET_NEG(y); - MPFR_RET(0); + MPFR_SET_INF (y); + MPFR_SET_NEG (y); + MPFR_RET (0); } - MPFR_SET_NAN(y); + MPFR_SET_NAN (y); MPFR_RET_NAN; } - MPFR_CLEAR_FLAGS(y); + MPFR_SAVE_EXPO_MARK (expo); /* General case */ { /* Declaration of the intermediary variable */ mpfr_t t; - /* Declaration of the size variable */ mp_prec_t Nx = MPFR_PREC(x); /* Precision of input variable */ - mp_prec_t Ny = MPFR_PREC(y); /* Precision of input variable */ - - mp_prec_t Nt; /* Precision of the intermediary variable */ - long int err; /* Precision of error */ - + mp_prec_t Ny = MPFR_PREC(y); /* Precision of input variable */ + mp_prec_t Nt; /* Precision of the intermediary variable */ + mp_exp_t err; /* Precision of error */ + MPFR_ZIV_DECL (loop); + /* compute the precision of intermediary variable */ - Nt = MAX(Nx,Ny); + Nt = MAX (Nx, Ny); /* the optimal number of bits : see algorithms.ps */ Nt = Nt + 5 + MPFR_INT_CEIL_LOG2 (Nt); /* initialise of intermediary variable */ - mpfr_init (t); + mpfr_init2 (t, Nt); /* First computation of cosh */ - do + MPFR_ZIV_INIT (loop, Nt); + for (;;) { - /* reactualisation of the precision */ - mpfr_set_prec (t, Nt); - /* compute log1p */ - mpfr_add_ui (t, x, 1, GMP_RNDN); /* 1+x */ + mpfr_add_ui (t, x, 1, GMP_RNDN); /* 1+x */ mpfr_log (t, t, GMP_RNDN); /* log(1+x)*/ /* estimation of the error */ /*err=Nt-(__gmpfr_ceil_log2(1+pow(2,1-MPFR_GET_EXP(t))));*/ err = Nt - (MAX (1 - MPFR_GET_EXP (t), 0) + 1); + if (MPFR_LIKELY (mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN)))) + break; + /* actualisation of the precision */ - Nt += 10; + MPFR_ZIV_NEXT (loop, Nt); + mpfr_set_prec (t, Nt); } - while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, - Ny + (rnd_mode == GMP_RNDN))); - + MPFR_ZIV_FREE (loop); inexact = mpfr_set (y, t, rnd_mode); mpfr_clear (t); } - return inexact; + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_check_range (y, inexact, rnd_mode); } @@ -1,6 +1,6 @@ /* mpfr_log2 -- log base 2 -Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -23,110 +23,109 @@ MA 02111-1307, USA. */ #include "mpfr-impl.h" /* The computation of r=log2(a) - - r=log2(a)=log(a)/log(2) - */ + r=log2(a)=log(a)/log(2) */ int mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) { - int inexact = 0; + int inexact; + MPFR_SAVE_EXPO_DECL (expo); - if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(a) )) + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a))) { /* If a is NaN, the result is NaN */ - if (MPFR_IS_NAN(a)) + if (MPFR_IS_NAN (a)) { - MPFR_SET_NAN(r); + MPFR_SET_NAN (r); MPFR_RET_NAN; } /* check for infinity before zero */ - else if (MPFR_IS_INF(a)) + else if (MPFR_IS_INF (a)) { - if (MPFR_IS_NEG(a)) + if (MPFR_IS_NEG (a)) /* log(-Inf) = NaN */ { - MPFR_SET_NAN(r); + MPFR_SET_NAN (r); MPFR_RET_NAN; } else /* log(+Inf) = +Inf */ { - MPFR_SET_INF(r); - MPFR_SET_POS(r); - MPFR_RET(0); + MPFR_SET_INF (r); + MPFR_SET_POS (r); + MPFR_RET (0); } } else /* a is zero */ { - MPFR_ASSERTD(MPFR_IS_ZERO(a)); - MPFR_SET_INF(r); - MPFR_SET_NEG(r); - MPFR_RET(0); /* log2(0) is an exact -infinity */ + MPFR_ASSERTD (MPFR_IS_ZERO (a)); + MPFR_SET_INF (r); + MPFR_SET_NEG (r); + MPFR_RET (0); /* log2(0) is an exact -infinity */ } } /* If a is negative, the result is NaN */ - if (MPFR_UNLIKELY(MPFR_IS_NEG(a))) + if (MPFR_UNLIKELY (MPFR_IS_NEG (a))) { - MPFR_SET_NAN(r); + MPFR_SET_NAN (r); MPFR_RET_NAN; } /* If a is 1, the result is 0 */ - if (mpfr_cmp_ui(a, 1) == 0) + if (MPFR_UNLIKELY (mpfr_cmp_ui (a, 1) == 0)) { - MPFR_SET_ZERO(r); - MPFR_SET_POS(r); - MPFR_RET(0); /* only "normal" case where the result is exact */ + MPFR_SET_ZERO (r); + MPFR_SET_POS (r); + MPFR_RET (0); /* only "normal" case where the result is exact */ } - /* If a is integer, log2(a) is exact*/ - if (mpfr_cmp_ui_2exp (a, 1, MPFR_GET_EXP (a) - 1) == 0) + /* If a is 2^N, log2(a) is exact*/ + if (MPFR_UNLIKELY (mpfr_cmp_ui_2exp (a, 1, MPFR_GET_EXP (a) - 1) == 0)) return mpfr_set_si(r, MPFR_GET_EXP (a) - 1, rnd_mode); + MPFR_SAVE_EXPO_MARK (expo); + /* General case */ { /* Declaration of the intermediary variable */ mpfr_t t, tt; - /* Declaration of the size variable */ mp_prec_t Nx = MPFR_PREC(a); /* Precision of input variable */ mp_prec_t Ny = MPFR_PREC(r); /* Precision of input variable */ - - mp_prec_t Nt; /* Precision of the intermediary variable */ - long int err; /* Precision of error */ - + mp_prec_t Nt; /* Precision of the intermediary variable */ + mp_exp_t err; /* Precision of error */ + MPFR_ZIV_DECL (loop); /* compute the precision of intermediary variable */ - Nt=MAX(Nx,Ny); + Nt = MAX (Nx, Ny); /* the optimal number of bits : see algorithms.ps */ - Nt=Nt + 3 + MPFR_INT_CEIL_LOG2 (Nt); + Nt = Nt + 3 + MPFR_INT_CEIL_LOG2 (Nt); /* initialise of intermediary variable */ - mpfr_init(t); - mpfr_init(tt); - + mpfr_init2 (t, Nt); + mpfr_init2 (tt, Nt); /* First computation of log2 */ - do + MPFR_ZIV_INIT (loop, Nt); + for (;;) { - /* reactualisation of the precision */ - mpfr_set_prec(t,Nt); - mpfr_set_prec(tt,Nt); - /* compute log2 */ mpfr_const_log2(t,GMP_RNDD); /* log(2) */ mpfr_log(tt,a,GMP_RNDN); /* log(a) */ mpfr_div(t,tt,t,GMP_RNDN); /* log(a)/log(2) */ - + /* estimation of the error */ - err=Nt-3; - + err = Nt-3; + if (MPFR_LIKELY (mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN)))) + break; + /* actualisation of the precision */ - Nt += 10; + MPFR_ZIV_NEXT (loop, Nt); + mpfr_set_prec (t, Nt); + mpfr_set_prec (tt, Nt); } - while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, - Ny + (rnd_mode == GMP_RNDN))); + MPFR_ZIV_FREE (loop); inexact = mpfr_set (r, t, rnd_mode); @@ -134,5 +133,6 @@ mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) mpfr_clear (tt); } - return inexact; + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_check_range (r, inexact, rnd_mode); } @@ -199,8 +199,9 @@ int mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) { int inexact = 1; + MPFR_SAVE_EXPO_DECL (expo); - if (MPFR_ARE_SINGULAR(x,y)) + if (MPFR_ARE_SINGULAR (x, y)) { /* pow(x, 0) returns 1 for any x, even a NaN. */ if (MPFR_UNLIKELY( MPFR_IS_ZERO(y) )) @@ -336,69 +337,63 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) MPFR_RET_NAN; } - MPFR_CLEAR_FLAGS(z); + MPFR_SAVE_EXPO_MARK (expo); /* General case */ { /* Declaration of the intermediary variable */ - mpfr_t t, te, ti; - int loop = 0, ok; - + mpfr_t t; + int loop = 0; /* Declaration of the size variable */ mp_prec_t Nx = MPFR_PREC(x); /* Precision of input variable */ mp_prec_t Ny = MPFR_PREC(y); /* Precision of input variable */ mp_prec_t Nz = MPFR_PREC(z); /* Precision of output variable */ - - mp_prec_t Nt; /* Precision of the intermediary variable */ - long int err; /* Precision of error */ + mp_prec_t Nt; /* Precision of the intermediary variable */ + mp_exp_t err, exp_te; /* Precision of error */ + MPFR_ZIV_DECL (ziv_loop); /* compute the precision of intermediary variable */ - Nt = MAX(Nx, Ny); - Nt = MAX(Nt, Nz); /* take account of the output precision too! */ + Nt = MAX (Nx, Ny); + Nt = MAX (Nt, Nz); /* take account of the output precision too! */ /* the optimal number of bits : see algorithms.ps */ Nt = Nt + 5 + MPFR_INT_CEIL_LOG2 (Nt); /* initialise of intermediary variable */ - mpfr_init (t); - mpfr_init (ti); - mpfr_init (te); + mpfr_init2 (t, Nt); - do + MPFR_ZIV_INIT (ziv_loop, Nt); + for (;;) { loop ++; - /* reactualisation of the precision */ - mpfr_set_prec (t, Nt); - mpfr_set_prec (ti, Nt); - mpfr_set_prec (te, Nt); - /* compute exp(y*ln(x)) */ - mpfr_log (ti, x, GMP_RNDU); /* ln(x) */ - mpfr_mul (te, y, ti, GMP_RNDU); /* y*ln(x) */ - mpfr_exp (t, te, GMP_RNDN); /* exp(y*ln(x))*/ - - MPFR_ASSERTN(MPFR_IS_FP(te)); - MPFR_ASSERTN(MPFR_NOTZERO(te)); - /* otherwise MPFR_EXP(te) below doesn't exist */ + mpfr_log (t, x, GMP_RNDU); /* ln(x) */ + mpfr_mul (t, y, t, GMP_RNDU); /* y*ln(x) */ + exp_te = MPFR_GET_EXP (t); /* FIXME: May overflow */ + mpfr_exp (t, t, GMP_RNDN); /* exp(y*ln(x))*/ + /* FIXME: May overflow */ /* estimate of the error -- see pow function in algorithms.ps */ - err = Nt - (MPFR_GET_EXP (te) + 3); - - /* actualisation of the precision */ - Nt += 10; - - ok = mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, - Nz + (rnd_mode == GMP_RNDN)); + err = Nt - (exp_te + 3); + if (MPFR_LIKELY (mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Nz + (rnd_mode == GMP_RNDN)))) + break; /* check exact power */ - if (ok == 0 && loop == 1) + if (loop == 1) { - ok = mpfr_pow_is_exact (x, y); - if (ok) - inexact = 0; - } + if (mpfr_pow_is_exact (x, y)) + { + inexact = 0; + break; + } + } + + /* reactualisation of the precision */ + MPFR_ZIV_NEXT (ziv_loop, Nt); + mpfr_set_prec (t, Nt); } - while (err < 0 || ok == 0); + MPFR_ZIV_FREE (ziv_loop); if (inexact) inexact = mpfr_set (z, t, rnd_mode); @@ -406,9 +401,8 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) mpfr_set (z, t, GMP_RNDN); mpfr_clear (t); - mpfr_clear (ti); - mpfr_clear (te); } - return inexact; + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_check_range (z, inexact, rnd_mode); } @@ -35,7 +35,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) return mpfr_pow_ui (y, x, n, rnd_mode); else { - if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) { if (MPFR_IS_NAN (x)) { @@ -53,7 +53,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) } else /* x is zero */ { - MPFR_ASSERTD (MPFR_IS_ZERO(x)); + MPFR_ASSERTD (MPFR_IS_ZERO (x)); MPFR_SET_INF(y); if (MPFR_IS_POS (x) || ((unsigned) n & 1) == 0) MPFR_SET_POS (y); @@ -62,7 +62,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) MPFR_RET(0); } } - MPFR_CLEAR_FLAGS(y); + MPFR_CLEAR_FLAGS (y); /* detect exact powers: x^(-n) is exact iff x is a power of 2 */ if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0) @@ -73,13 +73,13 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) expx --; MPFR_ASSERTD (n < 0); /* Warning n*expx may overflow! - Many systems abort with LONG_MIN / 1 or LONG_MIN/-1*/ + Some systems abort with LONG_MIN / 1 or LONG_MIN/-1*/ if (n != -1 && expx > 0 && -expx < MPFR_EXP_MIN / (-n)) MPFR_EXP (y) = MPFR_EMIN_MIN - 1; /* Underflow */ else if (n != -1 && expx < 0 && -expx > MPFR_EXP_MAX / (-n)) MPFR_EXP (y) = MPFR_EMAX_MAX + 1; /* Overflow */ else - MPFR_EXP(y) += n * expx; + MPFR_EXP (y) += n * expx; return mpfr_check_range (y, 0, rnd_mode); } @@ -89,18 +89,17 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) { /* Declaration of the intermediary variable */ mpfr_t t; - /* Declaration of the size variable */ - mp_prec_t Nx = MPFR_PREC(x); /* Precision of input variable */ - mp_prec_t Ny = MPFR_PREC(y); /* Precision of output variable */ - - mp_prec_t Nt; /* Precision of the intermediary variable */ - long int err; /* Precision of error */ + mp_prec_t Nx = MPFR_PREC (x); /* Precision of input variable */ + mp_prec_t Ny = MPFR_PREC (y); /* Precision of output variable */ + mp_prec_t Nt; /* Precision of the intermediary variable */ + mp_exp_t err; /* Precision of error */ int inexact; MPFR_SAVE_EXPO_DECL (expo); + MPFR_ZIV_DECL (loop); /* compute the precision of intermediary variable */ - Nt = MAX(Nx,Ny); + Nt = MAX (Nx, Ny); /* the optimal number of bits : see algorithms.ps */ Nt = Nt + 3 + MPFR_INT_CEIL_LOG2 (Nt); @@ -108,24 +107,28 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) /* initialise of intermediary variable */ mpfr_init2 (t, Nt); - - do { - /* reactualisation of the precision */ - mpfr_set_prec (t, Nt); - - /* compute 1/(x^n) n>0*/ - mpfr_pow_ui (t, x, (unsigned long int) n, GMP_RNDN); - inexact = MPFR_IS_ZERO (t) || MPFR_IS_INF (t); - mpfr_ui_div (t, 1, t, GMP_RNDN); - inexact = inexact || MPFR_IS_ZERO (t) || MPFR_IS_INF (t); - /* error estimate -- see pow function in algorithms.ps */ - err = Nt - 3; - - /* actualisation of the precision */ - Nt += BITS_PER_MP_LIMB; - } while (inexact == 0 && - !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, - Ny + (rnd_mode == GMP_RNDN))); + + MPFR_ZIV_INIT (loop, Nt); + for (;;) + { + /* compute 1/(x^n) n>0*/ + mpfr_pow_ui (t, x, (unsigned long int) n, GMP_RNDN); + inexact = MPFR_IS_ZERO (t) || MPFR_IS_INF (t); + mpfr_ui_div (t, 1, t, GMP_RNDN); + inexact = inexact || MPFR_IS_ZERO (t) || MPFR_IS_INF (t); + + /* error estimate -- see pow function in algorithms.ps */ + err = Nt - 3; + if (MPFR_LIKELY (inexact != 0 + || mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN)))) + break; + + /* actualisation of the precision */ + Nt += BITS_PER_MP_LIMB; + mpfr_set_prec (t, Nt); + } + MPFR_ZIV_FREE (loop); inexact = mpfr_set (y, t, rnd_mode); mpfr_clear (t); @@ -33,6 +33,7 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) int inexact; mp_rnd_t rnd1; MPFR_SAVE_EXPO_DECL (expo); + MPFR_ZIV_DECL (loop); if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (y))) { @@ -81,41 +82,48 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) /* y^2 = sqr(y) */ return mpfr_sqr (x, y, rnd); } - /* MPFR_CLEAR_FLAGS useless due to mpfr_set */ + /* Augment exponent range */ MPFR_SAVE_EXPO_MARK (expo); __gmpfr_emin -= 3; /* So that we can check for underflow properly */ - prec = MPFR_PREC (x); - mpfr_init2 (res, prec + 9); + /* setup initial precision */ + prec = MPFR_PREC (x) + 3 + BITS_PER_MP_LIMB; + mpfr_init2 (res, prec); rnd1 = MPFR_IS_POS (y) ? GMP_RNDU : GMP_RNDD; /* away */ - do { - int i; - prec += 3; - for (m = n, i = 0; m; i++, m >>= 1) - prec++; - /* now 2^(i-1) <= n < 2^i */ - mpfr_set_prec (res, prec); - err = prec <= (mpfr_prec_t) i ? 0 : prec - (mpfr_prec_t) i; - MPFR_ASSERTD (i >= 1); - mpfr_clear_overflow (); - mpfr_clear_underflow (); - /* First step: compute square from y */ - inexact = mpfr_sqr (res, y, GMP_RNDU); - if (n & (1UL << (i-2))) - inexact |= mpfr_mul (res, res, y, rnd1); - for (i -= 3; i >= 0 && !mpfr_overflow_p () && !mpfr_underflow_p (); i--) - { - inexact |= mpfr_sqr (res, res, GMP_RNDU); - if (n & (1UL << i)) - inexact |= mpfr_mul (res, res, y, rnd1); - } - } - while (inexact && !mpfr_overflow_p () && !mpfr_underflow_p () - && !mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ, - MPFR_PREC(x) + (rnd == GMP_RNDN))); + MPFR_ZIV_INIT (loop, prec); + for (;;) + { + int i; + for (m = n, i = 0; m; i++, m >>= 1); + /* now 2^(i-1) <= n < 2^i */ + err = prec <= (mpfr_prec_t) i ? 0 : prec - (mpfr_prec_t) i; + MPFR_ASSERTD (i >= 1); + mpfr_clear_overflow (); + mpfr_clear_underflow (); + /* First step: compute square from y */ + inexact = mpfr_sqr (res, y, GMP_RNDU); + if (n & (1UL << (i-2))) + inexact |= mpfr_mul (res, res, y, rnd1); + for (i -= 3; i >= 0 && !mpfr_overflow_p () && !mpfr_underflow_p (); i--) + { + inexact |= mpfr_sqr (res, res, GMP_RNDU); + if (n & (1UL << i)) + inexact |= mpfr_mul (res, res, y, rnd1); + } + if (MPFR_LIKELY (inexact == 0 + || mpfr_overflow_p () || mpfr_underflow_p () + || mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ, + MPFR_PREC(x) + (rnd == GMP_RNDN)))) + break; + /* Actualisation of the precision */ + MPFR_ZIV_NEXT (loop, prec); + mpfr_set_prec (res, prec); + } + MPFR_ZIV_FREE (loop); + inexact = mpfr_set (x, res, rnd); mpfr_clear (res); @@ -30,31 +30,32 @@ mpfr_pow_pos_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd) int inexact; mp_rnd_t rnd1; mpz_t absz; + mp_size_t size_z; + MPFR_ZIV_DECL (loop); MPFR_ASSERTD (mpz_sgn (z) != 0); if (MPFR_UNLIKELY (mpz_cmpabs_ui (z, 1) == 0)) return mpfr_set (y, x, rnd); - prec = MPFR_PREC (x); - mpfr_init2 (res, prec + 9); - rnd1 = MPFR_IS_POS (x) ? GMP_RNDU : GMP_RNDD; /* away */ absz[0] = z[0]; SIZ (absz) = ABS(SIZ(absz)); /* Hack to get abs(z) */ + MPFR_MPZ_SIZEINBASE2 (size_z, z); - do { - mp_size_t i; - - MPFR_MPZ_SIZEINBASE2 (i, z); + prec = MPFR_PREC (x) + 3 + size_z; + mpfr_init2 (res, prec); + + MPFR_ZIV_INIT (loop, prec); + for (;;) { + mp_size_t i = size_z; /* now 2^(i-1) <= z < 2^i */ - prec += 3 + i; - mpfr_set_prec (res, prec); err = prec <= (mpfr_prec_t) i ? 0 : prec - (mpfr_prec_t) i; MPFR_ASSERTD (i >= 2); mpfr_clear_overflow (); mpfr_clear_underflow (); - /* First step: compute square from y */ + + /* First step: compute square from y */ inexact = mpfr_mul (res, x, x, GMP_RNDU); if (mpz_tstbit (absz, i-2)) inexact |= mpfr_mul (res, res, x, rnd1); @@ -63,11 +64,20 @@ mpfr_pow_pos_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd) inexact |= mpfr_sqr (res, res, GMP_RNDU); if (mpz_tstbit (absz, i)) inexact |= mpfr_mul (res, res, x, rnd1); - } /* for */ - } while (inexact && !mpfr_overflow_p() && !mpfr_underflow_p () - && !mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ, - MPFR_PREC (y) + (rnd == GMP_RNDN))); - + } + + if (MPFR_LIKELY (inexact == 0 + || mpfr_overflow_p () || mpfr_underflow_p () + || mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ, + MPFR_PREC (y) + (rnd == GMP_RNDN)))) + break; + + /* Actualisation of the precision */ + MPFR_ZIV_NEXT (loop, prec); + mpfr_set_prec (res, prec); + } + MPFR_ZIV_FREE (loop); + inexact = mpfr_set (y, res, rnd); mpfr_clear (res); @@ -189,7 +199,8 @@ mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd) /* Declaration of the intermediary variable */ mpfr_t t; mp_prec_t Nt; /* Precision of the intermediary variable */ - + MPFR_ZIV_DECL (loop); + /* compute the precision of intermediary variable */ Nt = MAX (MPFR_PREC (x), MPFR_PREC (y)); @@ -198,7 +209,8 @@ mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd) /* initialise of intermediary variable */ mpfr_init2 (t, Nt); - + + MPFR_ZIV_INIT (loop, Nt); for (;;) { /* compute 1/(x^n) n>0 */ mpfr_pow_pos_z (t, x, z, GMP_RNDN); @@ -206,15 +218,16 @@ mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd) mpfr_ui_div (t, 1, t, GMP_RNDN); inexact = inexact || MPFR_IS_ZERO (t) || MPFR_IS_INF (t); - if (inexact != 0 - || mpfr_can_round (t, Nt - 3, GMP_RNDN, GMP_RNDZ, - MPFR_PREC (y) + (rnd == GMP_RNDN))) + if (MPFR_LIKELY (inexact != 0 + || mpfr_can_round (t, Nt - 3, GMP_RNDN, GMP_RNDZ, + MPFR_PREC (y) + (rnd == GMP_RNDN)))) break; /* actualisation of the precision */ - Nt += BITS_PER_MP_LIMB; + MPFR_ZIV_NEXT (loop, Nt); mpfr_set_prec (t, Nt); } - + MPFR_ZIV_FREE (loop); + inexact = mpfr_set (y, t, rnd); mpfr_clear (t); } @@ -409,6 +409,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) size_t pstr_size; mp_size_t ysize, real_ysize; int res, err; + MPFR_ZIV_DECL (loop); TMP_DECL (marker); /* determine the minimal precision for the computation */ @@ -416,6 +417,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) /* compute y as long as rounding is not possible */ TMP_MARK(marker); + MPFR_ZIV_INIT (loop, prec); for (;;) { /* initialize y to the value of 0.mant_s[0]...mant_s[pr-1] */ @@ -633,8 +635,9 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) break; /* update the prec for next loop */ - prec += BITS_PER_MP_LIMB; + MPFR_ZIV_NEXT (loop, prec); } /* loop */ + MPFR_ZIV_FREE (loop); /* round y */ if (mpfr_round_raw (MPFR_MANT (x), result, @@ -193,13 +193,13 @@ static int mpfr_list_sum_once (mpfr_ptr ret, mpfr_srcptr const tab[], int mpfr_sum (mpfr_ptr ret, mpfr_ptr const tab[], unsigned long n, mp_rnd_t rnd) { - mp_prec_t initial_f, current_f; + mp_prec_t prec; unsigned long k; mpfr_srcptr *perm; - unsigned int guard_digits; - unsigned int initial_guard_digits; int error_trap; mpfr_t cur_sum; + MPFR_ZIV_DECL (loop); + MPFR_SAVE_EXPO_DECL (expo); TMP_DECL(marker); TMP_MARK(marker); @@ -209,30 +209,39 @@ int mpfr_sum (mpfr_ptr ret, mpfr_ptr const tab[], unsigned long n, return 0; } + /* Sort */ perm = (mpfr_srcptr *) TMP_ALLOC(n * sizeof(mpfr_srcptr)); - mpfr_count_sort (tab, n, perm); - initial_f = MAX (MPFR_PREC(tab[0]), MPFR_PREC(ret)); + /* Initial precision */ + prec = MAX (MPFR_PREC (tab[0]), MPFR_PREC (ret)); k = MPFR_INT_CEIL_LOG2 (n) + 1; - mpfr_init2 (cur_sum, initial_f); - initial_guard_digits = k + 2; - guard_digits = initial_guard_digits; - do - { - current_f = initial_f + guard_digits; - mpfr_set_prec (cur_sum, current_f); - error_trap = mpfr_list_sum_once (cur_sum, perm, n, - current_f + k); - guard_digits *= 2; + prec += k + 2; + mpfr_init2 (cur_sum, prec); + + MPFR_SAVE_EXPO_MARK (expo); + + /* Ziv Loop */ + MPFR_ZIV_INIT (loop, prec); + for (;;) { + error_trap = mpfr_list_sum_once (cur_sum, perm, n, prec + k); + if (MPFR_LIKELY (error_trap == 0 + || mpfr_can_round (cur_sum, + MPFR_GET_EXP (cur_sum) - prec + 2, + GMP_RNDN, rnd, MPFR_PREC (ret)))) + break; + MPFR_ZIV_NEXT (loop, prec); + mpfr_set_prec (cur_sum, prec); } - while ((error_trap != 0) && - !(mpfr_can_round (cur_sum, MPFR_GET_EXP(cur_sum) - current_f + 2, - GMP_RNDN, rnd, MPFR_PREC(ret)))); + MPFR_ZIV_FREE (loop); + error_trap |= mpfr_set (ret, cur_sum, rnd); mpfr_clear (cur_sum); TMP_FREE(marker); - return error_trap; + + MPFR_SAVE_EXPO_FREE (expo); + error_trap |= mpfr_check_range (ret, 0, rnd); + return error_trap; /* It doesn't return the ternary value */ } diff --git a/ui_pow_ui.c b/ui_pow_ui.c index 478d94f2d..1303a1884 100644 --- a/ui_pow_ui.c +++ b/ui_pow_ui.c @@ -1,6 +1,7 @@ /* mpfr_ui_pow_ui -- compute the power beetween two machine integer -Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005 + Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -25,15 +26,15 @@ int mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int y, unsigned long int n, mp_rnd_t rnd) { - long int err; + mp_exp_t err; unsigned long m; mpfr_t res; mp_prec_t prec; + int size_n; int inexact; + MPFR_ZIV_DECL (loop); MPFR_SAVE_EXPO_DECL (expo); - MPFR_CLEAR_FLAGS(x); - if (MPFR_UNLIKELY (n == 0)) /* y^0 = 1 for any y */ return mpfr_set_ui (x, 1, rnd); @@ -42,18 +43,17 @@ mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int y, unsigned long int n, /* 0^n = 0 for any n > 0 */ return mpfr_set_ui (x, 0, rnd); + for (size_n = 0, m = n; m; size_n++, m >>= 1); + MPFR_SAVE_EXPO_MARK (expo); - prec = MPFR_PREC (x); - mpfr_init2 (res, 2*prec); + prec = MPFR_PREC (x) + 3 + size_n; + mpfr_init2 (res, prec); - do + MPFR_ZIV_INIT (loop, prec); + for (;;) { - int i; + int i = size_n; - prec += 3; - for (i = 0, m = n; m; i++, m >>= 1) - prec++; - mpfr_set_prec (res, prec); inexact = mpfr_set_ui (res, y, GMP_RNDU); err = 1; /* now 2^(i-1) <= n < 2^i: i=1+floor(log2(n)) */ @@ -68,9 +68,17 @@ mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int y, unsigned long int n, we have err = 1+floor(log2(n)). Since prec >= MPFR_PREC(x) + 4 + floor(log2(n)), prec > err */ err = prec - err; + + if (MPFR_LIKELY (inexact == 0 + || mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ, + MPFR_PREC(x) + (rnd == GMP_RNDN)))) + break; + + /* Actualisation of the precision */ + MPFR_ZIV_NEXT (loop, prec); + mpfr_set_prec (res, prec); } - while (inexact && !mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ, - MPFR_PREC(x) + (rnd == GMP_RNDN))); + MPFR_ZIV_FREE (loop); inexact = mpfr_set (x, res, rnd); @@ -1,6 +1,6 @@ /* mpfr_zeta -- compute the Riemann Zeta function -Copyright 2003, 2004 Free Software Foundation. +Copyright 2003, 2004, 2005 Free Software Foundation. Contributed by Jean-Luc Re'my and the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -22,9 +22,7 @@ MA 02111-1307, USA. */ /* #define DEBUG */ -#include <stdlib.h> -#include <stdio.h> - +#define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" /* @@ -76,11 +74,7 @@ mpfr_zeta_part_b (mpfr_t b, mpfr_srcptr s, int n, int p, mpfr_t *tc) mpfr_neg (s1, s1, GMP_RNDN); mpfr_ui_pow (u, n, s1, GMP_RNDN); mpfr_mul (b, d, u, GMP_RNDN); -#ifdef DEBUG - printf ("\npart b="); - mpfr_out_str (stdout, 10, 0, b, GMP_RNDN); - printf ("\n"); -#endif + MPFR_TRACE (MPFR_DUMP (b)); mpfr_clear (s1); mpfr_clear (d); mpfr_clear (u); @@ -126,7 +120,7 @@ mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n) int i, preca; mpfr_t u, s1; - preca = mpfr_get_prec (sum); + preca = MPFR_PREC (sum); mpfr_init2 (u, preca); mpfr_init2 (s1, preca); mpfr_neg (s1, s, GMP_RNDN); @@ -139,11 +133,7 @@ mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n) mpfr_add (sum, sum, u, GMP_RNDN); } mpfr_add_ui (sum, sum, 1, GMP_RNDN); -#ifdef DEBUG - printf ("\npart a = "); - mpfr_out_str (stdout, 10, 0, sum, GMP_RNDN); - printf ("\n"); -#endif + MPFR_TRACE (MPFR_DUMP (sum)); mpfr_clear (s1); mpfr_clear (u); } @@ -156,15 +146,22 @@ mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n) static int mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) { - int p, n, l, add, can_round; + int p, n, l, add; double beta, sd, dnep; mpfr_t a, b, c, z_pre, f, g, s1; mpfr_t *tc1; mp_prec_t precz, precs, d, dint; int inex; + MPFR_ZIV_DECL (loop); + + MPFR_TRACE (MPFR_DUMP (s)); - precz = mpfr_get_prec (z); - precs = mpfr_get_prec (s); + precz = MPFR_PREC (z); + precs = MPFR_PREC (s); + + d = precz + 11; + /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */ + mpfr_init2 (s1, MAX (precs, ((mpfr_uexp_t) MPFR_EXP (s)))); mpfr_init2 (a, MPFR_PREC_MIN); mpfr_init2 (b, MPFR_PREC_MIN); @@ -172,22 +169,15 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) mpfr_init2 (z_pre, MPFR_PREC_MIN); mpfr_init2 (f, MPFR_PREC_MIN); mpfr_init2 (g, MPFR_PREC_MIN); -#ifdef DEBUG - printf ("Warning: mpfr_zeta assumes that s and Zeta(s) are not both representable in mpfr\n"); - printf ("s="); - mpfr_print_binary (s); - printf ("\n"); -#endif - d = precz + 11; - /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */ - mpfr_init2 (s1, (precs > (mp_exp_unsigned_t) MPFR_EXP(s)) ? - precs : (mp_exp_unsigned_t) MPFR_EXP(s)); - do + + MPFR_ZIV_INIT (loop, d); + for (;;) { /* Principal loop: we compute, in z_pre, an approximation of Zeta(s), that we send to mpfr_can_round */ mpfr_sub_ui (s1, s, 1, GMP_RNDN); MPFR_ASSERTN (MPFR_IS_FP (s1)); + if (MPFR_IS_ZERO (s1)) { mpfr_set_inf (z, 1); @@ -199,11 +189,8 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) uses the approximation Zeta(s)=1/(s-1)+gamma, where gamma is Euler's constant */ { - dint = MAX(d + 3, precs); -#ifdef DEBUG - printf ("branch 1\n"); - printf ("internal precision=%d\n", dint); -#endif + dint = MAX (d + 3, precs); + MPFR_TRACE (printf ("branch 1\ninternal precision=%d\n", dint)); mpfr_set_prec (z_pre, dint); mpfr_set_prec (g, dint); mpfr_ui_div (z_pre, 1, s1, GMP_RNDN); @@ -213,9 +200,8 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) else /* Branch 2 */ { size_t size; -#ifdef DEBUG - printf ("branch 2\n"); -#endif + + MPFR_TRACE (printf ("branch 2\n")); /* Computation of parameters n, p and working precision */ dnep = (double) d * LOG2; sd = mpfr_get_d (s, GMP_RNDN); @@ -235,20 +221,18 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) p = 1 + (int) beta / 2; n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832); } -#ifdef DEBUG - printf ("\nn=%d\np=%d\n",n,p); -#endif + MPFR_TRACE (printf ("\nn=%d\np=%d\n",n,p)); /* add = 4 + floor(1.5 * log(d) / log (2)). We should have add >= 10, which is always fulfilled since d = precz + 11 >= 12, thus ceil(log2(d)) >= 4 */ - add = 4 + (3 * __gmpfr_ceil_log2 ((double) d)) / 2; + add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2; MPFR_ASSERTD(add >= 10); dint = d + add; if (dint < precs) dint = precs; -#ifdef DEBUG - printf("internal precision=%d\n",dint); -#endif + + MPFR_TRACE (printf("internal precision=%d\n",dint)); + size = (p + 1) * sizeof(mpfr_t); tc1 = (mpfr_t*) (*__gmp_allocate_func) (size); for (l=1; l<=p; l++) @@ -258,9 +242,9 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) mpfr_set_prec (c, dint); mpfr_set_prec (z_pre, dint); mpfr_set_prec (f, dint); -#ifdef DEBUG - printf ("precision of z =%d\n", precz); -#endif + + MPFR_TRACE (printf ("precision of z =%d\n", precz)); + /* Computation of the coefficients c_k */ mpfr_zeta_c (p, tc1); /* Computation of the 3 parts of the fonction Zeta. */ @@ -271,11 +255,7 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) mpfr_ui_div (c, 1, s1, GMP_RNDN); mpfr_ui_pow (f, n, s1, GMP_RNDN); mpfr_div (c, c, f, GMP_RNDN); -#ifdef DEBUG - printf ("\npart c="); - mpfr_out_str (stdout, 10, 0, c, GMP_RNDN); - printf ("\n"); -#endif + MPFR_TRACE (MPFR_DUMP (c)); mpfr_add (z_pre, a, c, GMP_RNDN); mpfr_add (z_pre, z_pre, b, GMP_RNDN); for (l=1; l<=p; l++) @@ -283,34 +263,18 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) (*__gmp_free_func) (tc1, size); /* End branch 2 */ } -#ifdef DEBUG - printf ("\nZeta(s) before rounding="); - mpfr_print_binary (z_pre); -#endif - can_round = mpfr_can_round (z_pre, d - 3, GMP_RNDN, GMP_RNDZ, - precz + (rnd_mode == GMP_RNDN)); - if (can_round) - { -#ifdef DEBUG - printf ("\nwe can round"); -#endif - } - else - { -#ifdef DEBUG - printf ("\nwe cannot round"); -#endif - d = d + 5; - } + + MPFR_TRACE (MPFR_DUMP (z_pre)); + if (MPFR_LIKELY (mpfr_can_round (z_pre, d - 3, GMP_RNDN, GMP_RNDZ, + precz + (rnd_mode == GMP_RNDN)))) + break; + MPFR_ZIV_NEXT (loop, d); } - while (can_round == 0); + MPFR_ZIV_FREE (loop); inex = mpfr_set (z, z_pre, rnd_mode); -#ifdef DEBUG - printf ("\nZeta(s) after rounding="); - mpfr_print_binary (z); - printf ("\n"); -#endif + MPFR_TRACE (MPFR_DUMP (z)); + clear_and_return: mpfr_clear (a); mpfr_clear (b); @@ -327,81 +291,86 @@ int mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) { double sd, eps, m1, c; - int can_round; long add; mpfr_t z_pre, s1, s2, y, p; mp_prec_t precz, prec1, precs, precs1; int inex; + MPFR_SAVE_EXPO_DECL (expo); /* Zero, Nan or Inf ? */ - if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(s) )) + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s))) { - if (MPFR_IS_NAN(s)) + if (MPFR_IS_NAN (s)) { MPFR_SET_NAN (z); MPFR_RET_NAN; } - else if (MPFR_IS_INF(s)) + else if (MPFR_IS_INF (s)) { - if (MPFR_SIGN(s) > 0) + if (MPFR_IS_POS (s)) return mpfr_set_ui (z, 1, GMP_RNDN); /* Zeta(+Inf) = 1 */ MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */ MPFR_RET_NAN; } else /* s iz zero */ { - MPFR_ASSERTD(MPFR_IS_ZERO(s)); + MPFR_ASSERTD (MPFR_IS_ZERO (s)); mpfr_set_ui (z, 1, rnd_mode); mpfr_div_2ui (z, z, 1, rnd_mode); - MPFR_CHANGE_SIGN(z); - MPFR_RET(0); + MPFR_CHANGE_SIGN (z); + MPFR_RET (0); } } MPFR_CLEAR_FLAGS(z); /* s is neither Nan, nor Inf, nor Zero */ - mpfr_init2 (s2, mpfr_get_prec (s)); + mpfr_init2 (s2, MPFR_PREC (s)); mpfr_div_2ui (s2, s, 1, rnd_mode); - if (MPFR_IS_NEG(s) && mpfr_floor(s2, s2) == 0) /* Case s = -2n */ + if (MPFR_IS_NEG (s) && mpfr_floor (s2, s2) == 0) /* Case s = -2n */ { mpfr_clear (s2); return mpfr_set_ui (z, 0, rnd_mode); } mpfr_clear (s2); - precz = mpfr_get_prec (z); - precs = mpfr_get_prec (s); - /* Precision precs1 needed to represent 1 - s, and s + 2, - without any truncation */ - precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s)); - sd = mpfr_get_d (s, GMP_RNDN) - 1.0; - if (sd < 0.0) - sd = -sd; /* now sd = abs(s-1.0) */ - /* Precision prec1 is the precision on elementary computations; it ensures - a final precision prec1 - add for zeta(s) */ - /* eps = pow (2.0, - (double) precz - 14.0); */ - eps = __gmpfr_ceil_exp2 (- (double) precz - 14.0); - m1 = 1.0 + MAX(1.0 / eps, 2.0 * sd) * (1.0 + eps); - c = (1.0 + eps) * (1.0 + eps * MAX(8.0, m1)); - /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */ - add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1)); - prec1 = precz + add; /* Note that prec1 is still incremented by 10 at the first entry of loop below */ - prec1 = MAX(prec1, precs1); + + MPFR_SAVE_EXPO_MARK (expo); + + /* Compute Zeta */ if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */ inex = mpfr_zeta_pos (z, s, rnd_mode); else /* use reflection formula zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */ { - mpfr_init2 (z_pre, MPFR_PREC_MIN); - mpfr_init2 (s1, MPFR_PREC_MIN); - mpfr_init2 (y, MPFR_PREC_MIN); - mpfr_init2 (p, MPFR_PREC_MIN); - do + MPFR_ZIV_DECL (loop); + + precz = MPFR_PREC (z); + precs = MPFR_PREC (s); + + /* Precision precs1 needed to represent 1 - s, and s + 2, + without any truncation */ + precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s)); + sd = mpfr_get_d (s, GMP_RNDN) - 1.0; + if (sd < 0.0) + sd = -sd; /* now sd = abs(s-1.0) */ + /* Precision prec1 is the precision on elementary computations; + it ensures a final precision prec1 - add for zeta(s) */ + /* eps = pow (2.0, - (double) precz - 14.0); */ + eps = __gmpfr_ceil_exp2 (- (double) precz - 14.0); + m1 = 1.0 + MAX(1.0 / eps, 2.0 * sd) * (1.0 + eps); + c = (1.0 + eps) * (1.0 + eps * MAX(8.0, m1)); + /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */ + add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1)); + prec1 = precz + add; + prec1 = MAX (prec1, precs1) + 10; + + mpfr_init2 (z_pre, prec1); + mpfr_init2 (s1, prec1); + mpfr_init2 (y, prec1); + mpfr_init2 (p, prec1); + + MPFR_ZIV_INIT (loop, prec1); + for (;;) { - prec1 = prec1 + 10; - mpfr_set_prec (z_pre, prec1); - mpfr_set_prec (s1, prec1); - mpfr_set_prec (y, prec1); - mpfr_set_prec (p, prec1); mpfr_ui_sub (s1, 1, s, GMP_RNDN); /* s1 = 1-s */ mpfr_zeta_pos (z_pre, s1, GMP_RNDN); /* zeta(1-s) */ mpfr_gamma (y, s1, GMP_RNDN); /* gamma(1-s) */ @@ -416,15 +385,29 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) mpfr_pow (y, y, s1, GMP_RNDN); /* (2*Pi)^(s-1) */ mpfr_mul (z_pre, z_pre, y, GMP_RNDN); mpfr_mul_2ui (z_pre, z_pre, 1, GMP_RNDN); - can_round = mpfr_can_round (z_pre, prec1 - add, GMP_RNDN, GMP_RNDZ, - precz + (rnd_mode == GMP_RNDN)); + + if (MPFR_LIKELY (mpfr_can_round (z_pre, prec1 - add, GMP_RNDN, + GMP_RNDZ, + precz + (rnd_mode == GMP_RNDN)))) + break; + + /* Actualisation of the precision */ + MPFR_ZIV_NEXT (loop, prec1); + mpfr_set_prec (z_pre, prec1); + mpfr_set_prec (s1, prec1); + mpfr_set_prec (y, prec1); + mpfr_set_prec (p, prec1); } - while (can_round == 0); + MPFR_ZIV_FREE (loop); + inex = mpfr_set (z, z_pre, rnd_mode); + mpfr_clear(z_pre); mpfr_clear(s1); mpfr_clear(y); mpfr_clear(p); } - return inex; + + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_check_range (z, inex, rnd_mode); } |