1 /* mpc_acos -- arccosine of a complex number.
3 Copyright (C) 2009, 2010, 2011, 2012 INRIA
5 This file is part of GNU MPC.
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.
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
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/ .
21 #include <stdio.h> /* for MPC_ASSERT */
25 mpc_acos (mpc_ptr rop
, mpc_srcptr op
, mpc_rnd_t rnd
)
27 int inex_re
, inex_im
, inex
;
28 mpfr_prec_t p_re
, p_im
, p
;
39 if (mpfr_nan_p (mpc_realref (op
)) || mpfr_nan_p (mpc_imagref (op
)))
41 if (mpfr_inf_p (mpc_realref (op
)) || mpfr_inf_p (mpc_imagref (op
)))
43 mpfr_set_inf (mpc_imagref (rop
), mpfr_signbit (mpc_imagref (op
)) ? +1 : -1);
44 mpfr_set_nan (mpc_realref (rop
));
46 else if (mpfr_zero_p (mpc_realref (op
)))
48 inex_re
= set_pi_over_2 (mpc_realref (rop
), +1, MPC_RND_RE (rnd
));
49 mpfr_set_nan (mpc_imagref (rop
));
53 mpfr_set_nan (mpc_realref (rop
));
54 mpfr_set_nan (mpc_imagref (rop
));
57 return MPC_INEX (inex_re
, 0);
60 if (mpfr_inf_p (mpc_realref (op
)) || mpfr_inf_p (mpc_imagref (op
)))
62 if (mpfr_inf_p (mpc_realref (op
)))
64 if (mpfr_inf_p (mpc_imagref (op
)))
66 if (mpfr_sgn (mpc_realref (op
)) > 0)
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
);
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)
84 prec
= mpfr_get_prec (mpc_realref (rop
));
89 p
+= mpc_ceil_log2 (p
);
91 mpfr_const_pi (x
, GMP_RNDD
);
92 mpfr_mul_ui (x
, x
, 3, GMP_RNDD
);
94 mpfr_can_round (x
, p
- 1, GMP_RNDD
, MPC_RND_RE (rnd
),
95 prec
+(MPC_RND_RE (rnd
) == GMP_RNDN
));
99 mpfr_div_2ui (mpc_realref (rop
), x
, 2, MPC_RND_RE (rnd
));
105 if (mpfr_sgn (mpc_realref (op
)) > 0)
106 mpfr_set_ui (mpc_realref (rop
), 0, GMP_RNDN
);
108 inex_re
= mpfr_const_pi (mpc_realref (rop
), MPC_RND_RE (rnd
));
112 inex_re
= set_pi_over_2 (mpc_realref (rop
), +1, MPC_RND_RE (rnd
));
114 mpfr_set_inf (mpc_imagref (rop
), mpfr_signbit (mpc_imagref (op
)) ? +1 : -1);
116 return MPC_INEX (inex_re
, 0);
119 /* pure real argument */
120 if (mpfr_zero_p (mpc_imagref (op
)))
123 s_im
= mpfr_signbit (mpc_imagref (op
));
125 if (mpfr_cmp_ui (mpc_realref (op
), 1) > 0)
128 inex_im
= mpfr_acosh (mpc_imagref (rop
), mpc_realref (op
),
131 inex_im
= -mpfr_acosh (mpc_imagref (rop
), mpc_realref (op
),
132 INV_RND (MPC_RND_IM (rnd
)));
134 mpfr_set_ui (mpc_realref (rop
), 0, GMP_RNDN
);
136 else if (mpfr_cmp_si (mpc_realref (op
), -1) < 0)
139 minus_op_re
[0] = mpc_realref (op
)[0];
140 MPFR_CHANGE_SIGN (minus_op_re
);
143 inex_im
= mpfr_acosh (mpc_imagref (rop
), minus_op_re
,
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
));
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
));
157 mpc_conj (rop
, rop
, MPC_RNDNN
);
159 return MPC_INEX (inex_re
, inex_im
);
162 /* pure imaginary argument */
163 if (mpfr_zero_p (mpc_realref (op
)))
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
);
170 return MPC_INEX (inex_re
, inex_im
);
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
));
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
;
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
);
193 p
+= mpc_ceil_log2 (p
) + 3;
195 mpfr_set_prec (mpc_realref(z1
), p
);
196 mpfr_set_prec (pi_over_2
, p
);
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
)))
207 /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-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 */
217 if (mpfr_can_round (mpc_realref(z1
), p
- e1
, GMP_RNDN
, GMP_RNDZ
,
218 p_re
+ (MPC_RND_RE(rnd
) == GMP_RNDN
)))
222 inex
= mpc_set (rop
, z1
, rnd
);
223 inex_re
= MPC_INEX_RE(inex
);
225 mpfr_clear (pi_over_2
);
227 return MPC_INEX(inex_re
, inex_im
);