diff options
| author | zimmerma <zimmerma@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2009-06-23 13:23:27 +0000 |
|---|---|---|
| committer | zimmerma <zimmerma@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2009-06-23 13:23:27 +0000 |
| commit | 8ec9875d717be2483344a3054d56915f3b1fecdc (patch) | |
| tree | a80ad2e7852c5b6b9cc431e73e8a964ec96b7734 | |
| parent | 49981fe688dfd8575766a3838ffbae2b47646a57 (diff) | |
| download | mpc-8ec9875d717be2483344a3054d56915f3b1fecdc.tar.gz | |
[pow.c] improved sign of result for x^0
git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@614 211d60ee-9f03-0410-a15a-8952a2c7a4e4
| -rw-r--r-- | src/pow.c | 27 | ||||
| -rw-r--r-- | tests/pow.dat | 2 |
2 files changed, 24 insertions, 5 deletions
@@ -356,14 +356,31 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) if (y_real && mpfr_zero_p (MPC_RE (y))) /* case y zero */ { - if (x_real && mpfr_zero_p (MPC_RE (x))) /* 0^0 = NaN +i*NaN */ + if (x_real && mpfr_zero_p (MPC_RE(x))) /* 0^0 = NaN +i*NaN */ { - mpfr_set_nan (MPC_RE (z)); - mpfr_set_nan (MPC_IM (z)); + mpfr_set_nan (MPC_RE(z)); + mpfr_set_nan (MPC_IM(z)); return 0; } - else /* x^0 = 1 +i*0 even for NaN */ - return mpc_set_ui_ui (z, 1, 0, rnd); + else /* x^0 = 1 +/- i*0 even for x=NaN */ + { + int sa, sb, sc, sd; + mpc_init2 (u, 2); + mpc_log (u, x, MPC_RNDNN); + /* sa = sign(Re(log(x))), sb = sign(Im(log(x))) */ + sa = mpfr_signbit (MPC_RE(u)) ? -2 : 2; + if (mpfr_cmp_ui (MPC_RE(x), 1) == 0) + sa = 1; + sb = mpfr_signbit (MPC_IM(u)) ? -2 : 2; + /* sc = sign(Re(y)), sd = sign(Im(y)) */ + sc = mpfr_signbit (MPC_RE(y)) ? -2 : 2; + sd = mpfr_signbit (MPC_IM(y)) ? -2 : 2; + ret = mpc_set_ui_ui (z, 1, 0, rnd); + if (sa * sd + sb * sc < 0) + mpfr_neg (MPC_IM(z), MPC_IM(z), MPC_RND_IM(rnd)); + mpc_clear (u); + return ret; + } } if (x_real) /* case x real */ diff --git a/tests/pow.dat b/tests/pow.dat index b945121..cffc2b9 100644 --- a/tests/pow.dat +++ b/tests/pow.dat @@ -161,6 +161,8 @@ 0 0 53 +inf 53 nan 53 +0 53 +0 53 -inf 53 -1 N N 0 0 53 +1 53 +0 53 +0 53 +1 53 +0 53 +0 N N +# (1-0*I)^(+0+0*I) -> 1-0*I +0 0 53 +1 53 -0 53 +1 53 -0 53 +0 53 +0 N N # exact cases # (-4)^(1/4) = 1+i |