src/pow.c - master - niceSOFT/mpc

src/pow.c
977 строк32 KB
Paul Zimmermann
avoid another compiler warning
22 янв 2026, 13:51
22 янв 2026, 13:51c2fd42a
/* mpc_pow -- Raise a complex number to the power of another complex number.

Copyright (C) 2009, 2010, 2011, 2012, 2014, 2015, 2016, 2018, 2020, 2022, 2024, 2025 INRIA

This file is part of GNU MPC.

GNU MPC is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
more details.

You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see http://www.gnu.org/licenses/ .
*/

#include <stdio.h> /* for MPC_ASSERT */
#include "mpc-impl.h"

/* Return non-zero iff c+i*d is an exact square (a+i*b)^2,
   with a, b both of the form m*2^e with m, e integers.
   If so, returns in a+i*b the corresponding square root, with a >= 0.
   The variables a, b must not overlap with c, d.

   We have c = a^2 - b^2 and d = 2*a*b.

   If one of a, b is exact, then both are (see algorithms.tex).

   Case 1: a <> 0 and b <> 0.
   Let a = m*2^e and b = n*2^f with m, e, n, f integers, m and n odd
   (we will treat apart the case a = 0 or b = 0).
   Then 2*a*b = m*n*2^(e+f+1), thus necessarily e+f >= -1.
   Assume e < 0, then f >= 0, then a^2 - b^2 = m^2*2^(2e) - n^2*2^(2f) cannot
   be an integer, since n^2*2^(2f) is an integer, and m^2*2^(2e) is not.
   Similarly when f < 0 (and thus e >= 0).
   Thus we have e, f >= 0, and a, b are both integers.
   Let A = 2a^2, then eliminating b between c = a^2 - b^2 and d = 2*a*b
   gives A^2 - 2c*A - d^2 = 0, which has solutions c +/- sqrt(c^2+d^2).
   We thus need c^2+d^2 to be a square, and c + sqrt(c^2+d^2) --- the solution
   we are interested in --- to be two times a square. Then b = d/(2a) is
   necessarily an integer.

   Case 2: a = 0. Then d is necessarily zero, thus it suffices to check
   whether c = -b^2, i.e., if -c is a square.

   Case 3: b = 0. Then d is necessarily zero, thus it suffices to check
   whether c = a^2, i.e., if c is a square.

   *sc and *sd are the input signs of c and d, update them as output with
   the signs of a and b (in case of an exact square).
*/
static int
mpc_perfect_square_p (mpz_t a, mpz_t b, mpz_t c, mpz_t d, int *sc, int *sd)
{
  if (mpz_cmp_ui (d, 0) == 0) /* case a = 0 or b = 0 */
    {
      /* necessarily c < 0 here, since we have already considered the case
         where x is real non-negative and y is real */
      MPC_ASSERT (mpz_cmp_ui (c, 0) < 0);
      mpz_neg (b, c);
      if (mpz_perfect_square_p (b)) /* case 2 above */
        {
          mpz_sqrt (b, b);
          if (*sd < 0)
            mpz_neg (b, b); /* sqrt(-b^2 - 0*i) = +0 - i*b */
          mpz_set_ui (a, 0);
          *sc = 1; /* a = +0 */
          /* we don't set *sd, since it had the correct sign */
          return 1; /* c + i*d = (0 + i*b)^2 */
        }
    }
  else /* both a and b are non-zero */
    {
      if (mpz_divisible_2exp_p (d, 1) == 0)
        return 0; /* d must be even */
      mpz_mul (a, c, c);
      mpz_addmul (a, d, d); /* c^2 + d^2 */
      if (mpz_perfect_square_p (a))
        {
          mpz_sqrt (a, a);
          mpz_add (a, c, a); /* c + sqrt(c^2+d^2) */
          if (mpz_divisible_2exp_p (a, 1))
            {
              mpz_tdiv_q_2exp (a, a, 1);
              if (mpz_perfect_square_p (a))
                {
                  mpz_sqrt (a, a);
                  *sc = 1; /* a > 0 */
                  mpz_tdiv_q_2exp (b, d, 1); /* d/2 */
                  mpz_divexact (b, b, a); /* d/(2a) */
                  /* sign of b is that of d */
                  return 1;
                }
            }
        }
    }
  return 0; /* not a square */
}

