diff options
| -rw-r--r-- | agm.c | 8 | ||||
| -rw-r--r-- | exp2.c | 2 | ||||
| -rw-r--r-- | exp3.c | 6 | ||||
| -rw-r--r-- | get_str.c | 6 | ||||
| -rw-r--r-- | log2.c | 6 | ||||
| -rw-r--r-- | mpfr.h | 4 | ||||
| -rw-r--r-- | pi.c | 6 | ||||
| -rw-r--r-- | pow.c | 4 | ||||
| -rw-r--r-- | sin_cos.c | 6 |
9 files changed, 23 insertions, 25 deletions
@@ -25,13 +25,9 @@ MA 02111-1307, USA. */ #include "gmp-impl.h" #include "mpfr.h" - /*Memory gestion */ #define MON_INIT(xp, x, p, s) xp = (mp_ptr) TMP_ALLOC(s*BYTES_PER_MP_LIMB); x -> _mp_prec = p; x -> _mp_d = xp; x -> _mp_size = s; x -> _mp_exp = 0; - - - void #ifdef __STDC__ mpfr_agm(mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode) @@ -111,10 +107,10 @@ mpfr_agm(r, a, b, rnd_mode) eq=0; - err=ceil((3.0/2.0*log((double)p)/log(2.0)+1.0)*exp(-(double)p*log(2.0))+3.0*exp(-2.0*(double)p*uo*log(2.0)/(vo-uo))); + err=ceil((3.0/2.0*log((double)p)/LOG2+1.0)*exp(-(double)p*LOG2)+3.0*exp(-2.0*(double)p*uo*LOG2/(vo-uo))); if(p-err-3<=q) { p=q+err+4; - err=ceil((3.0/2.0*log((double)p)/log(2.0)+1.0)*exp(-(double)p*log(2.0))+3.0*exp(-2.0*(double)p*uo*log(2.0)/(vo-uo))); + err=ceil((3.0/2.0*log((double)p)/LOG2+1.0)*exp(-(double)p*LOG2)+3.0*exp(-2.0*(double)p*uo*LOG2/(vo-uo))); } /* Calculus of un and vn */ @@ -39,8 +39,6 @@ int mpfr_exp2_aux2 (mpz_t, mpfr_srcptr, int, int*); /* #define DEBUG */ -#define LOG2 0.69314718055994528622 /* log(2) rounded to zero on 53 bits */ - /* if k = the number of bits of z > q, divides z by 2^(k-q) and returns k-q. Otherwise do nothing and return 0. */ @@ -191,10 +191,10 @@ mp_rnd_t rnd_mode; ----- k bits -- */ prec_x = (int) ceil(log ((double) (MPFR_PREC(x)) / (double) BITS_PER_MP_LIMB) - /log(2.0)); + /LOG2); logn = (int) ceil(log ((double) prec_x+MPFR_PREC(y)) - /log(2.0)); + /LOG2); if (logn < 2) logn = 2; ttt = MPFR_EXP(x); mpfr_init2(x_copy,MPFR_PREC(x)); @@ -211,7 +211,7 @@ mp_rnd_t rnd_mode; Prec = realprec+shift+2+shift_x; k = (int) ceil(log ((double) (Prec) / (double) BITS_PER_MP_LIMB) - /log(2.0)); + /LOG2); /* Maintenant, il faut extraire : */ mpfr_init2(t, Prec); mpfr_init2(tmp, Prec); @@ -91,12 +91,12 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) /* the output exponent f is such that base^(f-1) <= |op| < base^f i.e. f = 1 + floor(log(|op|)/log(base)) = 1 + floor((log(|m|)+e*log(2))/log(base)) */ - f = 1 + (int) floor((log(d)+(double)e*log(2.0))/log((double)base)); + f = 1 + (int) floor((log(d)+(double)e*LOG2)/log((double)base)); if (n==0) { /* performs exact rounding, i.e. returns y such that for GMP_RNDU for example, we have: x*2^(e-p) <= y*base^(f-n) */ - n = (int) ((double)MPFR_PREC(op)*log(2.0)/log((double)base)); + n = (int) ((double)MPFR_PREC(op)*LOG2/log((double)base)); if (n==0) n=1; } #ifdef DEBUG @@ -105,7 +105,7 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) /* now the first n digits of the mantissa are obtained from rnd(op*base^(n-f)) */ if (pow2) prec = n*pow2; - else prec = (long) ceil((double)n*log((double)base)/log(2.0)); + else prec = (long) ceil((double)n*log((double)base)/LOG2); #ifdef DEBUG printf("prec=%d\n", prec); #endif @@ -69,12 +69,12 @@ mpfr_const_aux_log2(mylog, rnd_mode) mpfr_ptr mylog; mp_rnd_t rnd_mode; mpz_init(cst); logn = (int) ceil(log ((double) MPFR_PREC(mylog)) - /log(2.0)); + /LOG2); prec_x = prec_i_want + logn; while (!good){ prec = (int) ceil(log ((double) (prec_x)) - /log(2.0)); + /LOG2); mpfr_init2(tmp1, prec_x); mpfr_init2(result, prec_x); mpfr_init2(tmp2, prec_x); @@ -150,7 +150,7 @@ mpfr_const_log2(x, rnd_mode) mpfr_ptr x; mp_rnd_t rnd_mode; N=2; do { oldN = N; - N = precx + (int)ceil(log((double)N)/log(2.0)); + N = precx + (int)ceil(log((double)N)/LOG2); } while (N != oldN); mpz_init_set_ui(s,0); mpz_init(u); @@ -36,6 +36,10 @@ MA 02111-1307, USA. */ #define GMP_RNDU 2 #define GMP_RNDD 3 +/* Definition of constants */ + +#define LOG2 0.69314718055994528622 /* log(2) rounded to zero on 53 bits */ + /* Definitions of types and their semantics */ typedef unsigned long int mp_prec_t; /* easy to change if necessary */ @@ -55,13 +55,13 @@ mpfr_pi_machin3(mylog, rnd_mode) mpz_t cst; logn = (int) ceil(log ((double) MPFR_PREC(mylog)) - /log(2.0)); + /LOG2); prec_x = prec_i_want + logn + 5; mpz_init(cst); while (!good){ prec = (int) ceil(log ((double) (prec_x)) - /log(2.0)); + /LOG2); mpfr_init2(tmp1, prec_x); mpfr_init2(tmp2, prec_x); @@ -183,7 +183,7 @@ mpfr_const_pi(x, rnd_mode) N=1; do { oldN = N; - N = (prec+3)/4 + (int)ceil(log((double)N+1.0)/log(2.0)); + N = (prec+3)/4 + (int)ceil(log((double)N+1.0)/LOG2); } while (N != oldN); mpz_init(pi); mpz_init(num); mpz_init(den); mpz_init(d3); mpz_init(d2); mpz_init(tmp); @@ -56,7 +56,7 @@ mpfr_pow_ui (x, y, n, rnd) mpfr_mul(x, x, x, rnd); err = 2.0*err+1.0; if (n & (1<<i)) { mpfr_mul(x, x, y, rnd); err += 1.0; } } - return (int) ceil(log(err)/log(2.0)); + return (int) ceil(log(err)/LOG2); } /* sets x to y^n, and returns ceil(log2(max ulp error)) */ @@ -92,5 +92,5 @@ mpfr_ui_pow_ui (x, y, n, rnd) err = err + 1.0; } } - return (int) ceil(log(err)/log(2.0)); + return (int) ceil(log(err)/LOG2); } @@ -100,7 +100,7 @@ mp_rnd_t rnd_mode; prec_x = (int) ceil(log ((double) (MPFR_PREC(x)) / (double) BITS_PER_MP_LIMB) - /log(2.0)); + /LOG2); ttt = MPFR_EXP(x); mpfr_init2(x_copy,MPFR_PREC(x)); mpfr_set(x_copy,x,GMP_RNDD); @@ -114,7 +114,7 @@ mp_rnd_t rnd_mode; } target_prec = MPFR_PREC(sinus); if (MPFR_PREC(cosinus) > target_prec) target_prec = MPFR_PREC(cosinus); - logn = (int) ceil(log((double)target_prec)/log(2.0)); + logn = (int) ceil(log((double)target_prec)/LOG2); if (logn < 2) logn = 2; factor = logn; realprec = target_prec + logn; @@ -122,7 +122,7 @@ mp_rnd_t rnd_mode; Prec = realprec + 2*shift + 2 + shift_x + factor; k = (int) ceil(log ((double) (Prec) / (double) BITS_PER_MP_LIMB) - /log(2.0)); + /LOG2); /* Maintenant, il faut extraire : */ mpfr_init2(t_sin, Prec); mpfr_init2(t_cos, Prec); |