src/exp.c

   1 /* mpc_exp -- exponential of a complex number.
   2
   3 Copyright (C) 2002, 2009 Andreas Enge, Paul Zimmermann, Philippe Th\'eveny
   4
   5 This file is part of the MPC Library.
   6
   7 The MPC Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 2.1 of the License, or (at your
  10 option) any later version.
  11
  12 The MPC Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the MPC Library; see the file COPYING.LIB.  If not, write to
  19 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  20 MA 02111-1307, USA. */
  21
  22 #include "mpc-impl.h"
  23
  24 int
  25 mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
  26 {
  27   mpfr_t x, y, z;
  28   mp_prec_t prec;
  29   int ok = 0;
  30   int inex_re, inex_im;
  31
  32   /* let op = a + i*b, then exp(op) = exp(a)*[cos(b) + i*sin(b)]
  33                                     = exp(a)*cos(b) + i*exp(a)*sin(b).
  34
  35      We use the following algorithm (same for the imaginary part):
  36
  37      (1) x = o(exp(a)) rounded towards +infinity:
  38      (2) y = o(cos(b)) rounded to nearest
  39      (3) r = o(x*y)
  40      then the error on r for the real part is at most 4 ulps:
  41      |r - exp(a)*cos(b)| <= ulp(r) + |x*y - exp(a)*cos(b)|
  42                          <= ulp(r) + |x*y - exp(a)*y| + exp(a) * |y - cos(b)|
  43                          <= ulp(r) + |y| ulp(x) + 1/2 * x * ulp(y)
  44                          <= ulp(r) + 2 * ulp(x*y) + ulp(x*y) [Rule 4]
  45                          <= 4 * ulp(r) [Rule 8]
  46   */
  47
  48   /* special values */
  49   if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op)))
  50     /* NaNs
  51        exp(nan +i*y) = nan -i*0   if y = -0,
  52                        nan +i*0   if y = +0,
  53                        nan +i*nan otherwise
  54        exp(x+i*nan) =   +/-0 +/-i*0 if x=-inf,
  55                       +/-inf +i*nan if x=+inf,
  56                          nan +i*nan otherwise */
  57     {
  58       if (mpfr_zero_p (MPC_IM (op)))
  59         return mpc_set (rop, op, MPC_RNDNN);
  60
  61       if (mpfr_inf_p (MPC_RE (op)))
  62         {
  63           if (mpfr_signbit (MPC_RE (op)))
  64             return mpc_set_ui_ui (rop, 0, 0, MPC_RNDNN);
  65           else
  66             {
  67               mpfr_set_inf (MPC_RE (rop), +1);
  68               mpfr_set_nan (MPC_IM (rop));
  69               return MPC_INEX(0, 0); /* Inf/NaN are exact */
  70             }
  71         }
  72       mpfr_set_nan (MPC_RE (rop));
  73       mpfr_set_nan (MPC_IM (rop));
  74       return MPC_INEX(0, 0); /* NaN is exact */
  75     }
  76
  77
  78   if (mpfr_zero_p (MPC_IM(op)))
  79     /* special case when the input is real
  80        exp(x-i*0) = exp(x) -i*0, even if x is NaN
  81        exp(x+i*0) = exp(x) +i*0, even if x is NaN */
  82     {
  83       inex_re = mpfr_exp (MPC_RE(rop), MPC_RE(op), MPC_RND_RE(rnd));
  84       inex_im = mpfr_set (MPC_IM(rop), MPC_IM(op), MPC_RND_IM(rnd));
  85       return MPC_INEX(inex_re, inex_im);
  86     }
  87
  88   if (mpfr_zero_p (MPC_RE (op)))
  89     /* special case when the input is imaginary  */
  90     {
  91       inex_re = mpfr_cos (MPC_RE (rop), MPC_IM (op), MPC_RND_RE(rnd));
  92       inex_im = mpfr_sin (MPC_IM (rop), MPC_IM (op), MPC_RND_IM(rnd));
  93       return MPC_INEX(inex_re, inex_im);
  94     }
  95
  96
  97   if (mpfr_inf_p (MPC_RE (op)))
  98     /* real part is an infinity,
  99        exp(-inf +i*y) = 0*(cos y +i*sin y)
 100        exp(+inf +i*y) = +/-inf +i*nan         if y = +/-inf
 101                         +inf*(cos y +i*sin y) if 0 < |y| < inf */
 102     {
 103       mpfr_t n;
 104
 105       mpfr_init2 (n, 2);
 106       if (mpfr_signbit (MPC_RE (op)))
 107         mpfr_set_ui (n, 0, GMP_RNDN);
 108       else
 109         mpfr_set_inf (n, +1);
 110
 111       if (mpfr_inf_p (MPC_IM (op)))
 112         {
 113           inex_re = mpfr_set (MPC_RE (rop), n, GMP_RNDN);
 114           if (mpfr_signbit (MPC_RE (op)))
 115             inex_im = mpfr_set (MPC_IM (rop), n, GMP_RNDN);
 116           else
 117             {
 118               mpfr_set_nan (MPC_IM (rop));
 119               inex_im = 0; /* NaN is exact */
 120             }
 121         }
 122       else
 123         {
 124           mpfr_t c, s;
 125           mpfr_init2 (c, 2);
 126           mpfr_init2 (s, 2);
 127
 128           mpfr_sin_cos (s, c, MPC_IM (op), GMP_RNDN);
 129           inex_re = mpfr_copysign (MPC_RE (rop), n, c, GMP_RNDN);
 130           inex_im = mpfr_copysign (MPC_IM (rop), n, s, GMP_RNDN);
 131
 132           mpfr_clear (s);
 133           mpfr_clear (c);
 134         }
 135
 136       mpfr_clear (n);
 137       return MPC_INEX(inex_re, inex_im);
 138     }
 139
 140   if (mpfr_inf_p (MPC_IM (op)))
 141     /* real part is finite non-zero number, imaginary part is an infinity */
 142     {
 143       mpfr_set_nan (MPC_RE (rop));
 144       mpfr_set_nan (MPC_IM (rop));
 145       return MPC_INEX(0, 0); /* NaN is exact */
 146     }
 147
 148
 149   /* from now on, both parts of op are regular numbers */
 150
 151   prec = MPC_MAX_PREC(rop);
 152
 153   mpfr_init2 (x, 2);
 154   mpfr_init2 (y, 2);
 155   mpfr_init2 (z, 2);
 156
 157   do
 158     {
 159       prec += mpc_ceil_log2 (prec) + 5;
 160
 161       mpfr_set_prec (x, prec);
 162       mpfr_set_prec (y, prec);
 163       mpfr_set_prec (z, prec);
 164
 165       /* FIXME: x may overflow so x.y does overflow too, while Re(exp(op))
 166          could be represented in the precision of rop. */
 167       mpfr_exp (x, MPC_RE(op), GMP_RNDN);
 168       mpfr_sin_cos (z, y, MPC_IM(op), GMP_RNDN);
 169       mpfr_mul (y, y, x, GMP_RNDN);
 170
 171       ok = mpfr_inf_p (y) || mpfr_zero_p (x)
 172         || mpfr_can_round (y, prec - 2, GMP_RNDN, GMP_RNDZ,
 173                        MPFR_PREC(MPC_RE(rop)) + (MPC_RND_RE(rnd) == GMP_RNDN));
 174       if (ok) /* compute imaginary part */
 175         {
 176           mpfr_mul (z, z, x, GMP_RNDN);
 177           ok = mpfr_inf_p (z) || mpfr_zero_p (x)
 178             || mpfr_can_round (z, prec - 2, GMP_RNDN, GMP_RNDZ,
 179                        MPFR_PREC(MPC_IM(rop)) + (MPC_RND_IM(rnd) == GMP_RNDN));
 180         }
 181     }
 182   while (ok == 0);
 183
 184   inex_re = mpfr_set (MPC_RE(rop), y, MPC_RND_RE(rnd));
 185   inex_im = mpfr_set (MPC_IM(rop), z, MPC_RND_IM(rnd));
 186
 187   mpfr_clear (x);
 188   mpfr_clear (y);
 189   mpfr_clear (z);
 190
 191   return MPC_INEX(inex_re, inex_im);
 192 }