src/asin.c

   1 /* mpc_asin -- arcsine of a complex number.
   2
   3 Copyright (C) 2009 Philippe Th\'eveny
   4
   5 This file is part of the MPC Library.
   6
   7 The MPC Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 2.1 of the License, or (at your
  10 option) any later version.
  11
  12 The MPC Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the MPC Library; see the file COPYING.LIB.  If not, write to
  19 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  20 MA 02111-1307, USA. */
  21
  22 #include "mpc-impl.h"
  23
  24 extern int set_pi_over_2 (mpfr_ptr rop, int s, mpfr_rnd_t rnd);
  25
  26 int
  27 mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
  28 {
  29   mp_prec_t p, p_re, p_im;
  30   mp_rnd_t rnd_re, rnd_im;
  31   mpc_t z1;
  32   int inex;
  33
  34   /* special values */
  35   if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op)))
  36     {
  37       if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op)))
  38         {
  39           mpfr_set_nan (MPC_RE (rop));
  40           mpfr_set_inf (MPC_IM (rop), mpfr_signbit (MPC_IM (op)) ? -1 : +1);
  41         }
  42       else if (mpfr_zero_p (MPC_RE (op)))
  43         {
  44           mpfr_set (MPC_RE (rop), MPC_RE (op), GMP_RNDN);
  45           mpfr_set_nan (MPC_IM (rop));
  46         }
  47       else
  48         {
  49           mpfr_set_nan (MPC_RE (rop));
  50           mpfr_set_nan (MPC_IM (rop));
  51         }
  52
  53       return 0;
  54     }
  55
  56   if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op)))
  57     {
  58       int inex_re;
  59       if (mpfr_inf_p (MPC_RE (op)))
  60         {
  61           inex_re = set_pi_over_2 (MPC_RE (rop), -mpfr_signbit (MPC_RE (op)),
  62                                    MPC_RND_RE (rnd));
  63           mpfr_set_inf (MPC_IM (rop), -mpfr_signbit (MPC_IM (op)));
  64
  65           if (mpfr_inf_p (MPC_IM (op)))
  66             mpfr_div_2ui (MPC_RE (rop), MPC_RE (rop), 1, GMP_RNDN);
  67         }
  68       else
  69         {
  70           int s;
  71           s = mpfr_signbit (MPC_RE (op));
  72           inex_re = mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN);
  73           if (s)
  74             mpfr_neg (MPC_RE (rop), MPC_RE (rop), GMP_RNDN);
  75           mpfr_set_inf (MPC_IM (rop), -mpfr_signbit (MPC_IM (op)));
  76         }
  77
  78       return MPC_INEX (inex_re, 0);
  79     }
  80
  81   /* pure real argument */
  82   if (mpfr_zero_p (MPC_IM (op)))
  83     {
  84       int inex_re;
  85       int inex_im;
  86       int s_im;
  87       s_im = mpfr_signbit (MPC_IM (op));
  88
  89       if (mpfr_cmp_ui (MPC_RE (op), 1) > 0)
  90         {
  91           if (s_im)
  92             inex_im = -mpfr_acosh (MPC_IM (rop), MPC_RE (op),
  93                                    INV_RND (MPC_RND_IM (rnd)));
  94           else
  95             inex_im = mpfr_acosh (MPC_IM (rop), MPC_RE (op),
  96                                   MPC_RND_IM (rnd));
  97           inex_re = set_pi_over_2 (MPC_RE (rop), -mpfr_signbit (MPC_RE (op)),
  98                                    MPC_RND_RE (rnd));
  99           if (s_im)
 100             mpc_conj (rop, rop, MPC_RNDNN);
 101         }
 102       else if (mpfr_cmp_si (MPC_RE (op), -1) < 0)
 103         {
 104           mpfr_t minus_op_re;
 105           minus_op_re[0] = MPC_RE (op)[0];
 106           MPFR_CHANGE_SIGN (minus_op_re);
 107
 108           if (s_im)
 109             inex_im = -mpfr_acosh (MPC_IM (rop), minus_op_re,
 110                                    INV_RND (MPC_RND_IM (rnd)));
 111           else
 112             inex_im = mpfr_acosh (MPC_IM (rop), minus_op_re,
 113                                   MPC_RND_IM (rnd));
 114           inex_re = set_pi_over_2 (MPC_RE (rop), -mpfr_signbit (MPC_RE (op)),
 115                                    MPC_RND_RE (rnd));
 116           if (s_im)
 117             mpc_conj (rop, rop, MPC_RNDNN);
 118         }
 119       else
 120         {
 121           inex_im = mpfr_set_ui (MPC_IM (rop), 0, MPC_RND_IM (rnd));
 122           if (s_im)
 123             mpfr_neg (MPC_IM (rop), MPC_IM (rop), GMP_RNDN);
 124           inex_re = mpfr_asin (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd));
 125         }
 126
 127       return MPC_INEX (inex_re, inex_im);
 128     }
 129
 130   /* pure imaginary argument */
 131   if (mpfr_zero_p (MPC_RE (op)))
 132     {
 133       int inex_im;
 134       int s;
 135       s = mpfr_signbit (MPC_RE (op));
 136       mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN);
 137       if (s)
 138         mpfr_neg (MPC_RE (rop), MPC_RE (rop), GMP_RNDN);
 139       inex_im = mpfr_asinh (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd));
 140
 141       return MPC_INEX (0, inex_im);
 142     }
 143
 144   /* regular complex: asin(z) = -i*log(i*z+sqrt(1-z^2)) */
 145   p_re = mpfr_get_prec (MPC_RE(rop));
 146   p_im = mpfr_get_prec (MPC_IM(rop));
 147   rnd_re = MPC_RND_RE(rnd);
 148   rnd_im = MPC_RND_IM(rnd);
 149   p = p_re >= p_im ? p_re : p_im;
 150   mpc_init2 (z1, p);
 151   while (1)
 152   {
 153     mp_exp_t ex, ey, err;
 154     p += mpc_ceil_log2 (p) + 3;
 155     mpfr_set_prec (MPC_RE(z1), p);
 156     mpfr_set_prec (MPC_IM(z1), p);
 157
 158     /* z1 <- z^2 */
 159     mpc_sqr (z1, op, MPC_RNDNN);
 160     /* err(x) <= 1/2 ulp(x), err(y) <= 1/2 ulp(y) */
 161     /* z1 <- 1-z1 */
 162     ex = mpfr_get_exp (MPC_RE(z1));
 163     mpfr_ui_sub (MPC_RE(z1), 1, MPC_RE(z1), GMP_RNDN);
 164     mpfr_neg (MPC_IM(z1), MPC_IM(z1), GMP_RNDN);
 165     ex = ex - mpfr_get_exp (MPC_RE(z1));
 166     ex = (ex <= 0) ? 0 : ex;
 167     /* err(x) <= 2^ex * ulp(x) */
 168     ex = ex + mpfr_get_exp (MPC_RE(z1)) - p;
 169     /* err(x) <= 2^ex */
 170     ey = mpfr_get_exp (MPC_IM(z1)) - p - 1;
 171     /* err(y) <= 2^ey */
 172     ex = (ex >= ey) ? ex : ey; /* err(x), err(y) <= 2^ex, i.e., the norm
 173                                   of the error is bounded by |h|<=2^(ex+1/2) */
 174     /* z1 <- sqrt(z1): if z1 = z + h, then sqrt(z1) = sqrt(z) + h/2/sqrt(t) */
 175     ey = mpfr_get_exp (MPC_RE(z1)) >= mpfr_get_exp (MPC_IM(z1))
 176       ? mpfr_get_exp (MPC_RE(z1)) : mpfr_get_exp (MPC_IM(z1));
 177     /* we have |z1| >= 2^(ey-1) thus 1/|z1| <= 2^(1-ey) */
 178     mpc_sqrt (z1, z1, MPC_RNDNN);
 179     ex = (2 * ex + 1) - 2 - (ey - 1); /* |h^2/4/|t| <= 2^ex */
 180     ex = (ex + 1) / 2; /* ceil(ex/2) */
 181     /* express ex in terms of ulp(z1) */
 182     ey = mpfr_get_exp (MPC_RE(z1)) <= mpfr_get_exp (MPC_IM(z1))
 183       ? mpfr_get_exp (MPC_RE(z1)) : mpfr_get_exp (MPC_IM(z1));
 184     ex = ex - ey + p;
 185     /* take into account the rounding error in the mpc_sqrt call */
 186     err = (ex <= 0) ? 1 : ex + 1;
 187     /* err(x) <= 2^err * ulp(x), err(y) <= 2^err * ulp(y) */
 188     /* z1 <- i*z + z1 */
 189     ex = mpfr_get_exp (MPC_RE(z1));
 190     ey = mpfr_get_exp (MPC_IM(z1));
 191     mpfr_sub (MPC_RE(z1), MPC_RE(z1), MPC_IM(op), GMP_RNDN);
 192     mpfr_add (MPC_IM(z1), MPC_IM(z1), MPC_RE(op), GMP_RNDN);
 193     if (mpfr_cmp_ui (MPC_RE(z1), 0) == 0 || mpfr_cmp_ui (MPC_IM(z1), 0) == 0)
 194       continue;
 195     ex -= mpfr_get_exp (MPC_RE(z1)); /* cancellation in x */
 196     ey -= mpfr_get_exp (MPC_IM(z1)); /* cancellation in y */
 197     ex = (ex >= ey) ? ex : ey; /* maximum cancellation */
 198     err += ex;
 199     err = (err <= 0) ? 1 : err + 1; /* rounding error in sub/add */
 200     /* z1 <- log(z1): if z1 = z + h, then log(z1) = log(z) + h/t with
 201        |t| >= min(|z1|,|z|) */
 202     ex = mpfr_get_exp (MPC_RE(z1));
 203     ey = mpfr_get_exp (MPC_IM(z1));
 204     ex = (ex >= ey) ? ex : ey;
 205     err += ex - p; /* revert to absolute error <= 2^err */
 206     mpc_log (z1, z1, GMP_RNDN);
 207     err -= ex - 1; /* 1/|t| <= 1/|z| <= 2^(1-ex) */
 208     /* express err in terms of ulp(z1) */
 209     ey = mpfr_get_exp (MPC_RE(z1)) <= mpfr_get_exp (MPC_IM(z1))
 210       ? mpfr_get_exp (MPC_RE(z1)) : mpfr_get_exp (MPC_IM(z1));
 211     err = err - ey + p;
 212     /* take into account the rounding error in the mpc_log call */
 213     err = (err <= 0) ? 1 : err + 1;
 214     /* z1 <- -i*z1 */
 215     mpfr_swap (MPC_RE(z1), MPC_IM(z1));
 216     mpfr_neg (MPC_IM(z1), MPC_IM(z1), GMP_RNDN);
 217     if (mpfr_can_round (MPC_RE(z1), p - err, GMP_RNDN, GMP_RNDZ,
 218                         p_re + (rnd_re == GMP_RNDN)) &&
 219         mpfr_can_round (MPC_IM(z1), p - err, GMP_RNDN, GMP_RNDZ,
 220                         p_im + (rnd_im == GMP_RNDN)))
 221       break;
 222   }
 223
 224   inex = mpc_set (rop, z1, rnd);
 225   mpc_clear (z1);
 226
 227   return inex;
 228 }