/* mpc_sqrt -- Take the square root of a complex number.

Copyright (C) 2002, 2008, 2009, 2010, 2011, 2012 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 "mpc-impl.h"

#if MPFR_VERSION_MAJOR < 3

#define mpfr_min_prec(x) \

   ( ((prec + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB) * BITS_PER_MP_LIMB \

     - mpn_scan1 (x->_mpfr_d, 0))

#endif

int

mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd)

{

  int ok_w, ok_t = 0;

  mpfr_t    w, t;

  mpfr_rnd_t  rnd_w, rnd_t;

  mpfr_prec_t prec_w, prec_t;

  /* the rounding mode and the precision required for w and t, which can */

  /* be either the real or the imaginary part of a */

  mpfr_prec_t prec;

  int inex_w, inex_t = 1, inex_re, inex_im, loops = 0;

  const int re_cmp = mpfr_cmp_ui (mpc_realref (b), 0),

            im_cmp = mpfr_cmp_ui (mpc_imagref (b), 0);

     /* comparison of the real/imaginary part of b with 0 */

  int repr_w, repr_t = 0 /* to avoid gcc warning */ ;

     /* flag indicating whether the computed value is already representable

        at the target precision */

  const int im_sgn = mpfr_signbit (mpc_imagref (b)) == 0 ? 0 : -1;

     /* we need to know the sign of Im(b) when it is +/-0 */

  const mpfr_rnd_t r = im_sgn ? MPFR_RNDD : MPFR_RNDU;

     /* rounding mode used when computing t */

  /* special values */

  if (!mpc_fin_p (b)) {

   /* sqrt(x +i*Inf) = +Inf +I*Inf, even if x = NaN */

   /* sqrt(x -i*Inf) = +Inf -I*Inf, even if x = NaN */

   if (mpfr_inf_p (mpc_imagref (b)))

      {

         mpfr_set_inf (mpc_realref (a), +1);

         mpfr_set_inf (mpc_imagref (a), im_sgn);

         return MPC_INEX (0, 0);

      }

   if (mpfr_inf_p (mpc_realref (b)))

      {

         if (mpfr_signbit (mpc_realref (b)))

         {

            if (mpfr_number_p (mpc_imagref (b)))

               {

               /* sqrt(-Inf +i*y) = +0 +i*Inf, when y positive */

               /* sqrt(-Inf +i*y) = +0 -i*Inf, when y positive */

               mpfr_set_ui (mpc_realref (a), 0, MPFR_RNDN);

               mpfr_set_inf (mpc_imagref (a), im_sgn);

               return MPC_INEX (0, 0);

               }

            else

               {

               /* sqrt(-Inf +i*NaN) = NaN +/-i*Inf */

               mpfr_set_nan (mpc_realref (a));

               mpfr_set_inf (mpc_imagref (a), im_sgn);

               return MPC_INEX (0, 0);

               }

         }

         else

         {

            if (mpfr_number_p (mpc_imagref (b)))

               {

               /* sqrt(+Inf +i*y) = +Inf +i*0, when y positive */

               /* sqrt(+Inf +i*y) = +Inf -i*0, when y positive */

               mpfr_set_inf (mpc_realref (a), +1);

               mpfr_set_ui (mpc_imagref (a), 0, MPFR_RNDN);

               if (im_sgn)

                  mpc_conj (a, a, MPC_RNDNN);

               return MPC_INEX (0, 0);

               }

            else

               {

               /* sqrt(+Inf -i*Inf) = +Inf -i*Inf */

               /* sqrt(+Inf +i*Inf) = +Inf +i*Inf */

               /* sqrt(+Inf +i*NaN) = +Inf +i*NaN */

               return mpc_set (a, b, rnd);

               }

         }

      }

   /* sqrt(x +i*NaN) = NaN +i*NaN, if x is not infinite */

   /* sqrt(NaN +i*y) = NaN +i*NaN, if y is not infinite */

   if (mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b)))

      {

         mpfr_set_nan (mpc_realref (a));

         mpfr_set_nan (mpc_imagref (a));

         return MPC_INEX (0, 0);

      }

  }

  /* purely real */

  if (im_cmp == 0)

    {

      if (re_cmp == 0)

        {

          mpc_set_ui_ui (a, 0, 0, MPC_RNDNN);

          if (im_sgn)

            mpc_conj (a, a, MPC_RNDNN);

          return MPC_INEX (0, 0);

        }

      else if (re_cmp > 0)

        {

          inex_w = mpfr_sqrt (mpc_realref (a), mpc_realref (b), MPC_RND_RE (rnd));

          mpfr_set_ui (mpc_imagref (a), 0, MPFR_RNDN);

          if (im_sgn)

            mpc_conj (a, a, MPC_RNDNN);

          return MPC_INEX (inex_w, 0);

        }

      else

        {

          mpfr_init2 (w, MPC_PREC_RE (b));

          mpfr_neg (w, mpc_realref (b), MPFR_RNDN);

          if (im_sgn)

            {

              inex_w = -mpfr_sqrt (mpc_imagref (a), w, INV_RND (MPC_RND_IM (rnd)));

              mpfr_neg (mpc_imagref (a), mpc_imagref (a), MPFR_RNDN);

            }

          else

            inex_w = mpfr_sqrt (mpc_imagref (a), w, MPC_RND_IM (rnd));

          mpfr_set_ui (mpc_realref (a), 0, MPFR_RNDN);

          mpfr_clear (w);

          return MPC_INEX (0, inex_w);

        }

    }

  /* purely imaginary */

  if (re_cmp == 0)

    {

      mpfr_t y;

      y[0] = mpc_imagref (b)[0];

      /* If y/2 underflows, so does sqrt(y/2) */

      mpfr_div_2ui (y, y, 1, MPFR_RNDN);

      if (im_cmp > 0)

        {

          inex_w = mpfr_sqrt (mpc_realref (a), y, MPC_RND_RE (rnd));

          inex_t = mpfr_sqrt (mpc_imagref (a), y, MPC_RND_IM (rnd));

        }

      else

        {

          mpfr_neg (y, y, MPFR_RNDN);

          inex_w = mpfr_sqrt (mpc_realref (a), y, MPC_RND_RE (rnd));

          inex_t = -mpfr_sqrt (mpc_imagref (a), y, INV_RND (MPC_RND_IM (rnd)));

          mpfr_neg (mpc_imagref (a), mpc_imagref (a), MPFR_RNDN);

        }

      return MPC_INEX (inex_w, inex_t);

    }

  prec = MPC_MAX_PREC(a);

  mpfr_init (w);

  mpfr_init (t);

   if (re_cmp > 0) {

      rnd_w = MPC_RND_RE (rnd);

      prec_w = MPC_PREC_RE (a);

      rnd_t = MPC_RND_IM(rnd);

      if (rnd_t == MPFR_RNDZ)

         /* force MPFR_RNDD or MPFR_RNDUP, using sign(t) = sign(y) */

         rnd_t = (im_cmp > 0 ? MPFR_RNDD : MPFR_RNDU);

      prec_t = MPC_PREC_IM (a);

   }

   else {

      prec_w = MPC_PREC_IM (a);

      prec_t = MPC_PREC_RE (a);

      if (im_cmp > 0) {

         rnd_w = MPC_RND_IM(rnd);

         rnd_t = MPC_RND_RE(rnd);

         if (rnd_t == MPFR_RNDZ)

            rnd_t = MPFR_RNDD;

      }

      else {

         rnd_w = INV_RND(MPC_RND_IM (rnd));

         rnd_t = INV_RND(MPC_RND_RE (rnd));

         if (rnd_t == MPFR_RNDZ)

            rnd_t = MPFR_RNDU;

      }

   }

  do

    {

      loops ++;

      prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2;

      mpfr_set_prec (w, prec);

      mpfr_set_prec (t, prec);

      /* let b = x + iy */

      /* w = sqrt ((|x| + sqrt (x^2 + y^2)) / 2), rounded down */

      /* total error bounded by 3 ulps */

      inex_w = mpc_abs (w, b, MPFR_RNDD);

      if (re_cmp < 0)

        inex_w |= mpfr_sub (w, w, mpc_realref (b), MPFR_RNDD);

      else

        inex_w |= mpfr_add (w, w, mpc_realref (b), MPFR_RNDD);

      inex_w |= mpfr_div_2ui (w, w, 1, MPFR_RNDD);

      inex_w |= mpfr_sqrt (w, w, MPFR_RNDD);

      repr_w = mpfr_min_prec (w) <= prec_w;

      if (!repr_w)

         /* use the usual trick for obtaining the ternary value */

         ok_w = mpfr_can_round (w, prec - 2, MPFR_RNDD, MPFR_RNDU,

                                prec_w + (rnd_w == MPFR_RNDN));

      else {

            /* w is representable in the target precision and thus cannot be

               rounded up */

         if (rnd_w == MPFR_RNDN)

            /* If w can be rounded to nearest, then actually no rounding

               occurs, and the ternary value is known from inex_w. */

            ok_w = mpfr_can_round (w, prec - 2, MPFR_RNDD, MPFR_RNDN, prec_w);

         else

            /* If w can be rounded down, then any direct rounding and the

               ternary flag can be determined from inex_w. */

            ok_w = mpfr_can_round (w, prec - 2, MPFR_RNDD, MPFR_RNDD, prec_w);

      }

      if (!inex_w || ok_w) {

         /* t = y / 2w, rounded away */

         /* total error bounded by 7 ulps */

         inex_t = mpfr_div (t, mpc_imagref (b), w, r);

         if (!inex_t && inex_w)

            /* The division was exact, but w was not. */

            inex_t = im_sgn ? -1 : 1;

         inex_t |= mpfr_div_2ui (t, t, 1, r);

         repr_t = mpfr_min_prec (t) <= prec_t;

         if (!repr_t)

             /* As for w; since t was rounded away, we check whether rounding to 0

                is possible. */

            ok_t = mpfr_can_round (t, prec - 3, r, MPFR_RNDZ,

                                   prec_t + (rnd_t == MPFR_RNDN));

         else {

            if (rnd_t == MPFR_RNDN)

               ok_t = mpfr_can_round (t, prec - 3, r, MPFR_RNDN, prec_t);

            else

               ok_t = mpfr_can_round (t, prec - 3, r, r, prec_t);

         }

      }

    }

    while ((inex_w && !ok_w) || (inex_t && !ok_t));

   if (re_cmp > 0) {

         inex_re = mpfr_set (mpc_realref (a), w, MPC_RND_RE(rnd));

         inex_im = mpfr_set (mpc_imagref (a), t, MPC_RND_IM(rnd));

   }

   else if (im_cmp > 0) {

      inex_re = mpfr_set (mpc_realref(a), t, MPC_RND_RE(rnd));

      inex_im = mpfr_set (mpc_imagref(a), w, MPC_RND_IM(rnd));

   }

   else {

      inex_re = mpfr_neg (mpc_realref (a), t, MPC_RND_RE(rnd));

      inex_im = mpfr_neg (mpc_imagref (a), w, MPC_RND_IM(rnd));

   }

   if (repr_w && inex_w) {

      if (rnd_w == MPFR_RNDN) {

         /* w has not been rounded with mpfr_set/mpfr_neg, determine ternary

            value from inex_w instead */

         if (re_cmp > 0)

            inex_re = inex_w;

         else if (im_cmp > 0)

            inex_im = inex_w;

         else

            inex_im = -inex_w;

      }

      else {

         /* determine ternary value, but also potentially add 1 ulp; can only

            be done now when we are in the target precision */

         if (re_cmp > 0) {

            if (rnd_w == MPFR_RNDU) {

               MPFR_ADD_ONE_ULP (mpc_realref (a));

               inex_re = +1;

            }

            else

               inex_re = -1;

         }

         else if (im_cmp > 0) {

            if (rnd_w == MPFR_RNDU) {

               MPFR_ADD_ONE_ULP (mpc_imagref (a));

               inex_im = +1;

            }

            else

               inex_im = -1;

         }

         else {

            if (rnd_w == MPFR_RNDU) {

               MPFR_ADD_ONE_ULP (mpc_imagref (a));

               inex_im = -1;

            }

            else

               inex_im = +1;

         }

      }

   }

   if (repr_t && inex_t) {

      if (rnd_t == MPFR_RNDN) {

         if (re_cmp > 0)

            inex_im = inex_t;

         else if (im_cmp > 0)

            inex_re = inex_t;

         else

            inex_re = -inex_t;

      }

      else {

         if (re_cmp > 0) {

            if (rnd_t == r)

               inex_im = inex_t;

            else {

               inex_im = -inex_t;

               /* im_cmp > 0 implies that Im(b) > 0, thus im_sgn = 0

                  and r = MPFR_RNDU.

                  im_cmp < 0 implies that Im(b) < 0, thus im_sgn = -1

                  and r = MPFR_RNDD. */

               MPFR_SUB_ONE_ULP (mpc_imagref (a));

            }

         }

         else if (im_cmp > 0) {

            if (rnd_t == r)

               inex_re = inex_t;

            else {

               inex_re = -inex_t;

               /* im_cmp > 0 implies r = MPFR_RNDU (see above) */

               MPFR_SUB_ONE_ULP (mpc_realref (a));

            }

         }

         else { /* im_cmp < 0 */

            if (rnd_t == r)

               inex_re = -inex_t;

            else {

               inex_re = inex_t;

               /* im_cmp < 0 implies r = MPFR_RNDD (see above) */

               MPFR_SUB_ONE_ULP (mpc_realref (a));

            }

         }

      }

   }

  mpfr_clear (w);

  mpfr_clear (t);

  return MPC_INEX (inex_re, inex_im);

}