contrib/mpc/src/asin.c

   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
  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 }