/* fix the sign of Re(z) or Im(z) in case it is zero,
   and Re(x) is zero.
   sign_eps is 0 if Re(x) = +0, 1 if Re(x) = -0
   sign_a is the sign bit of Im(x).
   Assume y is an integer (does nothing otherwise).
*/
static void
fix_sign (mpc_ptr z, int sign_eps, int sign_a, mpfr_srcptr y)
{
  int ymod4 = -1;
  mpfr_exp_t ey;
  mpz_t my;
  unsigned long int t;

  mpz_init (my);

  ey = mpfr_get_z_exp (my, y);
  /* normalize so that my is odd */
  t = mpz_scan1 (my, 0);
  ey += (mpfr_exp_t) t;
  mpz_tdiv_q_2exp (my, my, t);
  /* y = my*2^ey */

  /* compute y mod 4 (in case y is an integer) */
  if (ey >= 2)
    ymod4 = 0;
  else if (ey == 1)
    ymod4 = mpz_tstbit (my, 0) * 2; /* correct if my < 0 */
  else if (ey == 0)
    {
      ymod4 = mpz_tstbit (my, 1) * 2 + mpz_tstbit (my, 0);
      if (mpz_cmp_ui (my , 0) < 0)
        ymod4 = 4 - ymod4;
    }
  else /* y is not an integer */
    goto end;

  if (mpfr_zero_p (mpc_realref(z)))
    {
      /* we assume y is always integer in that case (FIXME: prove it):
         (eps+I*a)^y = +0 + I*a^y for y = 1 mod 4 and sign_eps = 0
         (eps+I*a)^y = -0 - I*a^y for y = 3 mod 4 and sign_eps = 0 */
      MPC_ASSERT (ymod4 == 1 || ymod4 == 3);
      if ((ymod4 == 3 && sign_eps == 0) ||
          (ymod4 == 1 && sign_eps == 1))
        mpfr_neg (mpc_realref(z), mpc_realref(z), MPFR_RNDZ);
    }
  else if (mpfr_zero_p (mpc_imagref(z)))
    {
      /* we assume y is always integer in that case (FIXME: prove it):
         (eps+I*a)^y =  a^y - 0*I for y = 0 mod 4 and sign_a = sign_eps
         (eps+I*a)^y =  -a^y +0*I for y = 2 mod 4 and sign_a = sign_eps */
      MPC_ASSERT (ymod4 == 0 || ymod4 == 2);
      if ((ymod4 == 0 && sign_a == sign_eps) ||
          (ymod4 == 2 && sign_a != sign_eps))
        mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPFR_RNDZ);
    }

 end:
  mpz_clear (my);
}

/* fix the sign of Im(z) when both x and y are real
   s1 = signbit(Re(x)), s2 = signbit(Im(x)),
   s3 = signbit(Re(y)), s4 = signbit(Im(y))
   cmp = mpfr_cmpabs (Re(x), 1)
   even = 1 iff Re(y) is even (only used when Re(x) < 0)
   See algorithms.tex */
static void
fix_sign_real (mpc_ptr z, int s1, int s2, int s3, int s4, int cmp, int even)
{
  /* Re(x) Im(x) Re(y) Im(y) Im(z)
  (a) +     +     +     +     +
  (b) +     +     +     -     +
  (c) +     +     -     +     - when |Re(x)| <= 1, + otherwise
  (d) +     +     -     -     + when |Re(x)| <  1, - otherwise
  (e) +     -     +     +     - when |Re(x)| <= 1, + otherwise
  (f) +     -     +     -     + when |Re(x)| <  1, - otherwise
  (g) +     -     -     +     +
  (h) +     -     -     -     +
  (i) -     +     +     +     - when |Re(x)| <= 1 and Re(y) even, + otherwise
  (j) -     +     +     -     - when |Re(x)| >= 1 and Re(y) even, + otherwise
  (k) -     +     -     +     - when |Re(x)| >= 1 and Re(y)  odd, + otherwise
  (l) -     +     -     -     - when |Re(x)| <= 1 and Re(y)  odd, + otherwise
  (m) -     -     +     +     - when |Re(x)| >= 1 and Re(y)  odd, + otherwise
  (n) -     -     +     -     - when |Re(x)| <= 1 and Re(y)  odd, + otherwise
  (o) -     -     -     +     - when |Re(x)| <= 1 and Re(y) even, + otherwise
  (p) -     -     -     -     - when |Re(x)| >= 1 and Re(y) even, + otherwise
  */
  int s; /* signbit of Im(z) (0 for sign(Im(z))=+1, 1 for sign(Im(z))=-1) */
  if (s1 == 0) {
    if (s2 == s3)     /* rules (a), (b), (g), (h) */
      s = 0;
    else if (s4 == 0) /* rules (c) and (e) */
      s = (cmp <= 0) ? 1 : 0;
    else              /* rules (d) and (f) */
      s = (cmp < 0) ? 0 : 1;
    }
  else {
    if (even) {
      if (s4 == 0)    /* rules (i) and (o) */
        s = (cmp <= 0) ? 1 : 0;
      else            /* rules (j) and (p) */
        s = (cmp >= 0) ? 1 : 0;
    }
    else {
      if (s4 == 0)    /* rules (k) and (m) */
        s = (cmp >= 0) ? 1 : 0;
      else            /* rules (l) and (n) */
        s = (cmp <= 0) ? 1 : 0;
    }
  }
  if (mpfr_signbit (mpc_imagref (z)) != s)
    mpfr_neg (mpc_imagref (z), mpc_imagref (z), MPFR_RNDN);
}

