contrib/mpc/src/exp.c

   1 /* mpc_exp -- exponential of a complex number.
   2
   3 Copyright (C) 2002, 2009, 2010, 2011 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 "mpc-impl.h"
  22
  23 int
  24 mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
  25 {
  26   mpfr_t x, y, z;
  27   mpfr_prec_t prec;
  28   int ok = 0;
  29   int inex_re, inex_im;
  30   int saved_underflow, saved_overflow;
  31
  32   /* special values */
  33   if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
  34     /* NaNs
  35        exp(nan +i*y) = nan -i*0   if y = -0,
  36                        nan +i*0   if y = +0,
  37                        nan +i*nan otherwise
  38        exp(x+i*nan) =   +/-0 +/-i*0 if x=-inf,
  39                       +/-inf +i*nan if x=+inf,
  40                          nan +i*nan otherwise */
  41     {
  42       if (mpfr_zero_p (mpc_imagref (op)))
  43         return mpc_set (rop, op, MPC_RNDNN);
  44
  45       if (mpfr_inf_p (mpc_realref (op)))
  46         {
  47           if (mpfr_signbit (mpc_realref (op)))
  48             return mpc_set_ui_ui (rop, 0, 0, MPC_RNDNN);
  49           else
  50             {
  51               mpfr_set_inf (mpc_realref (rop), +1);
  52               mpfr_set_nan (mpc_imagref (rop));
  53               return MPC_INEX(0, 0); /* Inf/NaN are exact */
  54             }
  55         }
  56       mpfr_set_nan (mpc_realref (rop));
  57       mpfr_set_nan (mpc_imagref (rop));
  58       return MPC_INEX(0, 0); /* NaN is exact */
  59     }
  60
  61
  62   if (mpfr_zero_p (mpc_imagref(op)))
  63     /* special case when the input is real
  64        exp(x-i*0) = exp(x) -i*0, even if x is NaN
  65        exp(x+i*0) = exp(x) +i*0, even if x is NaN */
  66     {
  67       inex_re = mpfr_exp (mpc_realref(rop), mpc_realref(op), MPC_RND_RE(rnd));
  68       inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref(op), MPC_RND_IM(rnd));
  69       return MPC_INEX(inex_re, inex_im);
  70     }
  71
  72   if (mpfr_zero_p (mpc_realref (op)))
  73     /* special case when the input is imaginary  */
  74     {
  75       inex_re = mpfr_cos (mpc_realref (rop), mpc_imagref (op), MPC_RND_RE(rnd));
  76       inex_im = mpfr_sin (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM(rnd));
  77       return MPC_INEX(inex_re, inex_im);
  78     }
  79
  80
  81   if (mpfr_inf_p (mpc_realref (op)))
  82     /* real part is an infinity,
  83        exp(-inf +i*y) = 0*(cos y +i*sin y)
  84        exp(+inf +i*y) = +/-inf +i*nan         if y = +/-inf
  85                         +inf*(cos y +i*sin y) if 0 < |y| < inf */
  86     {
  87       mpfr_t n;
  88
  89       mpfr_init2 (n, 2);
  90       if (mpfr_signbit (mpc_realref (op)))
  91         mpfr_set_ui (n, 0, GMP_RNDN);
  92       else
  93         mpfr_set_inf (n, +1);
  94
  95       if (mpfr_inf_p (mpc_imagref (op)))
  96         {
  97           inex_re = mpfr_set (mpc_realref (rop), n, GMP_RNDN);
  98           if (mpfr_signbit (mpc_realref (op)))
  99             inex_im = mpfr_set (mpc_imagref (rop), n, GMP_RNDN);
 100           else
 101             {
 102               mpfr_set_nan (mpc_imagref (rop));
 103               inex_im = 0; /* NaN is exact */
 104             }
 105         }
 106       else
 107         {
 108           mpfr_t c, s;
 109           mpfr_init2 (c, 2);
 110           mpfr_init2 (s, 2);
 111
 112           mpfr_sin_cos (s, c, mpc_imagref (op), GMP_RNDN);
 113           inex_re = mpfr_copysign (mpc_realref (rop), n, c, GMP_RNDN);
 114           inex_im = mpfr_copysign (mpc_imagref (rop), n, s, GMP_RNDN);
 115
 116           mpfr_clear (s);
 117           mpfr_clear (c);
 118         }
 119
 120       mpfr_clear (n);
 121       return MPC_INEX(inex_re, inex_im);
 122     }
 123
 124   if (mpfr_inf_p (mpc_imagref (op)))
 125     /* real part is finite non-zero number, imaginary part is an infinity */
 126     {
 127       mpfr_set_nan (mpc_realref (rop));
 128       mpfr_set_nan (mpc_imagref (rop));
 129       return MPC_INEX(0, 0); /* NaN is exact */
 130     }
 131
 132
 133   /* from now on, both parts of op are regular numbers */
 134
 135   prec = MPC_MAX_PREC(rop)
 136          + MPC_MAX (MPC_MAX (-mpfr_get_exp (mpc_realref (op)), 0),
 137                    -mpfr_get_exp (mpc_imagref (op)));
 138     /* When op is close to 0, then exp is close to 1+Re(op), while
 139        cos is close to 1-Im(op); to decide on the ternary value of exp*cos,
 140        we need a high enough precision so that none of exp or cos is
 141        computed as 1. */
 142   mpfr_init2 (x, 2);
 143   mpfr_init2 (y, 2);
 144   mpfr_init2 (z, 2);
 145
 146   /* save the underflow or overflow flags from MPFR */
 147   saved_underflow = mpfr_underflow_p ();
 148   saved_overflow = mpfr_overflow_p ();
 149
 150   do
 151     {
 152       prec += mpc_ceil_log2 (prec) + 5;
 153
 154       mpfr_set_prec (x, prec);
 155       mpfr_set_prec (y, prec);
 156       mpfr_set_prec (z, prec);
 157
 158       /* FIXME: x may overflow so x.y does overflow too, while Re(exp(op))
 159          could be represented in the precision of rop. */
 160       mpfr_clear_overflow ();
 161       mpfr_clear_underflow ();
 162       mpfr_exp (x, mpc_realref(op), GMP_RNDN); /* error <= 0.5ulp */
 163       mpfr_sin_cos (z, y, mpc_imagref(op), GMP_RNDN); /* errors <= 0.5ulp */
 164       mpfr_mul (y, y, x, GMP_RNDN); /* error <= 2ulp */
 165       ok = mpfr_overflow_p () || mpfr_zero_p (x)
 166         || mpfr_can_round (y, prec - 2, GMP_RNDN, GMP_RNDZ,
 167                        MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == GMP_RNDN));
 168       if (ok) /* compute imaginary part */
 169         {
 170           mpfr_mul (z, z, x, GMP_RNDN);
 171           ok = mpfr_overflow_p () || mpfr_zero_p (x)
 172             || mpfr_can_round (z, prec - 2, GMP_RNDN, GMP_RNDZ,
 173                        MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == GMP_RNDN));
 174         }
 175     }
 176   while (ok == 0);
 177
 178   inex_re = mpfr_set (mpc_realref(rop), y, MPC_RND_RE(rnd));
 179   inex_im = mpfr_set (mpc_imagref(rop), z, MPC_RND_IM(rnd));
 180   if (mpfr_overflow_p ()) {
 181     /* overflow in real exponential, inex is sign of infinite result */
 182     inex_re = mpfr_sgn (y);
 183     inex_im = mpfr_sgn (z);
 184   }
 185   else if (mpfr_underflow_p ()) {
 186     /* underflow in real exponential, inex is opposite of sign of 0 result */
 187     inex_re = (mpfr_signbit (y) ? +1 : -1);
 188     inex_im = (mpfr_signbit (z) ? +1 : -1);
 189   }
 190
 191   mpfr_clear (x);
 192   mpfr_clear (y);
 193   mpfr_clear (z);
 194
 195   /* restore underflow and overflow flags from MPFR */
 196   if (saved_underflow)
 197     mpfr_set_underflow ();
 198   if (saved_overflow)
 199     mpfr_set_overflow ();
 200
 201   return MPC_INEX(inex_re, inex_im);
 202 }