From 1e5272f459eaeace3e2391e6739b65b8d1c7855a Mon Sep 17 00:00:00 2001 From: enge Date: Mon, 1 Oct 2012 13:39:46 +0000 Subject: log10.c: rewrite to source out more to mpc_log git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@1280 211d60ee-9f03-0410-a15a-8952a2c7a4e4 --- src/log10.c | 291 ++++++++++++++++++++---------------------------------------- 1 file changed, 98 insertions(+), 193 deletions(-) diff --git a/src/log10.c b/src/log10.c index 0d35c76..a262c72 100644 --- a/src/log10.c +++ b/src/log10.c @@ -21,213 +21,118 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include /* for CHAR_BIT */ #include "mpc-impl.h" -/* Auxiliary function which implements Ziv's strategy for special cases. - Exact cases should be dealt with separately. */ -static int -mpc_log10_aux (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, int nb) -{ - mp_prec_t prec = (MPFR_PREC_MIN > 4) ? MPFR_PREC_MIN : 4; - mpc_t tmp; - mpfr_t log10; - int ok = 0, ret; - - prec = mpfr_get_prec (mpc_imagref (rop)); - prec += 10; - mpc_init2 (tmp, prec); - mpfr_init2 (log10, prec); - while (ok == 0) - { - mpfr_set_ui (log10, 10, MPFR_RNDN); /* exact since prec >= 4 */ - mpfr_log (log10, log10, MPFR_RNDN); - /* In each case we have two roundings, thus the final value is - x * (1+u)^2 where x is the exact value, and |u| <= 2^(-prec-1). - Thus the error is always less than 3 ulps. */ - if (nb == 0) - mpfr_atan2 (mpc_imagref (tmp), mpc_imagref (op), mpc_realref (op), - MPC_RND_IM (rnd)); - else - mpfr_const_pi (mpc_imagref (tmp), MPC_RND_IM (rnd)); - mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, MPFR_RNDN); - ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, MPFR_RNDN, - MPFR_RNDZ, MPC_PREC_IM(rop) + (MPC_RND_IM (rnd) == MPFR_RNDN)); - if (ok) - ret = mpfr_set (mpc_imagref (rop), mpc_imagref (tmp), - MPC_RND_IM (rnd)); - prec += prec / 2; - mpc_set_prec (tmp, prec); - mpfr_set_prec (log10, prec); - } - mpc_clear (tmp); - mpfr_clear (log10); - return ret; -} int mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { - int ok = 0, loops = 0, re_cmp, im_cmp, inex_re, inex_im, negative_zero; - mpfr_t w; - mpfr_prec_t prec; - mpfr_rnd_t rnd_im; - mpc_t ww; - mpc_rnd_t invrnd; - - /* special values: NaN and infinities: same as mpc_log */ - if (!mpc_fin_p (op)) { - if ( mpfr_nan_p (mpc_realref (op)) - || mpfr_nan_p (mpc_imagref (op)) - || ( mpfr_inf_p (mpc_realref (op)) - && mpfr_signbit (mpc_realref (op))== 0 - && mpfr_number_p (mpc_imagref (op)))) - return mpc_log (rop, op, rnd); + int ok = 0, loops = 0, special_re, special_im, inex, inex_re, inex_im; + mpfr_t w; + mpfr_prec_t prec; + mpc_t ww; + + mpfr_init2 (w, 2); + mpc_init2 (ww, 2); + prec = MPC_PREC_RE(rop); + /* compute log(op)/log(10) */ + while (ok == 0) { + loops ++; + prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2; + mpfr_set_prec (w, prec); + mpc_set_prec (ww, prec); + + inex = mpc_log (ww, op, rnd); /* error <= 1 ulp */ + + if (!mpfr_number_p (mpc_imagref (ww)) + || mpfr_zero_p (mpc_imagref (ww))) { + /* no need to divide by log(10) */ + special_im = 1; + ok = 1; + } else { - /* (xxx, Inf) -> (+Inf, atan2(Inf/xxx)) - (Inf, yyy) -> (+Inf, atan2(yyy/Inf)) */ - inex_im = mpc_log10_aux (rop, op, rnd, 0); - mpfr_set_inf (mpc_realref (rop), +1); - return MPC_INEX(0, inex_im); + special_im = 0; + mpfr_set_ui (w, 10, MPFR_RNDN); /* exact since prec >= 4 */ + mpfr_log (w, w, MPFR_RNDN); /* error <= 1/2 ulp */ + mpfr_div (mpc_imagref (ww), mpc_imagref (ww), w, MPFR_RNDN); + + ok = mpfr_can_round (mpc_imagref (ww), prec - 2, MPFR_RNDN, MPFR_RNDZ, + MPC_PREC_IM(rop) + (MPC_RND_IM (rnd) == MPFR_RNDN)); } - } - /* special cases: real and purely imaginary numbers */ - re_cmp = mpfr_cmp_ui (mpc_realref (op), 0); - im_cmp = mpfr_cmp_ui (mpc_imagref (op), 0); - if (im_cmp == 0) /* Im(op) = 0 */ - { - if (re_cmp == 0) /* Re(op) = 0 */ - { - if (mpfr_signbit (mpc_realref (op)) == 0) - inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), - mpc_realref (op), MPC_RND_IM (rnd)); - else - inex_im = mpc_log10_aux (rop, op, rnd, 0); - mpfr_set_inf (mpc_realref (rop), -1); - inex_re = 0; /* -Inf is exact */ - } - else if (re_cmp > 0) - { - inex_re = mpfr_log10 (mpc_realref (rop), mpc_realref (op), - MPC_RND_RE (rnd)); - inex_im = mpfr_set (mpc_imagref (rop), mpc_imagref (op), - MPC_RND_IM (rnd)); - } - else /* log10(x + 0*i) for negative x */ - { /* op = x + 0*i; let w = -x = |x| */ - negative_zero = mpfr_signbit (mpc_imagref (op)); - if (negative_zero) - rnd_im = INV_RND (MPC_RND_IM (rnd)); - else - rnd_im = MPC_RND_IM (rnd); - w [0] = *mpc_realref (op); - MPFR_CHANGE_SIGN (w); - inex_re = mpfr_log10 (mpc_realref (rop), w, MPC_RND_RE (rnd)); - inex_im = mpc_log10_aux (rop, op, MPC_RND (0,rnd_im), 2); - if (negative_zero) - { - mpc_conj (rop, rop, MPC_RNDNN); - inex_im = -inex_im; + if (ok) { + if (!mpfr_number_p (mpc_realref (ww)) + || mpfr_zero_p (mpc_realref (ww))) + special_re = 1; + else { + special_re = 0; + if (special_im) { + /* log10 not yet computed */ + mpfr_set_ui (w, 10, MPFR_RNDN); + mpfr_log (w, w, MPFR_RNDN); } - } - return MPC_INEX(inex_re, inex_im); - } - else if (re_cmp == 0) - { - if (im_cmp > 0) - { - inex_re = mpfr_log10 (mpc_realref (rop), mpc_imagref (op), MPC_RND_RE (rnd)); - inex_im = mpc_log10_aux (rop, op, rnd, 2); - /* division by 2 does not change the ternary flag */ - mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, MPFR_RNDN); - } - else - { - w [0] = *mpc_imagref (op); - MPFR_CHANGE_SIGN (w); - inex_re = mpfr_log10 (mpc_realref (rop), w, MPC_RND_RE (rnd)); - invrnd = MPC_RND (0, INV_RND (MPC_RND_IM (rnd))); - inex_im = mpc_log10_aux (rop, op, invrnd, 2); - /* division by 2 does not change the ternary flag */ - mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, MPFR_RNDN); - mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), MPFR_RNDN); - inex_im = -inex_im; /* negate the ternary flag */ + mpfr_div (mpc_realref (ww), mpc_realref (ww), w, MPFR_RNDN); + /* error <= 24/7 ulp < 4 ulp for prec >= 4, see algorithms.tex */ + + ok = mpfr_can_round (mpc_realref (ww), prec - 2, MPFR_RNDN, + MPFR_RNDZ, + MPC_PREC_RE(rop) + (MPC_RND_RE (rnd) == MPFR_RNDN)); } - return MPC_INEX(inex_re, inex_im); - } - - /* generic case: neither Re(op) nor Im(op) is NaN, Inf or zero */ - prec = MPC_PREC_RE(rop); - mpfr_init2 (w, prec); - mpc_init2 (ww, prec); - /* let op = x + iy; compute log(op)/log(10) */ - while (ok == 0) - { - loops ++; - prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2; - mpfr_set_prec (w, prec); - mpc_set_prec (ww, prec); - mpc_log (ww, op, MPC_RNDNN); - mpfr_set_ui (w, 10, MPFR_RNDN); /* exact since prec >= 4 */ - mpfr_log (w, w, MPFR_RNDN); - mpc_div_fr (ww, ww, w, MPC_RNDNN); - - ok = mpfr_can_round (mpc_realref (ww), prec - 2, MPFR_RNDN, MPFR_RNDZ, - MPC_PREC_RE(rop) + (MPC_RND_RE (rnd) == MPFR_RNDN)); - - /* Special code to deal with cases where the real part of log10(x+i*y) - is exact, like x=3 and y=1. Since Re(log10(x+i*y)) = log10(x^2+y^2)/2 - this happens whenever x^2+y^2 is a nonnegative power of 10. - Indeed x^2+y^2 cannot equal 10^(a/2^b) for a, b integers, a odd, b>0, - since x^2+y^2 is rational, and 10^(a/2^b) is irrational. - Similarly, for b=0, x^2+y^2 cannot equal 10^a for a < 0 since x^2+y^2 - is a rational with denominator a power of 2. - Now let x^2+y^2 = 10^s. Without loss of generality we can assume - x = u/2^e and y = v/2^e with u, v, e integers: u^2+v^2 = 10^s*2^(2e) - thus u^2+v^2 = 0 mod 2^(2e). By recurrence on e, necessarily - u = v = 0 mod 2^e, thus x and y are necessarily integers. - */ - if ((ok == 0) && (loops == 1) && mpfr_integer_p (mpc_realref (op)) && - mpfr_integer_p (mpc_imagref (op))) - { - mpz_t x, y; - unsigned long s, v; - - mpz_init (x); - mpz_init (y); - mpfr_get_z (x, mpc_realref (op), MPFR_RNDN); /* exact */ - mpfr_get_z (y, mpc_imagref (op), MPFR_RNDN); /* exact */ - mpz_mul (x, x, x); - mpz_mul (y, y, y); - mpz_add (x, x, y); /* x^2+y^2 */ - v = mpz_scan1 (x, 0); - /* if x = 10^s then necessarily s = v */ - s = mpz_sizeinbase (x, 10); - /* since s is either the number of digits of x or one more, - then x = 10^(s-1) or 10^(s-2) */ - if (s == v + 1 || s == v + 2) - { - mpz_div_2exp (x, x, v); - mpz_ui_pow_ui (y, 5, v); - if (mpz_cmp (y, x) == 0) /* Re(log10(x+i*y)) is exactly v/2 */ - { + /* Special code to deal with cases where the real part of log10(x+i*y) + is exact, like x=3 and y=1. Since Re(log10(x+i*y)) = log10(x^2+y^2)/2 + this happens whenever x^2+y^2 is a nonnegative power of 10. + Indeed x^2+y^2 cannot equal 10^(a/2^b) for a, b integers, a odd, b>0, + since x^2+y^2 is rational, and 10^(a/2^b) is irrational. + Similarly, for b=0, x^2+y^2 cannot equal 10^a for a < 0 since x^2+y^2 + is a rational with denominator a power of 2. + Now let x^2+y^2 = 10^s. Without loss of generality we can assume + x = u/2^e and y = v/2^e with u, v, e integers: u^2+v^2 = 10^s*2^(2e) + thus u^2+v^2 = 0 mod 2^(2e). By recurrence on e, necessarily + u = v = 0 mod 2^e, thus x and y are necessarily integers. + */ + if (!ok && (loops == 1) && mpfr_integer_p (mpc_realref (op)) && + mpfr_integer_p (mpc_imagref (op))) { + mpz_t x, y; + unsigned long s, v; + + mpz_init (x); + mpz_init (y); + mpfr_get_z (x, mpc_realref (op), MPFR_RNDN); /* exact */ + mpfr_get_z (y, mpc_imagref (op), MPFR_RNDN); /* exact */ + mpz_mul (x, x, x); + mpz_mul (y, y, y); + mpz_add (x, x, y); /* x^2+y^2 */ + v = mpz_scan1 (x, 0); + /* if x = 10^s then necessarily s = v */ + s = mpz_sizeinbase (x, 10); + /* since s is either the number of digits of x or one more, + then x = 10^(s-1) or 10^(s-2) */ + if (s == v + 1 || s == v + 2) { + mpz_div_2exp (x, x, v); + mpz_ui_pow_ui (y, 5, v); + if (mpz_cmp (y, x) == 0) { + /* Re(log10(x+i*y)) is exactly v/2 */ /* we reset the precision of Re(ww) so that v can be represented exactly */ mpfr_set_prec (mpc_realref (ww), sizeof(unsigned long)*CHAR_BIT); mpfr_set_ui_2exp (mpc_realref (ww), v, -1, MPFR_RNDN); /* exact */ ok = 1; - } + } } - mpz_clear (x); - mpz_clear (y); - } - - ok = ok && mpfr_can_round (mpc_imagref (ww), prec-2, MPFR_RNDN, MPFR_RNDZ, - MPC_PREC_IM(rop) + (MPC_RND_IM (rnd) == MPFR_RNDN)); - } - - inex_re = mpfr_set (mpc_realref(rop), mpc_realref (ww), MPC_RND_RE (rnd)); - inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref (ww), MPC_RND_IM (rnd)); - mpfr_clear (w); - mpc_clear (ww); - return MPC_INEX(inex_re, inex_im); + mpz_clear (x); + mpz_clear (y); + } + } + } + + inex_re = mpfr_set (mpc_realref(rop), mpc_realref (ww), MPC_RND_RE (rnd)); + if (special_re) + inex_re = MPC_INEX_RE (inex); + /* recover flag from call to mpc_log above */ + inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref (ww), MPC_RND_IM (rnd)); + if (special_im) + inex_im = MPC_INEX_IM (inex); + mpfr_clear (w); + mpc_clear (ww); + + return MPC_INEX(inex_re, inex_im); } -- cgit v1.2.1