1/* mpc_asin -- arcsine of a complex number.
2
3Copyright (C) 2009, 2010, 2011 INRIA
4
5This file is part of GNU MPC.
6
7GNU MPC is free software; you can redistribute it and/or modify it under
8the terms of the GNU Lesser General Public License as published by the
9Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
13WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along with this program. If not, see http://www.gnu.org/licenses/ .
19*/
20
21#include "mpc-impl.h"
22
23int
24mpc_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