contrib/mpc/src/pow.c

   1 /* mpc_pow -- Raise a complex number to the power of another complex number.
   2
   3 Copyright (C) 2009, 2010, 2011, 2012 INRIA
   4
   5 This file is part of GNU MPC.
   6
   7 GNU MPC is free software; you can redistribute it and/or modify it under
   8 the terms of the GNU Lesser General Public License as published by the
   9 Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
  13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
  15 more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with this program. If not, see http://www.gnu.org/licenses/ .
  19 */
  20
  21 #include <stdio.h> /* for MPC_ASSERT */
  22 #include "mpc-impl.h"
  23
  24 /* Return non-zero iff c+i*d is an exact square (a+i*b)^2,
  25    with a, b both of the form m*2^e with m, e integers.
  26    If so, returns in a+i*b the corresponding square root, with a >= 0.
  27    The variables a, b must not overlap with c, d.
  28
  29    We have c = a^2 - b^2 and d = 2*a*b.
  30
  31    If one of a, b is exact, then both are (see algorithms.tex).
  32
  33    Case 1: a <> 0 and b <> 0.
  34    Let a = m*2^e and b = n*2^f with m, e, n, f integers, m and n odd
  35    (we will treat apart the case a = 0 or b = 0).
  36    Then 2*a*b = m*n*2^(e+f+1), thus necessarily e+f >= -1.
  37    Assume e < 0, then f >= 0, then a^2 - b^2 = m^2*2^(2e) - n^2*2^(2f) cannot
  38    be an integer, since n^2*2^(2f) is an integer, and m^2*2^(2e) is not.
  39    Similarly when f < 0 (and thus e >= 0).
  40    Thus we have e, f >= 0, and a, b are both integers.
  41    Let A = 2a^2, then eliminating b between c = a^2 - b^2 and d = 2*a*b
  42    gives A^2 - 2c*A - d^2 = 0, which has solutions c +/- sqrt(c^2+d^2).
  43    We thus need c^2+d^2 to be a square, and c + sqrt(c^2+d^2) --- the solution
  44    we are interested in --- to be two times a square. Then b = d/(2a) is
  45    necessarily an integer.
  46
  47    Case 2: a = 0. Then d is necessarily zero, thus it suffices to check
  48    whether c = -b^2, i.e., if -c is a square.
  49
  50    Case 3: b = 0. Then d is necessarily zero, thus it suffices to check
  51    whether c = a^2, i.e., if c is a square.
  52 */
  53 static int
  54 mpc_perfect_square_p (mpz_t a, mpz_t b, mpz_t c, mpz_t d)
  55 {
  56   if (mpz_cmp_ui (d, 0) == 0) /* case a = 0 or b = 0 */
  57     {
  58       /* necessarily c < 0 here, since we have already considered the case
  59          where x is real non-negative and y is real */
  60       MPC_ASSERT (mpz_cmp_ui (c, 0) < 0);
  61       mpz_neg (b, c);
  62       if (mpz_perfect_square_p (b)) /* case 2 above */
  63         {
  64           mpz_sqrt (b, b);
  65           mpz_set_ui (a, 0);
  66           return 1; /* c + i*d = (0 + i*b)^2 */
  67         }
  68     }
  69   else /* both a and b are non-zero */
  70     {
  71       if (mpz_divisible_2exp_p (d, 1) == 0)
  72         return 0; /* d must be even */
  73       mpz_mul (a, c, c);
  74       mpz_addmul (a, d, d); /* c^2 + d^2 */
  75       if (mpz_perfect_square_p (a))
  76         {
  77           mpz_sqrt (a, a);
  78           mpz_add (a, c, a); /* c + sqrt(c^2+d^2) */
  79           if (mpz_divisible_2exp_p (a, 1))
  80             {
  81               mpz_tdiv_q_2exp (a, a, 1);
  82               if (mpz_perfect_square_p (a))
  83                 {
  84                   mpz_sqrt (a, a);
  85                   mpz_tdiv_q_2exp (b, d, 1); /* d/2 */
  86                   mpz_divexact (b, b, a); /* d/(2a) */
  87                   return 1;
  88                 }
  89             }
  90         }
  91     }
  92   return 0; /* not a square */
  93 }
  94
  95 /* fix the sign of Re(z) or Im(z) in case it is zero,
  96    and Re(x) is zero.
  97    sign_eps is 0 if Re(x) = +0, 1 if Re(x) = -0
  98    sign_a is the sign bit of Im(x).
  99    Assume y is an integer (does nothing otherwise).
 100 */
 101 static void
 102 fix_sign (mpc_ptr z, int sign_eps, int sign_a, mpfr_srcptr y)
 103 {
 104   int ymod4 = -1;
 105   mpfr_exp_t ey;
 106   mpz_t my;
 107   unsigned long int t;
 108
 109   mpz_init (my);
 110
 111   ey = mpfr_get_z_exp (my, y);
 112   /* normalize so that my is odd */
 113   t = mpz_scan1 (my, 0);
 114   ey += (mpfr_exp_t) t;
 115   mpz_tdiv_q_2exp (my, my, t);
 116   /* y = my*2^ey */
 117
 118   /* compute y mod 4 (in case y is an integer) */
 119   if (ey >= 2)
 120     ymod4 = 0;
 121   else if (ey == 1)
 122     ymod4 = mpz_tstbit (my, 0) * 2; /* correct if my < 0 */
 123   else if (ey == 0)
 124     {
 125       ymod4 = mpz_tstbit (my, 1) * 2 + mpz_tstbit (my, 0);
 126       if (mpz_cmp_ui (my , 0) < 0)
 127         ymod4 = 4 - ymod4;
 128     }
 129   else /* y is not an integer */
 130     goto end;
 131
 132   if (mpfr_zero_p (mpc_realref(z)))
 133     {
 134       /* we assume y is always integer in that case (FIXME: prove it):
 135          (eps+I*a)^y = +0 + I*a^y for y = 1 mod 4 and sign_eps = 0
 136          (eps+I*a)^y = -0 - I*a^y for y = 3 mod 4 and sign_eps = 0 */
 137       MPC_ASSERT (ymod4 == 1 || ymod4 == 3);
 138       if ((ymod4 == 3 && sign_eps == 0) ||
 139           (ymod4 == 1 && sign_eps == 1))
 140         mpfr_neg (mpc_realref(z), mpc_realref(z), GMP_RNDZ);
 141     }
 142   else if (mpfr_zero_p (mpc_imagref(z)))
 143     {
 144       /* we assume y is always integer in that case (FIXME: prove it):
 145          (eps+I*a)^y =  a^y - 0*I for y = 0 mod 4 and sign_a = sign_eps
 146          (eps+I*a)^y =  -a^y +0*I for y = 2 mod 4 and sign_a = sign_eps */
 147       MPC_ASSERT (ymod4 == 0 || ymod4 == 2);
 148       if ((ymod4 == 0 && sign_a == sign_eps) ||
 149           (ymod4 == 2 && sign_a != sign_eps))
 150         mpfr_neg (mpc_imagref(z), mpc_imagref(z), GMP_RNDZ);
 151     }
 152
 153  end:
 154   mpz_clear (my);
 155 }
 156
 157 /* If x^y is exactly representable (with maybe a larger precision than z),
 158    round it in z and return the (mpc) inexact flag in [0, 10].
 159
 160    If x^y is not exactly representable, return -1.
 161
 162    If intermediate computations lead to numbers of more than maxprec bits,
 163    then abort and return -2 (in that case, to avoid loops, mpc_pow_exact
 164    should be called again with a larger value of maxprec).
 165
 166    Assume one of Re(x) or Im(x) is non-zero, and y is non-zero (y is real).
 167
 168    Warning: z and x might be the same variable, same for Re(z) or Im(z) and y.
 169
 170    In case -1 or -2 is returned, z is not modified.
 171 */
 172 static int
 173 mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd,
 174                mpfr_prec_t maxprec)
 175 {
 176   mpfr_exp_t ec, ed, ey;
 177   mpz_t my, a, b, c, d, u;
 178   unsigned long int t;
 179   int ret = -2;
 180   int sign_rex = mpfr_signbit (mpc_realref(x));
 181   int sign_imx = mpfr_signbit (mpc_imagref(x));
 182   int x_imag = mpfr_zero_p (mpc_realref(x));
 183   int z_is_y = 0;
 184   mpfr_t copy_of_y;
 185
 186   if (mpc_realref (z) == y || mpc_imagref (z) == y)
 187     {
 188       z_is_y = 1;
 189       mpfr_init2 (copy_of_y, mpfr_get_prec (y));
 190       mpfr_set (copy_of_y, y, GMP_RNDN);
 191     }
 192
 193   mpz_init (my);
 194   mpz_init (a);
 195   mpz_init (b);
 196   mpz_init (c);
 197   mpz_init (d);
 198   mpz_init (u);
 199
 200   ey = mpfr_get_z_exp (my, y);
 201   /* normalize so that my is odd */
 202   t = mpz_scan1 (my, 0);
 203   ey += (mpfr_exp_t) t;
 204   mpz_tdiv_q_2exp (my, my, t);
 205   /* y = my*2^ey with my odd */
 206
 207   if (x_imag)
 208     {
 209       mpz_set_ui (c, 0);
 210       ec = 0;
 211     }
 212   else
 213     ec = mpfr_get_z_exp (c, mpc_realref(x));
 214   if (mpfr_zero_p (mpc_imagref(x)))
 215     {
 216       mpz_set_ui (d, 0);
 217       ed = ec;
 218     }
 219   else
 220     {
 221       ed = mpfr_get_z_exp (d, mpc_imagref(x));
 222       if (x_imag)
 223         ec = ed;
 224     }
 225   /* x = c*2^ec + I * d*2^ed */
 226   /* equalize the exponents of x */
 227   if (ec < ed)
 228     {
 229       mpz_mul_2exp (d, d, (unsigned long int) (ed - ec));
 230       if ((mpfr_prec_t) mpz_sizeinbase (d, 2) > maxprec)
 231         goto end;
 232     }
 233   else if (ed < ec)
 234     {
 235       mpz_mul_2exp (c, c, (unsigned long int) (ec - ed));
 236       if ((mpfr_prec_t) mpz_sizeinbase (c, 2) > maxprec)
 237         goto end;
 238       ec = ed;
 239     }
 240   /* now ec=ed and x = (c + I * d) * 2^ec */
 241
 242   /* divide by two if possible */
 243   if (mpz_cmp_ui (c, 0) == 0)
 244     {
 245       t = mpz_scan1 (d, 0);
 246       mpz_tdiv_q_2exp (d, d, t);
 247       ec += (mpfr_exp_t) t;
 248     }
 249   else if (mpz_cmp_ui (d, 0) == 0)
 250     {
 251       t = mpz_scan1 (c, 0);
 252       mpz_tdiv_q_2exp (c, c, t);
 253       ec += (mpfr_exp_t) t;
 254     }
 255   else /* neither c nor d is zero */
 256     {
 257       unsigned long v;
 258       t = mpz_scan1 (c, 0);
 259       v = mpz_scan1 (d, 0);
 260       if (v < t)
 261         t = v;
 262       mpz_tdiv_q_2exp (c, c, t);
 263       mpz_tdiv_q_2exp (d, d, t);
 264       ec += (mpfr_exp_t) t;
 265     }
 266
 267   /* now either one of c, d is odd */
 268
 269   while (ey < 0)
 270     {
 271       /* check if x is a square */
 272       if (ec & 1)
 273         {
 274           mpz_mul_2exp (c, c, 1);
 275           mpz_mul_2exp (d, d, 1);
 276           ec --;
 277         }
 278       /* now ec is even */
 279       if (mpc_perfect_square_p (a, b, c, d) == 0)
 280         break;
 281       mpz_swap (a, c);
 282       mpz_swap (b, d);
 283       ec /= 2;
 284       ey ++;
 285     }
 286
 287   if (ey < 0)
 288     {
 289       ret = -1; /* not representable */
 290       goto end;
 291     }
 292
 293   /* Now ey >= 0, it thus suffices to check that x^my is representable.
 294      If my > 0, this is always true. If my < 0, we first try to invert
 295      (c+I*d)*2^ec.
 296   */
 297   if (mpz_cmp_ui (my, 0) < 0)
 298     {
 299       /* If my < 0, 1 / (c + I*d) = (c - I*d)/(c^2 + d^2), thus a sufficient
 300          condition is that c^2 + d^2 is a power of two, assuming |c| <> |d|.
 301          Assume a prime p <> 2 divides c^2 + d^2,
 302          then if p does not divide c or d, 1 / (c + I*d) cannot be exact.
 303          If p divides both c and d, then we can write c = p*c', d = p*d',
 304          and 1 / (c + I*d) = 1/p * 1/(c' + I*d'). This shows that if 1/(c+I*d)
 305          is exact, then 1/(c' + I*d') is exact too, and we are back to the
 306          previous case. In conclusion, a necessary and sufficient condition
 307          is that c^2 + d^2 is a power of two.
 308       */
 309       /* FIXME: we could first compute c^2+d^2 mod a limb for example */
 310       mpz_mul (a, c, c);
 311       mpz_addmul (a, d, d);
 312       t = mpz_scan1 (a, 0);
 313       if (mpz_sizeinbase (a, 2) != 1 + t) /* a is not a power of two */
 314         {
 315           ret = -1; /* not representable */
 316           goto end;
 317         }
 318       /* replace (c,d) by (c/(c^2+d^2), -d/(c^2+d^2)) */
 319       mpz_neg (d, d);
 320       ec = -ec - (mpfr_exp_t) t;
 321       mpz_neg (my, my);
 322     }
 323
 324   /* now ey >= 0 and my >= 0, and we want to compute
 325      [(c + I * d) * 2^ec] ^ (my * 2^ey).
 326
 327      We first compute [(c + I * d) * 2^ec]^my, then square ey times. */
 328   t = mpz_sizeinbase (my, 2) - 1;
 329   mpz_set (a, c);
 330   mpz_set (b, d);
 331   ed = ec;
 332   /* invariant: (a + I*b) * 2^ed = ((c + I*d) * 2^ec)^trunc(my/2^t) */
 333   while (t-- > 0)
 334     {
 335       unsigned long int v, w;
 336       /* square a + I*b */
 337       mpz_mul (u, a, b);
 338       mpz_mul (a, a, a);
 339       mpz_submul (a, b, b);
 340       mpz_mul_2exp (b, u, 1);
 341       ed *= 2;
 342       if (mpz_tstbit (my, t)) /* multiply by c + I*d */
 343         {
 344           mpz_mul (u, a, c);
 345           mpz_submul (u, b, d); /* ac-bd */
 346           mpz_mul (b, b, c);
 347           mpz_addmul (b, a, d); /* bc+ad */
 348           mpz_swap (a, u);
 349           ed += ec;
 350         }
 351       /* remove powers of two in (a,b) */
 352       if (mpz_cmp_ui (a, 0) == 0)
 353         {
 354           w = mpz_scan1 (b, 0);
 355           mpz_tdiv_q_2exp (b, b, w);
 356           ed += (mpfr_exp_t) w;
 357         }
 358       else if (mpz_cmp_ui (b, 0) == 0)
 359         {
 360           w = mpz_scan1 (a, 0);
 361           mpz_tdiv_q_2exp (a, a, w);
 362           ed += (mpfr_exp_t) w;
 363         }
 364       else
 365         {
 366           w = mpz_scan1 (a, 0);
 367           v = mpz_scan1 (b, 0);
 368           if (v < w)
 369             w = v;
 370           mpz_tdiv_q_2exp (a, a, w);
 371           mpz_tdiv_q_2exp (b, b, w);
 372           ed += (mpfr_exp_t) w;
 373         }
 374       if (   (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
 375           || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
 376         goto end;
 377     }
 378   /* now a+I*b = (c+I*d)^my */
 379
 380   while (ey-- > 0)
 381     {
 382       unsigned long sa, sb;
 383
 384       /* square a + I*b */
 385       mpz_mul (u, a, b);
 386       mpz_mul (a, a, a);
 387       mpz_submul (a, b, b);
 388       mpz_mul_2exp (b, u, 1);
 389       ed *= 2;
 390
 391       /* divide by largest 2^n possible, to avoid many loops for e.g.,
 392          (2+2*I)^16777216 */
 393       sa = mpz_scan1 (a, 0);
 394       sb = mpz_scan1 (b, 0);
 395       sa = (sa <= sb) ? sa : sb;
 396       mpz_tdiv_q_2exp (a, a, sa);
 397       mpz_tdiv_q_2exp (b, b, sa);
 398       ed += (mpfr_exp_t) sa;
 399
 400       if (   (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
 401           || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
 402         goto end;
 403     }
 404
 405   ret = mpfr_set_z (mpc_realref(z), a, MPC_RND_RE(rnd));
 406   ret = MPC_INEX(ret, mpfr_set_z (mpc_imagref(z), b, MPC_RND_IM(rnd)));
 407   mpfr_mul_2si (mpc_realref(z), mpc_realref(z), ed, MPC_RND_RE(rnd));
 408   mpfr_mul_2si (mpc_imagref(z), mpc_imagref(z), ed, MPC_RND_IM(rnd));
 409
 410  end:
 411   mpz_clear (my);
 412   mpz_clear (a);
 413   mpz_clear (b);
 414   mpz_clear (c);
 415   mpz_clear (d);
 416   mpz_clear (u);
 417
 418   if (ret >= 0 && x_imag)
 419     fix_sign (z, sign_rex, sign_imx, (z_is_y) ? copy_of_y : y);
 420
 421   if (z_is_y)
 422     mpfr_clear (copy_of_y);
 423
 424   return ret;
 425 }
 426
 427 /* Return 1 if y*2^k is an odd integer, 0 otherwise.
 428    Adapted from MPFR, file pow.c.
 429
 430    Examples: with k=0, check if y is an odd integer,
 431              with k=1, check if y is half-an-integer,
 432              with k=-1, check if y/2 is an odd integer.
 433 */
 434 #define MPFR_LIMB_HIGHBIT ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1))
 435 static int
 436 is_odd (mpfr_srcptr y, mpfr_exp_t k)
 437 {
 438   mpfr_exp_t expo;
 439   mpfr_prec_t prec;
 440   mp_size_t yn;
 441   mp_limb_t *yp;
 442
 443   expo = mpfr_get_exp (y) + k;
 444   if (expo <= 0)
 445     return 0;  /* |y| < 1 and not 0 */
 446
 447   prec = mpfr_get_prec (y);
 448   if ((mpfr_prec_t) expo > prec)
 449     return 0;  /* y is a multiple of 2^(expo-prec), thus not odd */
 450
 451   /* 0 < expo <= prec:
 452      y = 1xxxxxxxxxt.zzzzzzzzzzzzzzzzzz[000]
 453           expo bits   (prec-expo) bits
 454
 455      We have to check that:
 456      (a) the bit 't' is set
 457      (b) all the 'z' bits are zero
 458   */
 459
 460   prec = ((prec - 1) / BITS_PER_MP_LIMB + 1) * BITS_PER_MP_LIMB - expo;
 461   /* number of z+0 bits */
 462
 463   yn = prec / BITS_PER_MP_LIMB;
 464   /* yn is the index of limb containing the 't' bit */
 465
 466   yp = y->_mpfr_d;
 467   /* if expo is a multiple of BITS_PER_MP_LIMB, t is bit 0 */
 468   if (expo % BITS_PER_MP_LIMB == 0 ? (yp[yn] & 1) == 0
 469       : yp[yn] << ((expo % BITS_PER_MP_LIMB) - 1) != MPFR_LIMB_HIGHBIT)
 470     return 0;
 471   while (--yn >= 0)
 472     if (yp[yn] != 0)
 473       return 0;
 474   return 1;
 475 }
 476
 477 /* Put in z the value of x^y, rounded according to 'rnd'.
 478    Return the inexact flag in [0, 10]. */
 479 int
 480 mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
 481 {
 482   int ret = -2, loop, x_real, x_imag, y_real, z_real = 0, z_imag = 0;
 483   mpc_t t, u;
 484   mpfr_prec_t p, pr, pi, maxprec;
 485   int saved_underflow, saved_overflow;
 486
 487   /* save the underflow or overflow flags from MPFR */
 488   saved_underflow = mpfr_underflow_p ();
 489   saved_overflow = mpfr_overflow_p ();
 490
 491   x_real = mpfr_zero_p (mpc_imagref(x));
 492   y_real = mpfr_zero_p (mpc_imagref(y));
 493
 494   if (y_real && mpfr_zero_p (mpc_realref(y))) /* case y zero */
 495     {
 496       if (x_real && mpfr_zero_p (mpc_realref(x)))
 497         {
 498           /* we define 0^0 to be (1, +0) since the real part is
 499              coherent with MPFR where 0^0 gives 1, and the sign of the
 500              imaginary part cannot be determined                       */
 501           mpc_set_ui_ui (z, 1, 0, MPC_RNDNN);
 502           return 0;
 503         }
 504       else /* x^0 = 1 +/- i*0 even for x=NaN see algorithms.tex for the
 505               sign of zero */
 506         {
 507           mpfr_t n;
 508           int inex, cx1;
 509           int sign_zi;
 510           /* cx1 < 0 if |x| < 1
 511              cx1 = 0 if |x| = 1
 512              cx1 > 0 if |x| > 1
 513           */
 514           mpfr_init (n);
 515           inex = mpc_norm (n, x, GMP_RNDN);
 516           cx1 = mpfr_cmp_ui (n, 1);
 517           if (cx1 == 0 && inex != 0)
 518             cx1 = -inex;
 519
 520           sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
 521             || (cx1 == 0
 522                 && mpfr_signbit (mpc_imagref (x)) != mpfr_signbit (mpc_realref (y)))
 523             || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
 524
 525           /* warning: mpc_set_ui_ui does not set Im(z) to -0 if Im(rnd)=RNDD */
 526           ret = mpc_set_ui_ui (z, 1, 0, rnd);
 527
 528           if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi)
 529             mpc_conj (z, z, MPC_RNDNN);
 530
 531           mpfr_clear (n);
 532           return ret;
 533         }
 534     }
 535
 536   if (!mpc_fin_p (x) || !mpc_fin_p (y))
 537     {
 538       /* special values: exp(y*log(x)) */
 539       mpc_init2 (u, 2);
 540       mpc_log (u, x, MPC_RNDNN);
 541       mpc_mul (u, u, y, MPC_RNDNN);
 542       ret = mpc_exp (z, u, rnd);
 543       mpc_clear (u);
 544       goto end;
 545     }
 546
 547   if (x_real) /* case x real */
 548     {
 549       if (mpfr_zero_p (mpc_realref(x))) /* x is zero */
 550         {
 551           /* special values: exp(y*log(x)) */
 552           mpc_init2 (u, 2);
 553           mpc_log (u, x, MPC_RNDNN);
 554           mpc_mul (u, u, y, MPC_RNDNN);
 555           ret = mpc_exp (z, u, rnd);
 556           mpc_clear (u);
 557           goto end;
 558         }
 559
 560       /* Special case 1^y = 1 */
 561       if (mpfr_cmp_ui (mpc_realref(x), 1) == 0)
 562         {
 563           int s1, s2;
 564           s1 = mpfr_signbit (mpc_realref (y));
 565           s2 = mpfr_signbit (mpc_imagref (x));
 566
 567           ret = mpc_set_ui (z, +1, rnd);
 568           /* the sign of the zero imaginary part is known in some cases (see
 569              algorithm.tex). In such cases we have
 570              (x +s*0i)^(y+/-0i) = x^y + s*sign(y)*0i
 571              where s = +/-1.  We extend here this rule to fix the sign of the
 572              zero part.
 573
 574              Note that the sign must also be set explicitly when rnd=RNDD
 575              because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0.
 576           */
 577           if (MPC_RND_IM (rnd) == GMP_RNDD || s1 != s2)
 578             mpc_conj (z, z, MPC_RNDNN);
 579           goto end;
 580         }
 581
 582       /* x^y is real when:
 583          (a) x is real and y is integer
 584          (b) x is real non-negative and y is real */
 585       if (y_real && (mpfr_integer_p (mpc_realref(y)) ||
 586                      mpfr_cmp_ui (mpc_realref(x), 0) >= 0))
 587         {
 588           int s1, s2;
 589           s1 = mpfr_signbit (mpc_realref (y));
 590           s2 = mpfr_signbit (mpc_imagref (x));
 591
 592           ret = mpfr_pow (mpc_realref(z), mpc_realref(x), mpc_realref(y), MPC_RND_RE(rnd));
 593           ret = MPC_INEX(ret, mpfr_set_ui (mpc_imagref(z), 0, MPC_RND_IM(rnd)));
 594
 595           /* the sign of the zero imaginary part is known in some cases
 596              (see algorithm.tex). In such cases we have (x +s*0i)^(y+/-0i)
 597              = x^y + s*sign(y)*0i where s = +/-1.
 598              We extend here this rule to fix the sign of the zero part.
 599
 600              Note that the sign must also be set explicitly when rnd=RNDD
 601              because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0.
 602           */
 603           if (MPC_RND_IM(rnd) == GMP_RNDD || s1 != s2)
 604             mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPC_RND_IM(rnd));
 605           goto end;
 606         }
 607
 608       /* (-1)^(n+I*t) is real for n integer and t real */
 609       if (mpfr_cmp_si (mpc_realref(x), -1) == 0 && mpfr_integer_p (mpc_realref(y)))
 610         z_real = 1;
 611
 612       /* for x real, x^y is imaginary when:
 613          (a) x is negative and y is half-an-integer
 614          (b) x = -1 and Re(y) is half-an-integer
 615       */
 616       if ((mpfr_cmp_ui (mpc_realref(x), 0) < 0) && is_odd (mpc_realref(y), 1)
 617          && (y_real || mpfr_cmp_si (mpc_realref(x), -1) == 0))
 618         z_imag = 1;
 619     }
 620   else /* x non real */
 621     /* I^(t*I) and (-I)^(t*I) are real for t real,
 622        I^(n+t*I) and (-I)^(n+t*I) are real for n even and t real, and
 623        I^(n+t*I) and (-I)^(n+t*I) are imaginary for n odd and t real
 624        (s*I)^n is real for n even and imaginary for n odd */
 625     if ((mpc_cmp_si_si (x, 0, 1) == 0 || mpc_cmp_si_si (x, 0, -1) == 0 ||
 626          (mpfr_cmp_ui (mpc_realref(x), 0) == 0 && y_real)) &&
 627         mpfr_integer_p (mpc_realref(y)))
 628       { /* x is I or -I, and Re(y) is an integer */
 629         if (is_odd (mpc_realref(y), 0))
 630           z_imag = 1; /* Re(y) odd: z is imaginary */
 631         else
 632           z_real = 1; /* Re(y) even: z is real */
 633       }
 634     else /* (t+/-t*I)^(2n) is imaginary for n odd and real for n even */
 635       if (mpfr_cmpabs (mpc_realref(x), mpc_imagref(x)) == 0 && y_real &&
 636           mpfr_integer_p (mpc_realref(y)) && is_odd (mpc_realref(y), 0) == 0)
 637         {
 638           if (is_odd (mpc_realref(y), -1)) /* y/2 is odd */
 639             z_imag = 1;
 640           else
 641             z_real = 1;
 642         }
 643
 644   pr = mpfr_get_prec (mpc_realref(z));
 645   pi = mpfr_get_prec (mpc_imagref(z));
 646   p = (pr > pi) ? pr : pi;
 647   p += 12; /* experimentally, seems to give less than 10% of failures in
 648               Ziv's strategy; probably wrong now since q is not computed */
 649   if (p < 64)
 650     p = 64;
 651   mpc_init2 (u, p);
 652   mpc_init2 (t, p);
 653   pr += MPC_RND_RE(rnd) == GMP_RNDN;
 654   pi += MPC_RND_IM(rnd) == GMP_RNDN;
 655   maxprec = MPC_MAX_PREC (z);
 656   x_imag = mpfr_zero_p (mpc_realref(x));
 657   for (loop = 0;; loop++)
 658     {
 659       int ret_exp;
 660       mpfr_exp_t dr, di;
 661       mpfr_prec_t q;
 662
 663       mpc_log (t, x, MPC_RNDNN);
 664       mpc_mul (t, t, y, MPC_RNDNN);
 665
 666       /* Compute q such that |Re (y log x)|, |Im (y log x)| < 2^q.
 667          We recompute it at each loop since we might get different
 668          bounds if the precision is not enough. */
 669       q = mpfr_get_exp (mpc_realref(t)) > 0 ? mpfr_get_exp (mpc_realref(t)) : 0;
 670       if (mpfr_get_exp (mpc_imagref(t)) > (mpfr_exp_t) q)
 671         q = mpfr_get_exp (mpc_imagref(t));
 672
 673       mpfr_clear_overflow ();
 674       mpfr_clear_underflow ();
 675       ret_exp = mpc_exp (u, t, MPC_RNDNN);
 676       if (mpfr_underflow_p () || mpfr_overflow_p ()) {
 677          /* under- and overflow flags are set by mpc_exp */
 678          mpc_set (z, u, MPC_RNDNN);
 679          ret = ret_exp;
 680          goto exact;
 681       }
 682
 683       /* Since the error bound is global, we have to take into account the
 684          exponent difference between the real and imaginary parts. We assume
 685          either the real or the imaginary part of u is not zero.
 686       */
 687       dr = mpfr_zero_p (mpc_realref(u)) ? mpfr_get_exp (mpc_imagref(u))
 688         : mpfr_get_exp (mpc_realref(u));
 689       di = mpfr_zero_p (mpc_imagref(u)) ? dr : mpfr_get_exp (mpc_imagref(u));
 690       if (dr > di)
 691         {
 692           di = dr - di;
 693           dr = 0;
 694         }
 695       else
 696         {
 697           dr = di - dr;
 698           di = 0;
 699         }
 700       /* the term -3 takes into account the factor 4 in the complex error
 701          (see algorithms.tex) plus one due to the exponent difference: if
 702          z = a + I*b, where the relative error on z is at most 2^(-p), and
 703          EXP(a) = EXP(b) + k, the relative error on b is at most 2^(k-p) */
 704       if ((z_imag || (p > q + 3 + dr && mpfr_can_round (mpc_realref(u), p - q - 3 - dr, GMP_RNDN, GMP_RNDZ, pr))) &&
 705           (z_real || (p > q + 3 + di && mpfr_can_round (mpc_imagref(u), p - q - 3 - di, GMP_RNDN, GMP_RNDZ, pi))))
 706         break;
 707
 708       /* if Re(u) is not known to be zero, assume it is a normal number, i.e.,
 709          neither zero, Inf or NaN, otherwise we might enter an infinite loop */
 710       MPC_ASSERT (z_imag || mpfr_number_p (mpc_realref(u)));
 711       /* idem for Im(u) */
 712       MPC_ASSERT (z_real || mpfr_number_p (mpc_imagref(u)));
 713
 714       if (ret == -2) /* we did not yet call mpc_pow_exact, or it aborted
 715                         because intermediate computations had > maxprec bits */
 716         {
 717           /* check exact cases (see algorithms.tex) */
 718           if (y_real)
 719             {
 720               maxprec *= 2;
 721               ret = mpc_pow_exact (z, x, mpc_realref(y), rnd, maxprec);
 722               if (ret != -1 && ret != -2)
 723                 goto exact;
 724             }
 725           p += dr + di + 64;
 726         }
 727       else
 728         p += p / 2;
 729       mpc_set_prec (t, p);
 730       mpc_set_prec (u, p);
 731     }
 732
 733   if (z_real)
 734     {
 735       /* When the result is real (see algorithm.tex for details),
 736          Im(x^y) =
 737          + sign(imag(y))*0i,               if |x| > 1
 738          + sign(imag(x))*sign(real(y))*0i, if |x| = 1
 739          - sign(imag(y))*0i,               if |x| < 1
 740       */
 741       mpfr_t n;
 742       int inex, cx1;
 743       int sign_zi, sign_rex, sign_imx;
 744       /* cx1 < 0 if |x| < 1
 745          cx1 = 0 if |x| = 1
 746          cx1 > 0 if |x| > 1
 747       */
 748
 749       sign_rex = mpfr_signbit (mpc_realref (x));
 750       sign_imx = mpfr_signbit (mpc_imagref (x));
 751       mpfr_init (n);
 752       inex = mpc_norm (n, x, GMP_RNDN);
 753       cx1 = mpfr_cmp_ui (n, 1);
 754       if (cx1 == 0 && inex != 0)
 755         cx1 = -inex;
 756
 757       sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
 758         || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y)))
 759         || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
 760
 761       /* copy RE(y) to n since if z==y we will destroy Re(y) below */
 762       mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y)));
 763       mpfr_set (n, mpc_realref (y), GMP_RNDN);
 764       ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd));
 765       if (y_real && (x_real || x_imag))
 766         {
 767           /* FIXME: with y_real we assume Im(y) is really 0, which is the case
 768              for example when y comes from pow_fr, but in case Im(y) is +0 or
 769              -0, we might get different results */
 770           mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
 771           fix_sign (z, sign_rex, sign_imx, n);
 772           ret = MPC_INEX(ret, 0); /* imaginary part is exact */
 773         }
 774       else
 775         {
 776           ret = MPC_INEX (ret, mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd)));
 777           /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */
 778           if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi)
 779             mpc_conj (z, z, MPC_RNDNN);
 780         }
 781
 782       mpfr_clear (n);
 783     }
 784   else if (z_imag)
 785     {
 786       ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd));
 787       /* if z is imaginary and y real, then x cannot be real */
 788       if (y_real && x_imag)
 789         {
 790           int sign_rex = mpfr_signbit (mpc_realref (x));
 791
 792           /* If z overlaps with y we set Re(z) before checking Re(y) below,
 793              but in that case y=0, which was dealt with above. */
 794           mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd));
 795           /* Note: fix_sign only does something when y is an integer,
 796              then necessarily y = 1 or 3 (mod 4), and in that case the
 797              sign of Im(x) is irrelevant. */
 798           fix_sign (z, sign_rex, 0, mpc_realref (y));
 799           ret = MPC_INEX(0, ret);
 800         }
 801       else
 802         ret = MPC_INEX(mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd)), ret);
 803     }
 804   else
 805     ret = mpc_set (z, u, rnd);
 806  exact:
 807   mpc_clear (t);
 808   mpc_clear (u);
 809
 810   /* restore underflow and overflow flags from MPFR */
 811   if (saved_underflow)
 812     mpfr_set_underflow ();
 813   if (saved_overflow)
 814     mpfr_set_overflow ();
 815
 816  end:
 817   return ret;
 818 }