/*
 * IBM Accurate Mathematical Library
 * written by International Business Machines Corp.
 * Copyright (C) 2001-2024 Free Software Foundation, Inc.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, see <https://www.gnu.org/licenses/>.
 */
/*********************************************************************/
/*  MODULE_NAME: utan.c                                              */
/*                                                                   */
/*  FUNCTIONS: utan                                                  */
/*                                                                   */
/*  FILES NEEDED:dla.h endian.h mydefs.h utan.h                      */
/*               branred.c                                           */
/*               utan.tbl                                            */
/*                                                                   */
/*********************************************************************/

#include <errno.h>
#include <float.h>
#include "endian.h"
#include <dla.h>
#include "mydefs.h"
#include <math.h>
#include <math_private.h>
#include <fenv_private.h>
#include <math-underflow.h>
#include <libm-alias-double.h>
#include <fenv.h>

#ifndef SECTION
# define SECTION
#endif

/* tan with max ULP of ~0.619 based on random sampling.  */
double
SECTION
__tan (double x)
{
#include "utan.h"
#include "utan.tbl"

  int ux, i, n;
  double a, da, a2, b, db, c, dc, fi, gi, pz,
	 s, sy, t, t1, t2, t3, t4, w, x2, xn, y, ya,
         yya, z0, z, z2;
  mynumber num, v;

  double retval;

  int __branred (double, double *, double *);

  SET_RESTORE_ROUND_53BIT (FE_TONEAREST);

  /* x=+-INF, x=NaN */
  num.d = x;
  ux = num.i[HIGH_HALF];
  if ((ux & 0x7ff00000) == 0x7ff00000)
    {
      if ((ux & 0x7fffffff) == 0x7ff00000)
	__set_errno (EDOM);
      retval = x - x;
      goto ret;
    }

  w = (x < 0.0) ? -x : x;

  /* (I) The case abs(x) <= 1.259e-8 */
  if (w <= g1.d)
    {
      math_check_force_underflow_nonneg (w);
      retval = x;
      goto ret;
    }

  /* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
  if (w <= g2.d)
    {
      x2 = x * x;

      t2 = d9.d + x2 * d11.d;
      t2 = d7.d + x2 * t2;
      t2 = d5.d + x2 * t2;
      t2 = d3.d + x2 * t2;
      t2 *= x * x2;

      y = x + t2;
      retval = y;
      /* Max ULP is 0.504.  */
      goto ret;
    }

  /* (III) The case 0.0608 < abs(x) <= 0.787 */
  if (w <= g3.d)
    {
      i = ((int) (mfftnhf.d + 256 * w));
      z = w - xfg[i][0].d;
      z2 = z * z;
      s = (x < 0.0) ? -1 : 1;
      pz = z + z * z2 * (e0.d + z2 * e1.d);
      fi = xfg[i][1].d;
      gi = xfg[i][2].d;
      t2 = pz * (gi + fi) / (gi - pz);
      y = fi + t2;
      retval = (s * y);
      /* Max ULP is 0.60.  */
      goto ret;
    }

  /* (---) The case 0.787 < abs(x) <= 25 */
  if (w <= g4.d)
    {
      /* Range reduction by algorithm i */
      t = (x * hpinv.d + toint.d);
      xn = t - toint.d;
      v.d = t;
      t1 = (x - xn * mp1.d) - xn * mp2.d;
      n = v.i[LOW_HALF] & 0x00000001;
      da = xn * mp3.d;
      a = t1 - da;
      da = (t1 - a) - da;
      if (a < 0.0)
	{
	  ya = -a;
	  yya = -da;
	  sy = -1;
	}
      else
	{
	  ya = a;
	  yya = da;
	  sy = 1;
	}

      /* (VI) The case 0.787 < abs(x) <= 25,    0 < abs(y) <= 0.0608 */
      if (ya <= gy2.d)
	{
	  a2 = a * a;
	  t2 = d9.d + a2 * d11.d;
	  t2 = d7.d + a2 * t2;
	  t2 = d5.d + a2 * t2;
	  t2 = d3.d + a2 * t2;
	  t2 = da + a * a2 * t2;

	  if (n)
	    {
	      /* -cot */
	      EADD (a, t2, b, db);
	      DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
	      y = c + dc;
	      retval = (-y);
	      /* Max ULP is 0.506.  */
	      goto ret;
	    }
	  else
	    {
	      /* tan */
	      y = a + t2;
	      retval = y;
	      /* Max ULP is 0.506.  */
	      goto ret;
	    }
	}

      /* (VII) The case 0.787 < abs(x) <= 25,    0.0608 < abs(y) <= 0.787 */

      i = ((int) (mfftnhf.d + 256 * ya));
      z = (z0 = (ya - xfg[i][0].d)) + yya;
      z2 = z * z;
      pz = z + z * z2 * (e0.d + z2 * e1.d);
      fi = xfg[i][1].d;
      gi = xfg[i][2].d;

      if (n)
	{
	  /* -cot */
	  t2 = pz * (fi + gi) / (fi + pz);
	  y = gi - t2;
	  retval = (-sy * y);
	  /* Max ULP is 0.62.  */
	  goto ret;
	}
      else
	{
	  /* tan */
	  t2 = pz * (gi + fi) / (gi - pz);
	  y = fi + t2;
	  retval = (sy * y);
	  /* Max ULP is 0.62.  */
	  goto ret;
	}
    }

  /* (---) The case 25 < abs(x) <= 1e8 */
  if (w <= g5.d)
    {
      /* Range reduction by algorithm ii */
      t = (x * hpinv.d + toint.d);
      xn = t - toint.d;
      v.d = t;
      t1 = (x - xn * mp1.d) - xn * mp2.d;
      n = v.i[LOW_HALF] & 0x00000001;
      da = xn * pp3.d;
      t = t1 - da;
      da = (t1 - t) - da;
      t1 = xn * pp4.d;
      a = t - t1;
      da = ((t - a) - t1) + da;
      EADD (a, da, t1, t2);
      a = t1;
      da = t2;
      if (a < 0.0)
	{
	  ya = -a;
	  yya = -da;
	  sy = -1;
	}
      else
	{
	  ya = a;
	  yya = da;
	  sy = 1;
	}

      /* (VIII) The case 25 < abs(x) <= 1e8,    0 < abs(y) <= 0.0608 */
      if (ya <= gy2.d)
	{
	  a2 = a * a;
	  t2 = d9.d + a2 * d11.d;
	  t2 = d7.d + a2 * t2;
	  t2 = d5.d + a2 * t2;
	  t2 = d3.d + a2 * t2;
	  t2 = da + a * a2 * t2;

	  if (n)
	    {
	      /* -cot */
	      EADD (a, t2, b, db);
	      DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
	      y = c + dc;
	      retval = (-y);
	      /* Max ULP is 0.506.  */
	      goto ret;
	    }
	  else
	    {
	      /* tan */
	      y = a + t2;
	      retval = y;
	      /* Max ULP is 0.506.  */
	      goto ret;
	    }
	}

      /* (IX) The case 25 < abs(x) <= 1e8,    0.0608 < abs(y) <= 0.787 */
      i = ((int) (mfftnhf.d + 256 * ya));
      z = (z0 = (ya - xfg[i][0].d)) + yya;
      z2 = z * z;
      pz = z + z * z2 * (e0.d + z2 * e1.d);
      fi = xfg[i][1].d;
      gi = xfg[i][2].d;

      if (n)
	{
	  /* -cot */
	  t2 = pz * (fi + gi) / (fi + pz);
	  y = gi - t2;
	  retval = (-sy * y);
	  /* Max ULP is 0.62.  */
	  goto ret;
	}
      else
	{
	  /* tan */
	  t2 = pz * (gi + fi) / (gi - pz);
	  y = fi + t2;
	  retval = (sy * y);
	  /* Max ULP is 0.62.  */
	  goto ret;
	}
    }

  /* (---) The case 1e8 < abs(x) < 2**1024 */
  /* Range reduction by algorithm iii */
  n = (__branred (x, &a, &da)) & 0x00000001;
  EADD (a, da, t1, t2);
  a = t1;
  da = t2;
  if (a < 0.0)
    {
      ya = -a;
      yya = -da;
      sy = -1;
    }
  else
    {
      ya = a;
      yya = da;
      sy = 1;
    }

  /* (X) The case 1e8 < abs(x) < 2**1024,    0 < abs(y) <= 0.0608 */
  if (ya <= gy2.d)
    {
      a2 = a * a;
      t2 = d9.d + a2 * d11.d;
      t2 = d7.d + a2 * t2;
      t2 = d5.d + a2 * t2;
      t2 = d3.d + a2 * t2;
      t2 = da + a * a2 * t2;
      if (n)
	{
	  /* -cot */
	  EADD (a, t2, b, db);
	  DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
	  y = c + dc;
	  retval = (-y);
	  /* Max ULP is 0.506.  */
	  goto ret;
	}
      else
	{
	  /* tan */
	  y = a + t2;
	  retval = y;
	  /* Max ULP is 0.507.  */
	  goto ret;
	}
    }

  /* (XI) The case 1e8 < abs(x) < 2**1024,    0.0608 < abs(y) <= 0.787 */
  i = ((int) (mfftnhf.d + 256 * ya));
  z = (z0 = (ya - xfg[i][0].d)) + yya;
  z2 = z * z;
  pz = z + z * z2 * (e0.d + z2 * e1.d);
  fi = xfg[i][1].d;
  gi = xfg[i][2].d;

  if (n)
    {
      /* -cot */
      t2 = pz * (fi + gi) / (fi + pz);
      y = gi - t2;
      retval = (-sy * y);
      /* Max ULP is 0.62.  */
      goto ret;
    }
  else
    {
      /* tan */
      t2 = pz * (gi + fi) / (gi - pz);
      y = fi + t2;
      retval = (sy * y);
      /* Max ULP is 0.62.  */
      goto ret;
    }

ret:
  return retval;
}

#ifndef __tan
libm_alias_double (__tan, tan)
#endif