[pow.c] fixed two more bugs

[pow.dat] added some exact cases (thanks Andreas)


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

2 files changed, 20 insertions, 6 deletions

diff --git a/src/pow.c b/src/pow.c
index 44e3c9e..ab9eae0 100644
--- a/src/pow.c
+++ b/src/pow.c
@@ -184,7 +184,9 @@ mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd)
     }
 
   /* Now ey >= 0, it thus suffices to check that x^my is representable.
-     If my > 0, this is always true. */
+     If my > 0, this is always true. If my < 0, we first try to invert
+     (c+I*d)*2^ec.
+  */
   if (mpz_cmp_ui (my, 0) < 0)
     {
       /* If my < 0, 1 / (c + I*d) =
@@ -214,7 +216,7 @@ mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd)
           mpz_div_2exp (c, c, t); /* now c = sign(c0) */
           mpz_div_2exp (d, d, t); /* now d = sign(d0) */
           mpz_neg (d, d);         /* now d = -sign(d0) */
-          ec -= t + 1;
+          ec = -ec - (t + 1);
         }
       else
         {
@@ -382,11 +384,10 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
         }
 
       /* x^y is real when:
-         (a) x is real and y is real non-negative or integer
+         (a) x is real and y is integer
          (b) x is real non-negative and y is real */
-      if (y_real && (mpfr_cmp_ui (MPC_RE (x), 0) >= 0 ||
-                     mpfr_cmp_ui (MPC_RE (y), 0) >= 0 ||
-                     mpfr_integer_p (MPC_RE (y))))
+      if (y_real && (mpfr_integer_p (MPC_RE (y)) ||
+                     mpfr_cmp_ui (MPC_RE (x), 0) >= 0))
         {
           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 8e8dd4e..9a49398 100644
--- a/tests/pow.dat
+++ b/tests/pow.dat
@@ -20,3 +20,16 @@
 # MA 02111-1307, USA.
 #
 # For explanations on the file format, see add.dat.
+
+# exact cases
+# (-4)^(1/4) = 1+i
+0 0 2 1 2 1 2 -4 2 0 2 0x1p-2 2 0 N N
+# for an odd positive integer n, a positive integer m and an integer e:
+# (-4 m^4 16^e)^(n/4) = (1+i)^n m^n 2^(e n)
+# m=3 e=5 n=7
+0 0 12 0x88Bp38 12 -0x88Bp38 7 -0x51p22 7 0 3 0x7p-2 3 0 N N
+# (-4 16^e)^(-n/4)    = (1-i)^n 2^(- (e+1) n)
+# e=3 n=5
+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