1/* mpc_acos -- arccosine of a complex number. 2 3Copyright (C) INRIA, 2009, 2010 4 5This file is part of the MPC Library. 6 7The MPC Library is free software; you can redistribute it and/or modify 8it under the terms of the GNU Lesser General Public License as published by 9the Free Software Foundation; either version 2.1 of the License, or (at your 10option) any later version. 11 12The MPC Library is distributed in the hope that it will be useful, but 13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 14or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 15License for more details. 16 17You should have received a copy of the GNU Lesser General Public License 18along with the MPC Library; see the file COPYING.LIB. If not, write to 19the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, 20MA 02111-1307, USA. */ 21 22#include <stdio.h> /* for MPC_ASSERT */ 23#include "mpc-impl.h" 24 25int 26mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) 27{ 28 int inex_re, inex_im, inex; 29 mpfr_prec_t p_re, p_im, p; 30 mpc_t z1; 31 mpfr_t pi_over_2; 32 mpfr_exp_t e1, e2; 33 mpfr_rnd_t rnd_im; 34 mpc_rnd_t rnd1; 35 36 inex_re = 0; 37 inex_im = 0; 38 39 /* special values */ 40 if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op))) 41 { 42 if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op))) 43 { 44 mpfr_set_inf (MPC_IM (rop), mpfr_signbit (MPC_IM (op)) ? +1 : -1); 45 mpfr_set_nan (MPC_RE (rop)); 46 } 47 else if (mpfr_zero_p (MPC_RE (op))) 48 { 49 inex_re = set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd)); 50 mpfr_set_nan (MPC_IM (rop)); 51 } 52 else 53 { 54 mpfr_set_nan (MPC_RE (rop)); 55 mpfr_set_nan (MPC_IM (rop)); 56 } 57 58 return MPC_INEX (inex_re, 0); 59 } 60 61 if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op))) 62 { 63 if (mpfr_inf_p (MPC_RE (op))) 64 { 65 if (mpfr_inf_p (MPC_IM (op))) 66 { 67 if (mpfr_sgn (MPC_RE (op)) > 0) 68 { 69 inex_re = 70 set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd)); 71 mpfr_div_2ui (MPC_RE (rop), MPC_RE (rop), 1, GMP_RNDN); 72 } 73 else 74 { 75 76 /* the real part of the result is 3*pi/4 77 a = o(pi) error(a) < 1 ulp(a) 78 b = o(3*a) error(b) < 2 ulp(b) 79 c = b/4 exact 80 thus 1 bit is lost */ 81 mpfr_t x; 82 mpfr_prec_t prec; 83 int ok; 84 mpfr_init (x); 85 prec = mpfr_get_prec (MPC_RE (rop)); 86 p = prec; 87 88 do 89 { 90 p += mpc_ceil_log2 (p); 91 mpfr_set_prec (x, p); 92 mpfr_const_pi (x, GMP_RNDD); 93 mpfr_mul_ui (x, x, 3, GMP_RNDD); 94 ok = 95 mpfr_can_round (x, p - 1, GMP_RNDD, MPC_RND_RE (rnd), 96 prec+(MPC_RND_RE (rnd) == GMP_RNDN)); 97 98 } while (ok == 0); 99 inex_re = 100 mpfr_div_2ui (MPC_RE (rop), x, 2, MPC_RND_RE (rnd)); 101 mpfr_clear (x); 102 } 103 } 104 else 105 { 106 if (mpfr_sgn (MPC_RE (op)) > 0) 107 mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN); 108 else 109 inex_re = mpfr_const_pi (MPC_RE (rop), MPC_RND_RE (rnd)); 110 } 111 } 112 else 113 inex_re = set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd)); 114 115 mpfr_set_inf (MPC_IM (rop), mpfr_signbit (MPC_IM (op)) ? +1 : -1); 116 117 return MPC_INEX (inex_re, 0); 118 } 119 120 /* pure real argument */ 121 if (mpfr_zero_p (MPC_IM (op))) 122 { 123 int s_im; 124 s_im = mpfr_signbit (MPC_IM (op)); 125 126 if (mpfr_cmp_ui (MPC_RE (op), 1) > 0) 127 { 128 if (s_im) 129 inex_im = mpfr_acosh (MPC_IM (rop), MPC_RE (op), 130 MPC_RND_IM (rnd)); 131 else 132 inex_im = -mpfr_acosh (MPC_IM (rop), MPC_RE (op), 133 INV_RND (MPC_RND_IM (rnd))); 134 135 mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN); 136 } 137 else if (mpfr_cmp_si (MPC_RE (op), -1) < 0) 138 { 139 mpfr_t minus_op_re; 140 minus_op_re[0] = MPC_RE (op)[0]; 141 MPFR_CHANGE_SIGN (minus_op_re); 142 143 if (s_im) 144 inex_im = mpfr_acosh (MPC_IM (rop), minus_op_re, 145 MPC_RND_IM (rnd)); 146 else 147 inex_im = -mpfr_acosh (MPC_IM (rop), minus_op_re, 148 INV_RND (MPC_RND_IM (rnd))); 149 inex_re = mpfr_const_pi (MPC_RE (rop), MPC_RND_RE (rnd)); 150 } 151 else 152 { 153 inex_re = mpfr_acos (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd)); 154 mpfr_set_ui (MPC_IM (rop), 0, MPC_RND_IM (rnd)); 155 } 156 157 if (!s_im) 158 mpc_conj (rop, rop, MPC_RNDNN); 159 160 return MPC_INEX (inex_re, inex_im); 161 } 162 163 /* pure imaginary argument */ 164 if (mpfr_zero_p (MPC_RE (op))) 165 { 166 inex_re = set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd)); 167 inex_im = -mpfr_asinh (MPC_IM (rop), MPC_IM (op), 168 INV_RND (MPC_RND_IM (rnd))); 169 mpc_conj (rop,rop, MPC_RNDNN); 170 171 return MPC_INEX (inex_re, inex_im); 172 } 173 174 /* regular complex argument: acos(z) = Pi/2 - asin(z) */ 175 p_re = mpfr_get_prec (MPC_RE(rop)); 176 p_im = mpfr_get_prec (MPC_IM(rop)); 177 p = p_re; 178 mpc_init3 (z1, p, p_im); /* we round directly the imaginary part to p_im, 179 with rounding mode opposite to rnd_im */ 180 rnd_im = MPC_RND_IM(rnd); 181 /* the imaginary part of asin(z) has the same sign as Im(z), thus if 182 Im(z) > 0 and rnd_im = RNDZ, we want to round the Im(asin(z)) to -Inf 183 so that -Im(asin(z)) is rounded to zero */ 184 if (rnd_im == GMP_RNDZ) 185 rnd_im = mpfr_sgn (MPC_IM(op)) > 0 ? GMP_RNDD : GMP_RNDU; 186 else 187 rnd_im = rnd_im == GMP_RNDU ? GMP_RNDD 188 : rnd_im == GMP_RNDD ? GMP_RNDU 189 : rnd_im; /* both RNDZ and RNDA map to themselves for -asin(z) */ 190 rnd1 = RNDC(GMP_RNDN, rnd_im); 191 mpfr_init2 (pi_over_2, p); 192 for (;;) 193 { 194 p += mpc_ceil_log2 (p) + 3; 195 196 mpfr_set_prec (MPC_RE(z1), p); 197 mpfr_set_prec (pi_over_2, p); 198 199 mpfr_const_pi (pi_over_2, GMP_RNDN); 200 mpfr_div_2exp (pi_over_2, pi_over_2, 1, GMP_RNDN); /* Pi/2 */ 201 e1 = 1; /* Exp(pi_over_2) */ 202 inex = mpc_asin (z1, op, rnd1); /* asin(z) */ 203 MPC_ASSERT (mpfr_sgn (MPC_IM(z1)) * mpfr_sgn (MPC_IM(op)) > 0); 204 inex_im = MPC_INEX_IM(inex); /* inex_im is in {-1, 0, 1} */ 205 e2 = mpfr_get_exp (MPC_RE(z1)); 206 mpfr_sub (MPC_RE(z1), pi_over_2, MPC_RE(z1), GMP_RNDN); 207 if (!mpfr_zero_p (MPC_RE(z1))) 208 { 209 /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-p-1) + 210 2^(e2-p-1) */ 211 e1 = e1 >= e2 ? e1 + 1 : e2 + 1; 212 /* the error on x is bounded by 1/2 ulp(x) + 2^(e1-p-1) */ 213 e1 -= mpfr_get_exp (MPC_RE(z1)); 214 /* the error on x is bounded by 1/2 ulp(x) [1 + 2^e1] */ 215 e1 = e1 <= 0 ? 0 : e1; 216 /* the error on x is bounded by 2^e1 * ulp(x) */ 217 mpfr_neg (MPC_IM(z1), MPC_IM(z1), GMP_RNDN); /* exact */ 218 inex_im = -inex_im; 219 if (mpfr_can_round (MPC_RE(z1), p - e1, GMP_RNDN, GMP_RNDZ, 220 p_re + (MPC_RND_RE(rnd) == GMP_RNDN))) 221 break; 222 } 223 } 224 inex = mpc_set (rop, z1, rnd); 225 inex_re = MPC_INEX_RE(inex); 226 mpc_clear (z1); 227 mpfr_clear (pi_over_2); 228 229 return MPC_INEX(inex_re, inex_im); 230} 231