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
mpc_asin(mpc_ptr rop,mpc_srcptr op,mpc_rnd_t rnd)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