contrib/mpc/src/acos.c

   1 /* mpc_acos -- arccosine of a complex number.
   2
   3 Copyright (C) 2009, 2010, 2011, 2012 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 <stdio.h>    /* for MPC_ASSERT */
  22 #include "mpc-impl.h"
  23
  24 int
  25 mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
  26 {
  27   int inex_re, inex_im, inex;
  28   mpfr_prec_t p_re, p_im, p;
  29   mpc_t z1;
  30   mpfr_t pi_over_2;
  31   mpfr_exp_t e1, e2;
  32   mpfr_rnd_t rnd_im;
  33   mpc_rnd_t rnd1;
  34
  35   inex_re = 0;
  36   inex_im = 0;
  37
  38   /* special values */
  39   if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
  40     {
  41       if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op)))
  42         {
  43           mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? +1 : -1);
  44           mpfr_set_nan (mpc_realref (rop));
  45         }
  46       else if (mpfr_zero_p (mpc_realref (op)))
  47         {
  48           inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd));
  49           mpfr_set_nan (mpc_imagref (rop));
  50         }
  51       else
  52         {
  53           mpfr_set_nan (mpc_realref (rop));
  54           mpfr_set_nan (mpc_imagref (rop));
  55         }
  56
  57       return MPC_INEX (inex_re, 0);
  58     }
  59
  60   if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op)))
  61     {
  62       if (mpfr_inf_p (mpc_realref (op)))
  63         {
  64           if (mpfr_inf_p (mpc_imagref (op)))
  65             {
  66               if (mpfr_sgn (mpc_realref (op)) > 0)
  67                 {
  68                   inex_re =
  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);
  71                 }
  72               else
  73                 {
  74
  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)
  78                      c = b/4    exact
  79                      thus 1 bit is lost */
  80                   mpfr_t x;
  81                   mpfr_prec_t prec;
  82                   int ok;
  83                   mpfr_init (x);
  84                   prec = mpfr_get_prec (mpc_realref (rop));
  85                   p = prec;
  86
  87                   do
  88                     {
  89                       p += mpc_ceil_log2 (p);
  90                       mpfr_set_prec (x, p);
  91                       mpfr_const_pi (x, GMP_RNDD);
  92                       mpfr_mul_ui (x, x, 3, GMP_RNDD);
  93                       ok =
  94                         mpfr_can_round (x, p - 1, GMP_RNDD, MPC_RND_RE (rnd),
  95                                         prec+(MPC_RND_RE (rnd) == GMP_RNDN));
  96
  97                     } while (ok == 0);
  98                   inex_re =
  99                     mpfr_div_2ui (mpc_realref (rop), x, 2, MPC_RND_RE (rnd));
 100                   mpfr_clear (x);
 101                 }
 102             }
 103           else
 104             {
 105               if (mpfr_sgn (mpc_realref (op)) > 0)
 106                 mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN);
 107               else
 108                 inex_re = mpfr_const_pi (mpc_realref (rop), MPC_RND_RE (rnd));
 109             }
 110         }
 111       else
 112         inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd));
 113
 114       mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? +1 : -1);
 115
 116       return MPC_INEX (inex_re, 0);
 117     }
 118
 119   /* pure real argument */
 120   if (mpfr_zero_p (mpc_imagref (op)))
 121     {
 122       int s_im;
 123       s_im = mpfr_signbit (mpc_imagref (op));
 124
 125       if (mpfr_cmp_ui (mpc_realref (op), 1) > 0)
 126         {
 127           if (s_im)
 128             inex_im = mpfr_acosh (mpc_imagref (rop), mpc_realref (op),
 129                                   MPC_RND_IM (rnd));
 130           else
 131             inex_im = -mpfr_acosh (mpc_imagref (rop), mpc_realref (op),
 132                                    INV_RND (MPC_RND_IM (rnd)));
 133
 134           mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN);
 135         }
 136       else if (mpfr_cmp_si (mpc_realref (op), -1) < 0)
 137         {
 138           mpfr_t minus_op_re;
 139           minus_op_re[0] = mpc_realref (op)[0];
 140           MPFR_CHANGE_SIGN (minus_op_re);
 141
 142           if (s_im)
 143             inex_im = mpfr_acosh (mpc_imagref (rop), minus_op_re,
 144                                   MPC_RND_IM (rnd));
 145           else
 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));
 149         }
 150       else
 151         {
 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));
 154         }
 155
 156       if (!s_im)
 157         mpc_conj (rop, rop, MPC_RNDNN);
 158
 159       return MPC_INEX (inex_re, inex_im);
 160     }
 161
 162   /* pure imaginary argument */
 163   if (mpfr_zero_p (mpc_realref (op)))
 164     {
 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);
 169
 170       return MPC_INEX (inex_re, inex_im);
 171     }
 172
 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));
 176   p = p_re;
 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;
 185   else
 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);
 191   for (;;)
 192     {
 193       p += mpc_ceil_log2 (p) + 3;
 194
 195       mpfr_set_prec (mpc_realref(z1), p);
 196       mpfr_set_prec (pi_over_2, p);
 197
 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)))
 206         {
 207           /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-p-1) +
 208              2^(e2-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 */
 216           inex_im = -inex_im;
 217           if (mpfr_can_round (mpc_realref(z1), p - e1, GMP_RNDN, GMP_RNDZ,
 218                               p_re + (MPC_RND_RE(rnd) == GMP_RNDN)))
 219             break;
 220         }
 221     }
 222   inex = mpc_set (rop, z1, rnd);
 223   inex_re = MPC_INEX_RE(inex);
 224   mpc_clear (z1);
 225   mpfr_clear (pi_over_2);
 226
 227   return MPC_INEX(inex_re, inex_im);
 228 }