1 /* mpc_exp -- exponential of a complex number.
2
3 Copyright (C) 2002, 2009, 2010, 2011, 2012, 2020 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
mpc_exp(mpc_ptr rop,mpc_srcptr op,mpc_rnd_t rnd)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 mpfr_exp_t saved_emin, saved_emax;
32
33 /* special values */
34 if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
35 /* NaNs
36 exp(nan +i*y) = nan -i*0 if y = -0,
37 nan +i*0 if y = +0,
38 nan +i*nan otherwise
39 exp(x+i*nan) = +/-0 +/-i*0 if x=-inf,
40 +/-inf +i*nan if x=+inf,
41 nan +i*nan otherwise */
42 {
43 if (mpfr_zero_p (mpc_imagref (op)))
44 return mpc_set (rop, op, MPC_RNDNN);
45
46 if (mpfr_inf_p (mpc_realref (op)))
47 {
48 if (mpfr_signbit (mpc_realref (op)))
49 return mpc_set_ui_ui (rop, 0, 0, MPC_RNDNN);
50 else
51 {
52 mpfr_set_inf (mpc_realref (rop), +1);
53 mpfr_set_nan (mpc_imagref (rop));
54 return MPC_INEX(0, 0); /* Inf/NaN are exact */
55 }
56 }
57 mpfr_set_nan (mpc_realref (rop));
58 mpfr_set_nan (mpc_imagref (rop));
59 return MPC_INEX(0, 0); /* NaN is exact */
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 if (mpfr_inf_p (mpc_realref (op)))
81 /* real part is an infinity,
82 exp(-inf +i*y) = 0*(cos y +i*sin y)
83 exp(+inf +i*y) = +/-inf +i*nan if y = +/-inf
84 +inf*(cos y +i*sin y) if 0 < |y| < inf */
85 {
86 mpfr_t n;
87
88 mpfr_init2 (n, 2);
89 if (mpfr_signbit (mpc_realref (op)))
90 mpfr_set_ui (n, 0, MPFR_RNDN);
91 else
92 mpfr_set_inf (n, +1);
93
94 if (mpfr_inf_p (mpc_imagref (op)))
95 {
96 int real_sign = mpfr_signbit (mpc_realref (op));
97 inex_re = mpfr_set (mpc_realref (rop), n, MPFR_RNDN);
98 if (real_sign)
99 inex_im = mpfr_set (mpc_imagref (rop), n, MPFR_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), MPFR_RNDN);
113 inex_re = mpfr_copysign (mpc_realref (rop), n, c, MPFR_RNDN);
114 inex_im = mpfr_copysign (mpc_imagref (rop), n, s, MPFR_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 saved_emin = mpfr_get_emin ();
133 saved_emax = mpfr_get_emax ();
134 mpfr_set_emin (mpfr_get_emin_min ());
135 mpfr_set_emax (mpfr_get_emax_max ());
136
137 /* from now on, both parts of op are regular numbers */
138
139 prec = MPC_MAX_PREC(rop)
140 + MPC_MAX (MPC_MAX (-mpfr_get_exp (mpc_realref (op)), 0),
141 -mpfr_get_exp (mpc_imagref (op)));
142 /* When op is close to 0, then exp is close to 1+Re(op), while
143 cos is close to 1-Im(op); to decide on the ternary value of exp*cos,
144 we need a high enough precision so that none of exp or cos is
145 computed as 1. */
146 mpfr_init2 (x, 2);
147 mpfr_init2 (y, 2);
148 mpfr_init2 (z, 2);
149
150 /* save the underflow or overflow flags from MPFR */
151 saved_underflow = mpfr_underflow_p ();
152 saved_overflow = mpfr_overflow_p ();
153
154 do
155 {
156 prec += prec / 2 + mpc_ceil_log2 (prec) + 5;
157
158 mpfr_set_prec (x, prec);
159 mpfr_set_prec (y, prec);
160 mpfr_set_prec (z, prec);
161
162 /* FIXME: x may overflow so x.y does overflow too, while Re(exp(op))
163 could be represented in the precision of rop. */
164 mpfr_clear_overflow ();
165 mpfr_clear_underflow ();
166 mpfr_exp (x, mpc_realref(op), MPFR_RNDN); /* error <= 0.5ulp */
167 mpfr_sin_cos (z, y, mpc_imagref(op), MPFR_RNDN); /* errors <= 0.5ulp */
168 mpfr_mul (y, y, x, MPFR_RNDN); /* error <= 2ulp */
169 ok = mpfr_overflow_p () || mpfr_zero_p (x)
170 || mpfr_can_round (y, prec - 2, MPFR_RNDN, MPFR_RNDZ,
171 MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == MPFR_RNDN));
172 if (ok) /* compute imaginary part */
173 {
174 mpfr_mul (z, z, x, MPFR_RNDN);
175 ok = mpfr_overflow_p () || mpfr_zero_p (x)
176 || mpfr_can_round (z, prec - 2, MPFR_RNDN, MPFR_RNDZ,
177 MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == MPFR_RNDN));
178 }
179 }
180 while (ok == 0);
181
182 inex_re = mpfr_set (mpc_realref(rop), y, MPC_RND_RE(rnd));
183 inex_im = mpfr_set (mpc_imagref(rop), z, MPC_RND_IM(rnd));
184 if (mpfr_overflow_p ())
185 {
186 inex_re = mpc_fix_inf (mpc_realref(rop), MPC_RND_RE(rnd));
187 inex_im = mpc_fix_inf (mpc_imagref(rop), MPC_RND_IM(rnd));
188 }
189 else if (mpfr_underflow_p ())
190 {
191 inex_re = mpc_fix_zero (mpc_realref(rop), MPC_RND_RE(rnd));
192 inex_im = mpc_fix_zero (mpc_imagref(rop), MPC_RND_IM(rnd));
193 }
194
195 mpfr_clear (x);
196 mpfr_clear (y);
197 mpfr_clear (z);
198
199 /* restore underflow and overflow flags from MPFR */
200 if (saved_underflow)
201 mpfr_set_underflow ();
202 if (saved_overflow)
203 mpfr_set_overflow ();
204
205 /* restore the exponent range, and check the range of results */
206 mpfr_set_emin (saved_emin);
207 mpfr_set_emax (saved_emax);
208 inex_re = mpfr_check_range (mpc_realref (rop), inex_re, MPC_RND_RE (rnd));
209 inex_im = mpfr_check_range (mpc_imagref (rop), inex_im, MPC_RND_IM (rnd));
210
211 return MPC_INEX(inex_re, inex_im);
212 }
213