diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makefile.am | 13 | ||||
| -rw-r--r-- | src/acos.c | 37 | ||||
| -rw-r--r-- | src/acosh.c | 6 | ||||
| -rw-r--r-- | src/asin.c | 32 | ||||
| -rw-r--r-- | src/asinh.c | 4 | ||||
| -rw-r--r-- | src/atan.c | 84 | ||||
| -rw-r--r-- | src/atanh.c | 4 | ||||
| -rw-r--r-- | src/div.c | 136 | ||||
| -rw-r--r-- | src/div_2si.c (renamed from src/div_2exp.c) | 10 | ||||
| -rw-r--r-- | src/div_2ui.c (renamed from src/mul_2exp.c) | 10 | ||||
| -rw-r--r-- | src/div_fr.c | 4 | ||||
| -rw-r--r-- | src/exp.c | 30 | ||||
| -rw-r--r-- | src/fma.c | 34 | ||||
| -rw-r--r-- | src/fr_div.c | 4 | ||||
| -rw-r--r-- | src/get_version.c | 25 | ||||
| -rw-r--r-- | src/log.c | 129 | ||||
| -rw-r--r-- | src/log10.c | 96 | ||||
| -rw-r--r-- | src/mpc-impl.h | 57 | ||||
| -rw-r--r-- | src/mpc.h | 49 | ||||
| -rw-r--r-- | src/mul.c | 60 | ||||
| -rw-r--r-- | src/mul_2si.c | 32 | ||||
| -rw-r--r-- | src/mul_2ui.c | 32 | ||||
| -rw-r--r-- | src/mul_fr.c | 4 | ||||
| -rw-r--r-- | src/mul_i.c | 6 | ||||
| -rw-r--r-- | src/norm.c | 30 | ||||
| -rw-r--r-- | src/pow.c | 37 | ||||
| -rw-r--r-- | src/pow_fr.c | 4 | ||||
| -rw-r--r-- | src/pow_ui.c | 8 | ||||
| -rw-r--r-- | src/sin_cos.c | 42 | ||||
| -rw-r--r-- | src/sinh.c | 4 | ||||
| -rw-r--r-- | src/sqr.c | 35 | ||||
| -rw-r--r-- | src/sqrt.c | 98 | ||||
| -rw-r--r-- | src/strtoc.c | 4 | ||||
| -rw-r--r-- | src/tan.c | 24 | ||||
| -rw-r--r-- | src/tanh.c | 4 |
35 files changed, 648 insertions, 540 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index a9bfa9c..c52adcc 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,6 @@ ## src/Makefile.am -- Process this file with automake to produce Makefile.in ## -## Copyright (C) 2008, 2009, 2010, 2011 INRIA +## Copyright (C) 2008, 2009, 2010, 2011, 2012 INRIA ## ## This file is part of GNU MPC. ## @@ -18,13 +18,14 @@ ## along with this program. If not, see http://www.gnu.org/licenses/ . lib_LTLIBRARIES = libmpc.la -libmpc_la_LDFLAGS = $(MPC_LDFLAGS) -version-info 2:0:0 +libmpc_la_LDFLAGS = $(MPC_LDFLAGS) -version-info 3:0:0 libmpc_la_SOURCES = mpc-impl.h abs.c acos.c acosh.c add.c add_fr.c \ add_si.c add_ui.c arg.c asin.c asinh.c atan.c atanh.c clear.c cmp.c \ - cmp_si_si.c conj.c cos.c cosh.c div_2exp.c div.c div_fr.c div_ui.c exp.c \ - fma.c fr_div.c fr_sub.c get_prec2.c get_prec.c get_version.c get_x.c \ - imag.c init2.c init3.c inp_str.c log.c log10.c mem.c mul_2exp.c mul.c \ - mul_fr.c mul_i.c mul_si.c mul_ui.c neg.c norm.c out_str.c pow.c pow_fr.c \ + cmp_si_si.c conj.c cos.c cosh.c div_2si.c div_2ui.c div.c div_fr.c \ + div_ui.c exp.c fma.c fr_div.c fr_sub.c get_prec2.c get_prec.c \ + get_version.c get_x.c imag.c init2.c init3.c inp_str.c log.c log10.c \ + mem.c mul_2si.c mul_2ui.c mul.c mul_fr.c mul_i.c mul_si.c mul_ui.c \ + neg.c norm.c out_str.c pow.c pow_fr.c \ pow_ld.c pow_d.c pow_si.c pow_ui.c pow_z.c proj.c real.c rootofunity.c \ urandom.c set.c \ set_prec.c set_str.c set_x.c set_x_x.c sin.c sin_cos.c sinh.c sqr.c \ @@ -1,6 +1,6 @@ /* mpc_acos -- arccosine of a complex number. -Copyright (C) 2009, 2010, 2011 INRIA +Copyright (C) 2009, 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -67,7 +67,7 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd)); - mpfr_div_2ui (mpc_realref (rop), mpc_realref (rop), 1, GMP_RNDN); + mpfr_div_2ui (mpc_realref (rop), mpc_realref (rop), 1, MPFR_RNDN); } else { @@ -88,11 +88,11 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { p += mpc_ceil_log2 (p); mpfr_set_prec (x, p); - mpfr_const_pi (x, GMP_RNDD); - mpfr_mul_ui (x, x, 3, GMP_RNDD); + mpfr_const_pi (x, MPFR_RNDD); + mpfr_mul_ui (x, x, 3, MPFR_RNDD); ok = - mpfr_can_round (x, p - 1, GMP_RNDD, MPC_RND_RE (rnd), - prec+(MPC_RND_RE (rnd) == GMP_RNDN)); + mpfr_can_round (x, p - 1, MPFR_RNDD, MPC_RND_RE (rnd), + prec+(MPC_RND_RE (rnd) == MPFR_RNDN)); } while (ok == 0); inex_re = @@ -103,7 +103,7 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) else { if (mpfr_sgn (mpc_realref (op)) > 0) - mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN); else inex_re = mpfr_const_pi (mpc_realref (rop), MPC_RND_RE (rnd)); } @@ -131,7 +131,7 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) inex_im = -mpfr_acosh (mpc_imagref (rop), mpc_realref (op), INV_RND (MPC_RND_IM (rnd))); - mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN); } else if (mpfr_cmp_si (mpc_realref (op), -1) < 0) { @@ -180,13 +180,13 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* the imaginary part of asin(z) has the same sign as Im(z), thus if Im(z) > 0 and rnd_im = RNDZ, we want to round the Im(asin(z)) to -Inf so that -Im(asin(z)) is rounded to zero */ - if (rnd_im == GMP_RNDZ) - rnd_im = mpfr_sgn (mpc_imagref(op)) > 0 ? GMP_RNDD : GMP_RNDU; + if (rnd_im == MPFR_RNDZ) + rnd_im = mpfr_sgn (mpc_imagref(op)) > 0 ? MPFR_RNDD : MPFR_RNDU; else - rnd_im = rnd_im == GMP_RNDU ? GMP_RNDD - : rnd_im == GMP_RNDD ? GMP_RNDU + rnd_im = rnd_im == MPFR_RNDU ? MPFR_RNDD + : rnd_im == MPFR_RNDD ? MPFR_RNDU : rnd_im; /* both RNDZ and RNDA map to themselves for -asin(z) */ - rnd1 = RNDC(GMP_RNDN, rnd_im); + rnd1 = MPC_RND (MPFR_RNDN, rnd_im); mpfr_init2 (pi_over_2, p); for (;;) { @@ -195,14 +195,13 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_set_prec (mpc_realref(z1), p); mpfr_set_prec (pi_over_2, p); - mpfr_const_pi (pi_over_2, GMP_RNDN); - mpfr_div_2exp (pi_over_2, pi_over_2, 1, GMP_RNDN); /* Pi/2 */ + set_pi_over_2 (pi_over_2, +1, MPFR_RNDN); e1 = 1; /* Exp(pi_over_2) */ inex = mpc_asin (z1, op, rnd1); /* asin(z) */ MPC_ASSERT (mpfr_sgn (mpc_imagref(z1)) * mpfr_sgn (mpc_imagref(op)) > 0); inex_im = MPC_INEX_IM(inex); /* inex_im is in {-1, 0, 1} */ e2 = mpfr_get_exp (mpc_realref(z1)); - mpfr_sub (mpc_realref(z1), pi_over_2, mpc_realref(z1), GMP_RNDN); + mpfr_sub (mpc_realref(z1), pi_over_2, mpc_realref(z1), MPFR_RNDN); if (!mpfr_zero_p (mpc_realref(z1))) { /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-p-1) + @@ -213,10 +212,10 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* the error on x is bounded by 1/2 ulp(x) [1 + 2^e1] */ e1 = e1 <= 0 ? 0 : e1; /* the error on x is bounded by 2^e1 * ulp(x) */ - mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN); /* exact */ + mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), MPFR_RNDN); /* exact */ inex_im = -inex_im; - if (mpfr_can_round (mpc_realref(z1), p - e1, GMP_RNDN, GMP_RNDZ, - p_re + (MPC_RND_RE(rnd) == GMP_RNDN))) + if (mpfr_can_round (mpc_realref(z1), p - e1, MPFR_RNDN, MPFR_RNDZ, + p_re + (MPC_RND_RE(rnd) == MPFR_RNDN))) break; } } diff --git a/src/acosh.c b/src/acosh.c index a4eef4c..782f555 100644 --- a/src/acosh.c +++ b/src/acosh.c @@ -1,6 +1,6 @@ /* mpc_acosh -- inverse hyperbolic cosine of a complex number. -Copyright (C) 2009, 2011 INRIA +Copyright (C) 2009, 2011, 2012 INRIA This file is part of GNU MPC. @@ -46,7 +46,7 @@ mpc_acosh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) if (mpfr_signbit (mpc_imagref (op))) { inex = mpc_acos (a, op, - RNDC (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd))); + MPC_RND (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd))); /* change a to -i*a, i.e., -y+i*x to x+i*y */ tmp[0] = mpc_realref (a)[0]; @@ -58,7 +58,7 @@ mpc_acosh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) else { inex = mpc_acos (a, op, - RNDC (MPC_RND_IM (rnd), INV_RND(MPC_RND_RE (rnd)))); + MPC_RND (MPC_RND_IM (rnd), INV_RND(MPC_RND_RE (rnd)))); /* change a to i*a, i.e., y-i*x to x+i*y */ tmp[0] = mpc_realref (a)[0]; @@ -1,6 +1,6 @@ /* mpc_asin -- arcsine of a complex number. -Copyright (C) 2009, 2010, 2011 INRIA +Copyright (C) 2009, 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -38,7 +38,7 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) } else if (mpfr_zero_p (mpc_realref (op))) { - mpfr_set (mpc_realref (rop), mpc_realref (op), GMP_RNDN); + mpfr_set (mpc_realref (rop), mpc_realref (op), MPFR_RNDN); mpfr_set_nan (mpc_imagref (rop)); } else @@ -62,7 +62,7 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_set_inf (mpc_imagref (rop), (mpfr_signbit (mpc_imagref (op)) ? -1 : 1)); if (inf_im) - mpfr_div_2ui (mpc_realref (rop), mpc_realref (rop), 1, GMP_RNDN); + mpfr_div_2ui (mpc_realref (rop), mpc_realref (rop), 1, MPFR_RNDN); } else { @@ -116,7 +116,7 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { inex_im = mpfr_set_ui (mpc_imagref (rop), 0, MPC_RND_IM (rnd)); if (s_im) - mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN); + mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), MPFR_RNDN); inex_re = mpfr_asin (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); } @@ -129,9 +129,9 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) int inex_im; int s; s = mpfr_signbit (mpc_realref (op)); - mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN); if (s) - mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); + mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN); inex_im = mpfr_asinh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); return MPC_INEX (0, inex_im); @@ -158,8 +158,8 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* err(x) <= 1/2 ulp(x), err(y) <= 1/2 ulp(y) */ /* z1 <- 1-z1 */ ex = mpfr_get_exp (mpc_realref(z1)); - mpfr_ui_sub (mpc_realref(z1), 1, mpc_realref(z1), GMP_RNDN); - mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN); + mpfr_ui_sub (mpc_realref(z1), 1, mpc_realref(z1), MPFR_RNDN); + mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), MPFR_RNDN); ex = ex - mpfr_get_exp (mpc_realref(z1)); ex = (ex <= 0) ? 0 : ex; /* err(x) <= 2^ex * ulp(x) */ @@ -186,8 +186,8 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* z1 <- i*z + z1 */ ex = mpfr_get_exp (mpc_realref(z1)); ey = mpfr_get_exp (mpc_imagref(z1)); - mpfr_sub (mpc_realref(z1), mpc_realref(z1), mpc_imagref(op), GMP_RNDN); - mpfr_add (mpc_imagref(z1), mpc_imagref(z1), mpc_realref(op), GMP_RNDN); + mpfr_sub (mpc_realref(z1), mpc_realref(z1), mpc_imagref(op), MPFR_RNDN); + mpfr_add (mpc_imagref(z1), mpc_imagref(z1), mpc_realref(op), MPFR_RNDN); if (mpfr_cmp_ui (mpc_realref(z1), 0) == 0 || mpfr_cmp_ui (mpc_imagref(z1), 0) == 0) continue; ex -= mpfr_get_exp (mpc_realref(z1)); /* cancellation in x */ @@ -201,7 +201,7 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) ey = mpfr_get_exp (mpc_imagref(z1)); ex = (ex >= ey) ? ex : ey; err += ex - p; /* revert to absolute error <= 2^err */ - mpc_log (z1, z1, GMP_RNDN); + mpc_log (z1, z1, MPFR_RNDN); err -= ex - 1; /* 1/|t| <= 1/|z| <= 2^(1-ex) */ /* express err in terms of ulp(z1) */ ey = mpfr_get_exp (mpc_realref(z1)) <= mpfr_get_exp (mpc_imagref(z1)) @@ -211,11 +211,11 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) err = (err <= 0) ? 1 : err + 1; /* z1 <- -i*z1 */ mpfr_swap (mpc_realref(z1), mpc_imagref(z1)); - mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN); - if (mpfr_can_round (mpc_realref(z1), p - err, GMP_RNDN, GMP_RNDZ, - p_re + (rnd_re == GMP_RNDN)) && - mpfr_can_round (mpc_imagref(z1), p - err, GMP_RNDN, GMP_RNDZ, - p_im + (rnd_im == GMP_RNDN))) + mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), MPFR_RNDN); + if (mpfr_can_round (mpc_realref(z1), p - err, MPFR_RNDN, MPFR_RNDZ, + p_re + (rnd_re == MPFR_RNDN)) && + mpfr_can_round (mpc_imagref(z1), p - err, MPFR_RNDN, MPFR_RNDZ, + p_im + (rnd_im == MPFR_RNDN))) break; } diff --git a/src/asinh.c b/src/asinh.c index 89477a6..2807a8b 100644 --- a/src/asinh.c +++ b/src/asinh.c @@ -1,6 +1,6 @@ /* mpc_asinh -- inverse hyperbolic sine of a complex number. -Copyright (C) 2009, 2011 INRIA +Copyright (C) 2009, 2011, 2012 INRIA This file is part of GNU MPC. @@ -37,7 +37,7 @@ mpc_asinh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpc_init3 (a, MPC_PREC_IM(rop), MPC_PREC_RE(rop)); inex = mpc_asin (a, z, - RNDC (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd))); + MPC_RND (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd))); /* if a = asin(i*op) = x+i*y, and we want y-i*x */ @@ -32,11 +32,11 @@ set_pi_over_2 (mpfr_ptr rop, int s, mpfr_rnd_t rnd) int inex; inex = mpfr_const_pi (rop, s < 0 ? INV_RND (rnd) : rnd); - mpfr_div_2ui (rop, rop, 1, GMP_RNDN); + mpfr_div_2ui (rop, rop, 1, MPFR_RNDN); if (s < 0) { inex = -inex; - mpfr_neg (rop, rop, GMP_RNDN); + mpfr_neg (rop, rop, MPFR_RNDN); } return inex; @@ -64,7 +64,7 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_set_nan (mpc_realref (rop)); if (mpfr_zero_p (mpc_imagref (op)) || mpfr_inf_p (mpc_imagref (op))) { - mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref (rop), 0, MPFR_RNDN); if (s_im) mpc_conj (rop, rop, MPC_RNDNN); } @@ -76,7 +76,7 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) if (mpfr_inf_p (mpc_realref (op))) { inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd)); - mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref (rop), 0, MPFR_RNDN); } else { @@ -91,9 +91,9 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd)); - mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref (rop), 0, MPFR_RNDN); if (s_im) - mpc_conj (rop, rop, GMP_RNDN); + mpc_conj (rop, rop, MPFR_RNDN); return MPC_INEX (inex_re, 0); } @@ -103,9 +103,9 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { inex_re = mpfr_atan (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); - mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref (rop), 0, MPFR_RNDN); if (s_im) - mpc_conj (rop, rop, GMP_RNDN); + mpc_conj (rop, rop, MPFR_RNDN); return MPC_INEX (inex_re, 0); } @@ -125,9 +125,9 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* atan(+0+iy) = +0 +i*atanh(y), if |y| < 1 atan(-0+iy) = -0 +i*atanh(y), if |y| < 1 */ - mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN); if (s_re) - mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); + mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN); inex_im = mpfr_atanh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); } @@ -164,7 +164,7 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { p += mpc_ceil_log2 (p) + 2; mpfr_set_prec (y, p); - rnd_away = s_im == 0 ? GMP_RNDU : GMP_RNDD; + rnd_away = s_im == 0 ? MPFR_RNDU : MPFR_RNDD; inex_im = mpfr_ui_div (y, 1, mpc_imagref (op), rnd_away); /* FIXME: should we consider the case with unreasonably huge precision prec(y)>3*exp_min, where atanh(1/Im(op)) could be @@ -176,8 +176,8 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) inex_im |= mpfr_atanh (y, y, rnd_away); ok = inex_im == 0 - || mpfr_can_round (y, p - 2, rnd_away, GMP_RNDZ, - p_im + (rnd_im == GMP_RNDN)); + || mpfr_can_round (y, p - 2, rnd_away, MPFR_RNDZ, + p_im + (rnd_im == MPFR_RNDN)); } while (ok == 0); inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd)); @@ -196,7 +196,6 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_t minus_op_re; mpfr_exp_t op_re_exp, op_im_exp; mpfr_rnd_t rnd1, rnd2; - int loops = 0; mpfr_inits2 (MPFR_PREC_MIN, a, b, x, y, (mpfr_ptr) 0); @@ -227,8 +226,8 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* p: working precision */ p = (op_im_exp > 0 || prec > SAFE_ABS (mpfr_prec_t, op_im_exp)) ? prec : (prec - op_im_exp); - rnd1 = mpfr_sgn (mpc_realref (op)) > 0 ? GMP_RNDD : GMP_RNDU; - rnd2 = mpfr_sgn (mpc_realref (op)) < 0 ? GMP_RNDU : GMP_RNDD; + rnd1 = mpfr_sgn (mpc_realref (op)) > 0 ? MPFR_RNDD : MPFR_RNDU; + rnd2 = mpfr_sgn (mpc_realref (op)) < 0 ? MPFR_RNDU : MPFR_RNDD; do { @@ -250,7 +249,7 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) } else err = mpfr_get_exp (a); /* err = Exp(a) with the notations above */ - mpfr_atan2 (x, mpc_realref (op), a, GMP_RNDU); + mpfr_atan2 (x, mpc_realref (op), a, MPFR_RNDU); /* b = lower bound for atan (-x/(1+y)): for x negative, we need a lower bound on -x/(1+y), i.e., an upper bound on 1+y */ @@ -265,21 +264,21 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) } else expo = mpfr_get_exp (a); /* expo = Exp(c) with the notations above */ - mpfr_atan2 (b, minus_op_re, a, GMP_RNDD); + mpfr_atan2 (b, minus_op_re, a, MPFR_RNDD); err = err < expo ? err : expo; /* err = min(Exp(a),Exp(c)) */ - mpfr_sub (x, x, b, GMP_RNDU); + mpfr_sub (x, x, b, MPFR_RNDU); err = 5 + op_re_exp - err - mpfr_get_exp (x); /* error is bounded by [1 + 2^err] ulp(e) */ err = err < 0 ? 1 : err + 1; - mpfr_div_2ui (x, x, 1, GMP_RNDU); + mpfr_div_2ui (x, x, 1, MPFR_RNDU); /* Note: using RND2=RNDD guarantees that if x is exactly representable on prec + ... bits, mpfr_can_round will return 0 */ - ok = mpfr_can_round (x, p - err, GMP_RNDU, GMP_RNDD, - prec + (MPC_RND_RE (rnd) == GMP_RNDN)); + ok = mpfr_can_round (x, p - err, MPFR_RNDU, MPFR_RNDD, + prec + (MPC_RND_RE (rnd) == MPFR_RNDN)); } while (ok == 0); /* Imaginary part @@ -305,7 +304,6 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) */ err = 2; p = prec; /* working precision */ - rnd1 = mpfr_cmp_si (mpc_imagref (op), -1) > 0 ? GMP_RNDU : GMP_RNDD; do { @@ -315,23 +313,27 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_set_prec (y, p); /* a = upper bound for log(x^2 + (1+y)^2) */ - mpfr_add_ui (a, mpc_imagref (op), 1, rnd1); - mpfr_sqr (a, a, GMP_RNDU); - mpfr_sqr (y, mpc_realref (op), GMP_RNDU); - mpfr_add (a, a, y, GMP_RNDU); - mpfr_log (a, a, GMP_RNDU); + mpfr_add_ui (a, mpc_imagref (op), 1, MPFR_RNDA); + mpfr_sqr (a, a, MPFR_RNDU); + mpfr_sqr (y, mpc_realref (op), MPFR_RNDU); + mpfr_add (a, a, y, MPFR_RNDU); + mpfr_log (a, a, MPFR_RNDU); /* b = lower bound for log(x^2 + (1-y)^2) */ - mpfr_ui_sub (b, 1, mpc_imagref (op), GMP_RNDZ); - mpfr_sqr (b, b, GMP_RNDU); - /* mpfr_sqr (y, mpc_realref (op), GMP_RNDZ); */ + mpfr_ui_sub (b, 1, mpc_imagref (op), MPFR_RNDZ); /* round to zero */ + mpfr_sqr (b, b, MPFR_RNDZ); + /* we could write mpfr_sqr (y, mpc_realref (op), MPFR_RNDZ) but it is + more efficient to reuse the value of y (x^2) above and subtract + one ulp */ mpfr_nextbelow (y); - mpfr_add (b, b, y, GMP_RNDZ); - mpfr_log (b, b, GMP_RNDZ); + mpfr_add (b, b, y, MPFR_RNDZ); + mpfr_log (b, b, MPFR_RNDZ); - mpfr_sub (y, a, b, GMP_RNDU); + mpfr_sub (y, a, b, MPFR_RNDU); if (mpfr_zero_p (y)) + /* FIXME: happens when x and y have very different magnitudes; + could be handled more efficiently */ ok = 0; else { @@ -344,12 +346,16 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) else err = (expo < 0) ? 1 : expo + 2; - mpfr_div_2ui (y, y, 2, GMP_RNDN); - if (mpfr_zero_p (y)) - err = 2; /* underflow */ + mpfr_div_2ui (y, y, 2, MPFR_RNDN); + MPC_ASSERT (!mpfr_zero_p (y)); + /* FIXME: underflow. Since the main term of the Taylor series + in y=0 is 1/(x^2+1) * y, this means that y is very small + and/or x very large; but then the mpfr_zero_p (y) above + should be true. This needs a proof, or better yet, + special code. */ - ok = mpfr_can_round (y, p - err, GMP_RNDU, GMP_RNDD, - prec + (MPC_RND_IM (rnd) == GMP_RNDN)); + ok = mpfr_can_round (y, p - err, MPFR_RNDU, MPFR_RNDD, + prec + (MPC_RND_IM (rnd) == MPFR_RNDN)); } } while (ok == 0); diff --git a/src/atanh.c b/src/atanh.c index 3e1c448..e5716cc 100644 --- a/src/atanh.c +++ b/src/atanh.c @@ -1,6 +1,6 @@ /* mpc_atanh -- inverse hyperbolic tangent of a complex number. -Copyright (C) 2009, 2011 INRIA +Copyright (C) 2009, 2011, 2012 INRIA This file is part of GNU MPC. @@ -36,7 +36,7 @@ mpc_atanh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpc_init3 (a, MPC_PREC_IM(rop), MPC_PREC_RE(rop)); inex = mpc_atan (a, z, - RNDC (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd))); + MPC_RND (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd))); /* change a to -i*a, i.e., x+i*y to y-i*x */ tmp[0] = mpc_realref (a)[0]; @@ -68,28 +68,28 @@ mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w) if (a == 1) if (b == 1) { - mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); x = MPC_MPFR_SIGN (sign); - mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); y = MPC_MPFR_SIGN (sign); } else { /* b == -1 */ - mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); x = MPC_MPFR_SIGN (sign); - mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); y = -MPC_MPFR_SIGN (sign); } else /* a == -1 */ if (b == 1) { - mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN); + mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), MPFR_RNDN); x = MPC_MPFR_SIGN (sign); - mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); y = MPC_MPFR_SIGN (sign); } else { /* b == -1 */ - mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); x = -MPC_MPFR_SIGN (sign); - mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN); + mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), MPFR_RNDN); y = MPC_MPFR_SIGN (sign); } mpfr_clear (sign); @@ -121,25 +121,25 @@ mpc_div_fin_inf (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w) mpfr_init2 (x, 2); mpfr_init2 (y, 2); mpfr_init2 (zero, 2); - mpfr_set_ui (zero, 0ul, GMP_RNDN); + mpfr_set_ui (zero, 0ul, MPFR_RNDN); mpfr_init2 (a, mpfr_get_prec (mpc_realref (z))); mpfr_init2 (b, mpfr_get_prec (mpc_imagref (z))); - mpfr_set_ui (c, (mpfr_inf_p (mpc_realref (w)) ? 1 : 0), GMP_RNDN); - MPFR_COPYSIGN (c, c, mpc_realref (w), GMP_RNDN); - mpfr_set_ui (d, (mpfr_inf_p (mpc_imagref (w)) ? 1 : 0), GMP_RNDN); - MPFR_COPYSIGN (d, d, mpc_imagref (w), GMP_RNDN); + mpfr_set_ui (c, (mpfr_inf_p (mpc_realref (w)) ? 1 : 0), MPFR_RNDN); + MPFR_COPYSIGN (c, c, mpc_realref (w), MPFR_RNDN); + mpfr_set_ui (d, (mpfr_inf_p (mpc_imagref (w)) ? 1 : 0), MPFR_RNDN); + MPFR_COPYSIGN (d, d, mpc_imagref (w), MPFR_RNDN); - mpfr_mul (a, mpc_realref (z), c, GMP_RNDN); /* exact */ - mpfr_mul (b, mpc_imagref (z), d, GMP_RNDN); - mpfr_add (x, a, b, GMP_RNDN); + mpfr_mul (a, mpc_realref (z), c, MPFR_RNDN); /* exact */ + mpfr_mul (b, mpc_imagref (z), d, MPFR_RNDN); + mpfr_add (x, a, b, MPFR_RNDN); - mpfr_mul (b, mpc_imagref (z), c, GMP_RNDN); - mpfr_mul (a, mpc_realref (z), d, GMP_RNDN); - mpfr_sub (y, b, a, GMP_RNDN); + mpfr_mul (b, mpc_imagref (z), c, MPFR_RNDN); + mpfr_mul (a, mpc_realref (z), d, MPFR_RNDN); + mpfr_sub (y, b, a, MPFR_RNDN); - MPFR_COPYSIGN (mpc_realref (rop), zero, x, GMP_RNDN); - MPFR_COPYSIGN (mpc_imagref (rop), zero, y, GMP_RNDN); + MPFR_COPYSIGN (mpc_realref (rop), zero, x, MPFR_RNDN); + MPFR_COPYSIGN (mpc_imagref (rop), zero, y, MPFR_RNDN); mpfr_clear (c); mpfr_clear (d); @@ -173,10 +173,10 @@ mpc_div_real (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) inexact flags */ if (mpfr_zero_p (mpc_realref (rop))) mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis), - GMP_RNDN); /* exact */ + MPFR_RNDN); /* exact */ if (mpfr_zero_p (mpc_imagref (rop))) mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis), - GMP_RNDN); + MPFR_RNDN); return MPC_INEX(inex_re, inex_im); } @@ -204,7 +204,7 @@ mpc_div_imag (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) wloc[0] = mpc_imagref(w)[0]; /* copies mpfr struct IM(w) into wloc */ inex_re = mpfr_div (mpc_realref(dest), mpc_imagref(z), wloc, MPC_RND_RE(rnd)); - mpfr_neg (wloc, wloc, GMP_RNDN); + mpfr_neg (wloc, wloc, MPFR_RNDN); /* changes the sign only in wloc, not in w; no need to correct later */ inex_im = mpfr_div (mpc_imagref(dest), mpc_realref(z), wloc, MPC_RND_IM(rnd)); @@ -221,10 +221,10 @@ mpc_div_imag (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) inexact flags */ if (mpfr_zero_p (mpc_realref (rop))) mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis), - GMP_RNDN); /* exact */ + MPFR_RNDN); /* exact */ if (imag_z) mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis), - GMP_RNDN); + MPFR_RNDN); return MPC_INEX(inex_re, inex_im); } @@ -292,11 +292,11 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) /* first compute norm(c) */ mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_norm = mpc_norm (q, c, GMP_RNDU); + inexact_norm = mpc_norm (q, c, MPFR_RNDU); underflow_norm = mpfr_underflow_p (); overflow_norm = mpfr_overflow_p (); if (underflow_norm) - mpfr_set_ui (q, 0ul, GMP_RNDN); + mpfr_set_ui (q, 0ul, MPFR_RNDN); /* to obtain divisions by 0 later on */ /* now compute b*conjugate(c) */ @@ -312,69 +312,77 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) hopefully, the side-effects of mpc_mul do indeed raise the mpfr exceptions */ if (overflow_prod) { + int isinf = 0; tmpsgn = mpfr_sgn (mpc_realref(res)); if (tmpsgn > 0) - mpfr_nextabove (mpc_realref(res)); + { + mpfr_nextabove (mpc_realref(res)); + isinf = mpfr_inf_p (mpc_realref(res)); + mpfr_nextbelow (mpc_realref(res)); + } else if (tmpsgn < 0) - mpfr_nextbelow (mpc_realref(res)); - if (mpfr_inf_p (mpc_realref(res))) + { + mpfr_nextbelow (mpc_realref(res)); + isinf = mpfr_inf_p (mpc_realref(res)); + mpfr_nextabove (mpc_realref(res)); + } + if (isinf) { mpfr_set_inf (mpc_realref(res), tmpsgn); overflow_re = 1; } - else - if (tmpsgn > 0) - mpfr_nextbelow (mpc_realref(res)); - else if (tmpsgn < 0) - mpfr_nextabove (mpc_realref(res)); tmpsgn = mpfr_sgn (mpc_imagref(res)); + isinf = 0; if (tmpsgn > 0) - mpfr_nextabove (mpc_imagref(res)); + { + mpfr_nextabove (mpc_imagref(res)); + isinf = mpfr_inf_p (mpc_imagref(res)); + mpfr_nextbelow (mpc_imagref(res)); + } else if (tmpsgn < 0) - mpfr_nextbelow (mpc_imagref(res)); - if (mpfr_inf_p (mpc_imagref(res))) + { + mpfr_nextbelow (mpc_imagref(res)); + isinf = mpfr_inf_p (mpc_imagref(res)); + mpfr_nextabove (mpc_imagref(res)); + } + if (isinf) { mpfr_set_inf (mpc_imagref(res), tmpsgn); overflow_im = 1; } - else - if (tmpsgn > 0) - mpfr_nextbelow (mpc_imagref(res)); - else if (tmpsgn < 0) - mpfr_nextabove (mpc_imagref(res)); mpc_set (a, res, rnd); goto end; } /* divide the product by the norm */ if (inexact_norm == 0 && (inexact_re == 0 || inexact_im == 0)) { - /* The division has good chances to be exact in at least one part. */ - /* Since this can cause problems when not rounding to the nearest, */ - /* we use the division code of mpfr, which handles the situation. */ + /* The division has good chances to be exact in at least one part. */ + /* Since this can cause problems when not rounding to the nearest, */ + /* we use the division code of mpfr, which handles the situation. */ mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_re |= mpfr_div (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ); + inexact_re |= mpfr_div (mpc_realref (res), mpc_realref (res), q, MPFR_RNDZ); underflow_re = mpfr_underflow_p (); overflow_re = mpfr_overflow_p (); ok_re = !inexact_re || underflow_re || overflow_re - || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN, - GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN)); + || mpfr_can_round (mpc_realref (res), prec - 4, MPFR_RNDN, + MPFR_RNDZ, MPC_PREC_RE(a) + (rnd_re == MPFR_RNDN)); if (ok_re) /* compute imaginary part */ { mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_im |= mpfr_div (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ); + inexact_im |= mpfr_div (mpc_imagref (res), mpc_imagref (res), q, MPFR_RNDZ); underflow_im = mpfr_underflow_p (); overflow_im = mpfr_overflow_p (); ok_im = !inexact_im || underflow_im || overflow_im - || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN, - GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN)); + || mpfr_can_round (mpc_imagref (res), prec - 4, MPFR_RNDN, + MPFR_RNDZ, MPC_PREC_IM(a) + (rnd_im == MPFR_RNDN)); } } else { /* The division is inexact, so for efficiency reasons we invert q */ /* only once and multiply by the inverse. */ - if (mpfr_ui_div (q, 1ul, q, GMP_RNDZ) || inexact_norm) { + if (mpfr_ui_div (q, 1ul, q, MPFR_RNDZ) || inexact_norm) { /* if 1/q is inexact, the approximations of the real and imaginary part below will be inexact, unless RE(res) or IM(res) is zero */ @@ -383,22 +391,22 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) } mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_re |= mpfr_mul (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ); + inexact_re |= mpfr_mul (mpc_realref (res), mpc_realref (res), q, MPFR_RNDZ); underflow_re = mpfr_underflow_p (); overflow_re = mpfr_overflow_p (); ok_re = !inexact_re || underflow_re || overflow_re - || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN, - GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN)); + || mpfr_can_round (mpc_realref (res), prec - 4, MPFR_RNDN, + MPFR_RNDZ, MPC_PREC_RE(a) + (rnd_re == MPFR_RNDN)); if (ok_re) /* compute imaginary part */ { mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_im |= mpfr_mul (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ); + inexact_im |= mpfr_mul (mpc_imagref (res), mpc_imagref (res), q, MPFR_RNDZ); underflow_im = mpfr_underflow_p (); overflow_im = mpfr_overflow_p (); ok_im = !inexact_im || underflow_im || overflow_im - || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN, - GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN)); + || mpfr_can_round (mpc_imagref (res), prec - 4, MPFR_RNDN, + MPFR_RNDZ, MPC_PREC_IM(a) + (rnd_im == MPFR_RNDN)); } } } while ((!ok_re || !ok_im) && !underflow_norm && !overflow_norm @@ -416,16 +424,16 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) inexact_re = mpfr_sgn (mpc_realref (res)); } else if (underflow_re || (overflow_norm && !overflow_prod)) { - mpfr_set_zero (mpc_realref (a), mpfr_sgn (mpc_realref (res))); - inexact_re = -mpfr_sgn (mpc_realref (res)); + inexact_re = mpfr_signbit (mpc_realref (res)) ? 1 : -1; + mpfr_set_zero (mpc_realref (a), -inexact_re); } if (overflow_im || (underflow_norm && !underflow_prod)) { mpfr_set_inf (mpc_imagref (a), mpfr_sgn (mpc_imagref (res))); inexact_im = mpfr_sgn (mpc_imagref (res)); } else if (underflow_im || (overflow_norm && !overflow_prod)) { - mpfr_set_zero (mpc_imagref (a), mpfr_sgn (mpc_imagref (res))); - inexact_im = -mpfr_sgn (mpc_imagref (res)); + inexact_im = mpfr_signbit (mpc_imagref (res)) ? 1 : -1; + mpfr_set_zero (mpc_imagref (a), -inexact_im); } mpc_clear (res); diff --git a/src/div_2exp.c b/src/div_2si.c index 22e67e4..511f2cb 100644 --- a/src/div_2exp.c +++ b/src/div_2si.c @@ -1,6 +1,6 @@ -/* mpc_div_2exp -- Divide a complex number by 2^e. +/* mpc_div_2si -- Divide a complex number by 2^e. -Copyright (C) 2002, 2009, 2011 INRIA +Copyright (C) 2012 INRIA This file is part of GNU MPC. @@ -21,12 +21,12 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include "mpc-impl.h" int -mpc_div_2exp (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd) +mpc_div_2si (mpc_ptr a, mpc_srcptr b, long int c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_div_2exp (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_div_2exp (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); + inex_re = mpfr_div_2si (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_div_2si (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } diff --git a/src/mul_2exp.c b/src/div_2ui.c index ff2efe2..cd53855 100644 --- a/src/mul_2exp.c +++ b/src/div_2ui.c @@ -1,6 +1,6 @@ -/* mpc_mul_2exp -- Multiply a complex number by 2^e. +/* mpc_div_2ui -- Divide a complex number by 2^e. -Copyright (C) 2002, 2009, 2011 INRIA +Copyright (C) 2002, 2009, 2011, 2012 INRIA This file is part of GNU MPC. @@ -21,12 +21,12 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include "mpc-impl.h" int -mpc_mul_2exp (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd) +mpc_div_2ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_mul_2exp (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_mul_2exp (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); + inex_re = mpfr_div_2ui (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_div_2ui (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } diff --git a/src/div_fr.c b/src/div_fr.c index d5ea240..447c078 100644 --- a/src/div_fr.c +++ b/src/div_fr.c @@ -1,6 +1,6 @@ /* mpc_div_fr -- Divide a complex number by a floating-point number. -Copyright (C) 2002, 2008, 2009, 2010, 2011 INRIA +Copyright (C) 2002, 2008, 2009, 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -31,7 +31,7 @@ mpc_div_fr (mpc_ptr a, mpc_srcptr b, mpfr_srcptr c, mpc_rnd_t rnd) inex_re = mpfr_div (real, mpc_realref(b), c, MPC_RND_RE(rnd)); inex_im = mpfr_div (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); - mpfr_set (mpc_realref (a), real, GMP_RNDN); + mpfr_set (mpc_realref (a), real, MPFR_RNDN); mpfr_clear (real); @@ -1,6 +1,6 @@ /* mpc_exp -- exponential of a complex number. -Copyright (C) 2002, 2009, 2010, 2011 INRIA +Copyright (C) 2002, 2009, 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -88,15 +88,15 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_init2 (n, 2); if (mpfr_signbit (mpc_realref (op))) - mpfr_set_ui (n, 0, GMP_RNDN); + mpfr_set_ui (n, 0, MPFR_RNDN); else mpfr_set_inf (n, +1); if (mpfr_inf_p (mpc_imagref (op))) { - inex_re = mpfr_set (mpc_realref (rop), n, GMP_RNDN); + inex_re = mpfr_set (mpc_realref (rop), n, MPFR_RNDN); if (mpfr_signbit (mpc_realref (op))) - inex_im = mpfr_set (mpc_imagref (rop), n, GMP_RNDN); + inex_im = mpfr_set (mpc_imagref (rop), n, MPFR_RNDN); else { mpfr_set_nan (mpc_imagref (rop)); @@ -109,9 +109,9 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_init2 (c, 2); mpfr_init2 (s, 2); - mpfr_sin_cos (s, c, mpc_imagref (op), GMP_RNDN); - inex_re = mpfr_copysign (mpc_realref (rop), n, c, GMP_RNDN); - inex_im = mpfr_copysign (mpc_imagref (rop), n, s, GMP_RNDN); + mpfr_sin_cos (s, c, mpc_imagref (op), MPFR_RNDN); + inex_re = mpfr_copysign (mpc_realref (rop), n, c, MPFR_RNDN); + inex_im = mpfr_copysign (mpc_imagref (rop), n, s, MPFR_RNDN); mpfr_clear (s); mpfr_clear (c); @@ -159,18 +159,18 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) could be represented in the precision of rop. */ mpfr_clear_overflow (); mpfr_clear_underflow (); - mpfr_exp (x, mpc_realref(op), GMP_RNDN); /* error <= 0.5ulp */ - mpfr_sin_cos (z, y, mpc_imagref(op), GMP_RNDN); /* errors <= 0.5ulp */ - mpfr_mul (y, y, x, GMP_RNDN); /* error <= 2ulp */ + mpfr_exp (x, mpc_realref(op), MPFR_RNDN); /* error <= 0.5ulp */ + mpfr_sin_cos (z, y, mpc_imagref(op), MPFR_RNDN); /* errors <= 0.5ulp */ + mpfr_mul (y, y, x, MPFR_RNDN); /* error <= 2ulp */ ok = mpfr_overflow_p () || mpfr_zero_p (x) - || mpfr_can_round (y, prec - 2, GMP_RNDN, GMP_RNDZ, - MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == GMP_RNDN)); + || mpfr_can_round (y, prec - 2, MPFR_RNDN, MPFR_RNDZ, + MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == MPFR_RNDN)); if (ok) /* compute imaginary part */ { - mpfr_mul (z, z, x, GMP_RNDN); + mpfr_mul (z, z, x, MPFR_RNDN); ok = mpfr_overflow_p () || mpfr_zero_p (x) - || mpfr_can_round (z, prec - 2, GMP_RNDN, GMP_RNDZ, - MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == GMP_RNDN)); + || mpfr_can_round (z, prec - 2, MPFR_RNDN, MPFR_RNDZ, + MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == MPFR_RNDN)); } } while (ok == 0); @@ -51,10 +51,10 @@ mpc_fma_naive (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rn mpfr_init2 (ima_reb, mpfr_get_prec (mpc_imagref(a)) + mpfr_get_prec (mpc_realref(b))); mpfr_init2 (ima_imb, mpfr_get_prec (mpc_imagref(a)) + mpfr_get_prec (mpc_imagref(b))); - mpfr_mul (rea_reb, mpc_realref(a), mpc_realref(b), GMP_RNDZ); /* exact */ - mpfr_mul (rea_imb, mpc_realref(a), mpc_imagref(b), GMP_RNDZ); /* exact */ - mpfr_mul (ima_reb, mpc_imagref(a), mpc_realref(b), GMP_RNDZ); /* exact */ - mpfr_mul (ima_imb, mpc_imagref(a), mpc_imagref(b), GMP_RNDZ); /* exact */ + mpfr_mul (rea_reb, mpc_realref(a), mpc_realref(b), MPFR_RNDZ); /* exact */ + mpfr_mul (rea_imb, mpc_realref(a), mpc_imagref(b), MPFR_RNDZ); /* exact */ + mpfr_mul (ima_reb, mpc_imagref(a), mpc_realref(b), MPFR_RNDZ); /* exact */ + mpfr_mul (ima_imb, mpc_imagref(a), mpc_imagref(b), MPFR_RNDZ); /* exact */ /* Re(r) <- rea_reb - ima_imb + Re(c) */ @@ -67,7 +67,7 @@ mpc_fma_naive (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rn if (pre12 <= pre13 && pre12 <= pre23) /* (rea_reb - ima_imb) + Re(c) */ { mpfr_init2 (tmp, pre12); - mpfr_sub (tmp, rea_reb, ima_imb, GMP_RNDZ); /* exact */ + mpfr_sub (tmp, rea_reb, ima_imb, MPFR_RNDZ); /* exact */ inex_re = mpfr_add (mpc_realref(r), tmp, mpc_realref(c), MPC_RND_RE(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the real part of both, it is ok */ @@ -75,7 +75,7 @@ mpc_fma_naive (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rn else if (pre13 <= pre23) /* (rea_reb + Re(c)) - ima_imb */ { mpfr_init2 (tmp, pre13); - mpfr_add (tmp, rea_reb, mpc_realref(c), GMP_RNDZ); /* exact */ + mpfr_add (tmp, rea_reb, mpc_realref(c), MPFR_RNDZ); /* exact */ inex_re = mpfr_sub (mpc_realref(r), tmp, ima_imb, MPC_RND_RE(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the real part of both, it is ok */ @@ -83,7 +83,7 @@ mpc_fma_naive (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rn else /* rea_reb + (Re(c) - ima_imb) */ { mpfr_init2 (tmp, pre23); - mpfr_sub (tmp, mpc_realref(c), ima_imb, GMP_RNDZ); /* exact */ + mpfr_sub (tmp, mpc_realref(c), ima_imb, MPFR_RNDZ); /* exact */ inex_re = mpfr_add (mpc_realref(r), tmp, rea_reb, MPC_RND_RE(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the real part of both, it is ok */ @@ -99,7 +99,7 @@ mpc_fma_naive (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rn if (pim12 <= pim13 && pim12 <= pim23) /* (rea_imb + ima_reb) + Im(c) */ { mpfr_set_prec (tmp, pim12); - mpfr_add (tmp, rea_imb, ima_reb, GMP_RNDZ); /* exact */ + mpfr_add (tmp, rea_imb, ima_reb, MPFR_RNDZ); /* exact */ inex_im = mpfr_add (mpc_imagref(r), tmp, mpc_imagref(c), MPC_RND_IM(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the imaginary part of both, it is ok */ @@ -107,7 +107,7 @@ mpc_fma_naive (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rn else if (pim13 <= pim23) /* (rea_imb + Im(c)) + ima_reb */ { mpfr_set_prec (tmp, pim13); - mpfr_add (tmp, rea_imb, mpc_imagref(c), GMP_RNDZ); /* exact */ + mpfr_add (tmp, rea_imb, mpc_imagref(c), MPFR_RNDZ); /* exact */ inex_im = mpfr_add (mpc_imagref(r), tmp, ima_reb, MPC_RND_IM(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the imaginary part of both, it is ok */ @@ -115,7 +115,7 @@ mpc_fma_naive (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rn else /* rea_imb + (Im(c) + ima_reb) */ { mpfr_set_prec (tmp, pre23); - mpfr_add (tmp, mpc_imagref(c), ima_reb, GMP_RNDZ); /* exact */ + mpfr_add (tmp, mpc_imagref(c), ima_reb, MPFR_RNDZ); /* exact */ inex_im = mpfr_add (mpc_imagref(r), tmp, rea_imb, MPC_RND_IM(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the imaginary part of both, it is ok */ @@ -159,19 +159,19 @@ mpc_fma (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) break; diffre = mpfr_get_exp (mpc_realref(ab)); diffim = mpfr_get_exp (mpc_imagref(ab)); + mpc_add (ab, ab, c, MPC_RNDZZ); if (mpfr_zero_p (mpc_realref(ab)) || mpfr_zero_p (mpc_imagref(ab))) break; - mpc_add (ab, ab, c, MPC_RNDZZ); diffre -= mpfr_get_exp (mpc_realref(ab)); diffim -= mpfr_get_exp (mpc_imagref(ab)); diffre = (diffre > 0 ? diffre + 1 : 1); diffim = (diffim > 0 ? diffim + 1 : 1); - okre = diffre > wpre ? 0 : mpfr_can_round (mpc_realref(ab), - wpre - diffre, GMP_RNDN, GMP_RNDZ, - pre + (MPC_RND_RE (rnd) == GMP_RNDN)); - okim = diffim > wpim ? 0 : mpfr_can_round (mpc_imagref(ab), - wpim - diffim, GMP_RNDN, GMP_RNDZ, - pim + (MPC_RND_IM (rnd) == GMP_RNDN)); + okre = diffre > (mpfr_exp_t) wpre ? 0 : mpfr_can_round (mpc_realref(ab), + wpre - diffre, MPFR_RNDN, MPFR_RNDZ, + pre + (MPC_RND_RE (rnd) == MPFR_RNDN)); + okim = diffim > (mpfr_exp_t) wpim ? 0 : mpfr_can_round (mpc_imagref(ab), + wpim - diffim, MPFR_RNDN, MPFR_RNDZ, + pim + (MPC_RND_IM (rnd) == MPFR_RNDN)); if (okre && okim) { inex = mpc_set (r, ab, rnd); diff --git a/src/fr_div.c b/src/fr_div.c index e57eced..4c0536f 100644 --- a/src/fr_div.c +++ b/src/fr_div.c @@ -1,6 +1,6 @@ /* mpc_fr_div -- Divide a floating-point number by a complex number. -Copyright (C) 2008, 2009, 2011 INRIA +Copyright (C) 2008, 2009, 2011, 2012 INRIA This file is part of GNU MPC. @@ -29,7 +29,7 @@ mpc_fr_div (mpc_ptr a, mpfr_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) mpc_realref (bc)[0] = b [0]; mpfr_init (mpc_imagref (bc)); /* we consider the operand b to have imaginary part +0 */ - mpfr_set_ui (mpc_imagref (bc), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref (bc), 0, MPFR_RNDN); inexact = mpc_div (a, bc, c, rnd); diff --git a/src/get_version.c b/src/get_version.c index 1c429c0..e3b95a1 100644 --- a/src/get_version.c +++ b/src/get_version.c @@ -1,6 +1,6 @@ /* mpc_get_version -- MPC version -Copyright (C) 2008, 2009, 2010, 2011 INRIA +Copyright (C) 2008, 2009, 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -20,29 +20,8 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include "mpc-impl.h" -#if MPFR_VERSION_MAJOR < 3 -/* The following are functions defined for compatibility with mpfr < 3; - logically, they should be defined in a separate file, but then gcc - complains about an empty translation unit with mpfr >= 3. */ - -void -mpfr_set_zero (mpfr_ptr z, int s) -{ - mpfr_set_ui (z, 0ul, GMP_RNDN); - if (s < 0) - mpfr_neg (z, z, GMP_RNDN); -} - -int -mpfr_regular_p (mpfr_srcptr z) -{ - return (mpfr_number_p (z) && !mpfr_zero_p (z)); -} -#endif /* mpfr < 3 */ - - const char * mpc_get_version (void) { - return "1.0.0dev"; + return "1.1dev"; } @@ -23,12 +23,16 @@ along with this program. If not, see http://www.gnu.org/licenses/ . int mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){ - int ok=0; - mpfr_t w; + int ok, underflow = 0; + mpfr_srcptr x, y; + mpfr_t v, w; mpfr_prec_t prec; - int loops = 0; + int loops; int re_cmp, im_cmp; int inex_re, inex_im; + int err; + mpfr_exp_t expw; + int sgnw; /* special values: NaN and infinities */ if (!mpc_fin_p (op)) { @@ -96,7 +100,7 @@ mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){ inex_re = mpfr_log (mpc_realref (rop), mpc_imagref (op), MPC_RND_RE (rnd)); inex_im = mpfr_const_pi (mpc_imagref (rop), MPC_RND_IM (rnd)); /* division by 2 does not change the ternary flag */ - mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); + mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, MPFR_RNDN); } else { w [0] = *mpc_imagref (op); @@ -104,45 +108,110 @@ mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){ inex_re = mpfr_log (mpc_realref (rop), w, MPC_RND_RE (rnd)); inex_im = mpfr_const_pi (mpc_imagref (rop), INV_RND (MPC_RND_IM (rnd))); /* division by 2 does not change the ternary flag */ - mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); - mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN); + 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 */ } return MPC_INEX(inex_re, inex_im); } prec = MPC_PREC_RE(rop); - mpfr_init2 (w, prec); - /* let op = x + iy; log = 1/2 log (x^2 + y^2) + i atan2 (y, x) */ - /* loop for the real part: log (x^2 + y^2) */ - do { - loops ++; - prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2; + mpfr_init2 (w, 2); + /* let op = x + iy; log = 1/2 log (x^2 + y^2) + i atan2 (y, x) */ + /* loop for the real part: 1/2 log (x^2 + y^2), fast, but unsafe */ + /* implementation */ + ok = 0; + for (loops = 1; !ok && loops <= 2; loops++) { + prec += mpc_ceil_log2 (prec) + 4; mpfr_set_prec (w, prec); - /* w is rounded down */ - mpc_norm (w, op, GMP_RNDD); - /* error 1 ulp */ - MPC_ASSERT (!mpfr_inf_p (w)); - /* FIXME: intermediate overflow; the logarithm may be representable */ - - mpfr_log (w, w, GMP_RNDD); - /* generic error of log: (2^(2 - exp(w)) + 1) ulp */ - - if (mpfr_get_exp (w) >= 2) - ok = mpfr_can_round (w, prec - 2, GMP_RNDD, - MPC_RND_RE(rnd), MPC_PREC_RE(rop)); - else - ok = mpfr_can_round (w, prec - 3 + mpfr_get_exp (w), GMP_RNDD, - MPC_RND_RE(rnd), MPC_PREC_RE(rop)); - } while (ok == 0); + mpc_abs (w, op, MPFR_RNDN); + /* error 0.5 ulp */ + if (mpfr_inf_p (w)) + /* intermediate overflow; the logarithm may be representable. + Intermediate underflow is impossible. */ + break; + + mpfr_log (w, w, MPFR_RNDN); + /* generic error of log: (2^(- exp(w)) + 0.5) ulp */ + + if (mpfr_zero_p (w)) + /* impossible to round, switch to second algorithm */ + break; + + err = MPC_MAX (-mpfr_get_exp (w), 0) + 1; + /* number of lost digits */ + ok = mpfr_can_round (w, prec - err, MPFR_RNDN, MPFR_RNDZ, + mpfr_get_prec (mpc_realref (rop)) + (MPC_RND_RE (rnd) == MPFR_RNDN)); + } + + if (!ok) { + prec = MPC_PREC_RE(rop); + mpfr_init2 (v, 2); + /* compute 1/2 log (x^2 + y^2) = log |x| + 1/2 * log (1 + (y/x)^2) + if |x| >= |y|; otherwise, exchange x and y */ + if (mpfr_cmpabs (mpc_realref (op), mpc_imagref (op)) >= 0) { + x = mpc_realref (op); + y = mpc_imagref (op); + } + else { + x = mpc_imagref (op); + y = mpc_realref (op); + } + + do { + prec += mpc_ceil_log2 (prec) + 4; + mpfr_set_prec (v, prec); + mpfr_set_prec (w, prec); + + mpfr_div (v, y, x, MPFR_RNDD); /* error 1 ulp */ + mpfr_sqr (v, v, MPFR_RNDD); + /* generic error of multiplication: + 1 + 2*1*(2+1*2^(1-prec)) <= 5.0625 since prec >= 6 */ + mpfr_log1p (v, v, MPFR_RNDD); + /* error 1 + 4*5.0625 = 21.25 , see algorithms.tex */ + mpfr_div_2ui (v, v, 1, MPFR_RNDD); + /* If the result is 0, then there has been an underflow somewhere. */ + + mpfr_abs (w, x, MPFR_RNDN); /* exact */ + mpfr_log (w, w, MPFR_RNDN); /* error 0.5 ulp */ + expw = mpfr_get_exp (w); + sgnw = mpfr_signbit (w); + + mpfr_add (w, w, v, MPFR_RNDN); + if (!sgnw) /* v is positive, so no cancellation; + error 22.25 ulp; error counts lost bits */ + err = 5; + else + err = MPC_MAX (5 + mpfr_get_exp (v), + /* 21.25 ulp (v) rewritten in ulp (result, now in w) */ + -1 + expw - mpfr_get_exp (w) + /* 0.5 ulp (previous w), rewritten in ulp (result) */ + ) + 2; + + /* handle one special case: |x|=1, and (y/x)^2 underflows; + then 1/2*log(x^2+y^2) \approx 1/2*y^2 also underflows. */ + if ( (mpfr_cmp_si (x, -1) == 0 || mpfr_cmp_ui (x, 1) == 0) + && mpfr_zero_p (w)) + underflow = 1; + + } while (!underflow && + !mpfr_can_round (w, prec - err, MPFR_RNDN, MPFR_RNDZ, + mpfr_get_prec (mpc_realref (rop)) + (MPC_RND_RE (rnd) == MPFR_RNDN))); + mpfr_clear (v); + } /* imaginary part */ inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op), MPC_RND_IM (rnd)); - /* set the real part; cannot be done before when rop==op */ - inex_re = mpfr_div_2ui (mpc_realref(rop), w, 1ul, MPC_RND_RE (rnd)); + /* set the real part; cannot be done before if rop==op */ + if (underflow) + /* create underflow in result */ + inex_re = mpfr_set_ui_2exp (mpc_realref (rop), 1, + mpfr_get_emin_min () - 2, MPC_RND_RE (rnd)); + else + inex_re = mpfr_set (mpc_realref (rop), w, MPC_RND_RE (rnd)); mpfr_clear (w); return MPC_INEX(inex_re, inex_im); } diff --git a/src/log10.c b/src/log10.c index e5972cc..6b7e687 100644 --- a/src/log10.c +++ b/src/log10.c @@ -39,8 +39,8 @@ mpc_log10_aux (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, int flag, int nb) mpfr_init2 (log10, prec); while (ok == 0) { - mpfr_set_ui (log10, 10, GMP_RNDN); /* exact since prec >= 4 */ - mpfr_log (log10, log10, GMP_RNDN); + 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. */ @@ -49,40 +49,40 @@ mpc_log10_aux (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, int flag, int nb) case 0: /* imag <- atan2(y/x) */ mpfr_atan2 (mpc_imagref (tmp), mpc_imagref (op), mpc_realref (op), MPC_RND_IM (rnd)); - mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, GMP_RNDN); - ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, GMP_RNDN, - GMP_RNDZ, MPC_PREC_IM(rop) + - (MPC_RND_IM (rnd) == GMP_RNDN)); + 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)); break; case 1: /* real <- log(x) */ mpfr_log (mpc_realref (tmp), mpc_realref (op), MPC_RND_RE (rnd)); - mpfr_div (mpc_realref (tmp), mpc_realref (tmp), log10, GMP_RNDN); - ok = mpfr_can_round (mpc_realref (tmp), prec - 2, GMP_RNDN, - GMP_RNDZ, MPC_PREC_RE(rop) + - (MPC_RND_RE (rnd) == GMP_RNDN)); + mpfr_div (mpc_realref (tmp), mpc_realref (tmp), log10, MPFR_RNDN); + ok = mpfr_can_round (mpc_realref (tmp), prec - 2, MPFR_RNDN, + MPFR_RNDZ, MPC_PREC_RE(rop) + + (MPC_RND_RE (rnd) == MPFR_RNDN)); if (ok) ret = mpfr_set (mpc_realref (rop), mpc_realref (tmp), MPC_RND_RE (rnd)); break; case 2: /* imag <- pi */ mpfr_const_pi (mpc_imagref (tmp), MPC_RND_IM (rnd)); - mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, GMP_RNDN); - ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, GMP_RNDN, - GMP_RNDZ, MPC_PREC_IM(rop) + - (MPC_RND_IM (rnd) == GMP_RNDN)); + 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)); break; case 3: /* real <- log(y) */ mpfr_log (mpc_realref (tmp), mpc_imagref (op), MPC_RND_RE (rnd)); - mpfr_div (mpc_realref (tmp), mpc_realref (tmp), log10, GMP_RNDN); - ok = mpfr_can_round (mpc_realref (tmp), prec - 2, GMP_RNDN, - GMP_RNDZ, MPC_PREC_RE(rop) + - (MPC_RND_RE (rnd) == GMP_RNDN)); + mpfr_div (mpc_realref (tmp), mpc_realref (tmp), log10, MPFR_RNDN); + ok = mpfr_can_round (mpc_realref (tmp), prec - 2, MPFR_RNDN, + MPFR_RNDZ, MPC_PREC_RE(rop) + + (MPC_RND_RE (rnd) == MPFR_RNDN)); if (ok) ret = mpfr_set (mpc_realref (rop), mpc_realref (tmp), MPC_RND_RE (rnd)); @@ -184,7 +184,7 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) inex_re = mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd)); else inex_re = mpc_log10_aux (rop, ww, rnd, 0, 1); - inex_im = mpc_log10_aux (rop, op, RNDC(0,rnd_im), 1, 2); + inex_im = mpc_log10_aux (rop, op, MPC_RND (0,rnd_im), 1, 2); if (negative_zero) { mpc_conj (rop, rop, MPC_RNDNN); @@ -200,7 +200,7 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) inex_re = mpc_log10_aux (rop, op, rnd, 0, 3); inex_im = mpc_log10_aux (rop, op, rnd, 1, 2); /* division by 2 does not change the ternary flag */ - mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); + mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, MPFR_RNDN); } else { @@ -208,11 +208,11 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) ww->im[0] = *mpc_imagref (op); MPFR_CHANGE_SIGN (ww->im); inex_re = mpc_log10_aux (rop, ww, rnd, 0, 3); - invrnd = RNDC(0, INV_RND (MPC_RND_IM (rnd))); + invrnd = MPC_RND (0, INV_RND (MPC_RND_IM (rnd))); inex_im = mpc_log10_aux (rop, op, invrnd, 1, 2); /* division by 2 does not change the ternary flag */ - mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); - mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN); + 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 */ } return MPC_INEX(inex_re, inex_im); @@ -231,12 +231,12 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpc_set_prec (ww, prec); mpc_log (ww, op, MPC_RNDNN); - mpfr_set_ui (w, 10, GMP_RNDN); /* exact since prec >= 4 */ - mpfr_log (w, w, GMP_RNDN); + 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, GMP_RNDN, GMP_RNDZ, - MPC_PREC_RE(rop) + (MPC_RND_RE (rnd) == GMP_RNDN)); + 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 @@ -254,39 +254,39 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_integer_p (mpc_imagref (op))) { mpz_t x, y; - unsigned long s; + unsigned long s, v; mpz_init (x); mpz_init (y); - mpfr_get_z (x, mpc_realref (op), GMP_RNDN); /* exact */ - mpfr_get_z (y, mpc_imagref (op), GMP_RNDN); /* exact */ + 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 */ - s = mpz_sizeinbase (x, 10); /* always exact or 1 too big */ - /* since s is the number of digits of x in base 10 (or one more), - we subtract 1 since we want to check whether x = 10^s */ - s --; - mpz_ui_pow_ui (y, 10, s); - if (mpz_cmp (y, x) > 0) /* might be 1 too big */ + 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_divexact_ui (y, y, 10); - s --; - } - if (mpz_cmp (y, x) == 0) /* Re(log10(x+i*y)) is exactly s/2 */ - { - /* we reset the precision of Re(ww) so that s can be - represented exactly */ - mpfr_set_prec (mpc_realref (ww), sizeof(unsigned long)*CHAR_BIT); - mpfr_set_ui_2exp (mpc_realref (ww), s, -1, GMP_RNDN); /* exact */ - ok = 1; + 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, GMP_RNDN, GMP_RNDZ, - MPC_PREC_IM(rop) + (MPC_RND_IM (rnd) == GMP_RNDN)); + 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)); diff --git a/src/mpc-impl.h b/src/mpc-impl.h index c3ca9af..5664955 100644 --- a/src/mpc-impl.h +++ b/src/mpc-impl.h @@ -55,46 +55,18 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #define MPFR_SIGNBIT(x) (mpfr_signbit (x) ? -1 : 1) #define MPC_MPFR_SIGN(x) (mpfr_zero_p (x) ? 0 : MPFR_SIGNBIT (x)) /* should be called MPFR_SIGN, but this is taken in mpfr.h */ -#define MPFR_CHANGE_SIGN(x) mpfr_neg(x,x,GMP_RNDN) +#define MPFR_CHANGE_SIGN(x) mpfr_neg(x,x,MPFR_RNDN) #define MPFR_COPYSIGN(x,y,z,rnd) (mpfr_nan_p (z) ? \ mpfr_setsign (x, y, 0, rnd) : \ mpfr_copysign (x, y, z, rnd)) /* work around spurious signs in nan */ -#define MPFR_ADD_ONE_ULP(x) mpfr_add_one_ulp (x, GMP_RNDN) -#define MPFR_SUB_ONE_ULP(x) mpfr_sub_one_ulp (x, GMP_RNDN) +#define MPFR_ADD_ONE_ULP(x) mpfr_add_one_ulp (x, MPFR_RNDN) +#define MPFR_SUB_ONE_ULP(x) mpfr_sub_one_ulp (x, MPFR_RNDN) /* drop unused rounding mode from macroes */ #define MPFR_SWAP(a,b) do { mpfr_srcptr tmp; tmp = a; a = b; b = tmp; } while (0) /* - * Macro implementing rounding away from zero, to ease compatibility with - * mpfr < 3. f is the complete function call with a rounding mode of - * MPFR_RNDA, rop the name of the variable containing the result; it is - * already contained in f, but needs to be repeated so that the macro can - * modify the variable. - * Usage: replace each call to a function such as - * mpfr_add (rop, a, b, MPFR_RNDA) - * by - * ROUND_AWAY (mpfr_add (rop, a, b, MPFR_RNDA), rop) -*/ -#if MPFR_VERSION_MAJOR < 3 - /* round towards zero, add 1 ulp if not exact */ -#define MPFR_RNDA GMP_RNDZ -#define ROUND_AWAY(f,rop) \ - ((f) ? MPFR_ADD_ONE_ULP (rop), MPFR_SIGNBIT (rop) : 0) -#else -#define ROUND_AWAY(f,rop) \ - (f) -#endif /* mpfr < 3 */ - -#if MPFR_VERSION_MAJOR < 3 -/* declare missing functions, defined in get_version.c */ -__MPC_DECLSPEC void mpfr_set_zero (mpfr_ptr, int); -__MPC_DECLSPEC int mpfr_regular_p (mpfr_srcptr); -#endif /* mpfr < 3 */ - - -/* * MPC macros */ @@ -103,7 +75,7 @@ __MPC_DECLSPEC int mpfr_regular_p (mpfr_srcptr); #define MPC_MAX_PREC(x) MPC_MAX(MPC_PREC_RE(x), MPC_PREC_IM(x)) #define INV_RND(r) \ - (((r) == GMP_RNDU) ? GMP_RNDD : (((r) == GMP_RNDD) ? GMP_RNDU : (r))) + (((r) == MPFR_RNDU) ? MPFR_RNDD : (((r) == MPFR_RNDD) ? MPFR_RNDU : (r))) #define mpc_inf_p(z) (mpfr_inf_p(mpc_realref(z))||mpfr_inf_p(mpc_imagref(z))) /* Convention in C99 (G.3): z is regarded as an infinity if at least one of @@ -137,6 +109,27 @@ __MPC_DECLSPEC int mpfr_regular_p (mpfr_srcptr); } while (0) #endif + +/* + * Debug macros + */ + +#define MPC_OUT(x) \ +do { \ + printf (#x "[%lu,%lu]=", (unsigned long int) MPC_PREC_RE (x), \ + (unsigned long int) MPC_PREC_IM (x)); \ + mpc_out_str (stdout, 2, 0, x, MPC_RNDNN); \ + printf ("\n"); \ +} while (0) + +#define MPFR_OUT(x) \ +do { \ + printf (#x "[%lu]=", (unsigned long int) mpfr_get_prec (x)); \ + mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); \ + printf ("\n"); \ +} while (0) + + /* * Constants */ @@ -24,16 +24,11 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include "gmp.h" #include "mpfr.h" -/* Backwards compatibility with mpfr<3.0.0 */ -#ifndef mpfr_exp_t -#define mpfr_exp_t mp_exp_t -#endif - /* Define MPC version number */ #define MPC_VERSION_MAJOR 1 -#define MPC_VERSION_MINOR 0 +#define MPC_VERSION_MINOR 1 #define MPC_VERSION_PATCHLEVEL 0 -#define MPC_VERSION_STRING "1.0.0dev" +#define MPC_VERSION_STRING "1.1dev" /* Macros dealing with MPC VERSION */ #define MPC_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c)) @@ -74,29 +69,29 @@ along with this program. If not, see http://www.gnu.org/licenses/ . we reserve four bits for a real rounding mode. */ typedef int mpc_rnd_t; -#define RNDC(r1,r2) (((int)(r1)) + ((int)(r2) << 4)) +#define MPC_RND(r1,r2) (((int)(r1)) + ((int)(r2) << 4)) #define MPC_RND_RE(x) ((mpfr_rnd_t)((x) & 0x0F)) #define MPC_RND_IM(x) ((mpfr_rnd_t)((x) >> 4)) -#define MPC_RNDNN RNDC(GMP_RNDN,GMP_RNDN) -#define MPC_RNDNZ RNDC(GMP_RNDN,GMP_RNDZ) -#define MPC_RNDNU RNDC(GMP_RNDN,GMP_RNDU) -#define MPC_RNDND RNDC(GMP_RNDN,GMP_RNDD) +#define MPC_RNDNN MPC_RND (MPFR_RNDN,MPFR_RNDN) +#define MPC_RNDNZ MPC_RND (MPFR_RNDN,MPFR_RNDZ) +#define MPC_RNDNU MPC_RND (MPFR_RNDN,MPFR_RNDU) +#define MPC_RNDND MPC_RND (MPFR_RNDN,MPFR_RNDD) -#define MPC_RNDZN RNDC(GMP_RNDZ,GMP_RNDN) -#define MPC_RNDZZ RNDC(GMP_RNDZ,GMP_RNDZ) -#define MPC_RNDZU RNDC(GMP_RNDZ,GMP_RNDU) -#define MPC_RNDZD RNDC(GMP_RNDZ,GMP_RNDD) +#define MPC_RNDZN MPC_RND (MPFR_RNDZ,MPFR_RNDN) +#define MPC_RNDZZ MPC_RND (MPFR_RNDZ,MPFR_RNDZ) +#define MPC_RNDZU MPC_RND (MPFR_RNDZ,MPFR_RNDU) +#define MPC_RNDZD MPC_RND (MPFR_RNDZ,MPFR_RNDD) -#define MPC_RNDUN RNDC(GMP_RNDU,GMP_RNDN) -#define MPC_RNDUZ RNDC(GMP_RNDU,GMP_RNDZ) -#define MPC_RNDUU RNDC(GMP_RNDU,GMP_RNDU) -#define MPC_RNDUD RNDC(GMP_RNDU,GMP_RNDD) +#define MPC_RNDUN MPC_RND (MPFR_RNDU,MPFR_RNDN) +#define MPC_RNDUZ MPC_RND (MPFR_RNDU,MPFR_RNDZ) +#define MPC_RNDUU MPC_RND (MPFR_RNDU,MPFR_RNDU) +#define MPC_RNDUD MPC_RND (MPFR_RNDU,MPFR_RNDD) -#define MPC_RNDDN RNDC(GMP_RNDD,GMP_RNDN) -#define MPC_RNDDZ RNDC(GMP_RNDD,GMP_RNDZ) -#define MPC_RNDDU RNDC(GMP_RNDD,GMP_RNDU) -#define MPC_RNDDD RNDC(GMP_RNDD,GMP_RNDD) +#define MPC_RNDDN MPC_RND (MPFR_RNDD,MPFR_RNDN) +#define MPC_RNDDZ MPC_RND (MPFR_RNDD,MPFR_RNDZ) +#define MPC_RNDDU MPC_RND (MPFR_RNDD,MPFR_RNDU) +#define MPC_RNDDD MPC_RND (MPFR_RNDD,MPFR_RNDD) /* Definitions of types and their semantics */ @@ -151,8 +146,10 @@ __MPC_DECLSPEC int mpc_div_fr (mpc_ptr, mpc_srcptr, mpfr_srcptr, mpc_rnd_t); __MPC_DECLSPEC int mpc_fr_div (mpc_ptr, mpfr_srcptr, mpc_srcptr, mpc_rnd_t); __MPC_DECLSPEC int mpc_div_ui (mpc_ptr, mpc_srcptr, unsigned long int, mpc_rnd_t); __MPC_DECLSPEC int mpc_ui_div (mpc_ptr, unsigned long int, mpc_srcptr, mpc_rnd_t); -__MPC_DECLSPEC int mpc_div_2exp (mpc_ptr, mpc_srcptr, unsigned long int, mpc_rnd_t); -__MPC_DECLSPEC int mpc_mul_2exp (mpc_ptr, mpc_srcptr, unsigned long int, mpc_rnd_t); +__MPC_DECLSPEC int mpc_div_2ui (mpc_ptr, mpc_srcptr, unsigned long int, mpc_rnd_t); +__MPC_DECLSPEC int mpc_mul_2ui (mpc_ptr, mpc_srcptr, unsigned long int, mpc_rnd_t); +__MPC_DECLSPEC int mpc_div_2si (mpc_ptr, mpc_srcptr, long int, mpc_rnd_t); +__MPC_DECLSPEC int mpc_mul_2si (mpc_ptr, mpc_srcptr, long int, mpc_rnd_t); __MPC_DECLSPEC int mpc_conj (mpc_ptr, mpc_srcptr, mpc_rnd_t); __MPC_DECLSPEC int mpc_neg (mpc_ptr, mpc_srcptr, mpc_rnd_t); __MPC_DECLSPEC int mpc_norm (mpfr_ptr, mpc_srcptr, mpfr_rnd_t); @@ -125,11 +125,11 @@ mul_real (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) /* Signs of zeroes may be wrong. Their correction does not change the inexact flag. */ if (mpfr_zero_p (mpc_realref (z))) - mpfr_setsign (mpc_realref (z), mpc_realref (z), MPC_RND_RE(rnd) == GMP_RNDD - || (xrs != yrs && xis == yis), GMP_RNDN); + mpfr_setsign (mpc_realref (z), mpc_realref (z), MPC_RND_RE(rnd) == MPFR_RNDD + || (xrs != yrs && xis == yis), MPFR_RNDN); if (mpfr_zero_p (mpc_imagref (z))) - mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == GMP_RNDD - || (xrs != yis && xis != yrs), GMP_RNDN); + mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == MPFR_RNDD + || (xrs != yis && xis != yrs), MPFR_RNDN); return inex; } @@ -154,15 +154,15 @@ mul_imag (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) inex_re = -mpfr_mul (mpc_realref (rop), mpc_imagref (x), mpc_imagref (y), INV_RND (MPC_RND_RE (rnd))); - mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); /* exact */ + mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN); /* exact */ inex_im = mpfr_mul (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y), MPC_RND_IM (rnd)); mpc_set (z, rop, MPC_RNDNN); /* Sign of zeroes may be wrong (note that Re(z) cannot be zero) */ if (mpfr_zero_p (mpc_imagref (z))) - mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == GMP_RNDD - || sign, GMP_RNDN); + mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == MPFR_RNDD + || sign, MPFR_RNDN); if (overlap) mpc_clear (rop); @@ -187,10 +187,10 @@ mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, /* u=a*b, v=sign*c*d exactly */ mpfr_init2 (u, mpfr_get_prec (a) + mpfr_get_prec (b)); mpfr_init2 (v, mpfr_get_prec (c) + mpfr_get_prec (d)); - mpfr_mul (u, a, b, GMP_RNDN); - mpfr_mul (v, c, d, GMP_RNDN); + mpfr_mul (u, a, b, MPFR_RNDN); + mpfr_mul (v, c, d, MPFR_RNDN); if (sign < 0) - mpfr_neg (v, v, GMP_RNDN); + mpfr_neg (v, v, MPFR_RNDN); /* tentatively compute z as u+v; here we need z to be distinct from a, b, c, d to not lose the latter */ @@ -198,7 +198,7 @@ mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, if (mpfr_inf_p (z)) { /* replace by "correctly rounded overflow" */ - mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), GMP_RNDN); + mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), MPFR_RNDN); inex = mpfr_mul_2ui (z, z, mpfr_get_emax (), rnd); } else if (mpfr_zero_p (u) && !mpfr_zero_p (v)) { @@ -238,13 +238,13 @@ mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, mpz_add_si (ev, ev, (long int) ed); /* recompute u and v and move exponents to eu and ev */ - mpfr_mul (u, a, b, GMP_RNDN); + mpfr_mul (u, a, b, MPFR_RNDN); /* exponent of u is non-positive */ mpz_sub_ui (eu, eu, (unsigned long int) (-mpfr_get_exp (u))); mpfr_set_exp (u, (mpfr_prec_t) 0); - mpfr_mul (v, c, d, GMP_RNDN); + mpfr_mul (v, c, d, MPFR_RNDN); if (sign < 0) - mpfr_neg (v, v, GMP_RNDN); + mpfr_neg (v, v, MPFR_RNDN); mpz_sub_ui (ev, ev, (unsigned long int) (-mpfr_get_exp (v))); mpfr_set_exp (v, (mpfr_prec_t) 0); @@ -426,23 +426,23 @@ mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd) mpfr_init2 (u, 2); mpfr_init2 (x, 2); - inexact = mpfr_mul (v, a, d, GMP_RNDN); + inexact = mpfr_mul (v, a, d, MPFR_RNDN); if (inexact) { /* over- or underflow */ ok = 0; goto clear; } if (mul_a == -1) - mpfr_neg (v, v, GMP_RNDN); + mpfr_neg (v, v, MPFR_RNDN); - inexact = mpfr_mul (w, b, c, GMP_RNDN); + inexact = mpfr_mul (w, b, c, MPFR_RNDN); if (inexact) { /* over- or underflow */ ok = 0; goto clear; } if (mul_c == -1) - mpfr_neg (w, w, GMP_RNDN); + mpfr_neg (w, w, MPFR_RNDN); /* compute sign(v-w) */ sign_x = mpfr_cmp_abs (v, w); @@ -477,21 +477,21 @@ mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd) /* first compute away(b +/- a) and store it in u */ inexact = (mul_a == -1 ? - ROUND_AWAY (mpfr_sub (u, b, a, MPFR_RNDA), u) : - ROUND_AWAY (mpfr_add (u, b, a, MPFR_RNDA), u)); + mpfr_sub (u, b, a, MPFR_RNDA) : + mpfr_add (u, b, a, MPFR_RNDA)); /* then compute away(+/-c - d) and store it in x */ inexact |= (mul_c == -1 ? - ROUND_AWAY (mpfr_add (x, c, d, MPFR_RNDA), x) : - ROUND_AWAY (mpfr_sub (x, c, d, MPFR_RNDA), x)); + mpfr_add (x, c, d, MPFR_RNDA) : + mpfr_sub (x, c, d, MPFR_RNDA)); if (mul_c == -1) - mpfr_neg (x, x, GMP_RNDN); + mpfr_neg (x, x, MPFR_RNDN); if (inexact == 0) - mpfr_prec_round (u, prec_u = 2 * prec, GMP_RNDN); + mpfr_prec_round (u, prec_u = 2 * prec, MPFR_RNDN); /* compute away(u*x) and store it in u */ - inexact |= ROUND_AWAY (mpfr_mul (u, u, x, MPFR_RNDA), u); + inexact |= mpfr_mul (u, u, x, MPFR_RNDA); /* (a+b)*(c-d) */ /* if all computations are exact up to here, it may be that @@ -508,20 +508,20 @@ mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd) if (prec_x > prec_u) prec_x = prec_u; if (prec_x > prec) - mpfr_prec_round (x, prec_x, GMP_RNDN); + mpfr_prec_round (x, prec_x, MPFR_RNDN); } - rnd_u = (sign_u > 0) ? GMP_RNDU : GMP_RNDD; + rnd_u = (sign_u > 0) ? MPFR_RNDU : MPFR_RNDD; inexact |= mpfr_sub (x, v, w, rnd_u); /* ad - bc */ /* in case u=0, ensure that rnd_u rounds x away from zero */ if (mpfr_sgn (u) == 0) - rnd_u = (mpfr_sgn (x) > 0) ? GMP_RNDU : GMP_RNDD; + rnd_u = (mpfr_sgn (x) > 0) ? MPFR_RNDU : MPFR_RNDD; inexact |= mpfr_add (u, u, x, rnd_u); /* ac - bd */ ok = inexact == 0 || - mpfr_can_round (u, prec_u - 3, rnd_u, GMP_RNDZ, - prec_re + (rnd_re == GMP_RNDN)); + mpfr_can_round (u, prec_u - 3, rnd_u, MPFR_RNDZ, + prec_re + (rnd_re == MPFR_RNDN)); /* this ensures both we can round correctly and determine the correct inexact flag (for rounding to nearest) */ } diff --git a/src/mul_2si.c b/src/mul_2si.c new file mode 100644 index 0000000..14d0ca2 --- /dev/null +++ b/src/mul_2si.c @@ -0,0 +1,32 @@ +/* mpc_mul_2si -- Multiply a complex number by 2^e. + +Copyright (C) 2012 INRIA + +This file is part of GNU MPC. + +GNU MPC is free software; you can redistribute it and/or modify it under +the terms of the GNU Lesser General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY +WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for +more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see http://www.gnu.org/licenses/ . +*/ + +#include "mpc-impl.h" + +int +mpc_mul_2si (mpc_ptr a, mpc_srcptr b, long int c, mpc_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_mul_2si (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_mul_2si (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/src/mul_2ui.c b/src/mul_2ui.c new file mode 100644 index 0000000..46aa788 --- /dev/null +++ b/src/mul_2ui.c @@ -0,0 +1,32 @@ +/* mpc_mul_2ui -- Multiply a complex number by 2^e. + +Copyright (C) 2002, 2009, 2011, 2012 INRIA + +This file is part of GNU MPC. + +GNU MPC is free software; you can redistribute it and/or modify it under +the terms of the GNU Lesser General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY +WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for +more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see http://www.gnu.org/licenses/ . +*/ + +#include "mpc-impl.h" + +int +mpc_mul_2ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_mul_2ui (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_mul_2ui (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/src/mul_fr.c b/src/mul_fr.c index bd3574d..437b3ca 100644 --- a/src/mul_fr.c +++ b/src/mul_fr.c @@ -1,6 +1,6 @@ /* mpc_mul_fr -- Multiply a complex number by a floating-point number. -Copyright (C) 2002, 2008, 2009, 2010, 2011 INRIA +Copyright (C) 2002, 2008, 2009, 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -34,7 +34,7 @@ mpc_mul_fr (mpc_ptr a, mpc_srcptr b, mpfr_srcptr c, mpc_rnd_t rnd) inex_re = mpfr_mul (real, mpc_realref(b), c, MPC_RND_RE(rnd)); inex_im = mpfr_mul (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); - mpfr_set (mpc_realref (a), real, GMP_RNDN); /* exact */ + mpfr_set (mpc_realref (a), real, MPFR_RNDN); /* exact */ if (c == mpc_realref (a)) mpfr_clear (real); diff --git a/src/mul_i.c b/src/mul_i.c index 591b0c6..511b051 100644 --- a/src/mul_i.c +++ b/src/mul_i.c @@ -1,6 +1,6 @@ /* mpc_mul_i -- Multiply a complex number by plus or minus i. -Copyright (C) 2005, 2009, 2010, 2011 INRIA +Copyright (C) 2005, 2009, 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -35,8 +35,8 @@ mpc_mul_i (mpc_ptr a, mpc_srcptr b, int sign, mpc_rnd_t rnd) mpfr_swap (mpc_realref (a), mpc_imagref (a)); else { - mpfr_set (mpc_realref (a), mpc_imagref (b), GMP_RNDN); - mpfr_set (mpc_imagref (a), mpc_realref (b), GMP_RNDN); + mpfr_set (mpc_realref (a), mpc_imagref (b), MPFR_RNDN); + mpfr_set (mpc_imagref (a), mpc_realref (b), MPFR_RNDN); } if (sign >= 0) MPFR_CHANGE_SIGN (mpc_realref (a)); @@ -1,6 +1,6 @@ /* mpc_norm -- Square of the norm of a complex number. -Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011 INRIA +Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -77,8 +77,8 @@ mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) mpfr_set_prec (u, prec_u); mpfr_set_prec (v, prec_v); - inexact = mpfr_sqr (u, mpc_realref(b), GMP_RNDD); /* err <= 1 ulp in prec */ - inexact |= mpfr_sqr (v, mpc_imagref(b), GMP_RNDD); /* err <= 1 ulp in prec */ + inexact = mpfr_sqr (u, mpc_realref(b), MPFR_RNDD); /* err <= 1 ulp in prec */ + inexact |= mpfr_sqr (v, mpc_imagref(b), MPFR_RNDD); /* err <= 1 ulp in prec */ /* If loops = max_loops, inexact should be 0 here, except in case of underflow or overflow. @@ -86,12 +86,12 @@ mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) while-loop since it only remains to add u and v into a. */ if (inexact) { mpfr_set_prec (res, prec); - mpfr_add (res, u, v, GMP_RNDD); /* err <= 3 ulp in prec */ + mpfr_add (res, u, v, MPFR_RNDD); /* err <= 3 ulp in prec */ } } while (loops < max_loops && inexact != 0 - && !mpfr_can_round (res, prec - 2, GMP_RNDD, GMP_RNDU, - mpfr_get_prec (a) + (rnd == GMP_RNDN))); + && !mpfr_can_round (res, prec - 2, MPFR_RNDD, MPFR_RNDU, + mpfr_get_prec (a) + (rnd == MPFR_RNDN))); if (!inexact) /* squarings were exact, neither underflow nor overflow */ @@ -100,7 +100,7 @@ mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) since the norm is larger, there is an overflow for the norm */ else if (mpfr_overflow_p ()) { /* replace by "correctly rounded overflow" */ - mpfr_set_ui (a, 1ul, GMP_RNDN); + mpfr_set_ui (a, 1ul, MPFR_RNDN); inexact = mpfr_mul_2ui (a, a, mpfr_get_emax (), rnd); } else if (mpfr_underflow_p ()) { @@ -123,14 +123,14 @@ mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) && mpfr_get_exp (u) - 2 * (mpfr_exp_t) prec_u > emin && mpfr_get_exp (u) > -10) { mpfr_set_prec (v, MPFR_PREC_MIN); - mpfr_set_ui_2exp (v, 1, emin - 1, GMP_RNDZ); + mpfr_set_ui_2exp (v, 1, emin - 1, MPFR_RNDZ); inexact = mpfr_add (a, u, v, rnd); } else if (!mpfr_zero_p (v) && mpfr_get_exp (v) - 2 * (mpfr_exp_t) prec_v > emin && mpfr_get_exp (v) > -10) { mpfr_set_prec (u, MPFR_PREC_MIN); - mpfr_set_ui_2exp (u, 1, emin - 1, GMP_RNDZ); + mpfr_set_ui_2exp (u, 1, emin - 1, MPFR_RNDZ); inexact = mpfr_add (a, u, v, rnd); } else { @@ -145,17 +145,17 @@ mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) integer overflow */ if (mpfr_zero_p (u)) { /* recompute the scaled value exactly */ - mpfr_mul_2ui (u, mpc_realref (b), scale, GMP_RNDN); - mpfr_sqr (u, u, GMP_RNDN); + mpfr_mul_2ui (u, mpc_realref (b), scale, MPFR_RNDN); + mpfr_sqr (u, u, MPFR_RNDN); } else /* just scale */ - mpfr_mul_2ui (u, u, 2*scale, GMP_RNDN); + mpfr_mul_2ui (u, u, 2*scale, MPFR_RNDN); if (mpfr_zero_p (v)) { - mpfr_mul_2ui (v, mpc_imagref (b), scale, GMP_RNDN); - mpfr_sqr (v, v, GMP_RNDN); + mpfr_mul_2ui (v, mpc_imagref (b), scale, MPFR_RNDN); + mpfr_sqr (v, v, MPFR_RNDN); } else - mpfr_mul_2ui (v, v, 2*scale, GMP_RNDN); + mpfr_mul_2ui (v, v, 2*scale, MPFR_RNDN); inexact = mpfr_add (a, u, v, rnd); mpfr_clear_underflow (); @@ -137,7 +137,7 @@ fix_sign (mpc_ptr z, int sign_eps, int sign_a, mpfr_srcptr y) MPC_ASSERT (ymod4 == 1 || ymod4 == 3); if ((ymod4 == 3 && sign_eps == 0) || (ymod4 == 1 && sign_eps == 1)) - mpfr_neg (mpc_realref(z), mpc_realref(z), GMP_RNDZ); + mpfr_neg (mpc_realref(z), mpc_realref(z), MPFR_RNDZ); } else if (mpfr_zero_p (mpc_imagref(z))) { @@ -147,7 +147,7 @@ fix_sign (mpc_ptr z, int sign_eps, int sign_a, mpfr_srcptr y) MPC_ASSERT (ymod4 == 0 || ymod4 == 2); if ((ymod4 == 0 && sign_a == sign_eps) || (ymod4 == 2 && sign_a != sign_eps)) - mpfr_neg (mpc_imagref(z), mpc_imagref(z), GMP_RNDZ); + mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPFR_RNDZ); } end: @@ -187,7 +187,7 @@ mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd, { z_is_y = 1; mpfr_init2 (copy_of_y, mpfr_get_prec (y)); - mpfr_set (copy_of_y, y, GMP_RNDN); + mpfr_set (copy_of_y, y, MPFR_RNDN); } mpz_init (my); @@ -202,7 +202,7 @@ mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd, t = mpz_scan1 (my, 0); ey += (mpfr_exp_t) t; mpz_tdiv_q_2exp (my, my, t); - /* y = my*2^ey */ + /* y = my*2^ey with my odd */ if (x_imag) { @@ -229,7 +229,6 @@ mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd, mpz_mul_2exp (d, d, (unsigned long int) (ed - ec)); if ((mpfr_prec_t) mpz_sizeinbase (d, 2) > maxprec) goto end; - ed = ec; } else if (ed < ec) { @@ -484,7 +483,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) mpc_t t, u; mpfr_prec_t p, pr, pi, maxprec; int saved_underflow, saved_overflow; - + /* save the underflow or overflow flags from MPFR */ saved_underflow = mpfr_underflow_p (); saved_overflow = mpfr_overflow_p (); @@ -513,7 +512,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) cx1 > 0 if |x| > 1 */ mpfr_init (n); - inex = mpc_norm (n, x, GMP_RNDN); + inex = mpc_norm (n, x, MPFR_RNDN); cx1 = mpfr_cmp_ui (n, 1); if (cx1 == 0 && inex != 0) cx1 = -inex; @@ -526,7 +525,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) /* warning: mpc_set_ui_ui does not set Im(z) to -0 if Im(rnd)=RNDD */ ret = mpc_set_ui_ui (z, 1, 0, rnd); - if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi) + if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi) mpc_conj (z, z, MPC_RNDNN); mpfr_clear (n); @@ -575,7 +574,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) Note that the sign must also be set explicitly when rnd=RNDD because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0. */ - if (MPC_RND_IM (rnd) == GMP_RNDD || s1 != s2) + if (MPC_RND_IM (rnd) == MPFR_RNDD || s1 != s2) mpc_conj (z, z, MPC_RNDNN); goto end; } @@ -601,7 +600,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) Note that the sign must also be set explicitly when rnd=RNDD because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0. */ - if (MPC_RND_IM(rnd) == GMP_RNDD || s1 != s2) + if (MPC_RND_IM(rnd) == MPFR_RNDD || s1 != s2) mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPC_RND_IM(rnd)); goto end; } @@ -614,8 +613,8 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) (a) x is negative and y is half-an-integer (b) x = -1 and Re(y) is half-an-integer */ - if (mpfr_cmp_ui (mpc_realref(x), 0) < 0 && is_odd (mpc_realref(y), 1) && - (y_real || mpfr_cmp_si (mpc_realref(x), -1) == 0)) + if ((mpfr_cmp_ui (mpc_realref(x), 0) < 0) && is_odd (mpc_realref(y), 1) + && (y_real || mpfr_cmp_si (mpc_realref(x), -1) == 0)) z_imag = 1; } else /* x non real */ @@ -651,8 +650,8 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) p = 64; mpc_init2 (u, p); mpc_init2 (t, p); - pr += MPC_RND_RE(rnd) == GMP_RNDN; - pi += MPC_RND_IM(rnd) == GMP_RNDN; + pr += MPC_RND_RE(rnd) == MPFR_RNDN; + pi += MPC_RND_IM(rnd) == MPFR_RNDN; maxprec = MPC_MAX_PREC (z); x_imag = mpfr_zero_p (mpc_realref(x)); for (loop = 0;; loop++) @@ -703,8 +702,8 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) (see algorithms.tex) plus one due to the exponent difference: if z = a + I*b, where the relative error on z is at most 2^(-p), and EXP(a) = EXP(b) + k, the relative error on b is at most 2^(k-p) */ - if ((z_imag || (p > q + 3 + dr && mpfr_can_round (mpc_realref(u), p - q - 3 - dr, GMP_RNDN, GMP_RNDZ, pr))) && - (z_real || (p > q + 3 + di && mpfr_can_round (mpc_imagref(u), p - q - 3 - di, GMP_RNDN, GMP_RNDZ, pi)))) + if ((z_imag || (p > q + 3 + dr && mpfr_can_round (mpc_realref(u), p - q - 3 - dr, MPFR_RNDN, MPFR_RNDZ, pr))) && + (z_real || (p > q + 3 + di && mpfr_can_round (mpc_imagref(u), p - q - 3 - di, MPFR_RNDN, MPFR_RNDZ, pi)))) break; /* if Re(u) is not known to be zero, assume it is a normal number, i.e., @@ -751,7 +750,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) sign_rex = mpfr_signbit (mpc_realref (x)); sign_imx = mpfr_signbit (mpc_imagref (x)); mpfr_init (n); - inex = mpc_norm (n, x, GMP_RNDN); + inex = mpc_norm (n, x, MPFR_RNDN); cx1 = mpfr_cmp_ui (n, 1); if (cx1 == 0 && inex != 0) cx1 = -inex; @@ -762,7 +761,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) /* copy RE(y) to n since if z==y we will destroy Re(y) below */ mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y))); - mpfr_set (n, mpc_realref (y), GMP_RNDN); + mpfr_set (n, mpc_realref (y), MPFR_RNDN); ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd)); if (y_real && (x_real || x_imag)) { @@ -777,7 +776,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) { ret = MPC_INEX (ret, mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd))); /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */ - if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi) + if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi) mpc_conj (z, z, MPC_RNDNN); } diff --git a/src/pow_fr.c b/src/pow_fr.c index 8c5d930..87e255b 100644 --- a/src/pow_fr.c +++ b/src/pow_fr.c @@ -1,6 +1,6 @@ /* mpc_pow_fr -- Raise a complex number to a floating-point power. -Copyright (C) 2009, 2011 INRIA +Copyright (C) 2009, 2011, 2012 INRIA This file is part of GNU MPC. @@ -29,7 +29,7 @@ mpc_pow_fr (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd) /* avoid copying the significand of y by copying only the struct */ mpc_realref(yy)[0] = y[0]; mpfr_init2 (mpc_imagref(yy), MPFR_PREC_MIN); - mpfr_set_ui (mpc_imagref(yy), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref(yy), 0, MPFR_RNDN); inex = mpc_pow (z, x, yy, rnd); mpfr_clear (mpc_imagref(yy)); return inex; diff --git a/src/pow_ui.c b/src/pow_ui.c index da82a94..fb59310 100644 --- a/src/pow_ui.c +++ b/src/pow_ui.c @@ -131,10 +131,10 @@ mpc_pow_usi (mpc_ptr z, mpc_srcptr x, unsigned long y, int sign, /* the factor on the imaginary part is 2+2^(diff+2) <= 4 for diff <= -1 and < 2^(diff+3) for diff >= 0 */ ei = (diff <= -1) ? l0 + 3 : l0 + diff + 3; - if (mpfr_can_round (mpc_realref(t), p - er, GMP_RNDN, GMP_RNDZ, - MPC_PREC_RE(z) + (MPC_RND_RE(rnd) == GMP_RNDN)) - && mpfr_can_round (mpc_imagref(t), p - ei, GMP_RNDN, GMP_RNDZ, - MPC_PREC_IM(z) + (MPC_RND_IM(rnd) == GMP_RNDN))) { + if (mpfr_can_round (mpc_realref(t), p - er, MPFR_RNDN, MPFR_RNDZ, + MPC_PREC_RE(z) + (MPC_RND_RE(rnd) == MPFR_RNDN)) + && mpfr_can_round (mpc_imagref(t), p - ei, MPFR_RNDN, MPFR_RNDZ, + MPC_PREC_IM(z) + (MPC_RND_IM(rnd) == MPFR_RNDN))) { inex = mpc_set (z, t, rnd); done = 1; } diff --git a/src/sin_cos.c b/src/sin_cos.c index 0cff45a..0bb00cf 100644 --- a/src/sin_cos.c +++ b/src/sin_cos.c @@ -1,6 +1,6 @@ /* mpc_sin_cos -- combined sine and cosine of a complex number. -Copyright (C) 2010, 2011 INRIA +Copyright (C) 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -81,7 +81,7 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_t s, c; mpfr_init2 (s, 2); mpfr_init2 (c, 2); - mpfr_sin_cos (s, c, mpc_realref (op_loc), GMP_RNDZ); + mpfr_sin_cos (s, c, mpc_realref (op_loc), MPFR_RNDZ); mpfr_set_inf (mpc_realref (rop_sin), MPFR_SIGN (s)); mpfr_set_inf (mpc_imagref (rop_sin), MPFR_SIGN (c)*MPFR_SIGN (mpc_imagref (op_loc))); mpfr_clear (s); @@ -157,7 +157,7 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_t s, c; mpfr_init2 (c, 2); mpfr_init2 (s, 2); - mpfr_sin_cos (s, c, mpc_realref (op_loc), GMP_RNDN); + mpfr_sin_cos (s, c, mpc_realref (op_loc), MPFR_RNDN); mpfr_set_inf (mpc_realref (rop_cos), mpfr_sgn (c)); mpfr_set_inf (mpc_imagref (rop_cos), (mpfr_sgn (mpc_imagref (op_loc)) == mpfr_sgn (s) ? -1 : +1)); @@ -203,7 +203,7 @@ mpc_sin_cos_real (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, arbitrary rounding mode will work. */ if (rop_sin != NULL) { - mpfr_set (mpc_realref (rop_sin), s, GMP_RNDN); /* exact */ + mpfr_set (mpc_realref (rop_sin), s, MPFR_RNDN); /* exact */ inex_sin_re = inex_s; mpfr_set_zero (mpc_imagref (rop_sin), ( ( sign_im && !mpfr_signbit(c)) @@ -211,7 +211,7 @@ mpc_sin_cos_real (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, } if (rop_cos != NULL) { - mpfr_set (mpc_realref (rop_cos), c, GMP_RNDN); /* exact */ + mpfr_set (mpc_realref (rop_cos), c, MPFR_RNDN); /* exact */ inex_cos_re = inex_c; mpfr_set_zero (mpc_imagref (rop_cos), ( ( sign_im && mpfr_signbit(s)) @@ -245,7 +245,7 @@ mpc_sin_cos_imag (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, if (rop_sin != NULL) { /* sin(+-O +i*y) = +-0 +i*sinh(y) */ - mpfr_set (mpc_realref(rop_sin), mpc_realref(op_loc), GMP_RNDN); + mpfr_set (mpc_realref(rop_sin), mpc_realref(op_loc), MPFR_RNDN); inex_sin_im = mpfr_sinh (mpc_imagref(rop_sin), mpc_imagref(op_loc), MPC_RND_IM(rnd_sin)); } @@ -325,43 +325,43 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_set_prec (sch, prec); mpfr_set_prec (csh, prec); - mpfr_sin_cos (s, c, mpc_realref(op), GMP_RNDN); - mpfr_sinh_cosh (sh, ch, mpc_imagref(op), GMP_RNDN); + mpfr_sin_cos (s, c, mpc_realref(op), MPFR_RNDN); + mpfr_sinh_cosh (sh, ch, mpc_imagref(op), MPFR_RNDN); if (rop_sin != NULL) { /* real part of sine */ - mpfr_mul (sch, s, ch, GMP_RNDN); + mpfr_mul (sch, s, ch, MPFR_RNDN); ok = (!mpfr_number_p (sch)) - || mpfr_can_round (sch, prec - 2, GMP_RNDN, GMP_RNDZ, + || mpfr_can_round (sch, prec - 2, MPFR_RNDN, MPFR_RNDZ, MPC_PREC_RE (rop_sin) - + (MPC_RND_RE (rnd_sin) == GMP_RNDN)); + + (MPC_RND_RE (rnd_sin) == MPFR_RNDN)); if (ok) { /* imaginary part of sine */ - mpfr_mul (csh, c, sh, GMP_RNDN); + mpfr_mul (csh, c, sh, MPFR_RNDN); ok = (!mpfr_number_p (csh)) - || mpfr_can_round (csh, prec - 2, GMP_RNDN, GMP_RNDZ, + || mpfr_can_round (csh, prec - 2, MPFR_RNDN, MPFR_RNDZ, MPC_PREC_IM (rop_sin) - + (MPC_RND_IM (rnd_sin) == GMP_RNDN)); + + (MPC_RND_IM (rnd_sin) == MPFR_RNDN)); } } if (rop_cos != NULL && ok) { /* real part of cosine */ - mpfr_mul (c, c, ch, GMP_RNDN); + mpfr_mul (c, c, ch, MPFR_RNDN); ok = (!mpfr_number_p (c)) - || mpfr_can_round (c, prec - 2, GMP_RNDN, GMP_RNDZ, + || mpfr_can_round (c, prec - 2, MPFR_RNDN, MPFR_RNDZ, MPC_PREC_RE (rop_cos) - + (MPC_RND_RE (rnd_cos) == GMP_RNDN)); + + (MPC_RND_RE (rnd_cos) == MPFR_RNDN)); if (ok) { /* imaginary part of cosine */ - mpfr_mul (s, s, sh, GMP_RNDN); - mpfr_neg (s, s, GMP_RNDN); + mpfr_mul (s, s, sh, MPFR_RNDN); + mpfr_neg (s, s, MPFR_RNDN); ok = (!mpfr_number_p (s)) - || mpfr_can_round (s, prec - 2, GMP_RNDN, GMP_RNDZ, + || mpfr_can_round (s, prec - 2, MPFR_RNDN, MPFR_RNDZ, MPC_PREC_IM (rop_cos) - + (MPC_RND_IM (rnd_cos) == GMP_RNDN)); + + (MPC_RND_IM (rnd_cos) == MPFR_RNDN)); } } } while (ok == 0); @@ -1,6 +1,6 @@ /* mpc_sinh -- hyperbolic sine of a complex number. -Copyright (C)2008, 2009, 2011 INRIA +Copyright (C)2008, 2009, 2011, 2012 INRIA This file is part of GNU MPC. @@ -36,7 +36,7 @@ mpc_sinh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpc_realref (sin_z)[0] = mpc_imagref (rop)[0]; mpc_imagref (sin_z)[0] = mpc_realref (rop)[0]; - inex = mpc_sin (sin_z, z, RNDC (MPC_RND_IM (rnd), MPC_RND_RE (rnd))); + inex = mpc_sin (sin_z, z, MPC_RND (MPC_RND_IM (rnd), MPC_RND_RE (rnd))); /* sin_z and rop parts share the same significands, copy the rest now. */ mpc_realref (rop)[0] = mpc_imagref (sin_z)[0]; @@ -36,8 +36,8 @@ mpfr_fsss (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr c, mpfr_rnd_t rnd) /* u=a^2, v=c^2 exactly */ mpfr_init2 (u, 2*mpfr_get_prec (a)); mpfr_init2 (v, 2*mpfr_get_prec (c)); - mpfr_sqr (u, a, GMP_RNDN); - mpfr_sqr (v, c, GMP_RNDN); + mpfr_sqr (u, a, MPFR_RNDN); + mpfr_sqr (v, c, MPFR_RNDN); /* tentatively compute z as u-v; here we need z to be distinct from a and c to not lose the latter */ @@ -45,7 +45,7 @@ mpfr_fsss (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr c, mpfr_rnd_t rnd) if (mpfr_inf_p (z)) { /* replace by "correctly rounded overflow" */ - mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), GMP_RNDN); + mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), MPFR_RNDN); inex = mpfr_mul_2ui (z, z, mpfr_get_emax (), rnd); } else if (mpfr_zero_p (u) && !mpfr_zero_p (v)) { @@ -81,11 +81,11 @@ mpfr_fsss (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr c, mpfr_rnd_t rnd) mpz_mul_2exp (ev, ev, 1); /* recompute u and v and move exponents to eu and ev */ - mpfr_sqr (u, a, GMP_RNDN); + mpfr_sqr (u, a, MPFR_RNDN); /* exponent of u is non-positive */ mpz_sub_ui (eu, eu, (unsigned long int) (-mpfr_get_exp (u))); mpfr_set_exp (u, (mpfr_prec_t) 0); - mpfr_sqr (v, c, GMP_RNDN); + mpfr_sqr (v, c, MPFR_RNDN); mpz_sub_ui (ev, ev, (unsigned long int) (-mpfr_get_exp (v))); mpfr_set_exp (v, (mpfr_prec_t) 0); if (mpfr_nan_p (z)) { @@ -210,7 +210,7 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) if (mpfr_zero_p (mpc_imagref(op))) { int same_sign = mpfr_signbit (mpc_realref (op)) == mpfr_signbit (mpc_imagref (op)); inex_re = mpfr_sqr (mpc_realref(rop), mpc_realref(op), MPC_RND_RE(rnd)); - inex_im = mpfr_set_ui (mpc_imagref(rop), 0ul, GMP_RNDN); + inex_im = mpfr_set_ui (mpc_imagref(rop), 0ul, MPFR_RNDN); if (!same_sign) mpc_conj (rop, rop, MPC_RNDNN); return MPC_INEX(inex_re, inex_im); @@ -218,8 +218,8 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) if (mpfr_zero_p (mpc_realref(op))) { int same_sign = mpfr_signbit (mpc_realref (op)) == mpfr_signbit (mpc_imagref (op)); inex_re = -mpfr_sqr (mpc_realref(rop), mpc_imagref(op), INV_RND (MPC_RND_RE(rnd))); - mpfr_neg (mpc_realref(rop), mpc_realref(rop), GMP_RNDN); - inex_im = mpfr_set_ui (mpc_imagref(rop), 0ul, GMP_RNDN); + mpfr_neg (mpc_realref(rop), mpc_realref(rop), MPFR_RNDN); + inex_im = mpfr_set_ui (mpc_imagref(rop), 0ul, MPFR_RNDN); if (!same_sign) mpc_conj (rop, rop, MPC_RNDNN); return MPC_INEX(inex_re, inex_im); @@ -228,7 +228,7 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) if (rop == op) { mpfr_init2 (x, MPC_PREC_RE (op)); - mpfr_set (x, op->re, GMP_RNDN); + mpfr_set (x, op->re, MPFR_RNDN); } else x [0] = op->re [0]; @@ -266,23 +266,20 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* The error is bounded above by 1 ulp. */ /* We first let inexact be 1 if the real part is not computed */ /* exactly and determine the sign later. */ - inexact = ROUND_AWAY (mpfr_add (u, x, mpc_imagref (op), MPFR_RNDA), u) - | ROUND_AWAY (mpfr_sub (v, x, mpc_imagref (op), MPFR_RNDA), v); + inexact = mpfr_add (u, x, mpc_imagref (op), MPFR_RNDA) + | mpfr_sub (v, x, mpc_imagref (op), MPFR_RNDA); /* compute the real part as u*v, rounded away */ /* determine also the sign of inex_re */ if (mpfr_sgn (u) == 0 || mpfr_sgn (v) == 0) { /* as we have rounded away, the result is exact */ - mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN); inex_re = 0; ok = 1; } else { - mpfr_rnd_t rnd_away; - /* FIXME: can be replaced by MPFR_RNDA in mpfr >= 3 */ - rnd_away = (mpfr_sgn (u) * mpfr_sgn (v) > 0 ? GMP_RNDU : GMP_RNDD); - inexact |= ROUND_AWAY (mpfr_mul (u, u, v, MPFR_RNDA), u); /* error 5 */ + inexact |= mpfr_mul (u, u, v, MPFR_RNDA); /* error 5 */ if (mpfr_get_exp (u) == emin || mpfr_inf_p (u)) { /* under- or overflow */ inex_re = mpfr_fsss (rop->re, x, op->im, MPC_RND_RE (rnd)); @@ -290,8 +287,8 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) } else { ok = (!inexact) | mpfr_can_round (u, prec - 3, - rnd_away, GMP_RNDZ, - MPC_PREC_RE (rop) + (MPC_RND_RE (rnd) == GMP_RNDN)); + MPFR_RNDA, MPFR_RNDZ, + MPC_PREC_RE (rop) + (MPC_RND_RE (rnd) == MPFR_RNDN)); if (ok) { inex_re = mpfr_set (mpc_realref (rop), u, MPC_RND_RE (rnd)); if (inex_re == 0) @@ -311,7 +308,7 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_clear_underflow (); inex_im = mpfr_mul (rop->im, x, op->im, MPC_RND_IM (rnd)); if (!mpfr_underflow_p ()) - inex_im |= mpfr_mul_2exp (rop->im, rop->im, 1, MPC_RND_IM (rnd)); + inex_im |= mpfr_mul_2ui (rop->im, rop->im, 1, MPC_RND_IM (rnd)); /* We must not multiply by 2 if rop->im has been set to the smallest representable number. */ if (saved_underflow) @@ -1,6 +1,6 @@ /* mpc_sqrt -- Take the square root of a complex number. -Copyright (C) 2002, 2008, 2009, 2010, 2011 INRIA +Copyright (C) 2002, 2008, 2009, 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -46,7 +46,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) at the target precision */ const int im_sgn = mpfr_signbit (mpc_imagref (b)) == 0 ? 0 : -1; /* we need to know the sign of Im(b) when it is +/-0 */ - const mpfr_rnd_t r = im_sgn ? GMP_RNDD : GMP_RNDU; + const mpfr_rnd_t r = im_sgn ? MPFR_RNDD : MPFR_RNDU; /* rounding mode used when computing t */ /* special values */ @@ -68,7 +68,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) { /* sqrt(-Inf +i*y) = +0 +i*Inf, when y positive */ /* sqrt(-Inf +i*y) = +0 -i*Inf, when y positive */ - mpfr_set_ui (mpc_realref (a), 0, GMP_RNDN); + mpfr_set_ui (mpc_realref (a), 0, MPFR_RNDN); mpfr_set_inf (mpc_imagref (a), im_sgn); return MPC_INEX (0, 0); } @@ -87,7 +87,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) /* sqrt(+Inf +i*y) = +Inf +i*0, when y positive */ /* sqrt(+Inf +i*y) = +Inf -i*0, when y positive */ mpfr_set_inf (mpc_realref (a), +1); - mpfr_set_ui (mpc_imagref (a), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref (a), 0, MPFR_RNDN); if (im_sgn) mpc_conj (a, a, MPC_RNDNN); return MPC_INEX (0, 0); @@ -125,7 +125,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) else if (re_cmp > 0) { inex_w = mpfr_sqrt (mpc_realref (a), mpc_realref (b), MPC_RND_RE (rnd)); - mpfr_set_ui (mpc_imagref (a), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref (a), 0, MPFR_RNDN); if (im_sgn) mpc_conj (a, a, MPC_RNDNN); return MPC_INEX (inex_w, 0); @@ -133,16 +133,16 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) else { mpfr_init2 (w, MPC_PREC_RE (b)); - mpfr_neg (w, mpc_realref (b), GMP_RNDN); + mpfr_neg (w, mpc_realref (b), MPFR_RNDN); if (im_sgn) { inex_w = -mpfr_sqrt (mpc_imagref (a), w, INV_RND (MPC_RND_IM (rnd))); - mpfr_neg (mpc_imagref (a), mpc_imagref (a), GMP_RNDN); + mpfr_neg (mpc_imagref (a), mpc_imagref (a), MPFR_RNDN); } else inex_w = mpfr_sqrt (mpc_imagref (a), w, MPC_RND_IM (rnd)); - mpfr_set_ui (mpc_realref (a), 0, GMP_RNDN); + mpfr_set_ui (mpc_realref (a), 0, MPFR_RNDN); mpfr_clear (w); return MPC_INEX (0, inex_w); } @@ -155,7 +155,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) y[0] = mpc_imagref (b)[0]; /* If y/2 underflows, so does sqrt(y/2) */ - mpfr_div_2ui (y, y, 1, GMP_RNDN); + mpfr_div_2ui (y, y, 1, MPFR_RNDN); if (im_cmp > 0) { inex_w = mpfr_sqrt (mpc_realref (a), y, MPC_RND_RE (rnd)); @@ -163,10 +163,10 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) } else { - mpfr_neg (y, y, GMP_RNDN); + mpfr_neg (y, y, MPFR_RNDN); inex_w = mpfr_sqrt (mpc_realref (a), y, MPC_RND_RE (rnd)); inex_t = -mpfr_sqrt (mpc_imagref (a), y, INV_RND (MPC_RND_IM (rnd))); - mpfr_neg (mpc_imagref (a), mpc_imagref (a), GMP_RNDN); + mpfr_neg (mpc_imagref (a), mpc_imagref (a), MPFR_RNDN); } return MPC_INEX (inex_w, inex_t); } @@ -180,9 +180,9 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) rnd_w = MPC_RND_RE (rnd); prec_w = MPC_PREC_RE (a); rnd_t = MPC_RND_IM(rnd); - if (rnd_t == GMP_RNDZ) - /* force GMP_RNDD or GMP_RNDUP, using sign(t) = sign(y) */ - rnd_t = (im_cmp > 0 ? GMP_RNDD : GMP_RNDU); + if (rnd_t == MPFR_RNDZ) + /* force MPFR_RNDD or MPFR_RNDUP, using sign(t) = sign(y) */ + rnd_t = (im_cmp > 0 ? MPFR_RNDD : MPFR_RNDU); prec_t = MPC_PREC_IM (a); } else { @@ -191,14 +191,14 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) if (im_cmp > 0) { rnd_w = MPC_RND_IM(rnd); rnd_t = MPC_RND_RE(rnd); - if (rnd_t == GMP_RNDZ) - rnd_t = GMP_RNDD; + if (rnd_t == MPFR_RNDZ) + rnd_t = MPFR_RNDD; } else { rnd_w = INV_RND(MPC_RND_IM (rnd)); rnd_t = INV_RND(MPC_RND_RE (rnd)); - if (rnd_t == GMP_RNDZ) - rnd_t = GMP_RNDU; + if (rnd_t == MPFR_RNDZ) + rnd_t = MPFR_RNDU; } } @@ -211,30 +211,30 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) /* let b = x + iy */ /* w = sqrt ((|x| + sqrt (x^2 + y^2)) / 2), rounded down */ /* total error bounded by 3 ulps */ - inex_w = mpc_abs (w, b, GMP_RNDD); + inex_w = mpc_abs (w, b, MPFR_RNDD); if (re_cmp < 0) - inex_w |= mpfr_sub (w, w, mpc_realref (b), GMP_RNDD); + inex_w |= mpfr_sub (w, w, mpc_realref (b), MPFR_RNDD); else - inex_w |= mpfr_add (w, w, mpc_realref (b), GMP_RNDD); - inex_w |= mpfr_div_2ui (w, w, 1, GMP_RNDD); - inex_w |= mpfr_sqrt (w, w, GMP_RNDD); + inex_w |= mpfr_add (w, w, mpc_realref (b), MPFR_RNDD); + inex_w |= mpfr_div_2ui (w, w, 1, MPFR_RNDD); + inex_w |= mpfr_sqrt (w, w, MPFR_RNDD); repr_w = mpfr_min_prec (w) <= prec_w; if (!repr_w) /* use the usual trick for obtaining the ternary value */ - ok_w = mpfr_can_round (w, prec - 2, GMP_RNDD, GMP_RNDU, - prec_w + (rnd_w == GMP_RNDN)); + ok_w = mpfr_can_round (w, prec - 2, MPFR_RNDD, MPFR_RNDU, + prec_w + (rnd_w == MPFR_RNDN)); else { /* w is representable in the target precision and thus cannot be rounded up */ - if (rnd_w == GMP_RNDN) + if (rnd_w == MPFR_RNDN) /* If w can be rounded to nearest, then actually no rounding occurs, and the ternary value is known from inex_w. */ - ok_w = mpfr_can_round (w, prec - 2, GMP_RNDD, GMP_RNDN, prec_w); + ok_w = mpfr_can_round (w, prec - 2, MPFR_RNDD, MPFR_RNDN, prec_w); else /* If w can be rounded down, then any direct rounding and the ternary flag can be determined from inex_w. */ - ok_w = mpfr_can_round (w, prec - 2, GMP_RNDD, GMP_RNDD, prec_w); + ok_w = mpfr_can_round (w, prec - 2, MPFR_RNDD, MPFR_RNDD, prec_w); } if (!inex_w || ok_w) { @@ -249,11 +249,11 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) if (!repr_t) /* As for w; since t was rounded away, we check whether rounding to 0 is possible. */ - ok_t = mpfr_can_round (t, prec - 3, r, GMP_RNDZ, - prec_t + (rnd_t == GMP_RNDN)); + ok_t = mpfr_can_round (t, prec - 3, r, MPFR_RNDZ, + prec_t + (rnd_t == MPFR_RNDN)); else { - if (rnd_t == GMP_RNDN) - ok_t = mpfr_can_round (t, prec - 3, r, GMP_RNDN, prec_t); + if (rnd_t == MPFR_RNDN) + ok_t = mpfr_can_round (t, prec - 3, r, MPFR_RNDN, prec_t); else ok_t = mpfr_can_round (t, prec - 3, r, r, prec_t); } @@ -275,7 +275,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) } if (repr_w && inex_w) { - if (rnd_w == GMP_RNDN) { + if (rnd_w == MPFR_RNDN) { /* w has not been rounded with mpfr_set/mpfr_neg, determine ternary value from inex_w instead */ if (re_cmp > 0) @@ -289,7 +289,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) /* determine ternary value, but also potentially add 1 ulp; can only be done now when we are in the target precision */ if (re_cmp > 0) { - if (rnd_w == GMP_RNDU) { + if (rnd_w == MPFR_RNDU) { MPFR_ADD_ONE_ULP (mpc_realref (a)); inex_re = +1; } @@ -297,7 +297,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) inex_re = -1; } else if (im_cmp > 0) { - if (rnd_w == GMP_RNDU) { + if (rnd_w == MPFR_RNDU) { MPFR_ADD_ONE_ULP (mpc_imagref (a)); inex_im = +1; } @@ -305,7 +305,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) inex_im = -1; } else { - if (rnd_w == GMP_RNDU) { + if (rnd_w == MPFR_RNDU) { MPFR_ADD_ONE_ULP (mpc_imagref (a)); inex_im = -1; } @@ -315,7 +315,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) } } if (repr_t && inex_t) { - if (rnd_t == GMP_RNDN) { + if (rnd_t == MPFR_RNDN) { if (re_cmp > 0) inex_im = inex_t; else if (im_cmp > 0) @@ -329,11 +329,11 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) inex_im = inex_t; else { inex_im = -inex_t; - if ( (im_cmp > 0 && r == GMP_RNDD) - || (im_cmp < 0 && r == GMP_RNDU)) - MPFR_ADD_ONE_ULP (mpc_imagref (a)); - else - MPFR_SUB_ONE_ULP (mpc_imagref (a)); + /* im_cmp > 0 implies that Im(b) > 0, thus im_sgn = 0 + and r = MPFR_RNDU. + im_cmp < 0 implies that Im(b) < 0, thus im_sgn = -1 + and r = MPFR_RNDD. */ + MPFR_SUB_ONE_ULP (mpc_imagref (a)); } } else if (im_cmp > 0) { @@ -341,21 +341,17 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) inex_re = inex_t; else { inex_re = -inex_t; - if (r == GMP_RNDD) - MPFR_ADD_ONE_ULP (mpc_realref (a)); - else - MPFR_SUB_ONE_ULP (mpc_realref (a)); + /* im_cmp > 0 implies r = MPFR_RNDU (see above) */ + MPFR_SUB_ONE_ULP (mpc_realref (a)); } } - else { + else { /* im_cmp < 0 */ if (rnd_t == r) inex_re = -inex_t; else { inex_re = inex_t; - if (r == GMP_RNDD) - MPFR_SUB_ONE_ULP (mpc_realref (a)); - else - MPFR_ADD_ONE_ULP (mpc_realref (a)); + /* im_cmp < 0 implies r = MPFR_RNDD (see above) */ + MPFR_SUB_ONE_ULP (mpc_realref (a)); } } } diff --git a/src/strtoc.c b/src/strtoc.c index b96ccee..81702e2 100644 --- a/src/strtoc.c +++ b/src/strtoc.c @@ -1,6 +1,6 @@ /* mpc_strtoc -- Read a complex number from a string. -Copyright (C) 2009, 2010, 2011 INRIA +Copyright (C) 2009, 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -57,7 +57,7 @@ mpc_strtoc (mpc_ptr rop, const char *nptr, char **endptr, int base, mpc_rnd_t rn p = end; if (!bracketed) - inex_im = mpfr_set_ui (mpc_imagref (rop), 0ul, GMP_RNDN); + inex_im = mpfr_set_ui (mpc_imagref (rop), 0ul, MPFR_RNDN); else { if (!isspace ((unsigned char)*p)) goto error; @@ -1,6 +1,6 @@ /* mpc_tan -- tangent of a complex number. -Copyright (C) 2008, 2009, 2010, 2011 INRIA +Copyright (C) 2008, 2009, 2010, 2011, 2012 INRIA This file is part of GNU MPC. @@ -83,7 +83,7 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) int inex_im; mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd)); - mpfr_setsign (mpc_realref (rop), mpc_realref (rop), sign_re, GMP_RNDN); + mpfr_setsign (mpc_realref (rop), mpc_realref (rop), sign_re, MPFR_RNDN); /* exact, unless 1 is not in exponent range */ inex_im = mpfr_set_si (mpc_imagref (rop), @@ -111,10 +111,10 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_init (c); mpfr_init (s); - mpfr_sin_cos (s, c, mpc_realref (op), GMP_RNDN); + mpfr_sin_cos (s, c, mpc_realref (op), MPFR_RNDN); mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd)); mpfr_setsign (mpc_realref (rop), mpc_realref (rop), - mpfr_signbit (c) != mpfr_signbit (s), GMP_RNDN); + mpfr_signbit (c) != mpfr_signbit (s), MPFR_RNDN); /* exact, unless 1 is not in exponent range */ inex_im = mpfr_set_si (mpc_imagref (rop), (mpfr_signbit (mpc_imagref (op)) ? -1 : +1), @@ -207,22 +207,22 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) sign(tan(Re(op)))*0 + sign(Im(op))*I, where sign(tan(Re(op))) = sign(Re(x))*sign(Re(y)). */ int inex_re, inex_im; - mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN); if (mpfr_sgn (mpc_realref (x)) * mpfr_sgn (mpc_realref (y)) < 0) { - mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); + mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN); inex_re = 1; } else inex_re = -1; /* +0 is rounded down */ if (mpfr_sgn (mpc_imagref (op)) > 0) { - mpfr_set_ui (mpc_imagref (rop), 1, GMP_RNDN); + mpfr_set_ui (mpc_imagref (rop), 1, MPFR_RNDN); inex_im = 1; } else { - mpfr_set_si (mpc_imagref (rop), -1, GMP_RNDN); + mpfr_set_si (mpc_imagref (rop), -1, MPFR_RNDN); inex_im = -1; } inex = MPC_INEX(inex_re, inex_im); @@ -261,15 +261,15 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* Can the real part be rounded? */ ok = (!mpfr_number_p (mpc_realref (x))) - || mpfr_can_round (mpc_realref(x), prec - err, GMP_RNDN, GMP_RNDZ, - MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == GMP_RNDN)); + || mpfr_can_round (mpc_realref(x), prec - err, MPFR_RNDN, MPFR_RNDZ, + MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == MPFR_RNDN)); if (ok) { /* Can the imaginary part be rounded? */ ok = (!mpfr_number_p (mpc_imagref (x))) - || mpfr_can_round (mpc_imagref(x), prec - 6, GMP_RNDN, GMP_RNDZ, - MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == GMP_RNDN)); + || mpfr_can_round (mpc_imagref(x), prec - 6, MPFR_RNDN, MPFR_RNDZ, + MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == MPFR_RNDN)); } } while (ok == 0); @@ -1,6 +1,6 @@ /* mpc_tanh -- hyperbolic tangent of a complex number. -Copyright (C) 2008, 2009, 2011 INRIA +Copyright (C) 2008, 2009, 2011, 2012 INRIA This file is part of GNU MPC. @@ -36,7 +36,7 @@ mpc_tanh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpc_realref (tan_z)[0] = mpc_imagref (rop)[0]; mpc_imagref (tan_z)[0] = mpc_realref (rop)[0]; - inex = mpc_tan (tan_z, z, RNDC (MPC_RND_IM (rnd), MPC_RND_RE (rnd))); + inex = mpc_tan (tan_z, z, MPC_RND (MPC_RND_IM (rnd), MPC_RND_RE (rnd))); /* tan_z and rop parts share the same significands, copy the rest now. */ mpc_realref (rop)[0] = mpc_imagref (tan_z)[0]; |