1 /* mpc_asin -- arcsine of a complex number. 2 3 Copyright (C) 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_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) 25 { 26 mpfr_prec_t p, p_re, p_im, incr_p = 0; 27 mpfr_rnd_t rnd_re, rnd_im; 28 mpc_t z1; 29 int inex; 30 31 /* special values */ 32 if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) 33 { 34 if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op))) 35 { 36 mpfr_set_nan (mpc_realref (rop)); 37 mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? -1 : +1); 38 } 39 else if (mpfr_zero_p (mpc_realref (op))) 40 { 41 mpfr_set (mpc_realref (rop), mpc_realref (op), GMP_RNDN); 42 mpfr_set_nan (mpc_imagref (rop)); 43 } 44 else 45 { 46 mpfr_set_nan (mpc_realref (rop)); 47 mpfr_set_nan (mpc_imagref (rop)); 48 } 49 50 return 0; 51 } 52 53 if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op))) 54 { 55 int inex_re; 56 if (mpfr_inf_p (mpc_realref (op))) 57 { 58 int inf_im = mpfr_inf_p (mpc_imagref (op)); 59 60 inex_re = set_pi_over_2 (mpc_realref (rop), 61 (mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd)); 62 mpfr_set_inf (mpc_imagref (rop), (mpfr_signbit (mpc_imagref (op)) ? -1 : 1)); 63 64 if (inf_im) 65 mpfr_div_2ui (mpc_realref (rop), mpc_realref (rop), 1, GMP_RNDN); 66 } 67 else 68 { 69 mpfr_set_zero (mpc_realref (rop), (mpfr_signbit (mpc_realref (op)) ? -1 : 1)); 70 inex_re = 0; 71 mpfr_set_inf (mpc_imagref (rop), (mpfr_signbit (mpc_imagref (op)) ? -1 : 1)); 72 } 73 74 return MPC_INEX (inex_re, 0); 75 } 76 77 /* pure real argument */ 78 if (mpfr_zero_p (mpc_imagref (op))) 79 { 80 int inex_re; 81 int inex_im; 82 int s_im; 83 s_im = mpfr_signbit (mpc_imagref (op)); 84 85 if (mpfr_cmp_ui (mpc_realref (op), 1) > 0) 86 { 87 if (s_im) 88 inex_im = -mpfr_acosh (mpc_imagref (rop), mpc_realref (op), 89 INV_RND (MPC_RND_IM (rnd))); 90 else 91 inex_im = mpfr_acosh (mpc_imagref (rop), mpc_realref (op), 92 MPC_RND_IM (rnd)); 93 inex_re = set_pi_over_2 (mpc_realref (rop), 94 (mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd)); 95 if (s_im) 96 mpc_conj (rop, rop, MPC_RNDNN); 97 } 98 else if (mpfr_cmp_si (mpc_realref (op), -1) < 0) 99 { 100 mpfr_t minus_op_re; 101 minus_op_re[0] = mpc_realref (op)[0]; 102 MPFR_CHANGE_SIGN (minus_op_re); 103 104 if (s_im) 105 inex_im = -mpfr_acosh (mpc_imagref (rop), minus_op_re, 106 INV_RND (MPC_RND_IM (rnd))); 107 else 108 inex_im = mpfr_acosh (mpc_imagref (rop), minus_op_re, 109 MPC_RND_IM (rnd)); 110 inex_re = set_pi_over_2 (mpc_realref (rop), 111 (mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd)); 112 if (s_im) 113 mpc_conj (rop, rop, MPC_RNDNN); 114 } 115 else 116 { 117 inex_im = mpfr_set_ui (mpc_imagref (rop), 0, MPC_RND_IM (rnd)); 118 if (s_im) 119 mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN); 120 inex_re = mpfr_asin (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); 121 } 122 123 return MPC_INEX (inex_re, inex_im); 124 } 125 126 /* pure imaginary argument */ 127 if (mpfr_zero_p (mpc_realref (op))) 128 { 129 int inex_im; 130 int s; 131 s = mpfr_signbit (mpc_realref (op)); 132 mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); 133 if (s) 134 mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); 135 inex_im = mpfr_asinh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); 136 137 return MPC_INEX (0, inex_im); 138 } 139 140 /* regular complex: asin(z) = -i*log(i*z+sqrt(1-z^2)) */ 141 p_re = mpfr_get_prec (mpc_realref(rop)); 142 p_im = mpfr_get_prec (mpc_imagref(rop)); 143 rnd_re = MPC_RND_RE(rnd); 144 rnd_im = MPC_RND_IM(rnd); 145 p = p_re >= p_im ? p_re : p_im; 146 mpc_init2 (z1, p); 147 while (1) 148 { 149 mpfr_exp_t ex, ey, err; 150 151 p += mpc_ceil_log2 (p) + 3 + incr_p; /* incr_p is zero initially */ 152 incr_p = p / 2; 153 mpfr_set_prec (mpc_realref(z1), p); 154 mpfr_set_prec (mpc_imagref(z1), p); 155 156 /* z1 <- z^2 */ 157 mpc_sqr (z1, op, MPC_RNDNN); 158 /* err(x) <= 1/2 ulp(x), err(y) <= 1/2 ulp(y) */ 159 /* z1 <- 1-z1 */ 160 ex = mpfr_get_exp (mpc_realref(z1)); 161 mpfr_ui_sub (mpc_realref(z1), 1, mpc_realref(z1), GMP_RNDN); 162 mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN); 163 ex = ex - mpfr_get_exp (mpc_realref(z1)); 164 ex = (ex <= 0) ? 0 : ex; 165 /* err(x) <= 2^ex * ulp(x) */ 166 ex = ex + mpfr_get_exp (mpc_realref(z1)) - p; 167 /* err(x) <= 2^ex */ 168 ey = mpfr_get_exp (mpc_imagref(z1)) - p - 1; 169 /* err(y) <= 2^ey */ 170 ex = (ex >= ey) ? ex : ey; /* err(x), err(y) <= 2^ex, i.e., the norm 171 of the error is bounded by |h|<=2^(ex+1/2) */ 172 /* z1 <- sqrt(z1): if z1 = z + h, then sqrt(z1) = sqrt(z) + h/2/sqrt(t) */ 173 ey = mpfr_get_exp (mpc_realref(z1)) >= mpfr_get_exp (mpc_imagref(z1)) 174 ? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1)); 175 /* we have |z1| >= 2^(ey-1) thus 1/|z1| <= 2^(1-ey) */ 176 mpc_sqrt (z1, z1, MPC_RNDNN); 177 ex = (2 * ex + 1) - 2 - (ey - 1); /* |h^2/4/|t| <= 2^ex */ 178 ex = (ex + 1) / 2; /* ceil(ex/2) */ 179 /* express ex in terms of ulp(z1) */ 180 ey = mpfr_get_exp (mpc_realref(z1)) <= mpfr_get_exp (mpc_imagref(z1)) 181 ? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1)); 182 ex = ex - ey + p; 183 /* take into account the rounding error in the mpc_sqrt call */ 184 err = (ex <= 0) ? 1 : ex + 1; 185 /* err(x) <= 2^err * ulp(x), err(y) <= 2^err * ulp(y) */ 186 /* z1 <- i*z + z1 */ 187 ex = mpfr_get_exp (mpc_realref(z1)); 188 ey = mpfr_get_exp (mpc_imagref(z1)); 189 mpfr_sub (mpc_realref(z1), mpc_realref(z1), mpc_imagref(op), GMP_RNDN); 190 mpfr_add (mpc_imagref(z1), mpc_imagref(z1), mpc_realref(op), GMP_RNDN); 191 if (mpfr_cmp_ui (mpc_realref(z1), 0) == 0 || mpfr_cmp_ui (mpc_imagref(z1), 0) == 0) 192 continue; 193 ex -= mpfr_get_exp (mpc_realref(z1)); /* cancellation in x */ 194 ey -= mpfr_get_exp (mpc_imagref(z1)); /* cancellation in y */ 195 ex = (ex >= ey) ? ex : ey; /* maximum cancellation */ 196 err += ex; 197 err = (err <= 0) ? 1 : err + 1; /* rounding error in sub/add */ 198 /* z1 <- log(z1): if z1 = z + h, then log(z1) = log(z) + h/t with 199 |t| >= min(|z1|,|z|) */ 200 ex = mpfr_get_exp (mpc_realref(z1)); 201 ey = mpfr_get_exp (mpc_imagref(z1)); 202 ex = (ex >= ey) ? ex : ey; 203 err += ex - p; /* revert to absolute error <= 2^err */ 204 mpc_log (z1, z1, GMP_RNDN); 205 err -= ex - 1; /* 1/|t| <= 1/|z| <= 2^(1-ex) */ 206 /* express err in terms of ulp(z1) */ 207 ey = mpfr_get_exp (mpc_realref(z1)) <= mpfr_get_exp (mpc_imagref(z1)) 208 ? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1)); 209 err = err - ey + p; 210 /* take into account the rounding error in the mpc_log call */ 211 err = (err <= 0) ? 1 : err + 1; 212 /* z1 <- -i*z1 */ 213 mpfr_swap (mpc_realref(z1), mpc_imagref(z1)); 214 mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN); 215 if (mpfr_can_round (mpc_realref(z1), p - err, GMP_RNDN, GMP_RNDZ, 216 p_re + (rnd_re == GMP_RNDN)) && 217 mpfr_can_round (mpc_imagref(z1), p - err, GMP_RNDN, GMP_RNDZ, 218 p_im + (rnd_im == GMP_RNDN))) 219 break; 220 } 221 222 inex = mpc_set (rop, z1, rnd); 223 mpc_clear (z1); 224 225 return inex; 226 } 227