/* mpc_log -- Take the logarithm of a complex number.

Copyright (C) 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 <stdio.h> /* for MPC_ASSERT */

#include "mpc-impl.h"

int

mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){

   int ok, underflow = 0;

   mpfr_srcptr x, y;

   mpfr_t v, w;

   mpfr_prec_t prec;

   int loops;

   int re_cmp, im_cmp;

   int inex_re, inex_im;

   int err;

   mpfr_exp_t expw;

   int sgnw;

   /* special values: NaN and infinities */

   if (!mpc_fin_p (op)) {

      if (mpfr_nan_p (mpc_realref (op))) {

         if (mpfr_inf_p (mpc_imagref (op)))

            mpfr_set_inf (mpc_realref (rop), +1);

         else

            mpfr_set_nan (mpc_realref (rop));

         mpfr_set_nan (mpc_imagref (rop));

         inex_im = 0; /* Inf/NaN is exact */

      }

      else if (mpfr_nan_p (mpc_imagref (op))) {

         if (mpfr_inf_p (mpc_realref (op)))

            mpfr_set_inf (mpc_realref (rop), +1);

         else

            mpfr_set_nan (mpc_realref (rop));

         mpfr_set_nan (mpc_imagref (rop));

         inex_im = 0; /* Inf/NaN is exact */

      }

      else /* We have an infinity in at least one part. */ {

         inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op),

                               MPC_RND_IM (rnd));

         mpfr_set_inf (mpc_realref (rop), +1);

      }

      return MPC_INEX(0, inex_im);

   }

   /* special cases: real and purely imaginary numbers */

   re_cmp = mpfr_cmp_ui (mpc_realref (op), 0);

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

   if (im_cmp == 0) {

      if (re_cmp == 0) {

         inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op),

                               MPC_RND_IM (rnd));

         mpfr_set_inf (mpc_realref (rop), -1);

         inex_re = 0; /* -Inf is exact */

      }

      else if (re_cmp > 0) {

         inex_re = mpfr_log (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));

         inex_im = mpfr_set (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd));

      }

      else {

         /* op = x + 0*y; let w = -x = |x| */

         int negative_zero;

         mpfr_rnd_t rnd_im;

         negative_zero = mpfr_signbit (mpc_imagref (op));

         if (negative_zero)

            rnd_im = INV_RND (MPC_RND_IM (rnd));

         else

            rnd_im = MPC_RND_IM (rnd);

         w [0] = *mpc_realref (op);

         MPFR_CHANGE_SIGN (w);

         inex_re = mpfr_log (mpc_realref (rop), w, MPC_RND_RE (rnd));

         inex_im = mpfr_const_pi (mpc_imagref (rop), rnd_im);

         if (negative_zero) {

            mpc_conj (rop, rop, MPC_RNDNN);

            inex_im = -inex_im;

         }

      }

      return MPC_INEX(inex_re, inex_im);

   }

   else if (re_cmp == 0) {

      if (im_cmp > 0) {

         inex_re = mpfr_log (mpc_realref (rop), mpc_imagref (op), MPC_RND_RE (rnd));

         inex_im = mpfr_const_pi (mpc_imagref (rop), MPC_RND_IM (rnd));

         /* division by 2 does not change the ternary flag */

         mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, MPFR_RNDN);

      }

      else {

         w [0] = *mpc_imagref (op);

         MPFR_CHANGE_SIGN (w);

         inex_re = mpfr_log (mpc_realref (rop), w, MPC_RND_RE (rnd));

         inex_im = mpfr_const_pi (mpc_imagref (rop), INV_RND (MPC_RND_IM (rnd)));

         /* division by 2 does not change the ternary flag */

         mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, MPFR_RNDN);

         mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), MPFR_RNDN);

         inex_im = -inex_im; /* negate the ternary flag */

      }

      return MPC_INEX(inex_re, inex_im);

   }

   prec = MPC_PREC_RE(rop);

   mpfr_init2 (w, 2);

   /* let op = x + iy; log = 1/2 log (x^2 + y^2) + i atan2 (y, x)   */

   /* loop for the real part: 1/2 log (x^2 + y^2), fast, but unsafe */

   /* implementation                                                */

   ok = 0;

   for (loops = 1; !ok && loops <= 2; loops++) {

      prec += mpc_ceil_log2 (prec) + 4;

      mpfr_set_prec (w, prec);

      mpc_abs (w, op, MPFR_RNDN);

         /* error 0.5 ulp */

      if (mpfr_inf_p (w))

         /* intermediate overflow; the logarithm may be representable.

            Intermediate underflow is impossible.                      */

         break;

      mpfr_log (w, w, MPFR_RNDN);

         /* generic error of log: (2^(- exp(w)) + 0.5) ulp */

      if (mpfr_zero_p (w))

         /* impossible to round, switch to second algorithm */

         break;

      err = MPC_MAX (-mpfr_get_exp (w), 0) + 1;

         /* number of lost digits */

      ok = mpfr_can_round (w, prec - err, MPFR_RNDN, MPFR_RNDZ,

         mpfr_get_prec (mpc_realref (rop)) + (MPC_RND_RE (rnd) == MPFR_RNDN));

   }

   if (!ok) {

      prec = MPC_PREC_RE(rop);

      mpfr_init2 (v, 2);

      /* compute 1/2 log (x^2 + y^2) = log |x| + 1/2 * log (1 + (y/x)^2)

            if |x| >= |y|; otherwise, exchange x and y                   */

      if (mpfr_cmpabs (mpc_realref (op), mpc_imagref (op)) >= 0) {

         x = mpc_realref (op);

         y = mpc_imagref (op);

      }

      else {

         x = mpc_imagref (op);

         y = mpc_realref (op);

      }

      do {

         prec += mpc_ceil_log2 (prec) + 4;

         mpfr_set_prec (v, prec);

         mpfr_set_prec (w, prec);

         mpfr_div (v, y, x, MPFR_RNDD); /* error 1 ulp */

         mpfr_sqr (v, v, MPFR_RNDD);

            /* generic error of multiplication:

               1 + 2*1*(2+1*2^(1-prec)) <= 5.0625 since prec >= 6 */

         mpfr_log1p (v, v, MPFR_RNDD);

            /* error 1 + 4*5.0625 = 21.25 , see algorithms.tex */

         mpfr_div_2ui (v, v, 1, MPFR_RNDD);

            /* If the result is 0, then there has been an underflow somewhere. */

         mpfr_abs (w, x, MPFR_RNDN); /* exact */

         mpfr_log (w, w, MPFR_RNDN); /* error 0.5 ulp */

         expw = mpfr_get_exp (w);

         sgnw = mpfr_signbit (w);

         mpfr_add (w, w, v, MPFR_RNDN);

         if (!sgnw) /* v is positive, so no cancellation;

                       error 22.25 ulp; error counts lost bits */

            err = 5;

         else

            err =   MPC_MAX (5 + mpfr_get_exp (v),

                  /* 21.25 ulp (v) rewritten in ulp (result, now in w) */

                           -1 + expw             - mpfr_get_exp (w)

                  /* 0.5 ulp (previous w), rewritten in ulp (result) */

                  ) + 2;

         /* handle one special case: |x|=1, and (y/x)^2 underflows;

            then 1/2*log(x^2+y^2) \approx 1/2*y^2 also underflows.  */

         if (   (mpfr_cmp_si (x, -1) == 0 || mpfr_cmp_ui (x, 1) == 0)

             && mpfr_zero_p (w))

            underflow = 1;

      } while (!underflow &&

               !mpfr_can_round (w, prec - err, MPFR_RNDN, MPFR_RNDZ,

               mpfr_get_prec (mpc_realref (rop)) + (MPC_RND_RE (rnd) == MPFR_RNDN)));

      mpfr_clear (v);

   }

   /* imaginary part */

   inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op),

                         MPC_RND_IM (rnd));

   /* set the real part; cannot be done before if rop==op */

   if (underflow)

      /* create underflow in result */

      inex_re = mpfr_set_ui_2exp (mpc_realref (rop), 1,

                                  mpfr_get_emin_min () - 2, MPC_RND_RE (rnd));

   else

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

   mpfr_clear (w);

   return MPC_INEX(inex_re, inex_im);

}