contrib/mpfr/src/atan2.c

   1 /* mpfr_atan2 -- arc-tan 2 of a floating-point number
   2
   3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   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 3 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.LESSER.  If not, see
  20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 static int
  27 pi_div_2ui (mpfr_ptr dest, int i, int neg, mpfr_rnd_t rnd_mode)
  28 {
  29   int inexact;
  30   MPFR_SAVE_EXPO_DECL (expo);
  31
  32   MPFR_SAVE_EXPO_MARK (expo);
  33   if (neg) /* -PI/2^i */
  34     {
  35       inexact = - mpfr_const_pi (dest, MPFR_INVERT_RND (rnd_mode));
  36       MPFR_CHANGE_SIGN (dest);
  37     }
  38   else /* PI/2^i */
  39     {
  40       inexact = mpfr_const_pi (dest, rnd_mode);
  41     }
  42   mpfr_div_2ui (dest, dest, i, rnd_mode);  /* exact */
  43   MPFR_SAVE_EXPO_FREE (expo);
  44   return mpfr_check_range (dest, inexact, rnd_mode);
  45 }
  46
  47 int
  48 mpfr_atan2 (mpfr_ptr dest, mpfr_srcptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
  49 {
  50   mpfr_t tmp, pi;
  51   int inexact;
  52   mpfr_prec_t prec;
  53   mpfr_exp_t e;
  54   MPFR_SAVE_EXPO_DECL (expo);
  55   MPFR_ZIV_DECL (loop);
  56
  57   MPFR_LOG_FUNC
  58     (("y[%Pu]=%.*Rg x[%Pu]=%.*Rg rnd=%d",
  59       mpfr_get_prec (y), mpfr_log_prec, y,
  60       mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
  61      ("atan[%Pu]=%.*Rg inexact=%d",
  62       mpfr_get_prec (dest), mpfr_log_prec, dest, inexact));
  63
  64   /* Special cases */
  65   if (MPFR_ARE_SINGULAR (x, y))
  66     {
  67       /* atan2(0, 0) does not raise the "invalid" floating-point
  68          exception, nor does atan2(y, 0) raise the "divide-by-zero"
  69          floating-point exception.
  70          -- atan2(±0, -0) returns ±pi.313)
  71          -- atan2(±0, +0) returns ±0.
  72          -- atan2(±0, x) returns ±pi, for x < 0.
  73          -- atan2(±0, x) returns ±0, for x > 0.
  74          -- atan2(y, ±0) returns -pi/2 for y < 0.
  75          -- atan2(y, ±0) returns pi/2 for y > 0.
  76          -- atan2(±oo, -oo) returns ±3pi/4.
  77          -- atan2(±oo, +oo) returns ±pi/4.
  78          -- atan2(±oo, x) returns ±pi/2, for finite x.
  79          -- atan2(±y, -oo) returns ±pi, for finite y > 0.
  80          -- atan2(±y, +oo) returns ±0, for finite y > 0.
  81       */
  82       if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y))
  83         {
  84           MPFR_SET_NAN (dest);
  85           MPFR_RET_NAN;
  86         }
  87       if (MPFR_IS_ZERO (y))
  88         {
  89           if (MPFR_IS_NEG (x)) /* +/- PI */
  90             {
  91             set_pi:
  92               if (MPFR_IS_NEG (y))
  93                 {
  94                   inexact =  mpfr_const_pi (dest, MPFR_INVERT_RND (rnd_mode));
  95                   MPFR_CHANGE_SIGN (dest);
  96                   return -inexact;
  97                 }
  98               else
  99                 return mpfr_const_pi (dest, rnd_mode);
 100             }
 101           else /* +/- 0 */
 102             {
 103             set_zero:
 104               MPFR_SET_ZERO (dest);
 105               MPFR_SET_SAME_SIGN (dest, y);
 106               return 0;
 107             }
 108         }
 109       if (MPFR_IS_ZERO (x))
 110         {
 111           return pi_div_2ui (dest, 1, MPFR_IS_NEG (y), rnd_mode);
 112         }
 113       if (MPFR_IS_INF (y))
 114         {
 115           if (!MPFR_IS_INF (x)) /* +/- PI/2 */
 116             return pi_div_2ui (dest, 1, MPFR_IS_NEG (y), rnd_mode);
 117           else if (MPFR_IS_POS (x)) /* +/- PI/4 */
 118             return pi_div_2ui (dest, 2, MPFR_IS_NEG (y), rnd_mode);
 119           else /* +/- 3*PI/4: Ugly since we have to round properly */
 120             {
 121               mpfr_t tmp2;
 122               MPFR_ZIV_DECL (loop2);
 123               mpfr_prec_t prec2 = MPFR_PREC (dest) + 10;
 124
 125               MPFR_SAVE_EXPO_MARK (expo);
 126               mpfr_init2 (tmp2, prec2);
 127               MPFR_ZIV_INIT (loop2, prec2);
 128               for (;;)
 129                 {
 130                   mpfr_const_pi (tmp2, MPFR_RNDN);
 131                   mpfr_mul_ui (tmp2, tmp2, 3, MPFR_RNDN); /* Error <= 2  */
 132                   mpfr_div_2ui (tmp2, tmp2, 2, MPFR_RNDN);
 133                   if (mpfr_round_p (MPFR_MANT (tmp2), MPFR_LIMB_SIZE (tmp2),
 134                                     MPFR_PREC (tmp2) - 2,
 135                                     MPFR_PREC (dest) + (rnd_mode == MPFR_RNDN)))
 136                     break;
 137                   MPFR_ZIV_NEXT (loop2, prec2);
 138                   mpfr_set_prec (tmp2, prec2);
 139                 }
 140               MPFR_ZIV_FREE (loop2);
 141               if (MPFR_IS_NEG (y))
 142                 MPFR_CHANGE_SIGN (tmp2);
 143               inexact = mpfr_set (dest, tmp2, rnd_mode);
 144               mpfr_clear (tmp2);
 145               MPFR_SAVE_EXPO_FREE (expo);
 146               return mpfr_check_range (dest, inexact, rnd_mode);
 147             }
 148         }
 149       MPFR_ASSERTD (MPFR_IS_INF (x));
 150       if (MPFR_IS_NEG (x))
 151         goto set_pi;
 152       else
 153         goto set_zero;
 154     }
 155
 156   /* When x is a power of two, we call directly atan(y/x) since y/x is
 157      exact. */
 158   if (MPFR_UNLIKELY (MPFR_IS_POWER_OF_2 (x)))
 159     {
 160       int r;
 161       mpfr_t yoverx;
 162       unsigned int saved_flags = __gmpfr_flags;
 163
 164       mpfr_init2 (yoverx, MPFR_PREC (y));
 165       if (MPFR_LIKELY (mpfr_div_2si (yoverx, y, MPFR_GET_EXP (x) - 1,
 166                                      MPFR_RNDN) == 0))
 167         {
 168           /* Here the flags have not changed due to mpfr_div_2si. */
 169           r = mpfr_atan (dest, yoverx, rnd_mode);
 170           mpfr_clear (yoverx);
 171           return r;
 172         }
 173       else
 174         {
 175           /* Division is inexact because of a small exponent range */
 176           mpfr_clear (yoverx);
 177           __gmpfr_flags = saved_flags;
 178         }
 179     }
 180
 181   MPFR_SAVE_EXPO_MARK (expo);
 182
 183   /* Set up initial prec */
 184   prec = MPFR_PREC (dest) + 3 + MPFR_INT_CEIL_LOG2 (MPFR_PREC (dest));
 185   mpfr_init2 (tmp, prec);
 186
 187   MPFR_ZIV_INIT (loop, prec);
 188   if (MPFR_IS_POS (x))
 189     /* use atan2(y,x) = atan(y/x) */
 190     for (;;)
 191       {
 192         int div_inex;
 193         MPFR_BLOCK_DECL (flags);
 194
 195         MPFR_BLOCK (flags, div_inex = mpfr_div (tmp, y, x, MPFR_RNDN));
 196         if (div_inex == 0)
 197           {
 198             /* Result is exact. */
 199             inexact = mpfr_atan (dest, tmp, rnd_mode);
 200             goto end;
 201           }
 202
 203         /* Error <= ulp (tmp) except in case of underflow or overflow. */
 204
 205         /* If the division underflowed, since |atan(z)/z| < 1, we have
 206            an underflow. */
 207         if (MPFR_UNDERFLOW (flags))
 208           {
 209             int sign;
 210
 211             /* In the case MPFR_RNDN with 2^(emin-2) < |y/x| < 2^(emin-1):
 212                The smallest significand value S > 1 of |y/x| is:
 213                  * 1 / (1 - 2^(-px))                        if py <= px,
 214                  * (1 - 2^(-px) + 2^(-py)) / (1 - 2^(-px))  if py >= px.
 215                Therefore S - 1 > 2^(-pz), where pz = max(px,py). We have:
 216                atan(|y/x|) > atan(z), where z = 2^(emin-2) * (1 + 2^(-pz)).
 217                            > z - z^3 / 3.
 218                            > 2^(emin-2) * (1 + 2^(-pz) - 2^(2 emin - 5))
 219                Assuming pz <= -2 emin + 5, we can round away from zero
 220                (this is what mpfr_underflow always does on MPFR_RNDN).
 221                In the case MPFR_RNDN with |y/x| <= 2^(emin-2), we round
 222                toward zero, as |atan(z)/z| < 1. */
 223             MPFR_ASSERTN (MPFR_PREC_MAX <=
 224                           2 * (mpfr_uexp_t) - MPFR_EMIN_MIN + 5);
 225             if (rnd_mode == MPFR_RNDN && MPFR_IS_ZERO (tmp))
 226               rnd_mode = MPFR_RNDZ;
 227             sign = MPFR_SIGN (tmp);
 228             mpfr_clear (tmp);
 229             MPFR_SAVE_EXPO_FREE (expo);
 230             return mpfr_underflow (dest, rnd_mode, sign);
 231           }
 232
 233         mpfr_atan (tmp, tmp, MPFR_RNDN);   /* Error <= 2*ulp (tmp) since
 234                                              abs(D(arctan)) <= 1 */
 235         /* TODO: check that the error bound is correct in case of overflow. */
 236         /* FIXME: Error <= ulp(tmp) ? */
 237         if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - 2, MPFR_PREC (dest),
 238                                          rnd_mode)))
 239           break;
 240         MPFR_ZIV_NEXT (loop, prec);
 241         mpfr_set_prec (tmp, prec);
 242       }
 243   else /* x < 0 */
 244     /*  Use sign(y)*(PI - atan (|y/x|)) */
 245     {
 246       mpfr_init2 (pi, prec);
 247       for (;;)
 248         {
 249           mpfr_div (tmp, y, x, MPFR_RNDN);   /* Error <= ulp (tmp) */
 250           /* If tmp is 0, we have |y/x| <= 2^(-emin-2), thus
 251              atan|y/x| < 2^(-emin-2). */
 252           MPFR_SET_POS (tmp);               /* no error */
 253           mpfr_atan (tmp, tmp, MPFR_RNDN);   /* Error <= 2*ulp (tmp) since
 254                                                abs(D(arctan)) <= 1 */
 255           mpfr_const_pi (pi, MPFR_RNDN);     /* Error <= ulp(pi) /2 */
 256           e = MPFR_NOTZERO(tmp) ? MPFR_GET_EXP (tmp) : __gmpfr_emin - 1;
 257           mpfr_sub (tmp, pi, tmp, MPFR_RNDN);          /* see above */
 258           if (MPFR_IS_NEG (y))
 259             MPFR_CHANGE_SIGN (tmp);
 260           /* Error(tmp) <= (1/2+2^(EXP(pi)-EXP(tmp)-1)+2^(e-EXP(tmp)+1))*ulp
 261                         <= 2^(MAX (MAX (EXP(PI)-EXP(tmp)-1, e-EXP(tmp)+1),
 262                                         -1)+2)*ulp(tmp) */
 263           e = MAX (MAX (MPFR_GET_EXP (pi)-MPFR_GET_EXP (tmp) - 1,
 264                         e - MPFR_GET_EXP (tmp) + 1), -1) + 2;
 265           if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - e, MPFR_PREC (dest),
 266                                            rnd_mode)))
 267             break;
 268           MPFR_ZIV_NEXT (loop, prec);
 269           mpfr_set_prec (tmp, prec);
 270           mpfr_set_prec (pi, prec);
 271         }
 272       mpfr_clear (pi);
 273     }
 274   inexact = mpfr_set (dest, tmp, rnd_mode);
 275
 276  end:
 277   MPFR_ZIV_FREE (loop);
 278   mpfr_clear (tmp);
 279   MPFR_SAVE_EXPO_FREE (expo);
 280   return mpfr_check_range (dest, inexact, rnd_mode);
 281 }