sysdeps/ieee754/dbl-64/s_tan.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2024 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: utan.c                                              */
  21 /*                                                                   */
  22 /*  FUNCTIONS: utan                                                  */
  23 /*                                                                   */
  24 /*  FILES NEEDED:dla.h endian.h mydefs.h utan.h                      */
  25 /*               branred.c                                           */
  26 /*               utan.tbl                                            */
  27 /*                                                                   */
  28 /*********************************************************************/
  29
  30 #include <errno.h>
  31 #include <float.h>
  32 #include "endian.h"
  33 #include <dla.h>
  34 #include "mydefs.h"
  35 #include <math.h>
  36 #include <math_private.h>
  37 #include <fenv_private.h>
  38 #include <math-underflow.h>
  39 #include <libm-alias-double.h>
  40 #include <fenv.h>
  41
  42 #ifndef SECTION
  43 # define SECTION
  44 #endif
  45
  46 /* tan with max ULP of ~0.619 based on random sampling.  */
  47 double
  48 SECTION
  49 __tan (double x)
  50 {
  51 #include "utan.h"
  52 #include "utan.tbl"
  53
  54   int ux, i, n;
  55   double a, da, a2, b, db, c, dc, fi, gi, pz,
  56          s, sy, t, t1, t2, t3, t4, w, x2, xn, y, ya,
  57          yya, z0, z, z2;
  58   mynumber num, v;
  59
  60   double retval;
  61
  62   int __branred (double, double *, double *);
  63
  64   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
  65
  66   /* x=+-INF, x=NaN */
  67   num.d = x;
  68   ux = num.i[HIGH_HALF];
  69   if ((ux & 0x7ff00000) == 0x7ff00000)
  70     {
  71       if ((ux & 0x7fffffff) == 0x7ff00000)
  72         __set_errno (EDOM);
  73       retval = x - x;
  74       goto ret;
  75     }
  76
  77   w = (x < 0.0) ? -x : x;
  78
  79   /* (I) The case abs(x) <= 1.259e-8 */
  80   if (w <= g1.d)
  81     {
  82       math_check_force_underflow_nonneg (w);
  83       retval = x;
  84       goto ret;
  85     }
  86
  87   /* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
  88   if (w <= g2.d)
  89     {
  90       x2 = x * x;
  91
  92       t2 = d9.d + x2 * d11.d;
  93       t2 = d7.d + x2 * t2;
  94       t2 = d5.d + x2 * t2;
  95       t2 = d3.d + x2 * t2;
  96       t2 *= x * x2;
  97
  98       y = x + t2;
  99       retval = y;
 100       /* Max ULP is 0.504.  */
 101       goto ret;
 102     }
 103
 104   /* (III) The case 0.0608 < abs(x) <= 0.787 */
 105   if (w <= g3.d)
 106     {
 107       i = ((int) (mfftnhf.d + 256 * w));
 108       z = w - xfg[i][0].d;
 109       z2 = z * z;
 110       s = (x < 0.0) ? -1 : 1;
 111       pz = z + z * z2 * (e0.d + z2 * e1.d);
 112       fi = xfg[i][1].d;
 113       gi = xfg[i][2].d;
 114       t2 = pz * (gi + fi) / (gi - pz);
 115       y = fi + t2;
 116       retval = (s * y);
 117       /* Max ULP is 0.60.  */
 118       goto ret;
 119     }
 120
 121   /* (---) The case 0.787 < abs(x) <= 25 */
 122   if (w <= g4.d)
 123     {
 124       /* Range reduction by algorithm i */
 125       t = (x * hpinv.d + toint.d);
 126       xn = t - toint.d;
 127       v.d = t;
 128       t1 = (x - xn * mp1.d) - xn * mp2.d;
 129       n = v.i[LOW_HALF] & 0x00000001;
 130       da = xn * mp3.d;
 131       a = t1 - da;
 132       da = (t1 - a) - da;
 133       if (a < 0.0)
 134         {
 135           ya = -a;
 136           yya = -da;
 137           sy = -1;
 138         }
 139       else
 140         {
 141           ya = a;
 142           yya = da;
 143           sy = 1;
 144         }
 145
 146       /* (VI) The case 0.787 < abs(x) <= 25,    0 < abs(y) <= 0.0608 */
 147       if (ya <= gy2.d)
 148         {
 149           a2 = a * a;
 150           t2 = d9.d + a2 * d11.d;
 151           t2 = d7.d + a2 * t2;
 152           t2 = d5.d + a2 * t2;
 153           t2 = d3.d + a2 * t2;
 154           t2 = da + a * a2 * t2;
 155
 156           if (n)
 157             {
 158               /* -cot */
 159               EADD (a, t2, b, db);
 160               DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
 161               y = c + dc;
 162               retval = (-y);
 163               /* Max ULP is 0.506.  */
 164               goto ret;
 165             }
 166           else
 167             {
 168               /* tan */
 169               y = a + t2;
 170               retval = y;
 171               /* Max ULP is 0.506.  */
 172               goto ret;
 173             }
 174         }
 175
 176       /* (VII) The case 0.787 < abs(x) <= 25,    0.0608 < abs(y) <= 0.787 */
 177
 178       i = ((int) (mfftnhf.d + 256 * ya));
 179       z = (z0 = (ya - xfg[i][0].d)) + yya;
 180       z2 = z * z;
 181       pz = z + z * z2 * (e0.d + z2 * e1.d);
 182       fi = xfg[i][1].d;
 183       gi = xfg[i][2].d;
 184
 185       if (n)
 186         {
 187           /* -cot */
 188           t2 = pz * (fi + gi) / (fi + pz);
 189           y = gi - t2;
 190           retval = (-sy * y);
 191           /* Max ULP is 0.62.  */
 192           goto ret;
 193         }
 194       else
 195         {
 196           /* tan */
 197           t2 = pz * (gi + fi) / (gi - pz);
 198           y = fi + t2;
 199           retval = (sy * y);
 200           /* Max ULP is 0.62.  */
 201           goto ret;
 202         }
 203     }
 204
 205   /* (---) The case 25 < abs(x) <= 1e8 */
 206   if (w <= g5.d)
 207     {
 208       /* Range reduction by algorithm ii */
 209       t = (x * hpinv.d + toint.d);
 210       xn = t - toint.d;
 211       v.d = t;
 212       t1 = (x - xn * mp1.d) - xn * mp2.d;
 213       n = v.i[LOW_HALF] & 0x00000001;
 214       da = xn * pp3.d;
 215       t = t1 - da;
 216       da = (t1 - t) - da;
 217       t1 = xn * pp4.d;
 218       a = t - t1;
 219       da = ((t - a) - t1) + da;
 220       EADD (a, da, t1, t2);
 221       a = t1;
 222       da = t2;
 223       if (a < 0.0)
 224         {
 225           ya = -a;
 226           yya = -da;
 227           sy = -1;
 228         }
 229       else
 230         {
 231           ya = a;
 232           yya = da;
 233           sy = 1;
 234         }
 235
 236       /* (VIII) The case 25 < abs(x) <= 1e8,    0 < abs(y) <= 0.0608 */
 237       if (ya <= gy2.d)
 238         {
 239           a2 = a * a;
 240           t2 = d9.d + a2 * d11.d;
 241           t2 = d7.d + a2 * t2;
 242           t2 = d5.d + a2 * t2;
 243           t2 = d3.d + a2 * t2;
 244           t2 = da + a * a2 * t2;
 245
 246           if (n)
 247             {
 248               /* -cot */
 249               EADD (a, t2, b, db);
 250               DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
 251               y = c + dc;
 252               retval = (-y);
 253               /* Max ULP is 0.506.  */
 254               goto ret;
 255             }
 256           else
 257             {
 258               /* tan */
 259               y = a + t2;
 260               retval = y;
 261               /* Max ULP is 0.506.  */
 262               goto ret;
 263             }
 264         }
 265
 266       /* (IX) The case 25 < abs(x) <= 1e8,    0.0608 < abs(y) <= 0.787 */
 267       i = ((int) (mfftnhf.d + 256 * ya));
 268       z = (z0 = (ya - xfg[i][0].d)) + yya;
 269       z2 = z * z;
 270       pz = z + z * z2 * (e0.d + z2 * e1.d);
 271       fi = xfg[i][1].d;
 272       gi = xfg[i][2].d;
 273
 274       if (n)
 275         {
 276           /* -cot */
 277           t2 = pz * (fi + gi) / (fi + pz);
 278           y = gi - t2;
 279           retval = (-sy * y);
 280           /* Max ULP is 0.62.  */
 281           goto ret;
 282         }
 283       else
 284         {
 285           /* tan */
 286           t2 = pz * (gi + fi) / (gi - pz);
 287           y = fi + t2;
 288           retval = (sy * y);
 289           /* Max ULP is 0.62.  */
 290           goto ret;
 291         }
 292     }
 293
 294   /* (---) The case 1e8 < abs(x) < 2**1024 */
 295   /* Range reduction by algorithm iii */
 296   n = (__branred (x, &a, &da)) & 0x00000001;
 297   EADD (a, da, t1, t2);
 298   a = t1;
 299   da = t2;
 300   if (a < 0.0)
 301     {
 302       ya = -a;
 303       yya = -da;
 304       sy = -1;
 305     }
 306   else
 307     {
 308       ya = a;
 309       yya = da;
 310       sy = 1;
 311     }
 312
 313   /* (X) The case 1e8 < abs(x) < 2**1024,    0 < abs(y) <= 0.0608 */
 314   if (ya <= gy2.d)
 315     {
 316       a2 = a * a;
 317       t2 = d9.d + a2 * d11.d;
 318       t2 = d7.d + a2 * t2;
 319       t2 = d5.d + a2 * t2;
 320       t2 = d3.d + a2 * t2;
 321       t2 = da + a * a2 * t2;
 322       if (n)
 323         {
 324           /* -cot */
 325           EADD (a, t2, b, db);
 326           DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
 327           y = c + dc;
 328           retval = (-y);
 329           /* Max ULP is 0.506.  */
 330           goto ret;
 331         }
 332       else
 333         {
 334           /* tan */
 335           y = a + t2;
 336           retval = y;
 337           /* Max ULP is 0.507.  */
 338           goto ret;
 339         }
 340     }
 341
 342   /* (XI) The case 1e8 < abs(x) < 2**1024,    0.0608 < abs(y) <= 0.787 */
 343   i = ((int) (mfftnhf.d + 256 * ya));
 344   z = (z0 = (ya - xfg[i][0].d)) + yya;
 345   z2 = z * z;
 346   pz = z + z * z2 * (e0.d + z2 * e1.d);
 347   fi = xfg[i][1].d;
 348   gi = xfg[i][2].d;
 349
 350   if (n)
 351     {
 352       /* -cot */
 353       t2 = pz * (fi + gi) / (fi + pz);
 354       y = gi - t2;
 355       retval = (-sy * y);
 356       /* Max ULP is 0.62.  */
 357       goto ret;
 358     }
 359   else
 360     {
 361       /* tan */
 362       t2 = pz * (gi + fi) / (gi - pz);
 363       y = fi + t2;
 364       retval = (sy * y);
 365       /* Max ULP is 0.62.  */
 366       goto ret;
 367     }
 368
 369 ret:
 370   return retval;
 371 }
 372
 373 #ifndef __tan
 374 libm_alias_double (__tan, tan)
 375 #endif