1 /* mpc_acos -- arccosine of a complex number. 2 3 Copyright (C) 2009, 2010, 2011, 2012 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 <stdio.h> /* for MPC_ASSERT */ 22 #include "mpc-impl.h" 23 24 int 25 mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) 26 { 27 int inex_re, inex_im, inex; 28 mpfr_prec_t p_re, p_im, p; 29 mpc_t z1; 30 mpfr_t pi_over_2; 31 mpfr_exp_t e1, e2; 32 mpfr_rnd_t rnd_im; 33 mpc_rnd_t rnd1; 34 35 inex_re = 0; 36 inex_im = 0; 37 38 /* special values */ 39 if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) 40 { 41 if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op))) 42 { 43 mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? +1 : -1); 44 mpfr_set_nan (mpc_realref (rop)); 45 } 46 else if (mpfr_zero_p (mpc_realref (op))) 47 { 48 inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd)); 49 mpfr_set_nan (mpc_imagref (rop)); 50 } 51 else 52 { 53 mpfr_set_nan (mpc_realref (rop)); 54 mpfr_set_nan (mpc_imagref (rop)); 55 } 56 57 return MPC_INEX (inex_re, 0); 58 } 59 60 if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op))) 61 { 62 if (mpfr_inf_p (mpc_realref (op))) 63 { 64 if (mpfr_inf_p (mpc_imagref (op))) 65 { 66 if (mpfr_sgn (mpc_realref (op)) > 0) 67 { 68 inex_re = 69 set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd)); 70 mpfr_div_2ui (mpc_realref (rop), mpc_realref (rop), 1, GMP_RNDN); 71 } 72 else 73 { 74 75 /* the real part of the result is 3*pi/4 76 a = o(pi) error(a) < 1 ulp(a) 77 b = o(3*a) error(b) < 2 ulp(b) 78 c = b/4 exact 79 thus 1 bit is lost */ 80 mpfr_t x; 81 mpfr_prec_t prec; 82 int ok; 83 mpfr_init (x); 84 prec = mpfr_get_prec (mpc_realref (rop)); 85 p = prec; 86 87 do 88 { 89 p += mpc_ceil_log2 (p); 90 mpfr_set_prec (x, p); 91 mpfr_const_pi (x, GMP_RNDD); 92 mpfr_mul_ui (x, x, 3, GMP_RNDD); 93 ok = 94 mpfr_can_round (x, p - 1, GMP_RNDD, MPC_RND_RE (rnd), 95 prec+(MPC_RND_RE (rnd) == GMP_RNDN)); 96 97 } while (ok == 0); 98 inex_re = 99 mpfr_div_2ui (mpc_realref (rop), x, 2, MPC_RND_RE (rnd)); 100 mpfr_clear (x); 101 } 102 } 103 else 104 { 105 if (mpfr_sgn (mpc_realref (op)) > 0) 106 mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); 107 else 108 inex_re = mpfr_const_pi (mpc_realref (rop), MPC_RND_RE (rnd)); 109 } 110 } 111 else 112 inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd)); 113 114 mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? +1 : -1); 115 116 return MPC_INEX (inex_re, 0); 117 } 118 119 /* pure real argument */ 120 if (mpfr_zero_p (mpc_imagref (op))) 121 { 122 int s_im; 123 s_im = mpfr_signbit (mpc_imagref (op)); 124 125 if (mpfr_cmp_ui (mpc_realref (op), 1) > 0) 126 { 127 if (s_im) 128 inex_im = mpfr_acosh (mpc_imagref (rop), mpc_realref (op), 129 MPC_RND_IM (rnd)); 130 else 131 inex_im = -mpfr_acosh (mpc_imagref (rop), mpc_realref (op), 132 INV_RND (MPC_RND_IM (rnd))); 133 134 mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); 135 } 136 else if (mpfr_cmp_si (mpc_realref (op), -1) < 0) 137 { 138 mpfr_t minus_op_re; 139 minus_op_re[0] = mpc_realref (op)[0]; 140 MPFR_CHANGE_SIGN (minus_op_re); 141 142 if (s_im) 143 inex_im = mpfr_acosh (mpc_imagref (rop), minus_op_re, 144 MPC_RND_IM (rnd)); 145 else 146 inex_im = -mpfr_acosh (mpc_imagref (rop), minus_op_re, 147 INV_RND (MPC_RND_IM (rnd))); 148 inex_re = mpfr_const_pi (mpc_realref (rop), MPC_RND_RE (rnd)); 149 } 150 else 151 { 152 inex_re = mpfr_acos (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); 153 mpfr_set_ui (mpc_imagref (rop), 0, MPC_RND_IM (rnd)); 154 } 155 156 if (!s_im) 157 mpc_conj (rop, rop, MPC_RNDNN); 158 159 return MPC_INEX (inex_re, inex_im); 160 } 161 162 /* pure imaginary argument */ 163 if (mpfr_zero_p (mpc_realref (op))) 164 { 165 inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd)); 166 inex_im = -mpfr_asinh (mpc_imagref (rop), mpc_imagref (op), 167 INV_RND (MPC_RND_IM (rnd))); 168 mpc_conj (rop,rop, MPC_RNDNN); 169 170 return MPC_INEX (inex_re, inex_im); 171 } 172 173 /* regular complex argument: acos(z) = Pi/2 - asin(z) */ 174 p_re = mpfr_get_prec (mpc_realref(rop)); 175 p_im = mpfr_get_prec (mpc_imagref(rop)); 176 p = p_re; 177 mpc_init3 (z1, p, p_im); /* we round directly the imaginary part to p_im, 178 with rounding mode opposite to rnd_im */ 179 rnd_im = MPC_RND_IM(rnd); 180 /* the imaginary part of asin(z) has the same sign as Im(z), thus if 181 Im(z) > 0 and rnd_im = RNDZ, we want to round the Im(asin(z)) to -Inf 182 so that -Im(asin(z)) is rounded to zero */ 183 if (rnd_im == GMP_RNDZ) 184 rnd_im = mpfr_sgn (mpc_imagref(op)) > 0 ? GMP_RNDD : GMP_RNDU; 185 else 186 rnd_im = rnd_im == GMP_RNDU ? GMP_RNDD 187 : rnd_im == GMP_RNDD ? GMP_RNDU 188 : rnd_im; /* both RNDZ and RNDA map to themselves for -asin(z) */ 189 rnd1 = MPC_RND (GMP_RNDN, rnd_im); 190 mpfr_init2 (pi_over_2, p); 191 for (;;) 192 { 193 p += mpc_ceil_log2 (p) + 3; 194 195 mpfr_set_prec (mpc_realref(z1), p); 196 mpfr_set_prec (pi_over_2, p); 197 198 set_pi_over_2 (pi_over_2, +1, GMP_RNDN); 199 e1 = 1; /* Exp(pi_over_2) */ 200 inex = mpc_asin (z1, op, rnd1); /* asin(z) */ 201 MPC_ASSERT (mpfr_sgn (mpc_imagref(z1)) * mpfr_sgn (mpc_imagref(op)) > 0); 202 inex_im = MPC_INEX_IM(inex); /* inex_im is in {-1, 0, 1} */ 203 e2 = mpfr_get_exp (mpc_realref(z1)); 204 mpfr_sub (mpc_realref(z1), pi_over_2, mpc_realref(z1), GMP_RNDN); 205 if (!mpfr_zero_p (mpc_realref(z1))) 206 { 207 /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-p-1) + 208 2^(e2-p-1) */ 209 e1 = e1 >= e2 ? e1 + 1 : e2 + 1; 210 /* the error on x is bounded by 1/2 ulp(x) + 2^(e1-p-1) */ 211 e1 -= mpfr_get_exp (mpc_realref(z1)); 212 /* the error on x is bounded by 1/2 ulp(x) [1 + 2^e1] */ 213 e1 = e1 <= 0 ? 0 : e1; 214 /* the error on x is bounded by 2^e1 * ulp(x) */ 215 mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN); /* exact */ 216 inex_im = -inex_im; 217 if (mpfr_can_round (mpc_realref(z1), p - e1, GMP_RNDN, GMP_RNDZ, 218 p_re + (MPC_RND_RE(rnd) == GMP_RNDN))) 219 break; 220 } 221 } 222 inex = mpc_set (rop, z1, rnd); 223 inex_re = MPC_INEX_RE(inex); 224 mpc_clear (z1); 225 mpfr_clear (pi_over_2); 226 227 return MPC_INEX(inex_re, inex_im); 228 } 229