diff options
| -rw-r--r-- | cbrt.c | 9 | ||||
| -rw-r--r-- | const_euler.c | 13 | ||||
| -rw-r--r-- | const_log2.c | 3 | ||||
| -rw-r--r-- | cos.c | 8 | ||||
| -rw-r--r-- | cosh.c | 19 | ||||
| -rw-r--r-- | div_ui.c | 7 | ||||
| -rw-r--r-- | eq.c | 21 | ||||
| -rw-r--r-- | exp.c | 7 | ||||
| -rw-r--r-- | exp_2.c | 23 | ||||
| -rw-r--r-- | get_ld.c | 3 | ||||
| -rw-r--r-- | hypot.c | 26 | ||||
| -rw-r--r-- | init2.c | 3 | ||||
| -rw-r--r-- | log.c | 27 | ||||
| -rw-r--r-- | log10.c | 17 | ||||
| -rw-r--r-- | log2.c | 27 | ||||
| -rw-r--r-- | reldiff.c | 14 | ||||
| -rw-r--r-- | save_expo.c | 11 | ||||
| -rw-r--r-- | set_f.c | 11 | ||||
| -rw-r--r-- | set_prc_raw.c | 14 | ||||
| -rw-r--r-- | set_prec.c | 5 | ||||
| -rw-r--r-- | set_str.c | 17 | ||||
| -rw-r--r-- | ui_div.c | 5 | ||||
| -rw-r--r-- | zeta.c | 3 |
23 files changed, 125 insertions, 168 deletions
@@ -1,6 +1,6 @@ /* mpfr_cbrt -- cube root function. -Copyright 2002, 2003 Free Software Foundation. +Copyright 2002, 2003, 2004 Free Software Foundation. Contributed by the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -20,7 +20,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include <stdlib.h> #include "gmp.h" #include "gmp-impl.h" @@ -46,7 +45,7 @@ MA 02111-1307, USA. */ */ int -mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) +mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) { mpz_t m; @@ -78,7 +77,7 @@ mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) else MPFR_ASSERTN(0); } - /* Useless due to mpz_init + /* Useless due to mpz_init MPFR_CLEAR_FLAGS(y);*/ mpz_init (m); @@ -95,7 +94,7 @@ mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) n = MPFR_PREC(y); if (rnd_mode == GMP_RNDN) n ++; - + /* we want 3*n-2 <= size_m + 3*sh + r <= 3*n i.e. 3*sh + size_m + r <= 3*n */ sh = (3 * n - size_m - r) / 3; diff --git a/const_euler.c b/const_euler.c index ee1fdb8c6..c62ec78ec 100644 --- a/const_euler.c +++ b/const_euler.c @@ -1,6 +1,6 @@ /* mpfr_const_euler -- Euler's constant -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2004 Free Software Foundation. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include <stdlib.h> #include "gmp.h" #include "gmp-impl.h" @@ -115,7 +114,7 @@ mpfr_const_euler_S (mpfr_t x, unsigned long n) mpz_clear (t); } -/* computes R(n) = exp(-n)/n * sum(k!/(-n)^k, k=0..n-2) +/* computes R(n) = exp(-n)/n * sum(k!/(-n)^k, k=0..n-2) with error at most 4*ulp(x). Assumes n>=2. Since x <= exp(-n)/n <= 1/8, then 4*ulp(x) <= ulp(1). */ @@ -137,7 +136,7 @@ mpfr_const_euler_R (mpfr_t x, unsigned long n) for (k = 1; k <= n; k++) { - mpz_mul_ui (a, a, k); + mpz_mul_ui (a, a, k); mpz_div_ui (a, a, n); /* the error e(k) on a is e(k) <= 1 + k/n*e(k-1) with e(0)=0, i.e. e(k) <= k */ @@ -148,11 +147,7 @@ mpfr_const_euler_R (mpfr_t x, unsigned long n) } /* the error on s is at most 1+2+...+n = n*(n+1)/2 */ mpz_div_ui (s, s, n); /* err <= 1 + (n+1)/2 */ - if (MPFR_PREC(x) < mpz_sizeinbase(s, 2)) - { - fprintf (stderr, "prec(x) is too small in mpfr_const_euler_R\n"); - exit (1); - } + MPFR_ASSERTN (MPFR_PREC(x) >= mpz_sizeinbase(s, 2)); mpfr_set_z (x, s, GMP_RNDD); /* exact */ mpfr_div_2ui (x, x, m, GMP_RNDD); /* now x = 1/n * sum(k!/(-n)^k, k=0..n-2) <= 1/n */ diff --git a/const_log2.c b/const_log2.c index 642173009..df00e9c98 100644 --- a/const_log2.c +++ b/const_log2.c @@ -1,6 +1,6 @@ /* mpfr_const_log2 -- compute natural logarithm of 2 -Copyright 1999, 2001, 2002, 2003 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" @@ -1,6 +1,6 @@ /* mpfr_cos -- cosine of a floating-point number -Copyright 2001, 2002, 2003 Free Software Foundation. +Copyright 2001, 2002, 2003, 2004 Free Software Foundation. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" @@ -27,8 +26,8 @@ MA 02111-1307, USA. */ static int mpfr_cos2_aux _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr)); -int -mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) +int +mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) { int K0, K, precy, m, k, l, inexact; mpfr_t r, s; @@ -138,4 +137,3 @@ mpfr_cos2_aux (mpfr_ptr s, mpfr_srcptr r) return l; } - @@ -1,6 +1,6 @@ /* mpfr_cosh -- hyperbolic cosine -Copyright 2001, 2002 Free Software Foundation. +Copyright 2001, 2002, 2004 Free Software Foundation. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" @@ -31,9 +30,9 @@ MA 02111-1307, USA. */ */ int -mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) +mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) { - + /****** Declaration ******/ mpfr_t x; mp_prec_t Nxt = MPFR_PREC(xt); @@ -41,13 +40,13 @@ mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(xt))) { - if (MPFR_IS_NAN(xt)) + if (MPFR_IS_NAN(xt)) { - MPFR_SET_NAN(y); + MPFR_SET_NAN(y); MPFR_RET_NAN; } else if (MPFR_IS_INF(xt)) - { + { MPFR_SET_INF(y); MPFR_SET_POS(y); MPFR_RET(0); @@ -65,7 +64,7 @@ mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) /* General case */ { /* Declaration of the intermediary variable */ - mpfr_t t, te,ti; + mpfr_t t, te,ti; /* Declaration of the size variable */ mp_prec_t Nx = Nxt; /* Precision of input variable */ @@ -73,7 +72,7 @@ mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) mp_prec_t Nt; /* Precision of the intermediary variable */ long int err; /* Precision of error */ - + /* compute the precision of intermediary variable */ Nt = MAX(Nx, Ny); /* the optimal number of bits : see algorithms.ps */ @@ -110,7 +109,7 @@ mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) } 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); @@ -1,6 +1,6 @@ /* mpfr_div_ui -- divide a floating-point number by a machine integer -Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" @@ -43,8 +42,8 @@ mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) MPFR_SET_NAN(y); MPFR_RET_NAN; } - else if (MPFR_IS_INF(x)) - { + else if (MPFR_IS_INF(x)) + { MPFR_SET_INF(y); MPFR_SET_SAME_SIGN(y, x); MPFR_RET(0); @@ -1,6 +1,6 @@ /* mpfr_eq -- Compare two floats up to a specified bit #. -Copyright 1999, 2001, 2003 Free Software Foundation, Inc. +Copyright 1999, 2001, 2003, 2004 Free Software Foundation, Inc. (Copied from GNU MP, file mpf_eq.) This file is part of the MPFR Library. @@ -20,7 +20,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" @@ -69,9 +68,9 @@ mpfr_eq (mpfr_srcptr u, mpfr_srcptr v, unsigned long int n_bits) { if ((unsigned long) vsize * BITS_PER_MP_LIMB < n_bits) { - k = usize - vsize - 1; - while (k >= 0 && !up[k]) --k; - if (k >= 0) + k = usize - vsize - 1; + while (k >= 0 && !up[k]) --k; + if (k >= 0) return 0; /* surely too different */ } size = vsize; @@ -80,9 +79,9 @@ mpfr_eq (mpfr_srcptr u, mpfr_srcptr v, unsigned long int n_bits) { if ((unsigned long) usize * BITS_PER_MP_LIMB < n_bits) { - k = vsize - usize - 1; - while (k >= 0 && !vp[k]) --k; - if (k >= 0) + k = vsize - usize - 1; + while (k >= 0 && !vp[k]) --k; + if (k >= 0) return 0; /* surely too different */ } size = usize; @@ -113,8 +112,8 @@ mpfr_eq (mpfr_srcptr u, mpfr_srcptr v, unsigned long int n_bits) return mpfr_cmp (u, v) == 0; if (n_bits & (BITS_PER_MP_LIMB - 1)) - return (up[i] >> (BITS_PER_MP_LIMB - (n_bits & (BITS_PER_MP_LIMB - 1))) == - vp[i] >> (BITS_PER_MP_LIMB - (n_bits & (BITS_PER_MP_LIMB - 1)))); + return (up[i] >> (BITS_PER_MP_LIMB - (n_bits & (BITS_PER_MP_LIMB - 1))) == + vp[i] >> (BITS_PER_MP_LIMB - (n_bits & (BITS_PER_MP_LIMB - 1)))); else - return (up[i] == vp[i]); + return (up[i] == vp[i]); } @@ -1,6 +1,6 @@ /* mpfr_exp -- exponential of a floating-point number -Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation. +Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation. Contributed by the Spaces project. This file is part of the MPFR Library. @@ -20,7 +20,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" @@ -32,8 +31,8 @@ MA 02111-1307, USA. */ where x = n*log(2)+(2^K)*r number of operations = O(K+prec(r)/K) */ -int -mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) +int +mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) { int expx, precy; double d; @@ -1,7 +1,7 @@ -/* mpfr_exp_2 -- exponential of a floating-point number +/* mpfr_exp_2 -- exponential of a floating-point number using Brent's algorithms in O(n^(1/2)*M(n)) and O(n^(1/3)*M(n)) -Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -20,7 +20,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" /* for count_leading_zeros */ @@ -105,8 +104,8 @@ mpz_normalize2 (mpz_t rop, mpz_t z, int expz, int target) together with Brent-Kung O(t^(1/2)) algorithm for the evaluation of power series. The resulting complexity is O(n^(1/3)*M(n)). */ -int -mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) +int +mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) { int n, K, precy, q, k, l, err, exps, inexact, error_r; mpfr_t r, s, t; mpz_t ss; @@ -162,7 +161,7 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) printf ("r="); mpfr_print_binary (r); putchar ('\n'); #endif mpfr_sub (r, x, r, GMP_RNDU); - /* possible cancellation here: the error on r is at most + /* possible cancellation here: the error on r is at most 3*2^(EXP(old_r)-EXP(new_r)) */ if (MPFR_SIGN(r) < 0) { /* initial approximation n was too large */ @@ -237,7 +236,7 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) using naive method with O(l) multiplications. Return the number of iterations l. The absolute error on s is less than 3*l*(l+1)*2^(-q). - Version using fixed-point arithmetic with mpz instead + Version using fixed-point arithmetic with mpz instead of mpfr for internal computations. s must have at least qn+1 limbs (qn should be enough, but currently fails since mpz_mul_2exp(s, s, q-1) reallocates qn+1 limbs) @@ -251,7 +250,7 @@ mpfr_exp2_aux (mpz_t s, mpfr_srcptr r, int q, int *exps) TMP_MARK(marker); qn = 1 + (q-1)/BITS_PER_MP_LIMB; - MY_INIT_MPZ(t, 2*qn+1); /* 2*qn+1 is needed since mpz_div_2exp may + MY_INIT_MPZ(t, 2*qn+1); /* 2*qn+1 is needed since mpz_div_2exp may need an extra limb */ MY_INIT_MPZ(rr, qn+1); mpz_set_ui(t, 1); expt=0; @@ -285,7 +284,7 @@ mpfr_exp2_aux (mpz_t s, mpfr_srcptr r, int q, int *exps) /* s <- 1 + r/1! + r^2/2! + ... + r^l/l! while MPFR_EXP(r^l/l!)+MPFR_EXPR(r)>-q using Brent/Kung method with O(sqrt(l)) multiplications. Return l. - Uses m multiplications of full size and 2l/m of decreasing size, + Uses m multiplications of full size and 2l/m of decreasing size, i.e. a total equivalent to about m+l/m full multiplications, i.e. 2*sqrt(l) for m=sqrt(l). Version using mpz. ss must have at least (sizer+1) limbs. @@ -343,9 +342,9 @@ mpfr_exp2_aux2 (mpz_t s, mpfr_srcptr r, int q, int *exps) expR[i] = mpz_normalize2(R[i], R[i], expR[i], 1-ql); } /* the absolute error on R[i]*rr is still 2*i-1 ulps */ - expt = mpz_normalize2(t, R[m-1], expR[m-1], 1-ql); + expt = mpz_normalize2(t, R[m-1], expR[m-1], 1-ql); /* err(t) <= 2*m-1 ulps */ - /* computes t = 1 + r/(l+1) + ... + r^(m-1)*l!/(l+m-1)! + /* computes t = 1 + r/(l+1) + ... + r^(m-1)*l!/(l+m-1)! using Horner's scheme */ for (i=m-2;i>=0;i--) { @@ -368,7 +367,7 @@ mpfr_exp2_aux2 (mpz_t s, mpfr_srcptr r, int q, int *exps) #endif mpz_add(s, s, t); /* no error here */ - /* updates rr, the multiplication of the factors l+i could be done + /* updates rr, the multiplication of the factors l+i could be done using binary splitting too, but it is not sure it would save much */ mpz_mul(t, rr, R[m]); /* err(t) <= err(rr) + 2m-1 */ expr += expR[m]; @@ -1,7 +1,7 @@ /* mpfr_get_ld -- convert a multiple precision floating-point number to a machine long double -Copyright 2002, 2003 Free Software Foundation, Inc. +Copyright 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -20,7 +20,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include <float.h> #include "gmp.h" @@ -1,6 +1,6 @@ /* mpfr_hypot -- Euclidean distance -Copyright 2001, 2002, 2003 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,8 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> -#include <stdlib.h> #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" @@ -32,7 +30,7 @@ MA 02111-1307, USA. */ */ int -mpfr_hypot (mpfr_ptr z, mpfr_srcptr x , mpfr_srcptr y , mp_rnd_t rnd_mode) +mpfr_hypot (mpfr_ptr z, mpfr_srcptr x , mpfr_srcptr y , mp_rnd_t rnd_mode) { int inexact; /* Flag exact computation */ @@ -58,7 +56,7 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x , mpfr_srcptr y , mp_rnd_t rnd_mode) MPFR_RET(0); } else if (MPFR_IS_ZERO(x)) - return mpfr_abs (z, y, rnd_mode); + return mpfr_abs (z, y, rnd_mode); else if (MPFR_IS_ZERO(y)) return mpfr_abs (z, x, rnd_mode); else @@ -84,7 +82,7 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x , mpfr_srcptr y , mp_rnd_t rnd_mode) /* we have x < 2^Ex thus x^2 < 2^(2*Ex), and ulp(x) = 2^(Ex-Nx) thus ulp(x^2) >= 2^(2*Ex-2*Nx). - y does not overlap with the result when + y does not overlap with the result when x^2+y^2 < (|x| + 1/2*ulp(x,Nz))^2 = x^2 + |x|*ulp(x,Nz) + 1/4*ulp(x,Nz)^2, i.e. a sufficient condition is y^2 < |x|*ulp(x,Nz), or 2^(2*Ey) <= 2^(2*Ex-1-Nz), i.e. 2*diff_exp > Nz @@ -109,7 +107,7 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x , mpfr_srcptr y , mp_rnd_t rnd_mode) Nx = MPFR_PREC(x); /* Precision of input variable */ Ny = MPFR_PREC(y); /* Precision of input variable */ - + /* compute the working precision -- see algorithms.ps */ Nt = MAX(MAX(MAX(Nx, Ny), Nz), 8); Nt = Nt - 8 + __gmpfr_ceil_log2 (Nt); @@ -125,13 +123,13 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x , mpfr_srcptr y , mp_rnd_t rnd_mode) do { - Nt += 10; + Nt += 10; not_exact = 0; /* reactualization 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); /* computations of hypot */ mpfr_div_2ui (te, x, sh, GMP_RNDZ); /* exact since Nt >= Nx */ @@ -140,14 +138,14 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x , mpfr_srcptr y , mp_rnd_t rnd_mode) mpfr_div_2ui (ti, y, sh, GMP_RNDZ); /* exact since Nt >= Ny */ if (mpfr_mul (ti, ti, ti, GMP_RNDZ)) /* y^2 */ - not_exact = 1; - + not_exact = 1; + if (mpfr_add (t, te, ti, GMP_RNDZ)) /* x^2+y^2 */ not_exact = 1; if (mpfr_sqrt (t, t, GMP_RNDZ)) /* sqrt(x^2+y^2)*/ not_exact = 1; - + } while (not_exact && !mpfr_can_round (t, Nt - 2, GMP_RNDN, GMP_RNDZ, Nz + (rnd_mode == GMP_RNDN))); @@ -1,6 +1,6 @@ /* mpfr_init2 -- initialize a floating-point number with given precision -Copyright 2001, 2002, 2003 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include <stdlib.h> #include "mpfr-impl.h" @@ -1,6 +1,6 @@ /* mpfr_log -- natural logarithm of a floating-point number -Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation. +Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" @@ -42,7 +41,7 @@ MA 02111-1307, USA. */ /* #define DEBUG */ int -mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) +mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) { int m, go_on, size, cancel, inexact = 0; mp_prec_t p, q; @@ -63,7 +62,7 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) /* check for infinity before zero */ else if (MPFR_IS_INF(a)) { - if (MPFR_IS_NEG(a)) + if (MPFR_IS_NEG(a)) /* log(-Inf) = NaN */ { MPFR_SET_NAN(r); @@ -75,7 +74,7 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) MPFR_SET_POS(r); MPFR_RET(0); } - } + } else if (MPFR_IS_ZERO(a)) { MPFR_SET_INF(r); @@ -85,7 +84,7 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) else MPFR_ASSERTN(0); } - + /* If a is negative, the result is NaN */ if (MPFR_UNLIKELY( MPFR_IS_NEG(a) )) { @@ -103,7 +102,7 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) MPFR_CLEAR_FLAGS(r); q=MPFR_PREC(r); - + ref = mpfr_get_d1 (a) - 1.0; if (ref<0) ref=-ref; @@ -127,11 +126,11 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) /* All the mpfr_t needed have a precision of p */ TMP_MARK(marker); size=(p-1)/BITS_PER_MP_LIMB+1; - MPFR_TMP_INIT(cstp, cst, p, size); + MPFR_TMP_INIT(cstp, cst, p, size); MPFR_TMP_INIT(rapportp, rapport, p, size); MPFR_TMP_INIT(agmp, agm, p, size); - MPFR_TMP_INIT(tmp1p, tmp1, p, size); - MPFR_TMP_INIT(tmp2p, tmp2, p, size); + MPFR_TMP_INIT(tmp1p, tmp1, p, size); + MPFR_TMP_INIT(tmp2p, tmp2, p, size); MPFR_TMP_INIT(sp, s, p, size); MPFR_TMP_INIT(mmp, mm, p, size); @@ -146,8 +145,8 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) mpfr_div (tmp2, cst, tmp1, GMP_RNDN); /* pi/2*AG(1,4/s), err<=5ulps */ mpfr_const_log2 (cst, GMP_RNDN); /* compute log(2), err<=1ulp */ mpfr_mul(tmp1,cst,mm,GMP_RNDN); /* I compute m*log(2), err<=2ulps */ - cancel = MPFR_GET_EXP (tmp2); - mpfr_sub(cst,tmp2,tmp1,GMP_RNDN); /* log(a), err<=7ulps+cancel */ + cancel = MPFR_GET_EXP (tmp2); + mpfr_sub(cst,tmp2,tmp1,GMP_RNDN); /* log(a), err<=7ulps+cancel */ cancel -= MPFR_GET_EXP (cst); #ifdef DEBUG printf("canceled bits=%d\n", cancel); @@ -173,9 +172,7 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) /* We clean */ TMP_FREE(marker); - + } return inexact; /* result is inexact */ } - - @@ -1,6 +1,6 @@ /* mpfr_log10 -- logarithm in base 10. -Copyright 2001, 2002, 2003 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" @@ -69,7 +68,7 @@ mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) else MPFR_ASSERTN(0); } - + /* If a is negative, the result is NaN */ if (MPFR_UNLIKELY( MPFR_IS_NEG(a) )) { @@ -97,10 +96,10 @@ mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) /* 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 output variable */ - + mp_prec_t Nt; /* Precision of the intermediary variable */ long int err; /* Precision of error */ - + /* compute the precision of intermediary variable */ Nt = MAX(Nx, Ny); /* the optimal number of bits : see algorithms.ps */ @@ -110,14 +109,14 @@ mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) mpfr_init (t); mpfr_init (tt); - + /* First computation of log10 */ do { /* reactualisation of the precision */ mpfr_set_prec (t, Nt); - mpfr_set_prec (tt, Nt); - + mpfr_set_prec (tt, Nt); + /* compute log10 */ mpfr_set_ui (t, 10, GMP_RNDN); /* 10 */ mpfr_log (t, t, GMP_RNDD); /* log(10) */ @@ -140,7 +139,7 @@ mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) /* actualisation of the precision */ Nt += 10; } while ((err < 0) || !ok); - + inexact = mpfr_set (r, t, rnd_mode); mpfr_clear (t); @@ -1,6 +1,6 @@ /* mpfr_log2 -- log base 2 -Copyright 2001, 2002, 2003 Free Software Foundation, Inc. +Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" @@ -39,7 +38,7 @@ mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) { /* If a is NaN, the result is NaN */ if (MPFR_IS_NAN(a)) - { + { MPFR_SET_NAN(r); MPFR_RET_NAN; } @@ -68,7 +67,7 @@ mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) else MPFR_ASSERTN(0); } - + /* If a is negative, the result is NaN */ if (MPFR_UNLIKELY(MPFR_IS_NEG(a))) { @@ -96,10 +95,10 @@ mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) /* 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 */ - + /* compute the precision of intermediary variable */ Nt=MAX(Nx,Ny); @@ -107,22 +106,22 @@ mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) Nt=Nt+3+__gmpfr_ceil_log2(Nt); /* initialise of intermediary variable */ - mpfr_init(t); - mpfr_init(tt); + mpfr_init(t); + mpfr_init(tt); + - /* First computation of log2 */ do { /* reactualisation of the precision */ - mpfr_set_prec(t,Nt); - mpfr_set_prec(tt,Nt); - + 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; @@ -131,7 +130,7 @@ mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) } while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, Ny + (rnd_mode == GMP_RNDN))); - + inexact = mpfr_set (r, t, rnd_mode); mpfr_clear (t); @@ -1,6 +1,6 @@ /* mpfr_reldiff -- compute relative difference of two floating-point numbers. -Copyright 2000, 2001 Free Software Foundation, Inc. +Copyright 2000, 2001, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,29 +19,28 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" #include "mpfr-impl.h" /* reldiff(b, c) = abs(b-c)/b */ -void +void mpfr_reldiff (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) { mpfr_t b_copy; if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) { MPFR_CLEAR_FLAGS(a); MPFR_SET_NAN(a); return; } - if (MPFR_IS_INF(b)) - { + if (MPFR_IS_INF(b)) + { if (MPFR_IS_INF(c) && (MPFR_SIGN(c) == MPFR_SIGN(b))) { MPFR_CLEAR_FLAGS(a); MPFR_SET_ZERO(a); return; } else { MPFR_CLEAR_FLAGS(a); MPFR_SET_NAN(a); return; } } - if (MPFR_IS_INF(c)) + if (MPFR_IS_INF(c)) { MPFR_SET_SAME_SIGN(a, b); MPFR_CLEAR_FLAGS(a); @@ -62,9 +61,8 @@ mpfr_reldiff (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) mpfr_sub (a, b, c, rnd_mode); mpfr_abs (a, a, rnd_mode); /* for compatibility with MPF */ mpfr_div (a, a, (a == b) ? b_copy : b, rnd_mode); - + if (a == b) mpfr_clear (b_copy); } } - diff --git a/save_expo.c b/save_expo.c index 5f2a83651..33c4d70d3 100644 --- a/save_expo.c +++ b/save_expo.c @@ -1,6 +1,6 @@ /* Save/restore the minimum and maximum exponents. -Copyright 2001, 2002 Free Software Foundation, Inc. +Copyright 2001, 2002, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,8 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> -#include <stdlib.h> #include "gmp.h" #include "mpfr.h" #include "mpfr-impl.h" @@ -41,12 +39,9 @@ mpfr_save_emin_emax (void) __gmpfr_emin = MPFR_EMIN_MIN; __gmpfr_emax = MPFR_EMAX_MAX; } - else if (save_ctr == 0) + else { - fprintf(stderr, - "Error: Too many consecutive calls to mpfr_save_emin_emax!\n" - "Probably a bug.\n"); - exit(1); + MPFR_ASSERTN (save_ctr != 0); } } @@ -1,6 +1,6 @@ /* mpfr_set_f -- set a MPFR number from a GNU MPF number -Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" @@ -50,7 +49,7 @@ mpfr_set_f (mpfr_ptr y, mpf_srcptr x, mp_rnd_t rnd_mode) return 0; /* 0 is exact */ } - count_leading_zeros(cnt, mx[sx - 1]); + count_leading_zeros(cnt, mx[sx - 1]); if (sy <= sx) /* we may have to round even when sy = sx */ { @@ -69,14 +68,14 @@ mpfr_set_f (mpfr_ptr y, mpf_srcptr x, mp_rnd_t rnd_mode) else { if (cnt) - mpn_lshift(my + sy - sx, mx, sx, cnt); + mpn_lshift(my + sy - sx, mx, sx, cnt); else - MPN_COPY(my + sy - sx, mx, sy); + MPN_COPY(my + sy - sx, mx, sy); MPN_ZERO(my, sy - sx); /* no rounding necessary, since y has a larger mantissa */ inexact = 0; } - + MPFR_SET_EXP(y, EXP(x) * BITS_PER_MP_LIMB - cnt); return inexact; diff --git a/set_prc_raw.c b/set_prc_raw.c index 2bc8e1d45..aa1eb5afe 100644 --- a/set_prc_raw.c +++ b/set_prc_raw.c @@ -1,6 +1,6 @@ /* mpfr_set_prec_raw -- reset the precision of a floating-point number -Copyright 2000, 2001 Free Software Foundation, Inc. +Copyright 2000, 2001, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,8 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> -#include <stdlib.h> #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" @@ -29,13 +27,7 @@ MA 02111-1307, USA. */ void mpfr_set_prec_raw (mpfr_ptr x, mpfr_prec_t p) { - MPFR_ASSERTN(p >= MPFR_PREC_MIN && p <= MPFR_PREC_MAX); - - if (p > (mpfr_prec_t) MPFR_GET_ALLOC_SIZE(x) * BITS_PER_MP_LIMB) - { - fprintf (stderr, "*** precision too large for allocated space\n"); - exit (1); - } - + MPFR_ASSERTN (p >= MPFR_PREC_MIN && p <= MPFR_PREC_MAX); + MPFR_ASSERTN (p <= (mpfr_prec_t) MPFR_GET_ALLOC_SIZE(x) * BITS_PER_MP_LIMB); MPFR_PREC(x) = p; } diff --git a/set_prec.c b/set_prec.c index 113da143b..7a01cbd8a 100644 --- a/set_prec.c +++ b/set_prec.c @@ -1,6 +1,6 @@ /* mpfr_set_prec -- reset the precision of a floating-point number -Copyright 1999, 2001, 2002 Free Software Foundation, Inc. +Copyright 1999, 2001, 2002, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" @@ -34,7 +33,7 @@ mpfr_set_prec (mpfr_ptr x, mpfr_prec_t p) /* first, check if p is correct */ MPFR_ASSERTN(p >= MPFR_PREC_MIN && p <= MPFR_PREC_MAX); - /* Calculate the new number of limbs */ + /* Calculate the new number of limbs */ xsize = (p - 1) / BITS_PER_MP_LIMB + 1; /* Realloc only if the new size is greater than the old */ @@ -1,6 +1,6 @@ /* mpfr_set_str -- set a floating-point number from a string -Copyright 2000, 2001, 2002, 2003 Free Software Foundation, Inc. +Copyright 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include <stdlib.h> #include <string.h> #include <ctype.h> @@ -144,7 +143,7 @@ mpfr_set_str (mpfr_t x, const char *str, int base, mp_rnd_t rnd) if (**endptr != '\0') { res = -1; /* invalid input: garbage after exponent */ - goto end; + goto end; } break; } @@ -160,18 +159,18 @@ mpfr_set_str (mpfr_t x, const char *str, int base, mp_rnd_t rnd) break; } else if (*str == '.') - { + { point = 1; continue; /* go to next character */ } else if ((c = digit_value_in_base (*str, base)) != -1) *mant_s = c; /* valid character */ - else + else { res = -1; /* invalid input */ goto end; } - if (!point) + if (!point) exp_s ++; mant_s ++; } @@ -209,7 +208,7 @@ mpfr_set_str (mpfr_t x, const char *str, int base, mp_rnd_t rnd) /* remember that y - n is allocated for n limbs */ /* required precision for str */ - pr = (mp_exp_t) __gmpfr_ceil + pr = (mp_exp_t) __gmpfr_ceil (((double) (n * BITS_PER_MP_LIMB) - 1.0) * log_b2[base-2]) + 1; /* check is there are enough digits in str */ @@ -278,7 +277,7 @@ mpfr_set_str (mpfr_t x, const char *str, int base, mp_rnd_t rnd) /* compute the exponent of y */ exp_y += exp_z + n * BITS_PER_MP_LIMB; - + /* normalize result */ if ((result[2 * n - 1] & MPFR_LIMB_HIGHBIT) == 0) { @@ -355,7 +354,7 @@ mpfr_set_str (mpfr_t x, const char *str, int base, mp_rnd_t rnd) negative, MPFR_PREC(x), rnd )) { /* overflaw when rounding y */ - MPFR_MANT(x)[MPFR_LIMB_SIZE(x) - 1] + MPFR_MANT(x)[MPFR_LIMB_SIZE(x) - 1] = MPFR_LIMB_HIGHBIT; exp_y ++; } @@ -1,6 +1,6 @@ /* mpfr_ui_div -- divide a machine integer by a floating-point number -Copyright 2000, 2001, 2002, 2003 Free Software Foundation, Inc. +Copyright 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" @@ -59,7 +58,7 @@ mpfr_ui_div (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mp_rnd_t rnd_mode) { /* 0 / 0 */ MPFR_SET_NAN(y); - MPFR_RET_NAN; + MPFR_RET_NAN; } } else @@ -1,6 +1,6 @@ /* mpfr_zeta -- compute the Riemann Zeta function -Copyright 2003 Free Software Foundation. +Copyright 2003, 2004 Free Software Foundation. Contributed by Jean-Luc Re'my and the Spaces project, INRIA Lorraine. This file is part of the MPFR Library. @@ -22,7 +22,6 @@ MA 02111-1307, USA. */ /* #define DEBUG */ -#include <stdio.h> #include <stdlib.h> #include "gmp.h" #include "gmp-impl.h" |