src/atan.c

   1 /* mpc_atan -- arctangent of a complex number.
   2
   3 Copyright (C) 2009 Philippe Th\'eveny, Paul Zimmermann
   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 /* 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_RE (op));
  57   s_im = mpfr_signbit (MPC_IM (op));
  58
  59   /* special values */
  60   if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op)))
  61     {
  62       if (mpfr_nan_p (MPC_RE (op)))
  63         {
  64           mpfr_set_nan (MPC_RE (rop));
  65           if (mpfr_zero_p (MPC_IM (op)) || mpfr_inf_p (MPC_IM (op)))
  66             {
  67               mpfr_set_ui (MPC_IM (rop), 0, GMP_RNDN);
  68               if (s_im)
  69                 mpc_conj (rop, rop, MPC_RNDNN);
  70             }
  71           else
  72             mpfr_set_nan (MPC_IM (rop));
  73         }
  74       else
  75         {
  76           if (mpfr_inf_p (MPC_RE (op)))
  77             {
  78               inex_re = set_pi_over_2 (MPC_RE (rop), -s_re, MPC_RND_RE (rnd));
  79               mpfr_set_ui (MPC_IM (rop), 0, GMP_RNDN);
  80             }
  81           else
  82             {
  83               mpfr_set_nan (MPC_RE (rop));
  84               mpfr_set_nan (MPC_IM (rop));
  85             }
  86         }
  87       return MPC_INEX (inex_re, 0);
  88     }
  89
  90   if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op)))
  91     {
  92       inex_re = set_pi_over_2 (MPC_RE (rop), -s_re, MPC_RND_RE (rnd));
  93
  94       mpfr_set_ui (MPC_IM (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_IM (op)))
 103     {
 104       inex_re = mpfr_atan (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd));
 105
 106       mpfr_set_ui (MPC_IM (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_RE (op)))
 115     {
 116       int cmp_1;
 117
 118       if (s_im)
 119         cmp_1 = -mpfr_cmp_si (MPC_IM (op), -1);
 120       else
 121         cmp_1 = mpfr_cmp_ui (MPC_IM (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_RE (rop), 0, GMP_RNDN);
 129           if (s_re)
 130             mpfr_neg (MPC_RE (rop), MPC_RE (rop), GMP_RNDN);
 131
 132           inex_im = mpfr_atanh (MPC_IM (rop), MPC_IM (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_RE (rop));
 139           mpfr_set_inf (MPC_IM (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           mp_rnd_t rnd_im, rnd_away;
 146           mpfr_t y;
 147           mp_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_IM (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_IM (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_RE (rop), -s_re, MPC_RND_RE (rnd));
 184           inex_im = mpfr_set (MPC_IM (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     mp_prec_t prec, p;
 194     mp_exp_t err, expo;
 195     int ok = 0;
 196     mpfr_t minus_op_re;
 197     mp_exp_t op_re_exp;
 198     mp_exp_t op_im_exp;
 199     mp_rnd_t rnd1, rnd2;
 200
 201     mpfr_inits (a, b, x, y, (mpfr_ptr) 0);
 202
 203     /* real part: Re(arctan(x+i*y)) = [arctan2(x,1-y) - arctan2(-x,1+y)]/2 */
 204     minus_op_re[0] = MPC_RE (op)[0];
 205     MPFR_CHANGE_SIGN (minus_op_re);
 206     op_re_exp = MPFR_EXP (MPC_RE (op));
 207     op_im_exp = MPFR_EXP (MPC_IM (op));
 208
 209     prec = mpfr_get_prec (MPC_RE (rop)); /* result precision */
 210
 211     /* a = o(1-y)         error(a) < 1 ulp(a)
 212        b = o(atan2(x,a))  error(b) < [1+2^{3+Exp(x)-Exp(a)-Exp(b)}] ulp(b)
 213                                      = kb ulp(b)
 214        c = o(1+y)         error(c) < 1 ulp(c)
 215        d = o(atan2(-x,c)) error(d) < [1+2^{3+Exp(x)-Exp(c)-Exp(d)}] ulp(d)
 216                                      = kd ulp(d)
 217        e = o(b - d)       error(e) < [1 + kb*2^{Exp(b}-Exp(e)}
 218                                         + kd*2^{Exp(d)-Exp(e)}] ulp(e)
 219                           error(e) < [1 + 2^{4+Exp(x)-Exp(a)-Exp(e)}
 220                                         + 2^{4+Exp(x)-Exp(c)-Exp(e)}] ulp(e)
 221                           because |atan(u)| < |u|
 222                                    < [1 + 2^{5+Exp(x)-min(Exp(a),Exp(c))
 223                                              -Exp(e)}] ulp(e)
 224        f = e/2            exact
 225     */
 226
 227     /* p: working precision */
 228     p = (op_im_exp > 0 || prec > SAFE_ABS (mp_prec_t, op_im_exp)) ? prec
 229       : (prec - op_im_exp);
 230     rnd1 = mpfr_sgn (MPC_RE (op)) > 0 ? GMP_RNDD : GMP_RNDU;
 231     rnd2 = mpfr_sgn (MPC_RE (op)) < 0 ? GMP_RNDU : GMP_RNDD;
 232
 233     do
 234       {
 235         p += mpc_ceil_log2 (p) + 2;
 236         mpfr_set_prec (a, p);
 237         mpfr_set_prec (b, p);
 238         mpfr_set_prec (x, p);
 239
 240         /* x = upper bound for atan (x/(1-y)). Since atan is increasing, we
 241            need an upper bound on x/(1-y), i.e., a lower bound on 1-y for
 242            x positive, and an upper bound on 1-y for x negative */
 243         mpfr_ui_sub (a, 1, MPC_IM (op), rnd1);
 244         if (mpfr_sgn (a) == 0) /* y is near 1, thus 1+y is near 2, and
 245                                   expo will be 1 or 2 below */
 246           {
 247             if (mpfr_cmp_ui (MPC_IM(op), 1) != 0)
 248               continue;
 249             err = 2; /* ensures err will be expo below */
 250           }
 251         else
 252           err = MPFR_EXP (a); /* err = Exp(a) with the notations above */
 253         mpfr_atan2 (x, MPC_RE (op), a, GMP_RNDU);
 254
 255         /* b = lower bound for atan (-x/(1+y)): for x negative, we need a
 256            lower bound on -x/(1+y), i.e., an upper bound on 1+y */
 257         mpfr_add_ui (a, MPC_IM(op), 1, rnd2);
 258         /* if a is zero but inexact, try again with a larger precision
 259            if a is exactly zero, i.e., Im(op) = -1, then the error on a is 0,
 260            and we can simply ignore the terms involving Exp(a) in the error */
 261         if (mpfr_sgn (a) == 0)
 262           {
 263             if (mpfr_cmp_si (MPC_IM(op), -1) != 0)
 264               continue;
 265             expo = err; /* will leave err unchanged below */
 266           }
 267         else
 268           expo = MPFR_EXP (a); /* expo = Exp(c) with the notations above */
 269         mpfr_atan2 (b, minus_op_re, a, GMP_RNDD);
 270
 271         err = err < expo ? err : expo; /* err = min(Exp(a),Exp(c)) */
 272         mpfr_sub (x, x, b, GMP_RNDU);
 273
 274         err = 5 + op_re_exp - err - MPFR_EXP (x);
 275         /* error is bounded by [1 + 2^err] ulp(e) */
 276         err = err < 0 ? 1 : err + 1;
 277
 278         mpfr_div_2ui (x, x, 1, GMP_RNDU);
 279
 280         /* Note: using RND2=RNDD guarantees that if x is exactly representable
 281            on prec + ... bits, mpfr_can_round will return 0 */
 282         ok = mpfr_can_round (x, p - err, GMP_RNDU, GMP_RNDD,
 283                              prec + (MPC_RND_RE (rnd) == GMP_RNDN));
 284       } while (ok == 0);
 285
 286     /* Imaginary part
 287        Im(atan(x+I*y)) = 1/4 * [log(x^2+(1+y)^2) - log (x^2 +(1-y)^2)] */
 288     prec = mpfr_get_prec (MPC_IM (rop)); /* result precision */
 289
 290     /* a = o(1+y)    error(a) < 1 ulp(a)
 291        b = o(a^2)    error(b) < 5 ulp(b)
 292        c = o(x^2)    error(c) < 1 ulp(c)
 293        d = o(b+c)    error(d) < 7 ulp(d)
 294        e = o(log(d)) error(e) < [1 + 7*2^{2-Exp(e)}] ulp(e) = ke ulp(e)
 295        f = o(1-y)    error(f) < 1 ulp(f)
 296        g = o(f^2)    error(g) < 5 ulp(g)
 297        h = o(c+f)    error(h) < 7 ulp(h)
 298        i = o(log(h)) error(i) < [1 + 7*2^{2-Exp(i)}] ulp(i) = ki ulp(i)
 299        j = o(e-i)    error(j) < [1 + ke*2^{Exp(e)-Exp(j)}
 300                                    + ki*2^{Exp(i)-Exp(j)}] ulp(j)
 301                      error(j) < [1 + 2^{Exp(e)-Exp(j)} + 2^{Exp(i)-Exp(j)}
 302                                    + 7*2^{3-Exp(j)}] ulp(j)
 303                               < [1 + 2^{max(Exp(e),Exp(i))-Exp(j)+1}
 304                                    + 7*2^{3-Exp(j)}] ulp(j)
 305        k = j/4       exact
 306     */
 307     err = 2;
 308     p = prec; /* working precision */
 309     rnd1 = mpfr_cmp_si (MPC_IM (op), -1) > 0 ? GMP_RNDU : GMP_RNDD;
 310
 311     do
 312       {
 313         p += mpc_ceil_log2 (p) + err;
 314         mpfr_set_prec (a, p);
 315         mpfr_set_prec (b, p);
 316         mpfr_set_prec (y, p);
 317
 318         /* a = upper bound for log(x^2 + (1+y)^2) */
 319         mpfr_add_ui (a, MPC_IM (op), 1, rnd1);
 320         mpfr_sqr (a, a, GMP_RNDU);
 321         mpfr_sqr (y, MPC_RE (op), GMP_RNDU);
 322         mpfr_add (a, a, y, GMP_RNDU);
 323         mpfr_log (a, a, GMP_RNDU);
 324
 325         /* b = lower bound for log(x^2 + (1-y)^2) */
 326         mpfr_ui_sub (b, 1, MPC_IM (op), GMP_RNDZ);
 327         mpfr_sqr (b, b, GMP_RNDU);
 328         /* mpfr_sqr (y, MPC_RE (op), GMP_RNDZ); */
 329         mpfr_nextbelow (y);
 330         mpfr_add (b, b, y, GMP_RNDZ);
 331         mpfr_log (b, b, GMP_RNDZ);
 332
 333         mpfr_sub (y, a, b, GMP_RNDU);
 334
 335         expo = MPFR_EXP (a) < MPFR_EXP (b) ? MPFR_EXP (b) : MPFR_EXP (a);
 336         expo = expo - MPFR_EXP (y) + 1;
 337         err = 3 - MPFR_EXP (y);
 338         /* error(j) <= [1 + 2^expo + 7*2^err] ulp(j) */
 339         if (expo <= err) /* error(j) <= [1 + 2^{err+1}] ulp(j) */
 340           err = (err < 0) ? 1 : err + 2;
 341         else
 342           err = (expo < 0) ? 1 : expo + 2;
 343
 344         mpfr_div_2ui (y, y, 2, GMP_RNDN);
 345         if (mpfr_zero_p (y))
 346           err = 2; /* underflow */
 347
 348         ok = mpfr_can_round (y, p - err, GMP_RNDU, GMP_RNDD,
 349                              prec + (MPC_RND_IM (rnd) == GMP_RNDN));
 350       } while (ok == 0);
 351
 352     inex = mpc_set_fr_fr (rop, x, y, rnd);
 353
 354     mpfr_clears (a, b, x, y, (mpfr_ptr) 0);
 355     return inex;
 356   }
 357 }