1/* mpc_exp -- exponential of a complex number. 2 3Copyright (C) 2002, 2009, 2010, 2011, 2012, 2020 INRIA 4 5This file is part of GNU MPC. 6 7GNU MPC is free software; you can redistribute it and/or modify it under 8the terms of the GNU Lesser General Public License as published by the 9Free Software Foundation; either version 3 of the License, or (at your 10option) any later version. 11 12GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY 13WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for 15more details. 16 17You should have received a copy of the GNU Lesser General Public License 18along with this program. If not, see http://www.gnu.org/licenses/ . 19*/ 20 21#include "mpc-impl.h" 22 23int 24mpc_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