1 /* mpc_exp -- exponential of a complex number.
3 Copyright (C) 2002, 2009, 2010, 2011 INRIA
5 This file is part of GNU MPC.
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.
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
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/ .
24 mpc_exp (mpc_ptr rop
, mpc_srcptr op
, mpc_rnd_t rnd
)
30 int saved_underflow
, saved_overflow
;
33 if (mpfr_nan_p (mpc_realref (op
)) || mpfr_nan_p (mpc_imagref (op
)))
35 exp(nan +i*y) = nan -i*0 if y = -0,
38 exp(x+i*nan) = +/-0 +/-i*0 if x=-inf,
39 +/-inf +i*nan if x=+inf,
40 nan +i*nan otherwise */
42 if (mpfr_zero_p (mpc_imagref (op
)))
43 return mpc_set (rop
, op
, MPC_RNDNN
);
45 if (mpfr_inf_p (mpc_realref (op
)))
47 if (mpfr_signbit (mpc_realref (op
)))
48 return mpc_set_ui_ui (rop
, 0, 0, MPC_RNDNN
);
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 */
56 mpfr_set_nan (mpc_realref (rop
));
57 mpfr_set_nan (mpc_imagref (rop
));
58 return MPC_INEX(0, 0); /* NaN is exact */
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 */
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
);
72 if (mpfr_zero_p (mpc_realref (op
)))
73 /* special case when the input is imaginary */
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
);
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 */
90 if (mpfr_signbit (mpc_realref (op
)))
91 mpfr_set_ui (n
, 0, GMP_RNDN
);
95 if (mpfr_inf_p (mpc_imagref (op
)))
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
);
102 mpfr_set_nan (mpc_imagref (rop
));
103 inex_im
= 0; /* NaN is exact */
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
);
121 return MPC_INEX(inex_re
, inex_im
);
124 if (mpfr_inf_p (mpc_imagref (op
)))
125 /* real part is finite non-zero number, imaginary part is an infinity */
127 mpfr_set_nan (mpc_realref (rop
));
128 mpfr_set_nan (mpc_imagref (rop
));
129 return MPC_INEX(0, 0); /* NaN is exact */
133 /* from now on, both parts of op are regular numbers */
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
146 /* save the underflow or overflow flags from MPFR */
147 saved_underflow
= mpfr_underflow_p ();
148 saved_overflow
= mpfr_overflow_p ();
152 prec
+= mpc_ceil_log2 (prec
) + 5;
154 mpfr_set_prec (x
, prec
);
155 mpfr_set_prec (y
, prec
);
156 mpfr_set_prec (z
, prec
);
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 */
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
));
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
);
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);
195 /* restore underflow and overflow flags from MPFR */
197 mpfr_set_underflow ();
199 mpfr_set_overflow ();
201 return MPC_INEX(inex_re
, inex_im
);