contrib/mpc/src/atan.c

   1 /* mpc_atan -- arctangent 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>
  22 #include "mpc-impl.h"
  23
  24 /* set rop to
  25    -pi/2 if s < 0
  26    +pi/2 else
  27    rounded in the direction rnd
  28 */
  29 int
  30 set_pi_over_2 (mpfr_ptr rop, int s, mpfr_rnd_t rnd)
  31 {
  32   int inex;
  33
  34   inex = mpfr_const_pi (rop, s < 0 ? INV_RND (rnd) : rnd);
  35   mpfr_div_2ui (rop, rop, 1, GMP_RNDN);
  36   if (s < 0)
  37     {
  38       inex = -inex;
  39       mpfr_neg (rop, rop, GMP_RNDN);
  40     }
  41
  42   return inex;
  43 }
  44
  45 int
  46 mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
  47 {
  48   int s_re;
  49   int s_im;
  50   int inex_re;
  51   int inex_im;
  52   int inex;
  53
  54   inex_re = 0;
  55   inex_im = 0;
  56   s_re = mpfr_signbit (mpc_realref (op));
  57   s_im = mpfr_signbit (mpc_imagref (op));
  58
  59   /* special values */
  60   if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
  61     {
  62       if (mpfr_nan_p (mpc_realref (op)))
  63         {
  64           mpfr_set_nan (mpc_realref (rop));
  65           if (mpfr_zero_p (mpc_imagref (op)) || mpfr_inf_p (mpc_imagref (op)))
  66             {
  67               mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN);
  68               if (s_im)
  69                 mpc_conj (rop, rop, MPC_RNDNN);
  70             }
  71           else
  72             mpfr_set_nan (mpc_imagref (rop));
  73         }
  74       else
  75         {
  76           if (mpfr_inf_p (mpc_realref (op)))
  77             {
  78               inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd));
  79               mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN);
  80             }
  81           else
  82             {
  83               mpfr_set_nan (mpc_realref (rop));
  84               mpfr_set_nan (mpc_imagref (rop));
  85             }
  86         }
  87       return MPC_INEX (inex_re, 0);
  88     }
  89
  90   if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op)))
  91     {
  92       inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd));
  93
  94       mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN);
  95       if (s_im)
  96         mpc_conj (rop, rop, GMP_RNDN);
  97
  98       return MPC_INEX (inex_re, 0);
  99     }
 100
 101   /* pure real argument */
 102   if (mpfr_zero_p (mpc_imagref (op)))
 103     {
 104       inex_re = mpfr_atan (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
 105
 106       mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN);
 107       if (s_im)
 108         mpc_conj (rop, rop, GMP_RNDN);
 109
 110       return MPC_INEX (inex_re, 0);
 111     }
 112
 113   /* pure imaginary argument */
 114   if (mpfr_zero_p (mpc_realref (op)))
 115     {
 116       int cmp_1;
 117
 118       if (s_im)
 119         cmp_1 = -mpfr_cmp_si (mpc_imagref (op), -1);
 120       else
 121         cmp_1 = mpfr_cmp_ui (mpc_imagref (op), +1);
 122
 123       if (cmp_1 < 0)
 124         {
 125           /* atan(+0+iy) = +0 +i*atanh(y), if |y| < 1
 126              atan(-0+iy) = -0 +i*atanh(y), if |y| < 1 */
 127
 128           mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN);
 129           if (s_re)
 130             mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN);
 131
 132           inex_im = mpfr_atanh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd));
 133         }
 134       else if (cmp_1 == 0)
 135         {
 136           /* atan(+/-0+i) = NaN +i*inf
 137              atan(+/-0-i) = NaN -i*inf */
 138           mpfr_set_nan (mpc_realref (rop));
 139           mpfr_set_inf (mpc_imagref (rop), s_im ? -1 : +1);
 140         }
 141       else
 142         {
 143           /* atan(+0+iy) = +pi/2 +i*atanh(1/y), if |y| > 1
 144              atan(-0+iy) = -pi/2 +i*atanh(1/y), if |y| > 1 */
 145           mpfr_rnd_t rnd_im, rnd_away;
 146           mpfr_t y;
 147           mpfr_prec_t p, p_im;
 148           int ok;
 149
 150           rnd_im = MPC_RND_IM (rnd);
 151           mpfr_init (y);
 152           p_im = mpfr_get_prec (mpc_imagref (rop));
 153           p = p_im;
 154
 155           /* a = o(1/y)      with error(a) < 1 ulp(a)
 156              b = o(atanh(a)) with error(b) < (1+2^{1+Exp(a)-Exp(b)}) ulp(b)
 157
 158              As |atanh (1/y)| > |1/y| we have Exp(a)-Exp(b) <=0 so, at most,
 159              2 bits of precision are lost.
 160
 161              We round atanh(1/y) away from 0.
 162           */
 163           do
 164             {
 165               p += mpc_ceil_log2 (p) + 2;
 166               mpfr_set_prec (y, p);
 167               rnd_away = s_im == 0 ? GMP_RNDU : GMP_RNDD;
 168               inex_im = mpfr_ui_div (y, 1, mpc_imagref (op), rnd_away);
 169               /* FIXME: should we consider the case with unreasonably huge
 170                  precision prec(y)>3*exp_min, where atanh(1/Im(op)) could be
 171                  representable while 1/Im(op) underflows ?
 172                  This corresponds to |y| = 0.5*2^emin, in which case the
 173                  result may be wrong. */
 174
 175               /* atanh cannot underflow: |atanh(x)| > |x| for |x| < 1 */
 176               inex_im |= mpfr_atanh (y, y, rnd_away);
 177
 178               ok = inex_im == 0
 179                 || mpfr_can_round (y, p - 2, rnd_away, GMP_RNDZ,
 180                                    p_im + (rnd_im == GMP_RNDN));
 181             } while (ok == 0);
 182
 183           inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd));
 184           inex_im = mpfr_set (mpc_imagref (rop), y, rnd_im);
 185           mpfr_clear (y);
 186         }
 187       return MPC_INEX (inex_re, inex_im);
 188     }
 189
 190   /* regular number argument */
 191   {
 192     mpfr_t a, b, x, y;
 193     mpfr_prec_t prec, p;
 194     mpfr_exp_t err, expo;
 195     int ok = 0;
 196     mpfr_t minus_op_re;
 197     mpfr_exp_t op_re_exp, op_im_exp;
 198     mpfr_rnd_t rnd1, rnd2;
 199
 200     mpfr_inits2 (MPFR_PREC_MIN, a, b, x, y, (mpfr_ptr) 0);
 201
 202     /* real part: Re(arctan(x+i*y)) = [arctan2(x,1-y) - arctan2(-x,1+y)]/2 */
 203     minus_op_re[0] = mpc_realref (op)[0];
 204     MPFR_CHANGE_SIGN (minus_op_re);
 205     op_re_exp = mpfr_get_exp (mpc_realref (op));
 206     op_im_exp = mpfr_get_exp (mpc_imagref (op));
 207
 208     prec = mpfr_get_prec (mpc_realref (rop)); /* result precision */
 209
 210     /* a = o(1-y)         error(a) < 1 ulp(a)
 211        b = o(atan2(x,a))  error(b) < [1+2^{3+Exp(x)-Exp(a)-Exp(b)}] ulp(b)
 212                                      = kb ulp(b)
 213        c = o(1+y)         error(c) < 1 ulp(c)
 214        d = o(atan2(-x,c)) error(d) < [1+2^{3+Exp(x)-Exp(c)-Exp(d)}] ulp(d)
 215                                      = kd ulp(d)
 216        e = o(b - d)       error(e) < [1 + kb*2^{Exp(b}-Exp(e)}
 217                                         + kd*2^{Exp(d)-Exp(e)}] ulp(e)
 218                           error(e) < [1 + 2^{4+Exp(x)-Exp(a)-Exp(e)}
 219                                         + 2^{4+Exp(x)-Exp(c)-Exp(e)}] ulp(e)
 220                           because |atan(u)| < |u|
 221                                    < [1 + 2^{5+Exp(x)-min(Exp(a),Exp(c))
 222                                              -Exp(e)}] ulp(e)
 223        f = e/2            exact
 224     */
 225
 226     /* p: working precision */
 227     p = (op_im_exp > 0 || prec > SAFE_ABS (mpfr_prec_t, op_im_exp)) ? prec
 228       : (prec - op_im_exp);
 229     rnd1 = mpfr_sgn (mpc_realref (op)) > 0 ? GMP_RNDD : GMP_RNDU;
 230     rnd2 = mpfr_sgn (mpc_realref (op)) < 0 ? GMP_RNDU : GMP_RNDD;
 231
 232     do
 233       {
 234         p += mpc_ceil_log2 (p) + 2;
 235         mpfr_set_prec (a, p);
 236         mpfr_set_prec (b, p);
 237         mpfr_set_prec (x, p);
 238
 239         /* x = upper bound for atan (x/(1-y)). Since atan is increasing, we
 240            need an upper bound on x/(1-y), i.e., a lower bound on 1-y for
 241            x positive, and an upper bound on 1-y for x negative */
 242         mpfr_ui_sub (a, 1, mpc_imagref (op), rnd1);
 243         if (mpfr_sgn (a) == 0) /* y is near 1, thus 1+y is near 2, and
 244                                   expo will be 1 or 2 below */
 245           {
 246             MPC_ASSERT (mpfr_cmp_ui (mpc_imagref(op), 1) == 0);
 247                /* check for intermediate underflow */
 248             err = 2; /* ensures err will be expo below */
 249           }
 250         else
 251           err = mpfr_get_exp (a); /* err = Exp(a) with the notations above */
 252         mpfr_atan2 (x, mpc_realref (op), a, GMP_RNDU);
 253
 254         /* b = lower bound for atan (-x/(1+y)): for x negative, we need a
 255            lower bound on -x/(1+y), i.e., an upper bound on 1+y */
 256         mpfr_add_ui (a, mpc_imagref(op), 1, rnd2);
 257         /* if a is exactly zero, i.e., Im(op) = -1, then the error on a is 0,
 258            and we can simply ignore the terms involving Exp(a) in the error */
 259         if (mpfr_sgn (a) == 0)
 260           {
 261             MPC_ASSERT (mpfr_cmp_si (mpc_imagref(op), -1) == 0);
 262                /* check for intermediate underflow */
 263             expo = err; /* will leave err unchanged below */
 264           }
 265         else
 266           expo = mpfr_get_exp (a); /* expo = Exp(c) with the notations above */
 267         mpfr_atan2 (b, minus_op_re, a, GMP_RNDD);
 268
 269         err = err < expo ? err : expo; /* err = min(Exp(a),Exp(c)) */
 270         mpfr_sub (x, x, b, GMP_RNDU);
 271
 272         err = 5 + op_re_exp - err - mpfr_get_exp (x);
 273         /* error is bounded by [1 + 2^err] ulp(e) */
 274         err = err < 0 ? 1 : err + 1;
 275
 276         mpfr_div_2ui (x, x, 1, GMP_RNDU);
 277
 278         /* Note: using RND2=RNDD guarantees that if x is exactly representable
 279            on prec + ... bits, mpfr_can_round will return 0 */
 280         ok = mpfr_can_round (x, p - err, GMP_RNDU, GMP_RNDD,
 281                              prec + (MPC_RND_RE (rnd) == GMP_RNDN));
 282       } while (ok == 0);
 283
 284     /* Imaginary part
 285        Im(atan(x+I*y)) = 1/4 * [log(x^2+(1+y)^2) - log (x^2 +(1-y)^2)] */
 286     prec = mpfr_get_prec (mpc_imagref (rop)); /* result precision */
 287
 288     /* a = o(1+y)    error(a) < 1 ulp(a)
 289        b = o(a^2)    error(b) < 5 ulp(b)
 290        c = o(x^2)    error(c) < 1 ulp(c)
 291        d = o(b+c)    error(d) < 7 ulp(d)
 292        e = o(log(d)) error(e) < [1 + 7*2^{2-Exp(e)}] ulp(e) = ke ulp(e)
 293        f = o(1-y)    error(f) < 1 ulp(f)
 294        g = o(f^2)    error(g) < 5 ulp(g)
 295        h = o(c+f)    error(h) < 7 ulp(h)
 296        i = o(log(h)) error(i) < [1 + 7*2^{2-Exp(i)}] ulp(i) = ki ulp(i)
 297        j = o(e-i)    error(j) < [1 + ke*2^{Exp(e)-Exp(j)}
 298                                    + ki*2^{Exp(i)-Exp(j)}] ulp(j)
 299                      error(j) < [1 + 2^{Exp(e)-Exp(j)} + 2^{Exp(i)-Exp(j)}
 300                                    + 7*2^{3-Exp(j)}] ulp(j)
 301                               < [1 + 2^{max(Exp(e),Exp(i))-Exp(j)+1}
 302                                    + 7*2^{3-Exp(j)}] ulp(j)
 303        k = j/4       exact
 304     */
 305     err = 2;
 306     p = prec; /* working precision */
 307
 308     do
 309       {
 310         p += mpc_ceil_log2 (p) + err;
 311         mpfr_set_prec (a, p);
 312         mpfr_set_prec (b, p);
 313         mpfr_set_prec (y, p);
 314
 315         /* a = upper bound for log(x^2 + (1+y)^2) */
 316         ROUND_AWAY (mpfr_add_ui (a, mpc_imagref (op), 1, MPFR_RNDA), a);
 317         mpfr_sqr (a, a, GMP_RNDU);
 318         mpfr_sqr (y, mpc_realref (op), GMP_RNDU);
 319         mpfr_add (a, a, y, GMP_RNDU);
 320         mpfr_log (a, a, GMP_RNDU);
 321
 322         /* b = lower bound for log(x^2 + (1-y)^2) */
 323         mpfr_ui_sub (b, 1, mpc_imagref (op), GMP_RNDZ); /* round to zero */
 324         mpfr_sqr (b, b, GMP_RNDZ);
 325         /* we could write mpfr_sqr (y, mpc_realref (op), GMP_RNDZ) but it is
 326            more efficient to reuse the value of y (x^2) above and subtract
 327            one ulp */
 328         mpfr_nextbelow (y);
 329         mpfr_add (b, b, y, GMP_RNDZ);
 330         mpfr_log (b, b, GMP_RNDZ);
 331
 332         mpfr_sub (y, a, b, GMP_RNDU);
 333
 334         if (mpfr_zero_p (y))
 335            /* FIXME: happens when x and y have very different magnitudes;
 336               could be handled more efficiently                           */
 337           ok = 0;
 338         else
 339           {
 340             expo = MPC_MAX (mpfr_get_exp (a), mpfr_get_exp (b));
 341             expo = expo - mpfr_get_exp (y) + 1;
 342             err = 3 - mpfr_get_exp (y);
 343             /* error(j) <= [1 + 2^expo + 7*2^err] ulp(j) */
 344             if (expo <= err) /* error(j) <= [1 + 2^{err+1}] ulp(j) */
 345               err = (err < 0) ? 1 : err + 2;
 346             else
 347               err = (expo < 0) ? 1 : expo + 2;
 348
 349             mpfr_div_2ui (y, y, 2, GMP_RNDN);
 350             MPC_ASSERT (!mpfr_zero_p (y));
 351                /* FIXME: underflow. Since the main term of the Taylor series
 352                   in y=0 is 1/(x^2+1) * y, this means that y is very small
 353                   and/or x very large; but then the mpfr_zero_p (y) above
 354                   should be true. This needs a proof, or better yet,
 355                   special code.                                              */
 356
 357             ok = mpfr_can_round (y, p - err, GMP_RNDU, GMP_RNDD,
 358                                  prec + (MPC_RND_IM (rnd) == GMP_RNDN));
 359           }
 360       } while (ok == 0);
 361
 362     inex = mpc_set_fr_fr (rop, x, y, rnd);
 363
 364     mpfr_clears (a, b, x, y, (mpfr_ptr) 0);
 365     return inex;
 366   }
 367 }