1 /* mpc_exp -- exponential of a complex number.
3 Copyright (C) 2002, 2009 Andreas Enge, Paul Zimmermann, Philippe Th\'eveny
5 This file is part of the MPC Library.
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.
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.
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. */
25 mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
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).
35 We use the following algorithm (same for the imaginary part):
37 (1) x = o(exp(a)) rounded towards +infinity:
38 (2) y = o(cos(b)) rounded to nearest
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]
49 if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op)))
51 exp(nan +i*y) = nan -i*0 if y = -0,
54 exp(x+i*nan) = +/-0 +/-i*0 if x=-inf,
55 +/-inf +i*nan if x=+inf,
56 nan +i*nan otherwise */
58 if (mpfr_zero_p (MPC_IM (op)))
59 return mpc_set (rop, op, MPC_RNDNN);
61 if (mpfr_inf_p (MPC_RE (op)))
63 if (mpfr_signbit (MPC_RE (op)))
64 return mpc_set_ui_ui (rop, 0, 0, MPC_RNDNN);
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 */
72 mpfr_set_nan (MPC_RE (rop));
73 mpfr_set_nan (MPC_IM (rop));
74 return MPC_INEX(0, 0); /* NaN is exact */
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 */
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);
88 if (mpfr_zero_p (MPC_RE (op)))
89 /* special case when the input is imaginary */
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);
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 */
106 if (mpfr_signbit (MPC_RE (op)))
107 mpfr_set_ui (n, 0, GMP_RNDN);
109 mpfr_set_inf (n, +1);
111 if (mpfr_inf_p (MPC_IM (op)))
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);
118 mpfr_set_nan (MPC_IM (rop));
119 inex_im = 0; /* NaN is exact */
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);
137 return MPC_INEX(inex_re, inex_im);
140 if (mpfr_inf_p (MPC_IM (op)))
141 /* real part is finite non-zero number, imaginary part is an infinity */
143 mpfr_set_nan (MPC_RE (rop));
144 mpfr_set_nan (MPC_IM (rop));
145 return MPC_INEX(0, 0); /* NaN is exact */
149 /* from now on, both parts of op are regular numbers */
151 prec = MPC_MAX_PREC(rop);
159 prec += mpc_ceil_log2 (prec) + 5;
161 mpfr_set_prec (x, prec);
162 mpfr_set_prec (y, prec);
163 mpfr_set_prec (z, prec);
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);
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 */
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));
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));
191 return MPC_INEX(inex_re, inex_im);