diff options
author | pelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4> | 2005-01-11 08:42:16 +0000 |
---|---|---|
committer | pelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4> | 2005-01-11 08:42:16 +0000 |
commit | 7e5c3aaba2e3f445c7ad7c0d7798446da297bf31 (patch) | |
tree | 34d80813914bef842c56e8b31a3c0da5c10435be /gmp_op.c | |
parent | 99b175a57ac0f15968aeb8627b2ec6845485def0 (diff) | |
download | mpfr-7e5c3aaba2e3f445c7ad7c0d7798446da297bf31.tar.gz |
Reduce size of code
Optimize a few too.
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@3178 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'gmp_op.c')
-rw-r--r-- | gmp_op.c | 231 |
1 files changed, 126 insertions, 105 deletions
@@ -1,6 +1,6 @@ /* Implementations of operations between mpfr and mpz/mpq data -Copyright 2001, 2003, 2004 Free Software Foundation, Inc. +Copyright 2001, 2003, 2004, 2005 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -20,65 +20,70 @@ the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include <stddef.h> + +#define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" - -int -mpfr_mul_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z,mp_rnd_t rnd_mode) + +/* Init and set a mpfr_t with enought precision to store a mpz */ +static void +init_set_z (mpfr_ptr t, mpz_srcptr z) { - mpfr_t t; - int res; - mpfr_init2(t, mpz_sizeinbase(z, 2) + MPFR_PREC_MIN ); - res = mpfr_set_z(t, z, rnd_mode); - MPFR_ASSERTD(res == 0); - res = mpfr_mul(y, x, t, rnd_mode); - mpfr_clear(t); - return(res); + mp_prec_t p; + int i; + + if (mpz_size (z) <= 1) + p = BITS_PER_MP_LIMB; + else + MPFR_MPZ_SIZEINBASE2 (p, z); + mpfr_init2 (t, p); + i = mpfr_set_z (t, z, GMP_RNDN); + MPFR_ASSERTD (i == 0); } -int -mpfr_div_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd_mode) +/* Init, set a mpfr_t with enought precision to store a mpz_t without round, + call the function, and clear the allocated mpfr_t */ +static int +foo (mpfr_ptr x, mpfr_srcptr y, mpz_srcptr z, mp_rnd_t r, + int (*f)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)) { mpfr_t t; - int res; - mpfr_init2(t, mpz_sizeinbase(z, 2) + MPFR_PREC_MIN ); - res = mpfr_set_z(t, z, rnd_mode); - MPFR_ASSERTD(res == 0); - res = mpfr_div(y, x, t, rnd_mode); - mpfr_clear(t); - return res; + int i; + init_set_z (t, z); + i = (*f) (x, y, t, r); + mpfr_clear (t); + return i; } int -mpfr_add_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd_mode) +mpfr_mul_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r) { - mpfr_t t; - int res; - - if (MPFR_UNLIKELY(mpz_sgn(z) == 0)) - return mpfr_set (y, x, rnd_mode); + return foo (y, x, z, r, mpfr_mul); +} - mpfr_init2 (t, mpz_sizeinbase(z, 2) + MPFR_PREC_MIN ); - res = mpfr_set_z (t, z, rnd_mode); - MPFR_ASSERTD (res == 0); - res = mpfr_add (y, x, t, rnd_mode); - mpfr_clear (t); - return res; +int +mpfr_div_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r) +{ + return foo (y, x, z, r, mpfr_div); } int -mpfr_sub_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z,mp_rnd_t rnd_mode) +mpfr_add_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r) { - mpfr_t t; - int res; + /* Mpz 0 is unsigned */ + if (MPFR_UNLIKELY (mpz_sgn (z) == 0)) + return mpfr_set (y, x, r); + else + return foo (y, x, z, r, mpfr_add); +} - if (MPFR_UNLIKELY(mpz_sgn(z) == 0)) - return mpfr_set (y, x, rnd_mode); - mpfr_init2 (t, mpz_sizeinbase(z, 2) + MPFR_PREC_MIN ); - res = mpfr_set_z (t, z, rnd_mode); - MPFR_ASSERTD (res == 0); - res=mpfr_sub (y, x, t, rnd_mode); - mpfr_clear (t); - return res; +int +mpfr_sub_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r) +{ + /* Mpz 0 is unsigned */ + if (MPFR_UNLIKELY (mpz_sgn (z) == 0)) + return mpfr_set (y, x, r); + else + return foo (y, x, z, r, mpfr_sub); } int @@ -86,26 +91,31 @@ mpfr_cmp_z (mpfr_srcptr x, mpz_srcptr z) { mpfr_t t; int res; - mpfr_init2 (t, mpz_sizeinbase(z, 2) + MPFR_PREC_MIN ); - res = mpfr_set_z (t, z, GMP_RNDN); - MPFR_ASSERTD(res == 0); + init_set_z (t, z); res = mpfr_cmp (x, t); mpfr_clear (t); return res; } int -mpfr_mul_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode) +mpfr_mul_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode) { mpfr_t tmp; int res; - - mpfr_init2 (tmp, MPFR_PREC(x) + mpz_sizeinbase(mpq_numref(z), 2) ); - res = mpfr_mul_z (tmp, x, mpq_numref(z), GMP_RNDN ); - MPFR_ASSERTD( res == 0 ); - res = mpfr_div_z (y, tmp, mpq_denref(z), rnd_mode); - mpfr_clear(tmp); - return(res); + mp_prec_t p; + + if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) + return mpfr_mul_ui (y, x, 0, rnd_mode); + else + { + MPFR_MPZ_SIZEINBASE2 (p, mpq_numref (z)); + mpfr_init2 (tmp, MPFR_PREC (x) + p); + res = mpfr_mul_z (tmp, x, mpq_numref(z), GMP_RNDN ); + MPFR_ASSERTD (res == 0); + res = mpfr_div_z (y, tmp, mpq_denref(z), rnd_mode); + mpfr_clear (tmp); + return res; + } } int @@ -113,13 +123,20 @@ mpfr_div_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode) { mpfr_t tmp; int res; + mp_prec_t p; - mpfr_init2 (tmp, MPFR_PREC(x) + mpz_sizeinbase(mpq_denref(z), 2) ); + if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) + return mpfr_div_ui (y, x, 0, rnd_mode); + else if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z)) == 0)) + p = 0; + else + MPFR_MPZ_SIZEINBASE2 (p, mpq_denref (z)); + mpfr_init2 (tmp, MPFR_PREC(x) + p); res = mpfr_mul_z (tmp, x, mpq_denref(z), GMP_RNDN ); MPFR_ASSERTD( res == 0 ); res = mpfr_div_z (y, tmp, mpq_numref(z), rnd_mode); - mpfr_clear(tmp); - return(res); + mpfr_clear (tmp); + return res; } int @@ -130,63 +147,65 @@ mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode) mp_exp_t err; int res; - 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; } - else if (MPFR_IS_INF(x)) + else if (MPFR_IS_INF (x)) { + MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0); MPFR_SET_INF (y); MPFR_SET_SAME_SIGN (y, x); MPFR_RET (0); } else { - MPFR_ASSERTD( MPFR_IS_ZERO(x)); - - if (MPFR_UNLIKELY( mpq_sgn(z) == 0)) + MPFR_ASSERTD (MPFR_IS_ZERO (x)); + if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */ else return mpfr_set_q (y, z, rnd_mode); } } - mpfr_inits(t, q, NULL); - for(;;) + mpfr_init2 (t, p); + mpfr_init2 (q, p); + for (;;) { - mpfr_set_prec(t, p); - mpfr_set_prec(q, p); - - res = mpfr_set_q(q, z, GMP_RNDN); /* Error <= 1/2 ulp(q) */ + res = mpfr_set_q (q, z, GMP_RNDN); /* Error <= 1/2 ulp(q) */ /* If z if @INF@ (1/0), res = 0, so it quits immediately */ - if (res == 0) /* Result is exact so we can add it directly!*/ + if (MPFR_UNLIKELY (res == 0)) + /* Result is exact so we can add it directly! */ { - res = mpfr_add(y, x, q, rnd_mode); + res = mpfr_add (y, x, q, rnd_mode); break; } - mpfr_add(t, x, q, GMP_RNDN); /* Error <= 1/2 ulp(t) */ + mpfr_add (t, x, q, GMP_RNDN); /* Error <= 1/2 ulp(t) */ /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t)) If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t))) <= 2^(EXP(q)-EXP(t)) If EXP(q)-EXP(t)<0, <= 2^0 */ /* We can get 0, but we can't round since q is inexact */ - if (MPFR_LIKELY( !MPFR_IS_ZERO (t))) + if (MPFR_LIKELY (!MPFR_IS_ZERO (t))) { - err = (mp_exp_t) p - 1 - MAX(MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0); + 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 (res != 0) /* We can round! */ + if (MPFR_LIKELY (res != 0)) /* We can round! */ { res = mpfr_set(y, t, rnd_mode); break; } } - p += BITS_PER_MP_LIMB; /* Next precision if we continue */ + p += BITS_PER_MP_LIMB; /* Next precision if we continue */ + mpfr_set_prec (t, p); + mpfr_set_prec (q, p); } - mpfr_clears(t, q, NULL); + mpfr_clear (t); + mpfr_clear (q); return res; } @@ -194,75 +213,74 @@ 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 = MPFR_PREC (y)+10; int res; mp_exp_t err; - 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_RET_NAN; } - else if (MPFR_IS_INF(x)) + else if (MPFR_IS_INF (x)) { + MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0); MPFR_SET_INF (y); MPFR_SET_SAME_SIGN (y, x); MPFR_RET (0); } else { - MPFR_ASSERTD( MPFR_IS_ZERO(x)); + MPFR_ASSERTD (MPFR_IS_ZERO (x)); - if (MPFR_UNLIKELY( mpq_sgn(z) == 0)) + if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */ else { - if (rnd_mode == GMP_RNDU) - rnd_mode = GMP_RNDD; - else if (rnd_mode == GMP_RNDD) - rnd_mode = GMP_RNDU; - res = mpfr_set_q (y, z, rnd_mode); + res = mpfr_set_q (y, z, MPFR_INVERT_RND (rnd_mode)); MPFR_CHANGE_SIGN (y); return -res; } } } - mpfr_inits(t, q, NULL); + mpfr_init2 (t, p); + mpfr_init2 (q, p); for(;;) { - mpfr_set_prec(t, p); - mpfr_set_prec(q, p); - res = mpfr_set_q(q, z, GMP_RNDN); /* Error <= 1/2 ulp(q) */ /* If z if @INF@ (1/0), res = 0, so it quits immediately */ - if (res == 0) /* Result is exact so we can add it directly!*/ + if (MPFR_UNLIKELY (res == 0)) + /* Result is exact so we can add it directly!*/ { - res = mpfr_sub(y, x, q, rnd_mode); + res = mpfr_sub (y, x, q, rnd_mode); break; } - mpfr_sub(t, x, q, GMP_RNDN); /* Error <= 1/2 ulp(t) */ + mpfr_sub (t, x, q, GMP_RNDN); /* Error <= 1/2 ulp(t) */ /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t)) If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t))) <= 2^(EXP(q)-EXP(t)) If EXP(q)-EXP(t)<0, <= 2^0 */ /* We can get 0, but we can't round since q is inexact */ - if (MPFR_LIKELY( !MPFR_IS_ZERO (t))) + if (MPFR_LIKELY (!MPFR_IS_ZERO (t))) { - err = (mp_exp_t) p - 1 - MAX(MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0); + 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 (res != 0) /* We can round! */ + if (MPFR_LIKELY (res != 0)) /* We can round! */ { - 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_set_prec (t, p); + mpfr_set_prec (q, p); } - mpfr_clears(t, q, NULL); + mpfr_clear (t); + mpfr_clear (q); return res; } @@ -271,10 +289,13 @@ mpfr_cmp_q (mpfr_srcptr x, mpq_srcptr z) { mpfr_t t; int res; + mp_prec_t p; /* x < a/b ? <=> x*b < a */ - mpfr_init2 (t, MPFR_PREC(x) + mpz_sizeinbase (mpq_denref(z), 2) ); + MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0); + MPFR_MPZ_SIZEINBASE2 (p, mpq_denref (z)); + mpfr_init2 (t, MPFR_PREC(x) + p); res = mpfr_mul_z (t, x, mpq_denref (z), GMP_RNDN ); - MPFR_ASSERTD( res == 0 ); + MPFR_ASSERTD (res == 0); res = mpfr_cmp_z (t, mpq_numref (z) ); mpfr_clear (t); return res; @@ -288,7 +309,7 @@ mpfr_cmp_f (mpfr_srcptr x, mpf_srcptr z) mpfr_init2 (t, MPFR_PREC_MIN + ABS(SIZ(z)) * BITS_PER_MP_LIMB ); res = mpfr_set_f (t, z, GMP_RNDN); - MPFR_ASSERTD( res == 0 ); + MPFR_ASSERTD (res == 0); res = mpfr_cmp (x, t); mpfr_clear (t); return res; |