contrib/mpfr/atan2.c

   1 /* mpfr_atan2 -- arc-tan 2 of a floating-point number
   2
   3 Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library, and was contributed by Mathieu Dutour.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 2.1 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 int
  27 mpfr_atan2 (mpfr_ptr dest, mpfr_srcptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
  28 {
  29   mpfr_t tmp, pi;
  30   int inexact;
  31   mp_prec_t prec;
  32   mp_exp_t e;
  33   MPFR_SAVE_EXPO_DECL (expo);
  34   MPFR_ZIV_DECL (loop);
  35
  36   MPFR_LOG_FUNC (("y[%#R]=%R x[%#R]=%R rnd=%d", y, y, x, x, rnd_mode),
  37                  ("atan[%#R]=%R inexact=%d", dest, dest, inexact));
  38
  39   /* Special cases */
  40   if (MPFR_ARE_SINGULAR (x, y))
  41     {
  42       /* atan2(0, 0) does not raise the "invalid" floating-point
  43          exception, nor does atan2(y, 0) raise the "divide-by-zero"
  44          floating-point exception.
  45          -- atan2(±0, -0) returns ±pi.313)
  46          -- atan2(±0, +0) returns ±0.
  47          -- atan2(±0, x) returns ±pi, for x < 0.
  48          -- atan2(±0, x) returns ±0, for x > 0.
  49          -- atan2(y, ±0) returns -pi/2 for y < 0.
  50          -- atan2(y, ±0) returns pi/2 for y > 0.
  51          -- atan2(±oo, -oo) returns ±3pi/4.
  52          -- atan2(±oo, +oo) returns ±pi/4.
  53          -- atan2(±oo, x) returns ±pi/2, for finite x.
  54          -- atan2(±y, -oo) returns ±pi, for finite y > 0.
  55          -- atan2(±y, +oo) returns ±0, for finite y > 0.
  56       */
  57       if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y))
  58         {
  59           MPFR_SET_NAN (dest);
  60           MPFR_RET_NAN;
  61         }
  62       if (MPFR_IS_ZERO (y))
  63         {
  64           if (MPFR_IS_NEG (x)) /* +/- PI */
  65             {
  66             set_pi:
  67               if (MPFR_IS_NEG (y))
  68                 {
  69                   inexact =  mpfr_const_pi (dest, MPFR_INVERT_RND (rnd_mode));
  70                   MPFR_CHANGE_SIGN (dest);
  71                   return -inexact;
  72                 }
  73               else
  74                 return mpfr_const_pi (dest, rnd_mode);
  75             }
  76           else /* +/- 0 */
  77             {
  78             set_zero:
  79               MPFR_SET_ZERO (dest);
  80               MPFR_SET_SAME_SIGN (dest, y);
  81               return 0;
  82             }
  83         }
  84       if (MPFR_IS_ZERO (x))
  85         {
  86         set_pi_2:
  87           if (MPFR_IS_NEG (y)) /* -PI/2 */
  88             {
  89               inexact = mpfr_const_pi (dest, MPFR_INVERT_RND(rnd_mode));
  90               MPFR_CHANGE_SIGN (dest);
  91               mpfr_div_2ui (dest, dest, 1, rnd_mode);
  92               return -inexact;
  93             }
  94           else /* PI/2 */
  95             {
  96               inexact = mpfr_const_pi (dest, rnd_mode);
  97               mpfr_div_2ui (dest, dest, 1, rnd_mode);
  98               return inexact;
  99             }
 100         }
 101       if (MPFR_IS_INF (y))
 102         {
 103           if (!MPFR_IS_INF (x)) /* +/- PI/2 */
 104             goto set_pi_2;
 105           else if (MPFR_IS_POS (x)) /* +/- PI/4 */
 106             {
 107               if (MPFR_IS_NEG (y))
 108                 {
 109                   rnd_mode = MPFR_INVERT_RND (rnd_mode);
 110                   inexact = mpfr_const_pi (dest, rnd_mode);
 111                   MPFR_CHANGE_SIGN (dest);
 112                   mpfr_div_2ui (dest, dest, 2, rnd_mode);
 113                   return -inexact;
 114                 }
 115               else
 116                 {
 117                   inexact = mpfr_const_pi (dest, rnd_mode);
 118                   mpfr_div_2ui (dest, dest, 2, rnd_mode);
 119                   return inexact;
 120                 }
 121             }
 122           else /* +/- 3*PI/4: Ugly since we have to round properly */
 123             {
 124               mpfr_t tmp2;
 125               MPFR_ZIV_DECL (loop2);
 126               mp_prec_t prec2 = MPFR_PREC (dest) + BITS_PER_MP_LIMB;
 127
 128               mpfr_init2 (tmp2, prec2);
 129               MPFR_ZIV_INIT (loop2, prec2);
 130               for (;;)
 131                 {
 132                   mpfr_const_pi (tmp2, GMP_RNDN);
 133                   mpfr_mul_ui (tmp2, tmp2, 3, GMP_RNDN); /* Error <= 2  */
 134                   mpfr_div_2ui (tmp2, tmp2, 2, GMP_RNDN);
 135                   if (mpfr_round_p (MPFR_MANT (tmp2), MPFR_LIMB_SIZE (tmp2),
 136                                     MPFR_PREC (tmp2) - 2,
 137                                     MPFR_PREC (dest) + (rnd_mode == GMP_RNDN)))
 138                     break;
 139                   MPFR_ZIV_NEXT (loop2, prec2);
 140                   mpfr_set_prec (tmp2, prec2);
 141                 }
 142               MPFR_ZIV_FREE (loop2);
 143               if (MPFR_IS_NEG (y))
 144                 MPFR_CHANGE_SIGN (tmp2);
 145               inexact = mpfr_set (dest, tmp2, rnd_mode);
 146               mpfr_clear (tmp2);
 147               return inexact;
 148             }
 149         }
 150       MPFR_ASSERTD (MPFR_IS_INF (x));
 151       if (MPFR_IS_NEG (x))
 152         goto set_pi;
 153       else
 154         goto set_zero;
 155     }
 156
 157   /* When x=1, atan2(y,x) = atan(y). FIXME: more generally, if x is a power
 158      of two, we could call directly atan(y/x) since y/x is exact. */
 159   if (mpfr_cmp_ui (x, 1) == 0)
 160     return mpfr_atan (dest, y, rnd_mode);
 161
 162   MPFR_SAVE_EXPO_MARK (expo);
 163
 164   /* Set up initial prec */
 165   prec = MPFR_PREC (dest) + 3 + MPFR_INT_CEIL_LOG2 (MPFR_PREC (dest));
 166   mpfr_init2 (tmp, prec);
 167
 168   MPFR_ZIV_INIT (loop, prec);
 169   if (MPFR_IS_POS (x))
 170     /* use atan2(y,x) = atan(y/x) */
 171     for (;;)
 172       {
 173         int div_inex;
 174         MPFR_BLOCK_DECL (flags);
 175
 176         MPFR_BLOCK (flags, div_inex = mpfr_div (tmp, y, x, GMP_RNDN));
 177         if (div_inex == 0)
 178           {
 179             /* Result is exact. */
 180             inexact = mpfr_atan (dest, tmp, rnd_mode);
 181             goto end;
 182           }
 183
 184         /* Error <= ulp (tmp) except in case of underflow or overflow. */
 185
 186         /* If the division underflowed, since |atan(z)/z| < 1, we have
 187            an underflow. */
 188         if (MPFR_UNDERFLOW (flags))
 189           {
 190             int sign;
 191
 192             /* In the case GMP_RNDN with 2^(emin-2) < |y/x| < 2^(emin-1):
 193                The smallest significand value S > 1 of |y/x| is:
 194                  * 1 / (1 - 2^(-px))                        if py <= px,
 195                  * (1 - 2^(-px) + 2^(-py)) / (1 - 2^(-px))  if py >= px.
 196                Therefore S - 1 > 2^(-pz), where pz = max(px,py). We have:
 197                atan(|y/x|) > atan(z), where z = 2^(emin-2) * (1 + 2^(-pz)).
 198                            > z - z^3 / 3.
 199                            > 2^(emin-2) * (1 + 2^(-pz) - 2^(2 emin - 5))
 200                Assuming pz <= -2 emin + 5, we can round away from zero
 201                (this is what mpfr_underflow always does on GMP_RNDN).
 202                In the case GMP_RNDN with |y/x| <= 2^(emin-2), we round
 203                towards zero, as |atan(z)/z| < 1. */
 204             MPFR_ASSERTN (MPFR_PREC_MAX <=
 205                           2 * (mpfr_uexp_t) - MPFR_EMIN_MIN + 5);
 206             if (rnd_mode == GMP_RNDN && MPFR_IS_ZERO (tmp))
 207               rnd_mode = GMP_RNDZ;
 208             sign = MPFR_SIGN (tmp);
 209             mpfr_clear (tmp);
 210             MPFR_SAVE_EXPO_FREE (expo);
 211             return mpfr_underflow (dest, rnd_mode, sign);
 212           }
 213
 214         mpfr_atan (tmp, tmp, GMP_RNDN);   /* Error <= 2*ulp (tmp) since
 215                                              abs(D(arctan)) <= 1 */
 216         /* TODO: check that the error bound is correct in case of overflow. */
 217         /* FIXME: Error <= ulp(tmp) ? */
 218         if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - 2, MPFR_PREC (dest),
 219                                          rnd_mode)))
 220           break;
 221         MPFR_ZIV_NEXT (loop, prec);
 222         mpfr_set_prec (tmp, prec);
 223       }
 224   else /* x < 0 */
 225     /*  Use sign(y)*(PI - atan (|y/x|)) */
 226     {
 227       mpfr_init2 (pi, prec);
 228       for (;;)
 229         {
 230           mpfr_div (tmp, y, x, GMP_RNDN);   /* Error <= ulp (tmp) */
 231           /* If tmp is 0, we have |y/x| <= 2^(-emin-2), thus
 232              atan|y/x| < 2^(-emin-2). */
 233           MPFR_SET_POS (tmp);               /* no error */
 234           mpfr_atan (tmp, tmp, GMP_RNDN);   /* Error <= 2*ulp (tmp) since
 235                                                abs(D(arctan)) <= 1 */
 236           mpfr_const_pi (pi, GMP_RNDN);     /* Error <= ulp(pi) /2 */
 237           e = MPFR_NOTZERO(tmp) ? MPFR_GET_EXP (tmp) : __gmpfr_emin - 1;
 238           mpfr_sub (tmp, pi, tmp, GMP_RNDN);          /* see above */
 239           if (MPFR_IS_NEG (y))
 240             MPFR_CHANGE_SIGN (tmp);
 241           /* Error(tmp) <= (1/2+2^(EXP(pi)-EXP(tmp)-1)+2^(e-EXP(tmp)+1))*ulp
 242                         <= 2^(MAX (MAX (EXP(PI)-EXP(tmp)-1, e-EXP(tmp)+1),
 243                                         -1)+2)*ulp(tmp) */
 244           e = MAX (MAX (MPFR_GET_EXP (pi)-MPFR_GET_EXP (tmp) - 1,
 245                         e - MPFR_GET_EXP (tmp) + 1), -1) + 2;
 246           if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - e, MPFR_PREC (dest),
 247                                            rnd_mode)))
 248             break;
 249           MPFR_ZIV_NEXT (loop, prec);
 250           mpfr_set_prec (tmp, prec);
 251           mpfr_set_prec (pi, prec);
 252         }
 253       mpfr_clear (pi);
 254     }
 255   inexact = mpfr_set (dest, tmp, rnd_mode);
 256
 257  end:
 258   MPFR_ZIV_FREE (loop);
 259   mpfr_clear (tmp);
 260   MPFR_SAVE_EXPO_FREE (expo);
 261   return mpfr_check_range (dest, inexact, rnd_mode);
 262 }