/* If x^y is exactly representable (with maybe a larger precision than z),
   round it in z and return the (mpc) inexact flag in [0, 10].

   If x^y is not exactly representable, return -1.

   If intermediate computations lead to numbers of more than maxprec bits,
   then abort and return -2 (in that case, to avoid loops, mpc_pow_exact
   should be called again with a larger value of maxprec).

   Assume one of Re(x) or Im(x) is non-zero, and y is non-zero (y is real).

   Warning: z and x might be the same variable, same for Re(z) or Im(z) and y.

   In case -1 or -2 is returned, z is not modified.
*/
static int
mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd,
               mpfr_prec_t maxprec)
{
  mpfr_exp_t ec, ed, ey;
  mpz_t my, a, b, c, d, u;
  int sa, sb, sc, sd;
  unsigned long int t;
  int ret = -2;
  int x_imag = mpfr_zero_p (mpc_realref(x));
  int z_is_y = 0;
  mpfr_t copy_of_y;
  int inex_im;

  if (mpc_realref (z) == y || mpc_imagref (z) == y)
    {
      z_is_y = 1;
      mpfr_init2 (copy_of_y, mpfr_get_prec (y));
      mpfr_set (copy_of_y, y, MPFR_RNDN);
    }

  mpz_init (my);
  mpz_init (a);
  mpz_init (b);
  mpz_init (c);
  mpz_init (d);
  mpz_init (u);

  ey = mpfr_get_z_exp (my, y);
  /* normalize so that my is odd */
  t = mpz_scan1 (my, 0);
  ey += (mpfr_exp_t) t;
  mpz_tdiv_q_2exp (my, my, t);
  /* y = my*2^ey with my odd */

  if (x_imag)
    {
      mpz_set_ui (c, 0);
      ec = 0;
    }
  else
    ec = mpfr_get_z_exp (c, mpc_realref(x));
  sc = mpfr_signbit (mpc_realref(x)) ? -1 : 1;
  if (mpfr_zero_p (mpc_imagref(x)))
    {
      mpz_set_ui (d, 0);
      ed = ec;
    }
  else
    {
      ed = mpfr_get_z_exp (d, mpc_imagref(x));
      if (x_imag)
        ec = ed;
    }
  sd = mpfr_signbit (mpc_imagref(x)) ? -1 : 1;

  /* x = c*2^ec + I * d*2^ed, where sc/sd are the signs of c, d
     (only valuable when c or d is zero, since we lost the sign bit
     in the conversion to integer) */

  /* equalize the exponents of x */
  if (ec < ed)
    {
      mpz_mul_2exp (d, d, (unsigned long int) (ed - ec));
      if ((mpfr_prec_t) mpz_sizeinbase (d, 2) > maxprec)
        goto end;
    }
  else if (ed < ec)
    {
      mpz_mul_2exp (c, c, (unsigned long int) (ec - ed));
      if ((mpfr_prec_t) mpz_sizeinbase (c, 2) > maxprec)
        goto end;
      ec = ed;
    }

  /* now ec=ed and x = (c + I * d) * 2^ec */

  /* divide by two if possible */
  if (mpz_cmp_ui (c, 0) == 0)
    {
      t = mpz_scan1 (d, 0);
      mpz_tdiv_q_2exp (d, d, t);
      ec += (mpfr_exp_t) t;
    }
  else if (mpz_cmp_ui (d, 0) == 0)
    {
      t = mpz_scan1 (c, 0);
      mpz_tdiv_q_2exp (c, c, t);
      ec += (mpfr_exp_t) t;
    }
  else /* neither c nor d is zero */
    {
      unsigned long v;
      t = mpz_scan1 (c, 0);
      v = mpz_scan1 (d, 0);
      if (v < t)
        t = v;
      mpz_tdiv_q_2exp (c, c, t);
      mpz_tdiv_q_2exp (d, d, t);
      ec += (mpfr_exp_t) t;
    }

  /* now either one of c, d is odd */

  while (ey < 0)
    {
      /* check if x is a square */
      if (ec & 1)
        {
          mpz_mul_2exp (c, c, 1);
          mpz_mul_2exp (d, d, 1);
          ec --;
        }
      /* now ec is even */
      if (mpc_perfect_square_p (a, b, c, d, &sc, &sd) == 0)
        break;
      mpz_swap (a, c);
      mpz_swap (b, d);
      ec /= 2;
      ey ++;
    }

  if (ey < 0)
    {
      ret = -1; /* not representable */
      goto end;
    }

  /* Now ey >= 0, it thus suffices to check that x^my is representable.
     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) = (c - I*d)/(c^2 + d^2), thus a sufficient
         condition is that c^2 + d^2 is a power of two, assuming |c| <> |d|.
         Assume a prime p <> 2 divides c^2 + d^2,
         then if p does not divide c or d, 1 / (c + I*d) cannot be exact.
         If p divides both c and d, then we can write c = p*c', d = p*d',
         and 1 / (c + I*d) = 1/p * 1/(c' + I*d'). This shows that if 1/(c+I*d)
         is exact, then 1/(c' + I*d') is exact too, and we are back to the
         previous case. In conclusion, a necessary and sufficient condition
         is that c^2 + d^2 is a power of two.
      */
      /* FIXME: we could first compute c^2+d^2 mod a limb for example */
      mpz_mul (a, c, c);
      mpz_addmul (a, d, d);
      t = mpz_scan1 (a, 0);
      if (mpz_sizeinbase (a, 2) != 1 + t) /* a is not a power of two */
        {
          ret = -1; /* not representable */
          goto end;
        }
      /* replace (c,d) by (c/(c^2+d^2), -d/(c^2+d^2)) */
      mpz_neg (d, d);
      ec = -ec - (mpfr_exp_t) t;
      sd = -sd; /* only the sign of d changes */
      mpz_neg (my, my);
    }

  /* now ey >= 0 and my >= 0, and we want to compute
     [(c + I * d) * 2^ec] ^ (my * 2^ey).

     We first compute [(c + I * d) * 2^ec]^my, then square ey times. */
  t = mpz_sizeinbase (my, 2) - 1;
  mpz_set (a, c);
  sa = sc;
  mpz_set (b, d);
  sb = sd;
  ed = ec;
  /* invariant: (a + I*b) * 2^ed = ((c + I*d) * 2^ec)^trunc(my/2^t) */
  while (t-- > 0)
    {
      unsigned long int v, w;
      /* square a + I*b */
      mpz_mul (u, a, b);
      mpz_mul (a, a, a);
      mpz_submul (a, b, b);
      /* The case a=0 cannot happen when the initial x is real or imaginary,
         which is the only case where the sign of zeros matters.
         Thus the sign of a (in case a=0) does not really matter. */
      mpz_mul_2exp (b, u, 1);
      sb = sa * sb;
      sa = mpz_cmp_ui (a, 0) >= 0 ? +1 : -1;
      ed *= 2;
      if (mpz_tstbit (my, t)) /* multiply by c + I*d */
        {
          int ta = sa * sc > 0 && sb * sd < 0;
          int tb = sb * sc > 0 && sa * sd < 0;
          mpz_mul (u, a, c);
          mpz_submul (u, b, d); /* ac-bd */
          mpz_mul (b, b, c);
          mpz_addmul (b, a, d); /* bc+ad */
          mpz_swap (a, u);
          if (mpz_cmp_ui (a, 0) != 0)
            sa = mpz_cmp_ui (a, 0) > 0 ? +1 : -1;
          else
            sa = (ta) ? +1 : -1;
          if (mpz_cmp_ui (b, 0) != 0)
            sb = mpz_cmp_ui (b, 0) > 0 ? +1 : -1;
          else
            sa = (tb) ? +1 : -1;
          ed += ec;
        }
      /* remove powers of two in (a,b) */
      if (mpz_cmp_ui (a, 0) == 0)
        {
          w = mpz_scan1 (b, 0);
          mpz_tdiv_q_2exp (b, b, w);
          ed += (mpfr_exp_t) w;
        }
      else if (mpz_cmp_ui (b, 0) == 0)
        {
          w = mpz_scan1 (a, 0);
          mpz_tdiv_q_2exp (a, a, w);
          ed += (mpfr_exp_t) w;
        }
      else
        {
          w = mpz_scan1 (a, 0);
          v = mpz_scan1 (b, 0);
          if (v < w)
            w = v;
          mpz_tdiv_q_2exp (a, a, w);
          mpz_tdiv_q_2exp (b, b, w);
          ed += (mpfr_exp_t) w;
        }
      if (   (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
          || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
        goto end;
    }
  /* now a+I*b = (c+I*d)^my */

  while (ey-- > 0)
    {
      unsigned long ta, tb;

      /* square a + I*b */
      mpz_mul (u, a, b);
      mpz_mul (a, a, a);
      mpz_submul (a, b, b);
      mpz_mul_2exp (b, u, 1);
      sb = sa * sb;
      /* a cannot be zero here when the initial x is real or imaginary,
         which is the only case where the initial signs of zero do matter,
         thus sa does not matter here when a=0. */
      sa = mpz_cmp_ui (a, 0) >= 0 ? +1 : -1;
      ed *= 2;

      /* divide by largest 2^n possible, to avoid many loops for e.g.,
         (2+2*I)^16777216 */
      ta = mpz_scan1 (a, 0);
      tb = mpz_scan1 (b, 0);
      ta = (ta <= tb) ? ta : tb;
      mpz_tdiv_q_2exp (a, a, ta);
      mpz_tdiv_q_2exp (b, b, ta);
      ed += (mpfr_exp_t) ta;

      if (   (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
          || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
        goto end;
    }

  ret = mpfr_set_z (mpc_realref(z), a, MPC_RND_RE(rnd));
  if (mpz_cmp_ui (a, 0) == 0 && sa < 0)
    mpfr_neg (mpc_realref(z), mpc_realref(z), MPC_RND_RE(rnd));
  inex_im = mpfr_set_z (mpc_imagref(z), b, MPC_RND_IM(rnd));
  if (mpz_cmp_ui (b, 0) == 0 && sb < 0)
    mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPC_RND_IM(rnd));
  ret = MPC_INEX(ret, inex_im);
  mpfr_mul_2si (mpc_realref(z), mpc_realref(z), ed, MPC_RND_RE(rnd));
  mpfr_mul_2si (mpc_imagref(z), mpc_imagref(z), ed, MPC_RND_IM(rnd));

 end:
  mpz_clear (my);
  mpz_clear (a);
  mpz_clear (b);
  mpz_clear (c);
  mpz_clear (d);
  mpz_clear (u);

  if (z_is_y)
    mpfr_clear (copy_of_y);

  return ret;
}

/* Return 1 if y*2^k is an odd integer, 0 otherwise.
   Adapted from MPFR, file pow.c.

   Examples: with k=0, check if y is an odd integer,
             with k=1, check if y is half-an-integer,
             with k=-1, check if y/2 is an odd integer.
*/
#define MPFR_LIMB_HIGHBIT ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1))
static int
is_odd (mpfr_srcptr y, mpfr_exp_t k)
{
  mpfr_exp_t expo;
  mpfr_prec_t prec;
  mp_size_t yn;
  mp_limb_t *yp;

  expo = mpfr_get_exp (y) + k;
  if (expo <= 0)
    return 0;  /* |y| < 1 and not 0 */

  prec = mpfr_get_prec (y);
  if ((mpfr_prec_t) expo > prec)
    return 0;  /* y is a multiple of 2^(expo-prec), thus not odd */

  /* 0 < expo <= prec:
     y = 1xxxxxxxxxt.zzzzzzzzzzzzzzzzzz[000]
          expo bits   (prec-expo) bits

     We have to check that:
     (a) the bit 't' is set
     (b) all the 'z' bits are zero
  */

  prec = ((prec - 1) / BITS_PER_MP_LIMB + 1) * BITS_PER_MP_LIMB - expo;
  /* number of z+0 bits */

  yn = prec / BITS_PER_MP_LIMB;
  /* yn is the index of limb containing the 't' bit */

  yp = y->_mpfr_d;
  /* if expo is a multiple of BITS_PER_MP_LIMB, t is bit 0 */
  if (expo % BITS_PER_MP_LIMB == 0 ? (yp[yn] & 1) == 0
      : yp[yn] << ((expo % BITS_PER_MP_LIMB) - 1) != MPFR_LIMB_HIGHBIT)
    return 0;
  while (--yn >= 0)
    if (yp[yn] != 0)
      return 0;
  return 1;
}

/* Put in z the value of x^y, rounded according to 'rnd'.
   Return the inexact flag in [0, 10]. */
int
mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
{
  int ret = -2, x_real, x_imag, y_real, z_real = 0, z_imag = 0,
     ramified = 0;
  mpc_t t, u;
  mpfr_prec_t p, pr, pi, maxprec;
  int saved_underflow, saved_overflow;
  int inex_re, inex_im;
  mpfr_exp_t saved_emin, saved_emax;

  /* save the underflow or overflow flags from MPFR */
  saved_underflow = mpfr_underflow_p ();
  saved_overflow = mpfr_overflow_p ();

  x_real = mpfr_zero_p (mpc_imagref(x));
  y_real = mpfr_zero_p (mpc_imagref(y));

  if (y_real && mpfr_zero_p (mpc_realref(y))) /* case y zero */
    {
      if (x_real && mpfr_zero_p (mpc_realref(x)))
        {
          /* we define 0^0 to be (1, +0) since the real part is
             coherent with MPFR where 0^0 gives 1, and the sign of the
             imaginary part cannot be determined                       */
          mpc_set_ui_ui (z, 1, 0, MPC_RNDNN);
          return 0;
        }
      else /* x^0 = 1 +/- i*0 even for x=NaN see algorithms.tex for the
              sign of zero */
        {
          mpfr_t n;
          int inex, cx1;
          int sign_zi;
          /* cx1 < 0 if |x| < 1
             cx1 = 0 if |x| = 1
             cx1 > 0 if |x| > 1
          */
          mpfr_init (n);
          inex = mpc_norm (n, x, MPFR_RNDN);
          cx1 = mpfr_cmp_ui (n, 1);
          if (cx1 == 0 && inex != 0)
            cx1 = -inex;

          sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
            || (cx1 == 0
                && mpfr_signbit (mpc_imagref (x)) != mpfr_signbit (mpc_realref (y)))
            || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));

          /* warning: mpc_set_ui_ui does not set Im(z) to -0 if Im(rnd)=RNDD */
          ret = mpc_set_ui_ui (z, 1, 0, rnd);

          if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi)
            mpc_conj (z, z, MPC_RNDNN);

          mpfr_clear (n);
          return ret;
        }
    }

  if (!mpc_fin_p (x) || !mpc_fin_p (y))
    {
      /* special values: exp(y*log(x)) */
      mpc_init2 (u, 2);
      mpc_log (u, x, MPC_RNDNN);
      mpc_mul (u, u, y, MPC_RNDNN);
      ret = mpc_exp (z, u, rnd);
      mpc_clear (u);
      goto end;
    }

  if (x_real) /* case x real */
    {
      int s1, s2;

      if (mpfr_zero_p (mpc_realref(x))) /* x is zero */
        {
          /* special values: exp(y*log(x)) */
          mpc_init2 (u, 2);
          mpc_log (u, x, MPC_RNDNN);
          mpc_mul (u, u, y, MPC_RNDNN);
          ret = mpc_exp (z, u, rnd);
          mpc_clear (u);
          goto end;
        }

      /* Special case 1^y = 1 */
      if (mpfr_cmp_ui (mpc_realref(x), 1) == 0)
        {
          s1 = mpfr_signbit (mpc_realref (y));
          s2 = mpfr_signbit (mpc_imagref (x));

          ret = mpc_set_ui (z, +1, rnd);
          /* the sign of the zero imaginary part is known in some cases (see
             algorithm.tex). In such cases we have
             (x +s*0i)^(y+/-0i) = x^y + s*sign(y)*0i
             where s = +/-1.  We extend here this rule to fix the sign of the
             zero part.

             Note that the sign must also be set explicitly when rnd=RNDD
             because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0.
          */
          if (MPC_RND_IM (rnd) == MPFR_RNDD || s1 != s2)
            mpc_conj (z, z, MPC_RNDNN);
          goto end;
        }

      /* x^y is real when:
         (a) x is real and y is integer
         (b) x is real non-negative and y is real
         (and we are in the case where x is real) */
      s1 = mpfr_signbit (mpc_realref (x));
      if (y_real && (mpfr_integer_p (mpc_realref(y)) || s1 == 0))
        {
          int s3, s4, cmp, even = 0;
          s2 = mpfr_signbit (mpc_imagref (x));
          s3 = mpfr_signbit (mpc_realref (y));
          s4 = mpfr_signbit (mpc_imagref (y));
          cmp = mpfr_cmpabs_ui (mpc_realref (x), 1);
          if (s1 != 0) { /* x is negative, so Re(y) is integer */
            mpz_t yint;
            mpz_init (yint);
            mpfr_get_z (yint, mpc_realref(y), MPFR_RNDN);
            even = mpz_even_p (yint);
            mpz_clear (yint);
          }
          ret = mpfr_pow (mpc_realref(z), mpc_realref(x), mpc_realref(y), MPC_RND_RE(rnd));
          inex_im = mpfr_set_ui (mpc_imagref(z), 0, MPC_RND_IM(rnd));
          ret = MPC_INEX(ret, inex_im);

          /* the choice of the sign of the zero imaginary part is detailed
             in algorithm.tex. */
          fix_sign_real (z, s1, s2, s3, s4, cmp, even);
          goto end;
        }

      /* (-1)^(n+I*t) is real for n integer and t real */
      if (mpfr_cmp_si (mpc_realref(x), -1) == 0 && mpfr_integer_p (mpc_realref(y)))
        z_real = 1;

      /* for x real, x^y is imaginary when:
         (a) x is negative and y is half-an-integer
         (b) x = -1 and Re(y) is half-an-integer
      */
      if ((mpfr_cmp_ui (mpc_realref(x), 0) < 0) && is_odd (mpc_realref(y), 1)
         && (y_real || mpfr_cmp_si (mpc_realref(x), -1) == 0))
        z_imag = 1;
    }
  else /* x non real */
    /* I^(t*I) and (-I)^(t*I) are real for t real,
       I^(n+t*I) and (-I)^(n+t*I) are real for n even and t real, and
       I^(n+t*I) and (-I)^(n+t*I) are imaginary for n odd and t real
       (s*I)^n is real for n even and imaginary for n odd */
    if ((mpc_cmp_si_si (x, 0, 1) == 0 || mpc_cmp_si_si (x, 0, -1) == 0 ||
         (mpfr_cmp_ui (mpc_realref(x), 0) == 0 && y_real)) &&
        mpfr_integer_p (mpc_realref(y)))
      { /* x is I or -I, and Re(y) is an integer */
        if (is_odd (mpc_realref(y), 0))
          z_imag = 1; /* Re(y) odd: z is imaginary */
        else
          z_real = 1; /* Re(y) even: z is real */
      }
    else /* (t+/-t*I)^(2n) is imaginary for n odd and real for n even */
      if (mpfr_cmpabs (mpc_realref(x), mpc_imagref(x)) == 0 && y_real &&
          mpfr_integer_p (mpc_realref(y)) && is_odd (mpc_realref(y), 0) == 0)
        {
          ramified = 1;
          if (is_odd (mpc_realref(y), -1)) /* y/2 is odd */
            z_imag = 1;
          else
            z_real = 1;
        }

  saved_emin = mpfr_get_emin ();
  saved_emax = mpfr_get_emax ();
  mpfr_set_emin (mpfr_get_emin_min ());
  mpfr_set_emax (mpfr_get_emax_max ());

  pr = mpfr_get_prec (mpc_realref(z));
  pi = mpfr_get_prec (mpc_imagref(z));
  p = (pr > pi) ? pr : pi;
  p += 12; /* experimentally, seems to give less than 10% of failures in
              Ziv's strategy; probably wrong now since q is not computed */
  if (p < 64)
    p = 64;
  mpc_init2 (u, p);
  mpc_init2 (t, p);
  pr += MPC_RND_RE(rnd) == MPFR_RNDN;
  pi += MPC_RND_IM(rnd) == MPFR_RNDN;
  maxprec = MPC_MAX_PREC (z);
  x_imag = mpfr_zero_p (mpc_realref(x));
  for (;;)
    {
      int ret_exp;
      mpfr_exp_t dr, di;
      mpfr_prec_t q;

      mpc_log (t, x, MPC_RNDNN);
      mpc_mul (t, t, y, MPC_RNDNN);

      /* Compute q such that |Re (y log x)|, |Im (y log x)| < 2^q.
         We recompute it at each loop since we might get different
         bounds if the precision is not enough. */
      q = mpfr_get_exp (mpc_realref(t)) > 0 ? mpfr_get_exp (mpc_realref(t)) : 0;
      if (mpfr_get_exp (mpc_imagref(t)) > (mpfr_exp_t) q)
        q = mpfr_get_exp (mpc_imagref(t));

      /* the signs of the real/imaginary parts of exp(t) are determined by the
         quadrant of exp(i*imag(t)), which depends on imag(t) mod (2pi).
         We ensure that p >= q + 64 to get enough precision, but this might
         be not enough in corner cases (FIXME). */
      if (p < q + 64)
        {
          p = q + 64;
          goto try_again;
        }

      mpfr_clear_overflow ();
      mpfr_clear_underflow ();
      ret_exp = mpc_exp (u, t, MPC_RNDNN);
      if (mpfr_underflow_p () || mpfr_overflow_p ()) {
         /* under- and overflow flags are set by mpc_exp */
         mpc_set (z, u, MPC_RNDNN);
         inex_re = MPC_INEX_RE(ret_exp);
         inex_im = MPC_INEX_IM(ret_exp);
         if (mpfr_inf_p (mpc_realref (z)))
           inex_re = mpc_fix_inf (mpc_realref (z), MPC_RND_RE(rnd));
         if (mpfr_inf_p (mpc_imagref (z)))
           {
             if (z_real)
               inex_im = mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM(rnd));
             else
               inex_im = mpc_fix_inf (mpc_imagref (z), MPC_RND_IM(rnd));
           }
         ret = MPC_INEX(inex_re,inex_im);
         goto exact;
      }

      /* Since the error bound is global, we have to take into account the
         exponent difference between the real and imaginary parts. We assume
         either the real or the imaginary part of u is not zero.
      */
      dr = mpfr_zero_p (mpc_realref(u)) ? mpfr_get_exp (mpc_imagref(u))
        : mpfr_get_exp (mpc_realref(u));
      di = mpfr_zero_p (mpc_imagref(u)) ? dr : mpfr_get_exp (mpc_imagref(u));
      if (dr > di)
        {
          di = dr - di;
          dr = 0;
        }
      else
        {
          dr = di - dr;
          di = 0;
        }
      /* the term -3 takes into account the factor 4 in the complex error
         (see algorithms.tex) plus one due to the exponent difference: if
         z = a + I*b, where the relative error on z is at most 2^(-p), and
         EXP(a) = EXP(b) + k, the relative error on b is at most 2^(k-p) */
      if ((z_imag || (p > q + 3 + dr && mpfr_can_round (mpc_realref(u), p - q - 3 - dr, MPFR_RNDN, MPFR_RNDZ, pr))) &&
          (z_real || (p > q + 3 + di && mpfr_can_round (mpc_imagref(u), p - q - 3 - di, MPFR_RNDN, MPFR_RNDZ, pi))))
        break;

      /* if Re(u) is not known to be zero, assume it is a normal number, i.e.,
         neither zero, Inf or NaN, otherwise we might enter an infinite loop */
      MPC_ASSERT (z_imag || mpfr_number_p (mpc_realref(u)));
      /* idem for Im(u) */
      MPC_ASSERT (z_real || mpfr_number_p (mpc_imagref(u)));

      if (ret == -2) /* we did not yet call mpc_pow_exact, or it aborted
                        because intermediate computations had > maxprec bits */
        {
          /* check exact cases (see algorithms.tex) */
          if (y_real)
            {
              maxprec *= 2;
              ret = mpc_pow_exact (z, x, mpc_realref(y), rnd, maxprec);
              if (ret != -1 && ret != -2)
                goto exact;
            }
          p += dr + di + 64;
        }
      else
        p += p / 2;
    try_again:
      mpc_set_prec (t, p);
      mpc_set_prec (u, p);
    }

  if (z_real)
    {
      /* When the result is real (see algorithm.tex for details) and
         x=x1 * (1 \pm i), y a positive integer divisible by 4, then
         Im(x^y) = 0i with a sign that cannot be determined (and is thus
         chosen as _1). Otherwise,
         Im(x^y) =
         + sign(imag(y))*0i,               if |x| > 1
         + sign(imag(x))*sign(real(y))*0i, if |x| = 1
         - sign(imag(y))*0i,               if |x| < 1
      */
      if (ramified)
         ret = MPC_INEX (
            mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd)),
            mpfr_set_ui (mpc_imagref (z), 0, MPFR_RNDN));
      else {
         mpfr_t n;
         int inex, cx1;
         int sign_zi, sign_rex, sign_imx;
         /* cx1 < 0 if |x| < 1
            cx1 = 0 if |x| = 1
            cx1 > 0 if |x| > 1
         */

         sign_rex = mpfr_signbit (mpc_realref (x));
         sign_imx = mpfr_signbit (mpc_imagref (x));
         mpfr_init (n);
         inex = mpc_norm (n, x, MPFR_RNDN);
         cx1 = mpfr_cmp_ui (n, 1);
         if (cx1 == 0 && inex != 0)
           cx1 = -inex;

         sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
           || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y)))
           || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));

         /* copy RE(y) to n since if z==y we will destroy Re(y) below */
         mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y)));
         mpfr_set (n, mpc_realref (y), MPFR_RNDN);
         ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd));
         if (y_real && (x_real || x_imag))
           {
             /* FIXME: with y_real we assume Im(y) is really 0, which is the case
                for example when y comes from pow_fr, but in case Im(y) is +0 or
                -0, we might get different results */
             mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
             fix_sign (z, sign_rex, sign_imx, n);
             ret = MPC_INEX(ret, 0); /* imaginary part is exact */
           }
         else
           {
             inex_im = mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
             ret = MPC_INEX (ret, inex_im);
             /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */
             if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi)
               mpc_conj (z, z, MPC_RNDNN);
           }

         mpfr_clear (n);
      }
    }
  else if (z_imag)
    {
      if (ramified)
         ret = MPC_INEX (
            mpfr_set_ui (mpc_realref (z), 0, MPFR_RNDN),
            mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd)));
      else
      {
         ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd));
         /* if z is imaginary and y real, then x cannot be real */
         if (y_real && x_imag)
           {
             int sign_rex = mpfr_signbit (mpc_realref (x));

             /* If z overlaps with y we set Re(z) before checking Re(y) below,
                but in that case y=0, which was dealt with above. */
             mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd));
             /* Note: fix_sign only does something when y is an integer,
                then necessarily y = 1 or 3 (mod 4), and in that case the
                sign of Im(x) is irrelevant. */
             fix_sign (z, sign_rex, 0, mpc_realref (y));
             ret = MPC_INEX(0, ret);
           }
         else
           {
             inex_re = mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd));
             ret = MPC_INEX(inex_re, ret);
           }
      }
    }
  else
    ret = mpc_set (z, u, rnd);
 exact:
  mpc_clear (t);
  mpc_clear (u);

  /* restore underflow and overflow flags from MPFR */
  if (saved_underflow)
    mpfr_set_underflow ();
  if (saved_overflow)
    mpfr_set_overflow ();

  /* restore the exponent range, and check the range of results */
  mpfr_set_emin (saved_emin);
  mpfr_set_emax (saved_emax);
  inex_re = mpfr_check_range (mpc_realref (z), MPC_INEX_RE(ret),
                              MPC_RND_RE (rnd));
  inex_im = mpfr_check_range (mpc_imagref (z), MPC_INEX_IM(ret),
                              MPC_RND_IM (rnd));
  ret = MPC_INEX(inex_re, inex_im);

 end:
  return ret;
}