[pow.c] more work on cases with sign of zero as input/output

git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@618 211d60ee-9f03-0410-a15a-8952a2c7a4e4

4 files changed, 66 insertions, 2 deletions

diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index 4dc9095..a1d09e2 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -1175,6 +1175,35 @@ A relevant reference is \cite{BrDiJeLeMeMuReStTo09}, especially Section 4.5
 which discusses complex floating-point numbers, and gives error bounds for
 multiplication, division and square root.
 
+\paragraph{Sign of zeroes.}
+When the output value has a zero real or imaginary part, its sign should be
+decided. When the inputs also have a zero real or imaginary part, we
+consider all possible limits, and if all those limits give the same sign,
+we take this as the sign of the zero part.
+Otherwise, we round to $+0$, except for rounding toward $-\infty$, where we
+round to $-0$.
+
+Example~: consider $x = 1 + 0 i$ and $y = -0 + 0 i$, we consider that $x$
+is the limit of $1 + \epsilon i$ for $\epsilon > 0$ tending to zero.
+Similarly, $y = -\delta + \gamma i$ with $\delta, \gamma > 0$.
+Now $\log x \approx \epsilon^2/2 + \epsilon i$, thus
+$y \log x \approx (-\epsilon^2 \delta/2 - \epsilon \gamma)
++ i (\epsilon^2 \gamma/2 - epsilon \delta)$.
+Thus if we neglect terms of order $3$ or more, 
+whatever the relative growth of $\epsilon, \delta, \gamma$, the real and
+imaginary parts of $y \log x$ are negative, thus we decide $x^y$ is rounded
+to $1 - 0 i$.
+
+For $x = 1 - 0 i$ and $y = 0 + 0 i$, we write
+$x = 1 - \epsilon i$ and $y = \delta + \gamma i$,
+which gives $\log x \approx \epsilon^2/2 - \epsilon i$
+and $y \log x \approx \epsilon \gamma - i \epsilon \delta + O(\epsilon^2)$,
+thus we also round $x^y$ to $1 - 0 i$.
+
+% (1 -0)^(-0 -0):
+% x = 1 - epsilon i, y = -delta-gamma*i
+% y log(x) = -epsilon gamma + i epsilon delta + O(epsilon^2)
+
 \bibliographystyle{acm}
 \bibliography{algorithms}
 
diff --git a/src/pow.c b/src/pow.c
index 2e9b17b..9849c46 100644
--- a/src/pow.c
+++ b/src/pow.c
@@ -473,8 +473,33 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
       if (y_real && (mpfr_integer_p (MPC_RE(y)) ||
                      mpfr_cmp_ui (MPC_RE(x), 0) >= 0))
         {
+          int s1, s2;
+          /* if x = x + eps*i and y = y + delta*i, then for x > 0:
+             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
+             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 */
+          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 */
+          if (s1 != s2)
+            {
+              if (MPC_RND_IM(rnd) == GMP_RNDD)
+                s1 = -1;
+              else
+                s1 = 1; /* take +0 as arbitrary sign */
+            }
+          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)));
+          if (s1 == -1)
+            mpfr_neg (MPC_IM(z), MPC_IM(z), MPC_RND_IM(rnd));
           goto end;
         }
 
diff --git a/tests/Makefile.am b/tests/Makefile.am
index b01ca1f..d9e8394 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -8,7 +8,7 @@ check_PROGRAMS = tabs tadd tadd_fr tadd_ui targ tconj tcos tcosh tdiv         \
 tdiv_2exp tdiv_fr tdiv_ui texp tfr_div tfr_sub tget_version timag tio_str     \
 tlog tmul tmul_2exp tmul_fr tmul_i tmul_si tmul_ui tneg tnorm tpow tprec      \
 tproj treal treimref tset tsin tsinh tsqr tsqrt tstrtoc tsub tsub_fr tsub_ui  \
-ttan ttanh tui_div tui_ui_sub 
+ttan ttanh tui_div tui_ui_sub
 
 check_LTLIBRARIES=libmpc-tests.la
 libmpc_tests_la_SOURCES=mpc-tests.h random.c tgeneric.c read_data.c           \
diff --git a/tests/pow.dat b/tests/pow.dat
index da4f19e..d58865a 100644
--- a/tests/pow.dat
+++ b/tests/pow.dat
@@ -161,8 +161,13 @@
 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
+# zero cases from algorithms.tex
+# (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
+# (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
+# (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
@@ -176,6 +181,11 @@
 0 0 2 -0x1p-18 2 0x1p-18 2 -0x1p14 2 0 3 -0x5p-2 3 0 N N
 # e=2 n=5
 0 0 2 -0x1p-13 2 0x1p-13 2 -0x1p10 2 0 3 -0x5p-2 3 0 N N
+# (+2 +0)^(-3 -0) -> (-1/8 -0)
+# x = 2 + epsilon*i, y = -3 - delta*i
+# 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
 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