diff options
Diffstat (limited to 'src/pow.c')
| -rw-r--r-- | src/pow.c | 37 |
1 files changed, 18 insertions, 19 deletions
@@ -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); } |