1 /* mpc_exp -- exponential of a complex number. 2 3 Copyright (C) 2002, 2009, 2010, 2011 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 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 32 /* special values */ 33 if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) 34 /* NaNs 35 exp(nan +i*y) = nan -i*0 if y = -0, 36 nan +i*0 if y = +0, 37 nan +i*nan otherwise 38 exp(x+i*nan) = +/-0 +/-i*0 if x=-inf, 39 +/-inf +i*nan if x=+inf, 40 nan +i*nan otherwise */ 41 { 42 if (mpfr_zero_p (mpc_imagref (op))) 43 return mpc_set (rop, op, MPC_RNDNN); 44 45 if (mpfr_inf_p (mpc_realref (op))) 46 { 47 if (mpfr_signbit (mpc_realref (op))) 48 return mpc_set_ui_ui (rop, 0, 0, MPC_RNDNN); 49 else 50 { 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 */ 54 } 55 } 56 mpfr_set_nan (mpc_realref (rop)); 57 mpfr_set_nan (mpc_imagref (rop)); 58 return MPC_INEX(0, 0); /* NaN is exact */ 59 } 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 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 */ 86 { 87 mpfr_t n; 88 89 mpfr_init2 (n, 2); 90 if (mpfr_signbit (mpc_realref (op))) 91 mpfr_set_ui (n, 0, GMP_RNDN); 92 else 93 mpfr_set_inf (n, +1); 94 95 if (mpfr_inf_p (mpc_imagref (op))) 96 { 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); 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), 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); 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 133 /* from now on, both parts of op are regular numbers */ 134 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 141 computed as 1. */ 142 mpfr_init2 (x, 2); 143 mpfr_init2 (y, 2); 144 mpfr_init2 (z, 2); 145 146 /* save the underflow or overflow flags from MPFR */ 147 saved_underflow = mpfr_underflow_p (); 148 saved_overflow = mpfr_overflow_p (); 149 150 do 151 { 152 prec += mpc_ceil_log2 (prec) + 5; 153 154 mpfr_set_prec (x, prec); 155 mpfr_set_prec (y, prec); 156 mpfr_set_prec (z, prec); 157 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 */ 169 { 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)); 174 } 175 } 176 while (ok == 0); 177 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); 184 } 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); 189 } 190 191 mpfr_clear (x); 192 mpfr_clear (y); 193 mpfr_clear (z); 194 195 /* restore underflow and overflow flags from MPFR */ 196 if (saved_underflow) 197 mpfr_set_underflow (); 198 if (saved_overflow) 199 mpfr_set_overflow (); 200 201 return MPC_INEX(inex_re, inex_im); 202 } 203