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