summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorpelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4>2004-08-24 14:00:42 +0000
committerpelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4>2004-08-24 14:00:42 +0000
commite154e3a84fa8802ae9303e634a1dbb758be4c53c (patch)
tree75e27e242276083fd3ce1e273d0d5aced92f1006
parent6f8448cc348b37bd97cd2ccdf945001b5ef20366 (diff)
downloadmpfr-e154e3a84fa8802ae9303e634a1dbb758be4c53c.tar.gz
Try to retype correctly the functions (replace int by the correct type).
Still some works to do. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@2956 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r--exp_2.c148
1 files changed, 84 insertions, 64 deletions
diff --git a/exp_2.c b/exp_2.c
index be9efc4f1..24448cc93 100644
--- a/exp_2.c
+++ b/exp_2.c
@@ -24,10 +24,14 @@ MA 02111-1307, USA. */
#define MPFR_NEED_LONGLONG_H /* for count_leading_zeros */
#include "mpfr-impl.h"
-static int mpfr_exp2_aux _MPFR_PROTO ((mpz_t, mpfr_srcptr, int, int*));
-static int mpfr_exp2_aux2 _MPFR_PROTO ((mpz_t, mpfr_srcptr, int, int*));
-static mp_exp_t mpz_normalize _MPFR_PROTO ((mpz_t, mpz_t, int));
-static int mpz_normalize2 _MPFR_PROTO ((mpz_t, mpz_t, int, int));
+static unsigned long
+mpfr_exp2_aux (mpz_t, mpfr_srcptr, mp_prec_t, mp_exp_t *);
+static unsigned long
+mpfr_exp2_aux2 (mpz_t, mpfr_srcptr, mp_prec_t, mp_exp_t *);
+static mp_exp_t
+mpz_normalize (mpz_t, mpz_t, mp_exp_t);
+static mp_exp_t
+mpz_normalize2 (mpz_t, mpz_t, mp_exp_t, mp_exp_t);
#define SWITCH 100 /* number of bits to switch from O(n^(1/2)*M(n)) method
to O(n^(1/3)*M(n)) method */
@@ -41,30 +45,33 @@ static int mpz_normalize2 _MPFR_PROTO ((mpz_t, mpz_t, int, int));
Otherwise do nothing and return 0.
*/
static mp_exp_t
-mpz_normalize (mpz_t rop, mpz_t z, int q)
+mpz_normalize (mpz_t rop, mpz_t z, mp_exp_t q)
{
- int k;
+ size_t k;
k = mpz_sizeinbase(z, 2);
- if (k > q) {
- mpz_div_2exp(rop, z, k-q);
- return k-q;
- }
- else {
- if (rop != z) mpz_set(rop, z);
- return 0;
- }
+ MPFR_ASSERTD (k == (mpfr_uexp_t) k);
+ if (q < 0 || (mpfr_uexp_t) k > (mpfr_uexp_t) q)
+ {
+ mpz_div_2exp(rop, z, (unsigned long) ((mpfr_uexp_t) k - q));
+ return (mp_exp_t) k - q;
+ }
+ if (MPFR_UNLIKELY(rop != z))
+ mpz_set(rop, z);
+ return 0;
}
/* if expz > target, shift z by (expz-target) bits to the left.
if expz < target, shift z by (target-expz) bits to the right.
Returns target.
*/
-static int
-mpz_normalize2 (mpz_t rop, mpz_t z, int expz, int target)
+static mp_exp_t
+mpz_normalize2 (mpz_t rop, mpz_t z, mp_exp_t expz, mp_exp_t target)
{
- if (target > expz) mpz_div_2exp(rop, z, target-expz);
- else mpz_mul_2exp(rop, z, expz-target);
+ if (target > expz)
+ mpz_div_2exp(rop, z, target-expz);
+ else
+ mpz_mul_2exp(rop, z, expz-target);
return target;
}
@@ -76,7 +83,11 @@ mpz_normalize2 (mpz_t rop, mpz_t z, int expz, int target)
int
mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
- int n, K, precy, q, k, l, err, exps, error_r;
+ long n;
+ unsigned long K, k, l, err; /* FIXME: Which type ? */
+ int error_r;
+ mp_exp_t exps;
+ mp_prec_t q, precy;
int inexact;
mpfr_t r, s, t;
mpz_t ss;
@@ -87,7 +98,7 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
MPFR_TRACE ( printf("Py=%d Px=%d", MPFR_PREC(y), MPFR_PREC(x)) );
MPFR_TRACE ( MPFR_DUMP (x) );
- n = (int) (mpfr_get_d1 (x) / LOG2);
+ n = (long) (mpfr_get_d1 (x) / LOG2);
/* error bounds the cancelled bits in x - n*log(2) */
if (MPFR_UNLIKELY(n == 0))
@@ -150,8 +161,9 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
MY_INIT_MPZ(ss, 3 + 2*((q-1)/BITS_PER_MP_LIMB));
exps = mpfr_get_z_exp (ss, s);
/* s <- 1 + r/1! + r^2/2! + ... + r^l/l! */
- l = (precy < SWITCH) ? mpfr_exp2_aux (ss, r, q, &exps) /* naive method */
- : mpfr_exp2_aux2 (ss, r, q, &exps); /* Brent/Kung method */
+ l = (precy < SWITCH) ?
+ mpfr_exp2_aux (ss, r, q, &exps) /* naive method */
+ : mpfr_exp2_aux2 (ss, r, q, &exps); /* Brent/Kung method */
MPFR_TRACE(printf("l=%d q=%d (K+l)*q^2=%1.3e\n", l, q, (K+l)*(double)q*q));
@@ -172,11 +184,10 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
mpfr_div_2ui(s, s, -n, GMP_RNDU);
/* error is at most 2^K*(3l*(l+1)) ulp for mpfr_exp2_aux */
- if (precy<SWITCH)
- l = 3*l*(l+1);
- else
- l = l*(l+4);
- k=0; while (l) { k++; l >>= 1; }
+ l = (precy < SWITCH) ? 3*l*(l+1) : l*(l+4) ;
+ k = MPFR_INT_CEIL_LOG2 (l);
+ /* k = 0; while (l) { k++; l >>= 1; } */
+
/* now k = ceil(log(error in ulps)/log(2)) */
K += k;
@@ -213,34 +224,39 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
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)
*/
-static int
-mpfr_exp2_aux (mpz_t s, mpfr_srcptr r, int q, int *exps)
+static unsigned long
+mpfr_exp2_aux (mpz_t s, mpfr_srcptr r, mp_prec_t q, mp_exp_t *exps)
{
- int l, dif, qn;
+ unsigned long l;
+ mp_exp_t dif;
+ mp_size_t qn;
mpz_t t, rr;
mp_exp_t expt, expr;
TMP_DECL(marker);
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
- need an extra limb */
+ expt = 0;
+ *exps = 1 - (mp_exp_t) q; /* s = 2^(q-1) */
+ MY_INIT_MPZ(t, 2*qn+1);
MY_INIT_MPZ(rr, qn+1);
- mpz_set_ui(t, 1); expt=0;
- mpz_set_ui(s, 1); mpz_mul_2exp(s, s, q-1); *exps = 1-q; /* s = 2^(q-1) */
- expr = mpfr_get_z_exp(rr, r); /* no error here */
+ mpz_set_ui(t, 1);
+ mpz_set_ui(s, 1);
+ mpz_mul_2exp(s, s, q-1);
+ expr = mpfr_get_z_exp(rr, r); /* no error here */
l = 0;
do {
l++;
- mpz_mul(t, t, rr); expt=expt+expr;
+ mpz_mul(t, t, rr);
+ expt += expr;
dif = *exps + mpz_sizeinbase(s, 2) - expt - mpz_sizeinbase(t, 2);
/* truncates the bits of t which are < ulp(s) = 2^(1-q) */
- expt += mpz_normalize(t, t, q-dif); /* error at most 2^(1-q) */
- mpz_div_ui(t, t, l); /* error at most 2^(1-q) */
+ expt += mpz_normalize(t, t, (mp_exp_t) q-dif); /* error at most 2^(1-q) */
+ mpz_div_ui(t, t, l); /* error at most 2^(1-q) */
/* the error wrt t^l/l! is here at most 3*l*ulp(s) */
MPFR_ASSERTD (expt == *exps);
- mpz_add(s, s, t); /* no error here: exact */
+ mpz_add(s, s, t); /* no error here: exact */
/* ensures rr has the same size as t: after several shifts, the error
on rr is still at most ulp(t)=ulp(s) */
expr += mpz_normalize(rr, rr, mpz_sizeinbase(t, 2));
@@ -259,45 +275,51 @@ mpfr_exp2_aux (mpz_t s, mpfr_srcptr r, int q, int *exps)
Version using mpz. ss must have at least (sizer+1) limbs.
The error is bounded by (l^2+4*l) ulps where l is the return value.
*/
-static int
-mpfr_exp2_aux2 (mpz_t s, mpfr_srcptr r, int q, int *exps)
+static unsigned long
+mpfr_exp2_aux2 (mpz_t s, mpfr_srcptr r, mp_prec_t q, mp_exp_t *exps)
{
- int expr, l, m, i, sizer, *expR, expt, ql;
+ mp_exp_t expr, *expR, expt;
+ mp_size_t sizer;
+ mp_prec_t ql;
+ unsigned long l, m, i;
mpz_t t, *R, rr, tmp;
TMP_DECL(marker);
/* estimate value of l */
+ MPFR_ASSERTD (MPFR_GET_EXP (r) < 0);
l = q / (- MPFR_GET_EXP (r));
- m = (int) __gmpfr_isqrt (l);
+ m = __gmpfr_isqrt (l);
/* we access R[2], thus we need m >= 2 */
if (m < 2)
m = 2;
+
TMP_MARK(marker);
- R = (mpz_t*) TMP_ALLOC((m+1)*sizeof(mpz_t)); /* R[i] stands for r^i */
- expR = (int*) TMP_ALLOC((m+1)*sizeof(int)); /* exponent for R[i] */
+ R = (mpz_t*) TMP_ALLOC((m+1)*sizeof(mpz_t)); /* R[i] is r^i */
+ expR = (mp_exp_t*) TMP_ALLOC((m+1)*sizeof(mp_exp_t)); /* exponent for R[i] */
sizer = 1 + (MPFR_PREC(r)-1)/BITS_PER_MP_LIMB;
mpz_init(tmp);
MY_INIT_MPZ(rr, sizer+2);
- MY_INIT_MPZ(t, 2*sizer); /* double size for products */
- mpz_set_ui(s, 0); *exps = 1-q; /* initialize s to zero, 1 ulp = 2^(1-q) */
- for (i=0;i<=m;i++)
+ MY_INIT_MPZ(t, 2*sizer); /* double size for products */
+ mpz_set_ui(s, 0);
+ *exps = 1-q; /* 1 ulp = 2^(1-q) */
+ for (i = 0 ; i <= m ; i++)
MY_INIT_MPZ(R[i], sizer+2);
expR[1] = mpfr_get_z_exp(R[1], r); /* exact operation: no error */
expR[1] = mpz_normalize2(R[1], R[1], expR[1], 1-q); /* error <= 1 ulp */
mpz_mul(t, R[1], R[1]); /* err(t) <= 2 ulps */
mpz_div_2exp(R[2], t, q-1); /* err(R[2]) <= 3 ulps */
expR[2] = 1-q;
- for (i=3;i<=m;i++)
+ for (i = 3 ; i <= m ; i++)
{
mpz_mul(t, R[i-1], R[1]); /* err(t) <= 2*i-2 */
mpz_div_2exp(R[i], t, q-1); /* err(R[i]) <= 2*i-1 ulps */
expR[i] = 1-q;
}
- mpz_set_ui(R[0], 1);
- mpz_mul_2exp(R[0], R[0], q-1);
- expR[0]=1-q; /* R[0]=1 */
- mpz_set_ui(rr, 1);
- expr=0; /* rr contains r^l/l! */
+ mpz_set_ui (R[0], 1);
+ mpz_mul_2exp (R[0], R[0], q-1);
+ expR[0] = 1-q; /* R[0]=1 */
+ mpz_set_ui (rr, 1);
+ expr = 0; /* rr contains r^l/l! */
/* by induction: err(rr) <= 2*l ulps */
l = 0;
@@ -305,17 +327,15 @@ mpfr_exp2_aux2 (mpz_t s, mpfr_srcptr r, int q, int *exps)
do
{
/* all R[i] must have exponent 1-ql */
- if (l)
- for (i=0;i<m;i++)
- {
- expR[i] = mpz_normalize2(R[i], R[i], expR[i], 1-ql);
- }
+ if (l != 0)
+ for (i = 0 ; i < m ; i++)
+ 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)!
using Horner's scheme */
- for (i=m-2;i>=0;i--)
+ for (i = m-1 ; i-- != 0 ; )
{
mpz_div_ui(t, t, l+i+1); /* err(t) += 1 ulp */
mpz_add(t, t, R[i]);
@@ -335,12 +355,12 @@ mpfr_exp2_aux2 (mpz_t s, mpfr_srcptr r, int q, int *exps)
mpz_mul(t, rr, R[m]); /* err(t) <= err(rr) + 2m-1 */
expr += expR[m];
mpz_set_ui (tmp, 1);
- for (i=1; i<=m; i++)
- mpz_mul_ui (tmp, tmp, l + i);
+ for (i = 1 ; i <= m ; i++)
+ mpz_mul_ui (tmp, tmp, l + i);
mpz_fdiv_q(t, t, tmp); /* err(t) <= err(rr) + 2m */
expr += mpz_normalize(rr, t, ql); /* err_rr(l+1) <= err_rr(l) + 2m+1 */
ql = q - *exps - mpz_sizeinbase(s, 2) + expr + mpz_sizeinbase(rr, 2);
- l+=m;
+ l += m;
}
while ((size_t) expr+mpz_sizeinbase(rr, 2) > (size_t)((int)-q));