[glibc.git] / sysdeps / ieee754 / dbl-64 / s_tan.c

/*
 * 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
Commit	Line	Data
	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