diff options
| author | zimmerma <zimmerma@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2009-06-27 08:42:49 +0000 |
|---|---|---|
| committer | zimmerma <zimmerma@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2009-06-27 08:42:49 +0000 |
| commit | 316fadea81084c3d5b06c386da6a3068746c5686 (patch) | |
| tree | 767a5cee11ad84cda19217ecafbb9695e0114aec | |
| parent | 552b9c61ef3414b4bccc12baa284edc300d71f44 (diff) | |
| download | mpc-316fadea81084c3d5b06c386da6a3068746c5686.tar.gz | |
[pow.c] more work on sign of zero
git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@619 211d60ee-9f03-0410-a15a-8952a2c7a4e4
| -rw-r--r-- | src/pow.c | 15 | ||||
| -rw-r--r-- | tests/pow.dat | 10 |
2 files changed, 18 insertions, 7 deletions
@@ -478,15 +478,16 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) log(x) ~ log(x) + eps/x*i thus y*log(x) ~ y*log(x) + [y/x*eps + log(x)*delta]*i For x < 0 and y integer: - log(x) ~ log|x| + [Pi + eps/x]*i thus - y*log(x) ~ y*log|x| + [y*Pi + y/x*eps + log|x|*delta]*i + log(x) ~ log|x| + sign(eps)*[Pi - |eps/x|]*i thus + y*log(x) ~ y*log|x| + [sign(eps)*y*Pi-y/|x|*eps+log|x|*delta]*i thus if y is even we get the same term as for x > 0 for the imaginary part, otherwise its opposite. */ - s1 = mpfr_signbit (MPC_RE(x)) ? -1 : 1; - s1 *= mpfr_signbit (MPC_RE(y)) ? -1 : 1; /* sign of y/x */ - s1 *= MPFR_SIGN(MPC_IM(x)); /* sign of y/x*eps */ + s1 = mpfr_signbit (MPC_RE(x)) ? 1 : -1; /* sign of x */ + s1 *= mpfr_signbit (MPC_RE(y)) ? 1 : -1; /* sign of y/x for x > 0 + or -y/|x| for x < 0 */ + s1 *= MPFR_SIGN(MPC_IM(x)); /* sign of -y*eps */ s2 = mpfr_get_exp (MPC_RE(x)) >= 1 ? 1 : -1; /* sign of log|x| */ - s2 *= MPFR_SIGN(MPC_IM(y)); /* sign of log(x)*delta */ + s2 *= MPFR_SIGN(MPC_IM(y)); /* sign of log|x|*delta */ if (s1 != s2) { if (MPC_RND_IM(rnd) == GMP_RNDD) @@ -494,7 +495,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) else s1 = 1; /* take +0 as arbitrary sign */ } - if (MPFR_SIGN(MPC_RE(x)) < 0 && is_odd (MPC_RE(y), 0)) + else if (MPFR_SIGN(MPC_RE(x)) < 0 && is_odd (MPC_RE(y), 0)) s1 = -s1; ret = mpfr_pow (MPC_RE(z), MPC_RE(x), MPC_RE(y), MPC_RND_RE(rnd)); ret = MPC_INEX(ret, mpfr_set_ui (MPC_IM(z), 0, MPC_RND_IM(rnd))); diff --git a/tests/pow.dat b/tests/pow.dat index d58865a..7baaad0 100644 --- a/tests/pow.dat +++ b/tests/pow.dat @@ -186,6 +186,16 @@ # log(x) = log(2) + epsilon/2*i + O(epsilon^2) # y*log(x) = [-3*log(2) + o(1)] + [-3*epsilon/2-delta*log(2)]*i 0 0 2 0x1p-3 2 -0 2 2 2 +0 2 -3 2 -0 N N +# (-2 -0)^(3 +0) -> (-8 -0) +# x = -2 - epsilon*i, y = 3 + delta*i +# log(x) = log(2) - [Pi-epsilon/2]*i + O(epsilon^2) +# y*log(x) ~ 3*log(2) + [-3*Pi+3*epsilon/2+delta*log(2)]*i +0 0 2 -8 2 -0 2 -2 2 -0 2 3 2 +0 N N +# (-2 -0)^(-3 -0) -> (-1/8 +0) +# x = -2 - epsilon*i, y = -3 - delta*i +# log(x) = log(2) - [Pi-epsilon/2]*i + O(epsilon^2) +# y*log(x) ~ -3*log(2) + [3*Pi-3*epsilon/2-delta*log(2)]*i +0 0 2 -0x1p-3 2 +0 2 -2 2 -0 2 -3 2 -0 N N 0 0 2 +0 2 -2 2 +0 2 0x1p-1 2 -1 2 -0 N N 0 - 2 +0 3 -5 2 +0 53 0xCCCCCCCCCCCCDp-54 2 -1 2 -0 N N - 0 5 -25 2 +0 2 +0 53 0xCCCCCCCCCCCCDp-54 2 -2 2 -0 N N |