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=0;
 662       /* to avoid warning message, real initialisation below */
 663
 664       mpc_log (t, x, MPC_RNDNN);
 665       mpc_mul (t, t, y, MPC_RNDNN);
 666
 667       if (loop == 0) {
 668          /* compute q such that |Re (y log x)|, |Im (y log x)| < 2^q */
 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
 674       mpfr_clear_overflow ();
 675       mpfr_clear_underflow ();
 676       ret_exp = mpc_exp (u, t, MPC_RNDNN);
 677       if (mpfr_underflow_p () || mpfr_overflow_p ()) {
 678          /* under- and overflow flags are set by mpc_exp */
 679          mpc_set (z, u, MPC_RNDNN);
 680          ret = ret_exp;
 681          goto exact;
 682       }
 683
 684       /* Since the error bound is global, we have to take into account the
 685          exponent difference between the real and imaginary parts. We assume
 686          either the real or the imaginary part of u is not zero.
 687       */
 688       dr = mpfr_zero_p (mpc_realref(u)) ? mpfr_get_exp (mpc_imagref(u))
 689         : mpfr_get_exp (mpc_realref(u));
 690       di = mpfr_zero_p (mpc_imagref(u)) ? dr : mpfr_get_exp (mpc_imagref(u));
 691       if (dr > di)
 692         {
 693           di = dr - di;
 694           dr = 0;
 695         }
 696       else
 697         {
 698           dr = di - dr;
 699           di = 0;
 700         }
 701       /* the term -3 takes into account the factor 4 in the complex error
 702          (see algorithms.tex) plus one due to the exponent difference: if
 703          z = a + I*b, where the relative error on z is at most 2^(-p), and
 704          EXP(a) = EXP(b) + k, the relative error on b is at most 2^(k-p) */
 705       if ((z_imag || (p > q + 3 + dr && mpfr_can_round (mpc_realref(u), p - q - 3 - dr, GMP_RNDN, GMP_RNDZ, pr))) &&
 706           (z_real || (p > q + 3 + di && mpfr_can_round (mpc_imagref(u), p - q - 3 - di, GMP_RNDN, GMP_RNDZ, pi))))
 707         break;
 708
 709       /* if Re(u) is not known to be zero, assume it is a normal number, i.e.,
 710          neither zero, Inf or NaN, otherwise we might enter an infinite loop */
 711       MPC_ASSERT (z_imag || mpfr_number_p (mpc_realref(u)));
 712       /* idem for Im(u) */
 713       MPC_ASSERT (z_real || mpfr_number_p (mpc_imagref(u)));
 714
 715       if (ret == -2) /* we did not yet call mpc_pow_exact, or it aborted
 716                         because intermediate computations had > maxprec bits */
 717         {
 718           /* check exact cases (see algorithms.tex) */
 719           if (y_real)
 720             {
 721               maxprec *= 2;
 722               ret = mpc_pow_exact (z, x, mpc_realref(y), rnd, maxprec);
 723               if (ret != -1 && ret != -2)
 724                 goto exact;
 725             }
 726           p += dr + di + 64;
 727         }
 728       else
 729         p += p / 2;
 730       mpc_set_prec (t, p);
 731       mpc_set_prec (u, p);
 732     }
 733
 734   if (z_real)
 735     {
 736       /* When the result is real (see algorithm.tex for details),
 737          Im(x^y) =
 738          + sign(imag(y))*0i,               if |x| > 1
 739          + sign(imag(x))*sign(real(y))*0i, if |x| = 1
 740          - sign(imag(y))*0i,               if |x| < 1
 741       */
 742       mpfr_t n;
 743       int inex, cx1;
 744       int sign_zi, sign_rex, sign_imx;
 745       /* cx1 < 0 if |x| < 1
 746          cx1 = 0 if |x| = 1
 747          cx1 > 0 if |x| > 1
 748       */
 749
 750       sign_rex = mpfr_signbit (mpc_realref (x));
 751       sign_imx = mpfr_signbit (mpc_imagref (x));
 752       mpfr_init (n);
 753       inex = mpc_norm (n, x, GMP_RNDN);
 754       cx1 = mpfr_cmp_ui (n, 1);
 755       if (cx1 == 0 && inex != 0)
 756         cx1 = -inex;
 757
 758       sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
 759         || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y)))
 760         || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
 761
 762       /* copy RE(y) to n since if z==y we will destroy Re(y) below */
 763       mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y)));
 764       mpfr_set (n, mpc_realref (y), GMP_RNDN);
 765       ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd));
 766       if (y_real && (x_real || x_imag))
 767         {
 768           /* FIXME: with y_real we assume Im(y) is really 0, which is the case
 769              for example when y comes from pow_fr, but in case Im(y) is +0 or
 770              -0, we might get different results */
 771           mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
 772           fix_sign (z, sign_rex, sign_imx, n);
 773           ret = MPC_INEX(ret, 0); /* imaginary part is exact */
 774         }
 775       else
 776         {
 777           ret = MPC_INEX (ret, mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd)));
 778           /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */
 779           if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi)
 780             mpc_conj (z, z, MPC_RNDNN);
 781         }
 782
 783       mpfr_clear (n);
 784     }
 785   else if (z_imag)
 786     {
 787       ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd));
 788       /* if z is imaginary and y real, then x cannot be real */
 789       if (y_real && x_imag)
 790         {
 791           int sign_rex = mpfr_signbit (mpc_realref (x));
 792
 793           /* If z overlaps with y we set Re(z) before checking Re(y) below,
 794              but in that case y=0, which was dealt with above. */
 795           mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd));
 796           /* Note: fix_sign only does something when y is an integer,
 797              then necessarily y = 1 or 3 (mod 4), and in that case the
 798              sign of Im(x) is irrelevant. */
 799           fix_sign (z, sign_rex, 0, mpc_realref (y));
 800           ret = MPC_INEX(0, ret);
 801         }
 802       else
 803         ret = MPC_INEX(mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd)), ret);
 804     }
 805   else
 806     ret = mpc_set (z, u, rnd);
 807  exact:
 808   mpc_clear (t);
 809   mpc_clear (u);
 810
 811   /* restore underflow and overflow flags from MPFR */
 812   if (saved_underflow)
 813     mpfr_set_underflow ();
 814   if (saved_overflow)
 815     mpfr_set_overflow ();
 816
 817  end:
 818   return ret;
 819 }