diff options
| -rw-r--r-- | acos.c | 3 | ||||
| -rw-r--r-- | acosh.c | 3 | ||||
| -rw-r--r-- | agm.c | 7 | ||||
| -rw-r--r-- | asin.c | 3 | ||||
| -rw-r--r-- | asinh.c | 3 | ||||
| -rw-r--r-- | atan.c | 6 | ||||
| -rw-r--r-- | atanh.c | 7 | ||||
| -rw-r--r-- | const_euler.c | 8 | ||||
| -rw-r--r-- | const_log2.c | 26 | ||||
| -rw-r--r-- | const_pi.c | 35 | ||||
| -rw-r--r-- | cos.c | 2 | ||||
| -rw-r--r-- | cosh.c | 56 | ||||
| -rw-r--r-- | erf.c | 13 | ||||
| -rw-r--r-- | exp2.c | 29 | ||||
| -rw-r--r-- | expm1.c | 9 | ||||
| -rw-r--r-- | factorial.c | 12 | ||||
| -rw-r--r-- | gamma.c | 10 | ||||
| -rw-r--r-- | hypot.c | 3 | ||||
| -rw-r--r-- | log.c | 13 | ||||
| -rw-r--r-- | log10.c | 3 | ||||
| -rw-r--r-- | log1p.c | 38 | ||||
| -rw-r--r-- | log2.c | 42 | ||||
| -rw-r--r-- | pow.c | 3 | ||||
| -rw-r--r-- | pow_si.c | 21 | ||||
| -rw-r--r-- | pow_ui.c | 25 | ||||
| -rw-r--r-- | sin.c | 3 | ||||
| -rw-r--r-- | sinh.c | 52 | ||||
| -rw-r--r-- | sqrt.c | 2 | ||||
| -rw-r--r-- | tan.c | 3 | ||||
| -rw-r--r-- | tanh.c | 59 | ||||
| -rw-r--r-- | ui_pow_ui.c | 4 | ||||
| -rw-r--r-- | zeta.c | 23 |
32 files changed, 298 insertions, 228 deletions
@@ -108,7 +108,8 @@ mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_div_2ui (tmp, tmp, 1, GMP_RNDN); mpfr_sub (arcc, tmp, arcc, GMP_RNDN); - if (mpfr_can_round (arcc, realprec, GMP_RNDN, rnd_mode, MPFR_PREC(acos))) + if (mpfr_can_round (arcc, realprec, GMP_RNDN, GMP_RNDZ, + MPFR_PREC(acos) + (rnd_mode == GMP_RNDN))) { inexact = mpfr_set (acos, arcc, rnd_mode); good = 1; @@ -108,7 +108,8 @@ mpfr_acosh (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode) /* actualisation of the precision */ Nt += 10; } - while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny)); + while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN))); inexact = mpfr_set (y, t, rnd_mode); @@ -145,13 +145,14 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode) while (mpfr_cmp2(u, v, &eq) != 0 && eq <= p - 2); /* Roundability of the result */ - can_round=mpfr_can_round(v,p-err-3,GMP_RNDN,rnd_mode,q); + can_round = mpfr_can_round (v, p - err - 3, GMP_RNDN, GMP_RNDZ, + q + (rnd_mode == GMP_RNDN)); if (can_round) - go_on=0; + go_on = 0; else { - go_on=1; + go_on = 1; p+=5; s=(p-1)/BITS_PER_MP_LIMB+1; MPFR_INIT(up, u, p, s); @@ -108,7 +108,8 @@ mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_sqrt (tmp, tmp, GMP_RNDN); mpfr_div (tmp, x, tmp, GMP_RNDN); mpfr_atan (arcs, tmp, GMP_RNDN); - if (mpfr_can_round (arcs, realprec, GMP_RNDN, rnd_mode, MPFR_PREC(asin))) + if (mpfr_can_round (arcs, realprec, GMP_RNDN, GMP_RNDZ, + MPFR_PREC(asin) + (rnd_mode == GMP_RNDN))) { inexact = mpfr_set (asin, arcs, rnd_mode); good = 1; @@ -102,7 +102,8 @@ mpfr_asinh (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) /* actualisation of the precision */ Nt += 10; } - while ((err < 0) || (!mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny) + while ((err < 0) || (!mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN)) || MPFR_IS_ZERO(t))); mpfr_restore_emin_emax (); @@ -222,10 +222,10 @@ mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_CHANGE_SIGN(arctgt); } - if (mpfr_can_round (arctgt, realprec, GMP_RNDN, rnd_mode, - MPFR_PREC (arctangent))) + if (mpfr_can_round (arctgt, realprec, GMP_RNDN, GMP_RNDZ, + MPFR_PREC (arctangent) + (rnd_mode == GMP_RNDN))) { - inexact = mpfr_set(arctangent, arctgt, rnd_mode); + inexact = mpfr_set (arctangent, arctgt, rnd_mode); good = 1; realprec += 1; } @@ -122,13 +122,14 @@ mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) Nt += 10; } - while ((err < 0) || - (!mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny) || MPFR_IS_ZERO(t))); + while ((err < 0) || (!mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN)) + || MPFR_IS_ZERO(t))); if(flag_neg) MPFR_CHANGE_SIGN(t); - inexact = mpfr_set(y,t,rnd_mode); + inexact = mpfr_set (y, t, rnd_mode); mpfr_clear(t); mpfr_clear(ti); diff --git a/const_euler.c b/const_euler.c index 0334f328b..97930751c 100644 --- a/const_euler.c +++ b/const_euler.c @@ -36,6 +36,7 @@ mpfr_const_euler (mpfr_t x, mp_rnd_t rnd) mp_prec_t prec = MPFR_PREC(x), m, log2m; mpfr_t y, z; unsigned long n; + int inexact; log2m = __gmpfr_ceil_log2 ((double) prec); m = prec + log2m; @@ -60,14 +61,15 @@ mpfr_const_euler (mpfr_t x, mp_rnd_t rnd) mpfr_const_euler_R (z, n); mpfr_sub (y, y, z, GMP_RNDN); } - while (!mpfr_can_round (y, m - 3, GMP_RNDN, rnd, prec)); + while (!mpfr_can_round (y, m - 3, GMP_RNDN, GMP_RNDZ, + prec + (rnd == GMP_RNDN))); - mpfr_set (x, y, rnd); + inexact = mpfr_set (x, y, rnd); mpfr_clear (y); mpfr_clear (z); - return 1; /* always inexact */ + return inexact; /* always inexact */ } /* computes S(n) = sum(n^k*(-1)^(k-1)/k!/k, k=1..ceil(4.319136566 * n)) diff --git a/const_log2.c b/const_log2.c index 61155b003..9e5761203 100644 --- a/const_log2.c +++ b/const_log2.c @@ -63,6 +63,7 @@ mpfr_const_aux_log2 (mpfr_ptr mylog, mp_rnd_t rnd_mode) int logn; mp_prec_t prec_i_want = MPFR_PREC(mylog); mp_prec_t prec_x; + int inexact; mpz_init(cst); logn = __gmpfr_ceil_log2 ((double) MPFR_PREC(mylog)); @@ -93,17 +94,20 @@ mpfr_const_aux_log2 (mpfr_ptr mylog, mp_rnd_t rnd_mode) mpfr_clear(tmp1); mpfr_clear(tmp2); mpfr_clear(tmp3); - if (mpfr_can_round(result, prec_x, GMP_RNDD, rnd_mode, prec_i_want)){ - mpfr_set(mylog, result, rnd_mode); - good = 1; - } else + if (mpfr_can_round (result, prec_x, GMP_RNDD, GMP_RNDZ, + prec_i_want + (rnd_mode == GMP_RNDN))) + { + inexact = mpfr_set (mylog, result, rnd_mode); + good = 1; + } + else { prec_x += logn; } - mpfr_clear(result); + mpfr_clear (result); } - mpz_clear(cst); - return 0; + mpz_clear (cst); + return inexact; } /* Cross-over point from nai"ve Taylor series to binary splitting, @@ -129,6 +133,7 @@ mpfr_const_log2 (mpfr_ptr x, mp_rnd_t rnd_mode) { mp_prec_t N, k, precx; mpz_t s, t, u; + int inexact; precx = MPFR_PREC(x); MPFR_CLEAR_FLAGS(x); @@ -138,7 +143,7 @@ mpfr_const_log2 (mpfr_ptr x, mp_rnd_t rnd_mode) { if ((rnd_mode == __mpfr_const_log2_rnd) || mpfr_can_round (__mpfr_const_log2, __gmpfr_const_log2_prec - 1, - __mpfr_const_log2_rnd, rnd_mode, precx)) + __mpfr_const_log2_rnd, GMP_RNDZ, precx + (rnd_mode == GMP_RNDN))) { return mpfr_set (x, __mpfr_const_log2, rnd_mode); } @@ -172,7 +177,7 @@ mpfr_const_log2 (mpfr_ptr x, mp_rnd_t rnd_mode) mpz_clear (u); } else /* use binary splitting method */ - mpfr_const_aux_log2(x, rnd_mode); + inexact = mpfr_const_aux_log2 (x, rnd_mode); /* store computed value */ if (__gmpfr_const_log2_prec == 0) @@ -183,5 +188,6 @@ mpfr_const_log2 (mpfr_ptr x, mp_rnd_t rnd_mode) mpfr_set (__mpfr_const_log2, x, rnd_mode); __gmpfr_const_log2_prec = precx; __mpfr_const_log2_rnd = rnd_mode; - return 1; + + return inexact; } diff --git a/const_pi.c b/const_pi.c index 62d5d72c8..e97bc02f6 100644 --- a/const_pi.c +++ b/const_pi.c @@ -51,6 +51,7 @@ mpfr_pi_machin3 (mpfr_ptr mylog, mp_rnd_t rnd_mode) int good = 0; mpfr_t tmp1, tmp2, result,tmp3,tmp4,tmp5,tmp6; mpz_t cst; + int inex; MPFR_CLEAR_FLAGS(mylog); logn = __gmpfr_ceil_log2 ((double) MPFR_PREC(mylog)); @@ -105,20 +106,20 @@ mpfr_pi_machin3 (mpfr_ptr mylog, mp_rnd_t rnd_mode) mpfr_clear(tmp4); mpfr_clear(tmp5); mpfr_clear(tmp6); - if (mpfr_can_round(result, prec_x - 5, GMP_RNDD, rnd_mode, prec_i_want)) + if (mpfr_can_round (result, prec_x - 5, GMP_RNDD, GMP_RNDZ, + prec_i_want + (rnd_mode == GMP_RNDN))) { - mpfr_set(mylog, result, rnd_mode); - mpfr_clear(result); + inex = mpfr_set (mylog, result, rnd_mode); good = 1; } else { - mpfr_clear(result); prec_x += logn; } + mpfr_clear (result); } - mpz_clear(cst); - return 0; + mpz_clear (cst); + return inex; } /* @@ -161,6 +162,7 @@ mpfr_const_pi (mpfr_ptr x, mp_rnd_t rnd_mode) int N, oldN, n, prec; mpz_t pi, num, den, d3, d2, tmp; mpfr_t y; + int inex; prec=MPFR_PREC(x); @@ -168,9 +170,10 @@ mpfr_const_pi (mpfr_ptr x, mp_rnd_t rnd_mode) if ((prec==__gmpfr_const_pi_prec && rnd_mode==__mpfr_const_pi_rnd) || (prec<=__gmpfr_const_pi_prec && mpfr_can_round(__mpfr_const_pi, __gmpfr_const_pi_prec, - __mpfr_const_pi_rnd, rnd_mode, prec))) + __mpfr_const_pi_rnd, GMP_RNDZ, + prec + (rnd_mode == GMP_RNDN)))) { - return mpfr_set(x, __mpfr_const_pi, rnd_mode); + return mpfr_set (x, __mpfr_const_pi, rnd_mode); } if (prec < 20000) @@ -215,11 +218,11 @@ mpfr_const_pi (mpfr_ptr x, mp_rnd_t rnd_mode) mpz_fdiv_q(tmp, tmp, den); mpz_add(pi, pi, tmp); } - mpfr_set_z(x, pi, rnd_mode); - mpfr_init2(y, mpfr_get_prec(x)); - mpz_add_ui(pi, pi, N+1); - mpfr_set_z(y, pi, rnd_mode); - if (mpfr_cmp(x, y) != 0) + mpfr_set_z (x, pi, rnd_mode); + mpfr_init2 (y, mpfr_get_prec(x)); + mpz_add_ui (pi, pi, N+1); + mpfr_set_z (y, pi, rnd_mode); + if (mpfr_cmp (x, y) != 0) { fprintf(stderr, "does not converge\n"); exit(1); @@ -234,7 +237,8 @@ mpfr_const_pi (mpfr_ptr x, mp_rnd_t rnd_mode) mpfr_clear(y); } else - mpfr_pi_machin3(x, rnd_mode); + inex = mpfr_pi_machin3 (x, rnd_mode); + /* store computed value */ if (__gmpfr_const_pi_prec==0) mpfr_init2(__mpfr_const_pi, prec); @@ -243,5 +247,6 @@ mpfr_const_pi (mpfr_ptr x, mp_rnd_t rnd_mode) mpfr_set(__mpfr_const_pi, x, rnd_mode); __gmpfr_const_pi_prec=prec; __mpfr_const_pi_rnd=rnd_mode; - return 1; + + return inex; } @@ -78,7 +78,7 @@ mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) /* now the error is bounded by 2^(k-m) = 2^(EXP(s)-err) */ l = mpfr_can_round (s, MPFR_GET_EXP (s) + m - k, - GMP_RNDN, rnd_mode, precy); + GMP_RNDN, GMP_RNDZ, precy + (rnd_mode == GMP_RNDN)); if (l == 0) { @@ -75,47 +75,51 @@ mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) long int 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+3+__gmpfr_ceil_log2(Nt); + Nt = Nt + 3 + __gmpfr_ceil_log2 (Nt); /* initialise of intermediary variable */ - mpfr_init(t); - mpfr_init(te); - mpfr_init(ti); + mpfr_init (t); + mpfr_init (te); + mpfr_init (ti); /* First computation of cosh */ - do { + do + { - /* reactualisation of the precision */ + /* reactualisation of the precision */ - mpfr_set_prec(t,Nt); - mpfr_set_prec(te,Nt); - mpfr_set_prec(ti,Nt); + mpfr_set_prec (t, Nt); + mpfr_set_prec (te, Nt); + mpfr_set_prec (ti, Nt); - /* compute cosh */ - mpfr_exp(te,x,GMP_RNDD); /* exp(x) */ - mpfr_ui_div(ti,1,te,GMP_RNDU); /* 1/exp(x) */ - mpfr_add(t,te,ti,GMP_RNDN); /* exp(x) + 1/exp(x)*/ - mpfr_div_2ui(t,t,1,GMP_RNDN); /* 1/2(exp(x) + 1/exp(x))*/ + /* compute cosh */ + mpfr_exp (te, x, GMP_RNDD); /* exp(x) */ + mpfr_ui_div (ti, 1, te, GMP_RNDU); /* 1/exp(x) */ + mpfr_add (t, te, ti, GMP_RNDN); /* exp(x) + 1/exp(x)*/ + mpfr_div_2ui (t, t, 1, GMP_RNDN); /* 1/2(exp(x) + 1/exp(x))*/ - /* estimation of the error */ - err=Nt-3; + /* estimation of the error */ + err = Nt - 3; - /* actualisation of the precision */ - Nt += 10; + /* actualisation of the precision */ + Nt += 10; - } while ((err <0) || !mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny)); + } + while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN))); - inexact = mpfr_set(y,t,rnd_mode); + inexact = mpfr_set (y, t, rnd_mode); - mpfr_clear(t); - mpfr_clear(ti); - mpfr_clear(te); + mpfr_clear (t); + mpfr_clear (ti); + mpfr_clear (te); } - mpfr_clear(x); - return inexact; + mpfr_clear (x); + + return inexact; } @@ -42,6 +42,7 @@ mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) int sign_x; mp_rnd_t rnd2; double n = (double) MPFR_PREC(y); + int inex; if (MPFR_IS_NAN(x)) { @@ -88,13 +89,13 @@ mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) } else /* use Taylor */ { - mpfr_erf_0 (y, x, rnd2); + inex = mpfr_erf_0 (y, x, rnd2); } if (sign_x < 0) MPFR_CHANGE_SIGN(y); - return 1; + return inex; } /* return x*2^e */ @@ -132,6 +133,7 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd_mode) unsigned int k; long log2tauk; int ok; + int inex; n = MPFR_PREC(res); /* target precision */ xf = mpfr_get_d (x, GMP_RNDN); @@ -189,7 +191,8 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd_mode) tauk = 4.0 * tauk + 11.0; /* final ulp-error on s */ log2tauk = __gmpfr_ceil_log2 (tauk); - ok = mpfr_can_round (s, m - log2tauk, GMP_RNDN, rnd_mode, n); + ok = mpfr_can_round (s, m - log2tauk, GMP_RNDN, GMP_RNDZ, + n + (rnd_mode == GMP_RNDN)); if (ok == 0) { @@ -199,14 +202,14 @@ mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, mp_rnd_t rnd_mode) } while (ok == 0); - mpfr_set (res, s, rnd_mode); + inex = mpfr_set (res, s, rnd_mode); mpfr_clear (y); mpfr_clear (t); mpfr_clear (u); mpfr_clear (s); - return 1; + return inex; } #if 0 @@ -113,24 +113,27 @@ mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_init (te); /* First computation */ - do { + do + { - /* reactualisation of the precision */ - mpfr_set_prec (t, Nt); - mpfr_set_prec (te, Nt); + /* reactualisation of the precision */ + mpfr_set_prec (t, Nt); + mpfr_set_prec (te, Nt); - /* compute exp(x*ln(2))*/ - mpfr_const_log2 (t, GMP_RNDU); /* ln(2) */ - mpfr_mul (te, x, t, GMP_RNDU); /* x*ln(2) */ - mpfr_exp (t, te, GMP_RNDN); /* exp(x*ln(2))*/ + /* compute exp(x*ln(2))*/ + mpfr_const_log2 (t, GMP_RNDU); /* ln(2) */ + mpfr_mul (te, x, t, GMP_RNDU); /* x*ln(2) */ + mpfr_exp (t, te, GMP_RNDN); /* exp(x*ln(2))*/ - /* estimate of the error -- see pow function in algorithms.ps*/ - err = Nt - (MPFR_GET_EXP (te) + 2); + /* estimate of the error -- see pow function in algorithms.ps*/ + err = Nt - (MPFR_GET_EXP (te) + 2); - /* actualisation of the precision */ - Nt += __gmpfr_isqrt (Nt) + 10; + /* actualisation of the precision */ + Nt += __gmpfr_isqrt (Nt) + 10; - } while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny)); + } + while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN))); inexact = mpfr_set (y, t, rnd_mode); @@ -104,12 +104,13 @@ mpfr_expm1 (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode) /* actualisation of the precision */ Nt += 10; } - while ((err <0) || !mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny)); + while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN))); - inexact = mpfr_set(y,t,rnd_mode); + inexact = mpfr_set (y, t, rnd_mode); - mpfr_clear(t); - mpfr_clear(te); + mpfr_clear (t); + mpfr_clear (te); } return inexact; diff --git a/factorial.c b/factorial.c index d4f1992b9..4ee38414b 100644 --- a/factorial.c +++ b/factorial.c @@ -76,26 +76,26 @@ mpfr_fac_ui (mpfr_ptr y, unsigned long int x , mp_rnd_t rnd_mode) err = Nt - 1 - (int) __gmpfr_ceil_log2 ((double) Nt); - round = !inexact || mpfr_can_round (t,err,GMP_RNDZ,rnd_mode,Ny); + round = !inexact || mpfr_can_round (t, err, GMP_RNDZ, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN)); if (round) { round = mpfr_set (y, t, rnd_mode); if (inexact == 0) inexact = round; - boucle=0; + boucle = 0; } else { - Nt=Nt+10; + Nt = Nt + 10; /*initialise of intermediary variable */ - mpfr_set_prec(t, Nt); + mpfr_set_prec (t, Nt); } } - mpfr_clear(t); + mpfr_clear (t); return inexact; - } } @@ -47,7 +47,6 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_t the_pi; mpfr_t GammaTrial; mpfr_t tmp, tmp2; - mp_prec_t Prec; mp_prec_t prec_gamma; mp_prec_t prec_nec; @@ -58,6 +57,7 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) int compared; int k; int sign; + int inex; /* Trivial cases */ if (MPFR_IS_NAN(x)) @@ -194,10 +194,10 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_out_str (stdout, 10, 0, GammaTrial, GMP_RNDD); printf ("\n"); #endif - if (mpfr_can_round (GammaTrial, realprec, GMP_RNDD, rnd_mode, - MPFR_PREC(gamma))) + if (mpfr_can_round (GammaTrial, realprec, GMP_RNDD, GMP_RNDZ, + MPFR_PREC(gamma) + (rnd_mode == GMP_RNDN))) { - mpfr_set (gamma, GammaTrial, rnd_mode); + inex = mpfr_set (gamma, GammaTrial, rnd_mode); good = 1; } else @@ -216,5 +216,5 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_clear (xp); - return 1; /* inexact result */ + return inex; /* inexact result */ } @@ -151,7 +151,8 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x , mpfr_srcptr y , mp_rnd_t rnd_mode) not_exact = 1; } - while (not_exact && mpfr_can_round (t, Nt - 2, GMP_RNDZ, rnd_mode, Nz) == 0); + while (not_exact && !mpfr_can_round (t, Nt - 2, GMP_RNDN, GMP_RNDZ, + Nz + (rnd_mode == GMP_RNDN))); inexact = mpfr_mul_2ui (z, t, sh, rnd_mode); @@ -159,13 +159,12 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) /* we have 7 ulps of error from the above roundings, 4 ulps from the 4/s^2 second order term, plus the canceled bits */ - if (mpfr_can_round (cst, p - cancel - 4, GMP_RNDN, rnd_mode, q) == 1) { - inexact = mpfr_set (r, cst, rnd_mode); -#ifdef DEBUG - printf("result="); mpfr_print_binary(r); putchar('\n'); -#endif - go_on=0; - } + if (mpfr_can_round (cst, p - cancel - 4, GMP_RNDN, GMP_RNDZ, + q + (rnd_mode == GMP_RNDN)) == 1) + { + inexact = mpfr_set (r, cst, rnd_mode); + go_on = 0; + } /* else we increase the precision */ else { p += BITS_PER_MP_LIMB + cancel; @@ -124,7 +124,8 @@ 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, rnd_mode, Ny); + 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)) @@ -95,35 +95,37 @@ mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) long int 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+5+__gmpfr_ceil_log2(Nt); + Nt = Nt + 5 + __gmpfr_ceil_log2 (Nt); /* initialise of intermediary variable */ - mpfr_init(t); + mpfr_init (t); /* First computation of cosh */ - do { - - /* reactualisation of the precision */ - mpfr_set_prec(t, Nt); + do + { + /* reactualisation of the precision */ + mpfr_set_prec (t, Nt); - /* compute log1p */ - 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); + /* compute log1p */ + mpfr_add_ui (t, x, 1, GMP_RNDN); /* 1+x */ + mpfr_log (t, t, GMP_RNDN); /* log(1+x)*/ - /* actualisation of the precision */ - Nt += 10; + /* 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); - } while ((err<0) || !mpfr_can_round(t, err, GMP_RNDN, rnd_mode, Ny)); + /* actualisation of the precision */ + Nt += 10; + } + while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN))); inexact = mpfr_set (y, t, rnd_mode); - mpfr_clear(t); + mpfr_clear (t); } + return inexact; } @@ -113,29 +113,31 @@ mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) /* First computation of log2 */ - do { - - /* reactualisation of the precision */ - mpfr_set_prec(t,Nt); - mpfr_set_prec(tt,Nt); + do + { + /* 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; - - /* actualisation of the precision */ - Nt += 10; - } while ((err<0) || !mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny)); + /* 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; + + /* actualisation of the precision */ + Nt += 10; + } + while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN))); - inexact = mpfr_set(r,t,rnd_mode); + inexact = mpfr_set (r, t, rnd_mode); - mpfr_clear(t); - mpfr_clear(tt); + mpfr_clear (t); + mpfr_clear (tt); } + return inexact; } @@ -349,7 +349,8 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) /* actualisation of the precision */ Nt += 10; - ok = mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Nz); + ok = mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Nz + (rnd_mode == GMP_RNDN)); /* check exact power */ if (ok == 0 && loop == 1) @@ -47,7 +47,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) MPFR_CLEAR_NAN(y); if (n == 0) - return mpfr_set_ui(y, 1, GMP_RNDN); + return mpfr_set_ui (y, 1, GMP_RNDN); if (MPFR_IS_INF(x)) { @@ -71,6 +71,14 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) MPFR_CLEAR_INF(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) + { + mpfr_set_si (y, (n % 2) ? MPFR_SIGN(x) : 1, rnd_mode); + MPFR_EXP(y) += n * (MPFR_EXP(x) - 1); + return 0; + } + n = -n; /* General case */ @@ -87,15 +95,15 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) int inexact; /* 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+__gmpfr_ceil_log2(Nt); + Nt = Nt + 3 + __gmpfr_ceil_log2 (Nt); mpfr_save_emin_emax (); /* initialise of intermediary variable */ - mpfr_init(t); - mpfr_init(ti); + mpfr_init (t); + mpfr_init (ti); do { @@ -113,7 +121,8 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) /* actualisation of the precision */ Nt += 10; } - while (err < 0 || !mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny)); + while (err < 0 || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN))); inexact = mpfr_set (y, t, rnd_mode); mpfr_clear (t); @@ -44,7 +44,7 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) MPFR_CLEAR_NAN(x); - if (n == 0) /* x^0 = 1 for any x */ + if (n == 0) /* y^0 = 1 for any y */ { /* The return mpfr_set_ui is important as 1 isn't necessarily in the exponent range. */ @@ -64,6 +64,12 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) MPFR_CLEAR_INF(x); + if (MPFR_IS_ZERO(y)) /* 0^n = 0 for any n */ + { + MPFR_SET_ZERO(x); + MPFR_RET(0); + } + mpfr_save_emin_emax (); mpfr_init (res); @@ -88,16 +94,23 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) if (mpfr_mul (res, res, y, rnd1)) inexact = 1; } + + /* check underflow */ + if (MPFR_EXP(res) <= (double) __gmpfr_emin) + { + mpfr_clear (res); + mpfr_restore_emin_emax (); + return mpfr_set_underflow (x, rnd, (n % 2) ? MPFR_SIGN(y) : 1); + } + err = prec - err; if (err < 0) err = 0; } - while (inexact && - mpfr_can_round (res, err, MPFR_SIGN(res) > 0 ? GMP_RNDU : GMP_RNDD, - rnd, MPFR_PREC(x)) == 0); + while (inexact && !mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ, + MPFR_PREC(x) + (rnd == GMP_RNDN))); - if (mpfr_set (x, res, rnd)) - inexact = 1; + inexact = mpfr_set (x, res, rnd); mpfr_clear (res); mpfr_restore_emin_emax (); return mpfr_check_range (x, inexact, rnd); @@ -124,7 +124,8 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) /* the absolute error on c is at most 2^(e-m) = 2^(EXP(c)-err) */ e = MPFR_GET_EXP (c) + m - e; - ok = (e >= 0) && mpfr_can_round (c, e, GMP_RNDN, rnd_mode, precy); + ok = (e >= 0) && mpfr_can_round (c, e, GMP_RNDN, GMP_RNDZ, + precy + (rnd_mode == GMP_RNDN)); if (ok == 0) { @@ -94,39 +94,41 @@ mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode) mpfr_init (ti); /* First computation of sinh */ - do { + do + { - /* reactualisation of the precision */ + /* reactualisation of the precision */ - mpfr_set_prec (t, Nt); - mpfr_set_prec (te, Nt); - mpfr_set_prec (ti, Nt); + mpfr_set_prec (t, Nt); + mpfr_set_prec (te, Nt); + mpfr_set_prec (ti, Nt); - /* compute sinh */ - mpfr_exp (te, x, GMP_RNDD); /* exp(x) */ - mpfr_ui_div (ti, 1, te, GMP_RNDU); /* 1/exp(x) */ - mpfr_sub (t, te, ti, GMP_RNDN); /* exp(x) - 1/exp(x) */ - mpfr_div_2ui (t, t, 1, GMP_RNDN); /* 1/2(exp(x) - 1/exp(x)) */ + /* compute sinh */ + mpfr_exp (te, x, GMP_RNDD); /* exp(x) */ + mpfr_ui_div (ti, 1, te, GMP_RNDU); /* 1/exp(x) */ + mpfr_sub (t, te, ti, GMP_RNDN); /* exp(x) - 1/exp(x) */ + mpfr_div_2ui (t, t, 1, GMP_RNDN); /* 1/2(exp(x) - 1/exp(x)) */ - /* it may be that t is zero (in fact, it can only occur when te=1, - and thus ti=1 too) */ + /* it may be that t is zero (in fact, it can only occur when te=1, + and thus ti=1 too) */ - if (MPFR_IS_ZERO(t)) - err = -1; - else - { - /* calculation of the error */ - d = MPFR_GET_EXP (te) - MPFR_GET_EXP (t) + 2; + if (MPFR_IS_ZERO(t)) + err = -1; + else + { + /* calculation of the error */ + d = MPFR_GET_EXP (te) - MPFR_GET_EXP (t) + 2; - /* error estimate */ - /* err = Nt-(__gmpfr_ceil_log2(1+pow(2,d)));*/ - err = Nt - (MAX(d,0) + 1); - } + /* error estimate */ + /* err = Nt-(__gmpfr_ceil_log2(1+pow(2,d)));*/ + err = Nt - (MAX(d,0) + 1); + } - /* actualisation of the precision */ - Nt += 10; + /* actualisation of the precision */ + Nt += 10; - } while ((err < 0) || !mpfr_can_round(t, err, GMP_RNDN, rnd_mode, Ny)); + } while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN))); if (flag_neg == 1) MPFR_CHANGE_SIGN(t); @@ -170,7 +170,7 @@ mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) printf ("(inexact = %lu)\n", q_limb); #endif - can_round = mpfr_can_round_raw(rp, rrsize, 1, err, + can_round = mpfr_can_round_raw (rp, rrsize, 1, err, GMP_RNDZ, rnd_mode, MPFR_PREC(r)); /* If we used all the limbs of both the dividend and the divisor, @@ -56,7 +56,8 @@ mpfr_tan (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) { mpfr_sin_cos (s, c, x, GMP_RNDN); /* err <= 1/2 ulp on s and c */ mpfr_div (c, s, c, GMP_RNDN); /* err <= 2 ulps */ - ok = mpfr_can_round (c, m - 1, GMP_RNDN, rnd_mode, precy); + ok = mpfr_can_round (c, m - 1, GMP_RNDN, GMP_RNDZ, + precy + (rnd_mode == GMP_RNDN)); if (ok == 0) { @@ -106,45 +106,48 @@ mpfr_tanh (y, xt, rnd_mode) /* First computation of cosh */ - do { + do + { - /* reactualisation of the precision */ - mpfr_set_prec(t,Nt); - mpfr_set_prec(te,Nt); - mpfr_set_prec(ta,Nt); - mpfr_set_prec(tb,Nt); + /* reactualisation of the precision */ + mpfr_set_prec (t, Nt); + mpfr_set_prec (te, Nt); + mpfr_set_prec (ta, Nt); + mpfr_set_prec (tb, Nt); - /* compute tanh */ - mpfr_mul_2ui(te,x,1,GMP_RNDN); /* 2x */ - mpfr_exp(te,te,GMP_RNDN); /* exp(2x) */ - mpfr_add_ui(ta,te,1,GMP_RNDD); /* exp(2x) + 1*/ - mpfr_sub_ui(tb,te,1,GMP_RNDU); /* exp(2x) - 1*/ - mpfr_div(t,tb,ta,GMP_RNDN); /* (exp(2x)-1)/(exp(2x)+1)*/ + /* compute tanh */ + mpfr_mul_2ui (te, x, 1, GMP_RNDN); /* 2x */ + mpfr_exp (te, te, GMP_RNDN); /* exp(2x) */ + mpfr_add_ui (ta, te, 1, GMP_RNDD); /* exp(2x) + 1*/ + mpfr_sub_ui (tb, te, 1, GMP_RNDU); /* exp(2x) - 1*/ + mpfr_div (t, tb, ta, GMP_RNDN); /* (exp(2x)-1)/(exp(2x)+1)*/ - /* calculation of the error*/ - d = MPFR_GET_EXP (te) - MPFR_GET_EXP (t); + /* calculation of the error*/ + d = MPFR_GET_EXP (te) - MPFR_GET_EXP (t); - /* estimation of the error */ - /*err = Nt-(__gmpfr_ceil_log2(7+pow(2,d+1)));*/ - err = Nt-(MAX(d+1,3)+1); + /* estimation of the error */ + /*err = Nt-(__gmpfr_ceil_log2(7+pow(2,d+1)));*/ + err = Nt - (MAX(d + 1, 3) + 1); - /* actualisation of the precision */ - Nt += 10; + /* actualisation of the precision */ + Nt += 10; - } while ((err <0)||!mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny)); + } + while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + Ny + (rnd_mode == GMP_RNDN))); if (flag_neg==1) - MPFR_CHANGE_SIGN(t); + MPFR_CHANGE_SIGN(t); - inexact = mpfr_set(y,t,rnd_mode); - mpfr_clear(t); - mpfr_clear(te); - mpfr_clear(ta); - mpfr_clear(tb); + inexact = mpfr_set (y, t, rnd_mode); + mpfr_clear (t); + mpfr_clear (te); + mpfr_clear (ta); + mpfr_clear (tb); } - mpfr_clear(x); - return inexact; + mpfr_clear (x); + return inexact; } diff --git a/ui_pow_ui.c b/ui_pow_ui.c index b8a512fb4..c8954e927 100644 --- a/ui_pow_ui.c +++ b/ui_pow_ui.c @@ -72,8 +72,8 @@ mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int y, unsigned long int n, if (err < 0) err = 0; } - while (inexact && (mpfr_can_round (res, err, - MPFR_SIGN(res) > 0 ? GMP_RNDU : GMP_RNDD, rnd, MPFR_PREC(x)) == 0)); + while (inexact && !mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ, + MPFR_PREC(x) + (rnd == GMP_RNDN))); if (mpfr_set (x, res, rnd)) inexact = 1; @@ -32,7 +32,6 @@ MA 02111-1307, USA. */ void mpfr_zeta_part_b _PROTO ((mpfr_t, mpfr_srcptr, int, int, mpfr_t *)); void mpfr_zeta_c _PROTO ((int, mpfr_t *)); void mpfr_zeta_part_a _PROTO ((mpfr_t, mpfr_srcptr, int)); -void mpfr_zeta_pos _PROTO ((mpfr_t, mpfr_srcptr, mp_rnd_t)); /* Parameters: @@ -157,7 +156,7 @@ mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n) Assumes s is neither NaN nor Infinite. Output: z - Zeta(s) rounded to the precision of z with direction rnd_mode */ -void +static int mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) { int p, n, l, dint, add, d, can_round; @@ -165,6 +164,7 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) mpfr_t a,b,c,z_pre,f,g,s1; mpfr_t *tc1; mp_prec_t precz, precs; + int inex; precz = mpfr_get_prec (z); precs = mpfr_get_prec (s); @@ -192,6 +192,7 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) if (MPFR_IS_ZERO (s1)) { mpfr_set_inf (z, 1); + inex = 0; goto clear_and_return; } else if (MPFR_GET_EXP (s1) <= -(d-3)/2) @@ -285,7 +286,8 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) printf ("\nZeta(s) before rounding="); mpfr_print_binary (z_pre); #endif - can_round = mpfr_can_round (z_pre, d - 3, GMP_RNDN, rnd_mode, precz); + can_round = mpfr_can_round (z_pre, d - 3, GMP_RNDN, GMP_RNDZ, + precz + (rnd_mode == GMP_RNDN)); if (can_round) { #ifdef DEBUG @@ -302,7 +304,7 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) } while (can_round == 0); - mpfr_set (z, z_pre, rnd_mode); + inex = mpfr_set (z, z_pre, rnd_mode); #ifdef DEBUG printf ("\nZeta(s) after rounding="); mpfr_print_binary (z); @@ -316,6 +318,8 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) mpfr_clear (f); mpfr_clear (g); mpfr_clear (s1); + + return inex; } int @@ -326,6 +330,7 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) long add; mpfr_t z_pre, s1, s2, y, sfrac, sint, p; mp_prec_t precz, prec1, precs, precs1; + int inex; if (mpfr_nan_p (s)) { @@ -377,7 +382,7 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) prec1 = precz + add; /* Note that prec1 is still incremented by 10 at the first entry of loop below */ prec1 = MAX(prec1, precs1); if (MPFR_SIGN (s) != -1 && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */ - mpfr_zeta_pos (z, s, rnd_mode); + inex = mpfr_zeta_pos (z, s, rnd_mode); else { mpfr_init2(z_pre, MPFR_PREC_MIN); @@ -422,11 +427,11 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) mpfr_pow (y, y, s1, GMP_RNDN); 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, - rnd_mode, precz); + can_round = mpfr_can_round (z_pre, prec1 - add, GMP_RNDN, GMP_RNDZ, + precz + (rnd_mode == GMP_RNDN)); } while (can_round == 0); - mpfr_set (z, z_pre, rnd_mode); + inex = mpfr_set (z, z_pre, rnd_mode); mpfr_clear(z_pre); mpfr_clear(s1); mpfr_clear(y); @@ -434,5 +439,5 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) mpfr_clear(sint); mpfr_clear(sfrac); } - return 1; + return inex; } |