sysdeps/ieee754/dbl-64/s_atan.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2023 Free Software Foundation, Inc.
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public License
  17  * along with this program; if not, see <https://www.gnu.org/licenses/>.
  18  */
  19 /************************************************************************/
  20 /*  MODULE_NAME: atnat.c                                                */
  21 /*                                                                      */
  22 /*  FUNCTIONS:  uatan                                                   */
  23 /*              signArctan                                              */
  24 /*                                                                      */
  25 /*  FILES NEEDED: dla.h endian.h mydefs.h atnat.h                       */
  26 /*                uatan.tbl                                             */
  27 /*                                                                      */
  28 /************************************************************************/
  29
  30 #include <dla.h>
  31 #include "mydefs.h"
  32 #include "uatan.tbl"
  33 #include "atnat.h"
  34 #include <fenv.h>
  35 #include <float.h>
  36 #include <libm-alias-double.h>
  37 #include <math.h>
  38 #include <fenv_private.h>
  39 #include <math-underflow.h>
  40
  41 #define  TWO52     0x1.0p52
  42
  43   /* Fix the sign of y and return */
  44 static double
  45 __signArctan (double x, double y)
  46 {
  47   return copysign (y, x);
  48 }
  49
  50 /* atan with max ULP of ~0.523 based on random sampling.  */
  51 double
  52 __atan (double x)
  53 {
  54   double cor, t1, t2, t3, u,
  55          v, w, ww, y, yy, z;
  56   int i, ux, dx;
  57   mynumber num;
  58
  59   num.d = x;
  60   ux = num.i[HIGH_HALF];
  61   dx = num.i[LOW_HALF];
  62
  63   /* x=NaN */
  64   if (((ux & 0x7ff00000) == 0x7ff00000)
  65       && (((ux & 0x000fffff) | dx) != 0x00000000))
  66     return x + x;
  67
  68   /* Regular values of x, including denormals +-0 and +-INF */
  69   SET_RESTORE_ROUND (FE_TONEAREST);
  70   u = (x < 0) ? -x : x;
  71   if (u < C)
  72     {
  73       if (u < B)
  74         {
  75           if (u < A)
  76             {
  77               math_check_force_underflow_nonneg (u);
  78               return x;
  79             }
  80           else
  81             {                   /* A <= u < B */
  82               v = x * x;
  83               yy = d11.d + v * d13.d;
  84               yy = d9.d + v * yy;
  85               yy = d7.d + v * yy;
  86               yy = d5.d + v * yy;
  87               yy = d3.d + v * yy;
  88               yy *= x * v;
  89
  90               y = x + yy;
  91               /* Max ULP is 0.511.  */
  92               return y;
  93             }
  94         }
  95       else
  96         {                       /* B <= u < C */
  97           i = (TWO52 + 256 * u) - TWO52;
  98           i -= 16;
  99           z = u - cij[i][0].d;
 100           yy = cij[i][5].d + z * cij[i][6].d;
 101           yy = cij[i][4].d + z * yy;
 102           yy = cij[i][3].d + z * yy;
 103           yy = cij[i][2].d + z * yy;
 104           yy *= z;
 105
 106           t1 = cij[i][1].d;
 107           y = t1 + yy;
 108           /* Max ULP is 0.56.  */
 109           return __signArctan (x, y);
 110         }
 111     }
 112   else
 113     {
 114       if (u < D)
 115         {                       /* C <= u < D */
 116           w = 1 / u;
 117           EMULV (w, u, t1, t2);
 118           ww = w * ((1 - t1) - t2);
 119           i = (TWO52 + 256 * w) - TWO52;
 120           i -= 16;
 121           z = (w - cij[i][0].d) + ww;
 122
 123           yy = cij[i][5].d + z * cij[i][6].d;
 124           yy = cij[i][4].d + z * yy;
 125           yy = cij[i][3].d + z * yy;
 126           yy = cij[i][2].d + z * yy;
 127           yy = HPI1 - z * yy;
 128
 129           t1 = HPI - cij[i][1].d;
 130           y = t1 + yy;
 131           /* Max ULP is 0.503.  */
 132           return __signArctan (x, y);
 133         }
 134       else
 135         {
 136           if (u < E)
 137             {                   /* D <= u < E */
 138               w = 1 / u;
 139               v = w * w;
 140               EMULV (w, u, t1, t2);
 141
 142               yy = d11.d + v * d13.d;
 143               yy = d9.d + v * yy;
 144               yy = d7.d + v * yy;
 145               yy = d5.d + v * yy;
 146               yy = d3.d + v * yy;
 147               yy *= w * v;
 148
 149               ww = w * ((1 - t1) - t2);
 150               ESUB (HPI, w, t3, cor);
 151               yy = ((HPI1 + cor) - ww) - yy;
 152               y = t3 + yy;
 153               /* Max ULP is 0.5003.  */
 154               return __signArctan (x, y);
 155             }
 156           else
 157             {
 158               /* u >= E */
 159               if (x > 0)
 160                 return HPI;
 161               else
 162                 return MHPI;
 163             }
 164         }
 165     }
 166 }
 167
 168 #ifndef __atan
 169 libm_alias_double (__atan, atan)
 170 #endif