[thirdparty/gcc.git] / libquadmath / math / tanq.c

/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */

/*
  Long double expansions are
  Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
  and are incorporated herein by permission of the author.  The author 
  reserves the right to distribute this material elsewhere under different
  copying permissions.  These modifications are distributed here under 
  the following terms:

    This library 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 library 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 library; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */

/* __quadmath_kernel_tanq( x, y, k )
 * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
 * Input x is assumed to be bounded by ~pi/4 in magnitude.
 * Input y is the tail of x.
 * Input k indicates whether tan (if k=1) or
 * -1/tan (if k= -1) is returned.
 *
 * Algorithm
 *	1. Since tan(-x) = -tan(x), we need only to consider positive x.
 *	2. if x < 2^-57, return x with inexact if x!=0.
 *	3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
 *          on [0,0.67433].
 *
 *	   Note: tan(x+y) = tan(x) + tan'(x)*y
 *		          ~ tan(x) + (1+x*x)*y
 *	   Therefore, for better accuracy in computing tan(x+y), let
 *		r = x^3 * R(x^2)
 *	   then
 *		tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
 *
 *      4. For x in [0.67433,pi/4],  let y = pi/4 - x, then
 *		tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
 *		       = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
 */

#include "quadmath-imp.h"


static const __float128
  one = 1.0Q,
  pio4hi = 7.8539816339744830961566084581987569936977E-1Q,
  pio4lo = 2.1679525325309452561992610065108379921906E-35Q,

  /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
     0 <= x <= 0.6743316650390625
     Peak relative error 8.0e-36  */
 TH =  3.333333333333333333333333333333333333333E-1Q,
 T0 = -1.813014711743583437742363284336855889393E7Q,
 T1 =  1.320767960008972224312740075083259247618E6Q,
 T2 = -2.626775478255838182468651821863299023956E4Q,
 T3 =  1.764573356488504935415411383687150199315E2Q,
 T4 = -3.333267763822178690794678978979803526092E-1Q,

 U0 = -1.359761033807687578306772463253710042010E8Q,
 U1 =  6.494370630656893175666729313065113194784E7Q,
 U2 = -4.180787672237927475505536849168729386782E6Q,
 U3 =  8.031643765106170040139966622980914621521E4Q,
 U4 = -5.323131271912475695157127875560667378597E2Q;
  /* 1.000000000000000000000000000000000000000E0 */


static __float128
__quadmath_kernel_tanq (__float128 x, __float128 y, int iy)
{
  __float128 z, r, v, w, s;
  int32_t ix, sign = 1;
  ieee854_float128 u, u1;

  u.value = x;
  ix = u.words32.w0 & 0x7fffffff;
  if (ix < 0x3fc60000)		/* x < 2**-57 */
    {
      if ((int) x == 0)
	{			/* generate inexact */
	  if ((ix | u.words32.w1 | u.words32.w2 | u.words32.w3
	       | (iy + 1)) == 0)
	    return one / fabsq (x);
	  else
	    return (iy == 1) ? x : -one / x;
	}
    }
  if (ix >= 0x3ffe5942) /* |x| >= 0.6743316650390625 */
    {
      if ((u.words32.w0 & 0x80000000) != 0)
	{
	  x = -x;
	  y = -y;
	  sign = -1;
	}
      else
	sign = 1;
      z = pio4hi - x;
      w = pio4lo - y;
      x = z + w;
      y = 0.0;
    }
  z = x * x;
  r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
  v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
  r = r / v;

  s = z * x;
  r = y + z * (s * r + y);
  r += TH * s;
  w = x + r;
  if (ix >= 0x3ffe5942)
    {
      v = (__float128) iy;
      w = (v - 2.0Q * (x - (w * w / (w + v) - r)));
      if (sign < 0)
	w = -w;
      return w;
    }
  if (iy == 1)
    return w;
  else
    {				/* if allow error up to 2 ulp,
				   simply return -1.0/(x+r) here */
      /*  compute -1.0/(x+r) accurately */
      u1.value = w;
      u1.words32.w2 = 0;
      u1.words32.w3 = 0;
      v = r - (u1.value - x);		/* u1+v = r+x */
      z = -1.0 / w;
      u.value = z;
      u.words32.w2 = 0;
      u.words32.w3 = 0;
      s = 1.0 + u.value * u1.value;
      return u.value + z * (s + u.value * v);
    }
}


\f


/* s_tanl.c -- long double version of s_tan.c.
 * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
 */
  
/* @(#)s_tan.c 5.1 93/09/24 */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */

/* tanl(x)
 * Return tangent function of x.
 *
 * kernel function:
 *	__kernel_tanq		... tangent function on [-pi/4,pi/4]
 *	__ieee754_rem_pio2q	... argument reduction routine
 *
 * Method.
 *      Let S,C and T denote the sin, cos and tan respectively on
 *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
 *	in [-pi/4 , +pi/4], and let n = k mod 4.
 *	We have
 *
 *          n        sin(x)      cos(x)        tan(x)
 *     ----------------------------------------------------------
 *	    0	       S	   C		 T
 *	    1	       C	  -S		-1/T
 *	    2	      -S	  -C		 T
 *	    3	      -C	   S		-1/T
 *     ----------------------------------------------------------
 *
 * Special cases:
 *      Let trig be any of sin, cos, or tan.
 *      trig(+-INF)  is NaN, with signals;
 *      trig(NaN)    is that NaN;
 *
 * Accuracy:
 *	TRIG(x) returns trig(x) nearly rounded
 */


__float128
tanq (__float128 x)
{
	__float128 y[2],z=0.0Q;
	int64_t n, ix;

    /* High word of x. */
	GET_FLT128_MSW64(ix,x);

    /* |x| ~< pi/4 */
	ix &= 0x7fffffffffffffffLL;
	if(ix <= 0x3ffe921fb54442d1LL) return __quadmath_kernel_tanq(x,z,1);

    /* tanl(Inf or NaN) is NaN */
	else if (ix>=0x7fff000000000000LL) {
	    if (ix == 0x7fff000000000000LL) {
		GET_FLT128_LSW64(n,x);
	    }
	    return x-x;		/* NaN */
	}

    /* argument reduction needed */
	else {
	    n = __quadmath_rem_pio2q(x,y);
					/*   1 -- n even, -1 -- n odd */
	    return __quadmath_kernel_tanq(y[0],y[1],1-((n&1)<<1));
	}
}
Commit	Line	Data
87969c8c	1	/*
	2	* ====================================================
	3	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
	4	*
	5	* Developed at SunPro, a Sun Microsystems, Inc. business.
	6	* Permission to use, copy, modify, and distribute this
	7	* software is freely granted, provided that this notice
	8	* is preserved.
	9	* ====================================================
	10	*/
	11
	12	/*
	13	Long double expansions are
	14	Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
	15	and are incorporated herein by permission of the author. The author
	16	reserves the right to distribute this material elsewhere under different
	17	copying permissions. These modifications are distributed here under
	18	the following terms:
	19
	20	This library is free software; you can redistribute it and/or
	21	modify it under the terms of the GNU Lesser General Public
	22	License as published by the Free Software Foundation; either
	23	version 2.1 of the License, or (at your option) any later version.
	24
	25	This library is distributed in the hope that it will be useful,
	26	but WITHOUT ANY WARRANTY; without even the implied warranty of
	27	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	28	Lesser General Public License for more details.
	29
	30	You should have received a copy of the GNU Lesser General Public
	31	License along with this library; if not, write to the Free Software
	32	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
	33
0aa903b3	34	/* __quadmath_kernel_tanq( x, y, k )
87969c8c	35	* kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
	36	* Input x is assumed to be bounded by ~pi/4 in magnitude.
	37	* Input y is the tail of x.
	38	* Input k indicates whether tan (if k=1) or
	39	* -1/tan (if k= -1) is returned.
	40	*
	41	* Algorithm
	42	* 1. Since tan(-x) = -tan(x), we need only to consider positive x.
	43	* 2. if x < 2^-57, return x with inexact if x!=0.
	44	* 3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
	45	* on [0,0.67433].
	46	*
	47	* Note: tan(x+y) = tan(x) + tan'(x)*y
	48	* ~ tan(x) + (1+xx)y
	49	* Therefore, for better accuracy in computing tan(x+y), let
	50	* r = x^3 * R(x^2)
	51	* then
	52	* tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
	53	*
	54	* 4. For x in [0.67433,pi/4], let y = pi/4 - x, then
	55	* tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
	56	* = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
	57	*/
	58
	59	#include "quadmath-imp.h"
	60
	61
	62
	63	static const __float128
	64	one = 1.0Q,
	65	pio4hi = 7.8539816339744830961566084581987569936977E-1Q,
	66	pio4lo = 2.1679525325309452561992610065108379921906E-35Q,
	67
	68	/* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
	69	0 <= x <= 0.6743316650390625
	70	Peak relative error 8.0e-36 */
	71	TH = 3.333333333333333333333333333333333333333E-1Q,
	72	T0 = -1.813014711743583437742363284336855889393E7Q,
	73	T1 = 1.320767960008972224312740075083259247618E6Q,
	74	T2 = -2.626775478255838182468651821863299023956E4Q,
	75	T3 = 1.764573356488504935415411383687150199315E2Q,
	76	T4 = -3.333267763822178690794678978979803526092E-1Q,
	77
	78	U0 = -1.359761033807687578306772463253710042010E8Q,
	79	U1 = 6.494370630656893175666729313065113194784E7Q,
	80	U2 = -4.180787672237927475505536849168729386782E6Q,
	81	U3 = 8.031643765106170040139966622980914621521E4Q,
	82	U4 = -5.323131271912475695157127875560667378597E2Q;
	83	/* 1.000000000000000000000000000000000000000E0 */
	84
	85
	86	static __float128
0aa903b3	87	__quadmath_kernel_tanq (__float128 x, __float128 y, int iy)
87969c8c	88	{
	89	__float128 z, r, v, w, s;
	90	int32_t ix, sign = 1;
	91	ieee854_float128 u, u1;
	92
	93	u.value = x;
	94	ix = u.words32.w0 & 0x7fffffff;
	95	if (ix < 0x3fc60000) /* x < 2*-57 /
	96	{
	97	if ((int) x == 0)
	98	{ /* generate inexact */
	99	if ((ix \| u.words32.w1 \| u.words32.w2 \| u.words32.w3
	100	\| (iy + 1)) == 0)
	101	return one / fabsq (x);
	102	else
	103	return (iy == 1) ? x : -one / x;
	104	}
	105	}
	106	if (ix >= 0x3ffe5942) /* \|x\| >= 0.6743316650390625 */
	107	{
	108	if ((u.words32.w0 & 0x80000000) != 0)
	109	{
	110	x = -x;
	111	y = -y;
	112	sign = -1;
	113	}
	114	else
	115	sign = 1;
	116	z = pio4hi - x;
	117	w = pio4lo - y;
	118	x = z + w;
	119	y = 0.0;
	120	}
	121	z = x * x;
	122	r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
	123	v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
	124	r = r / v;
	125
	126	s = z * x;
	127	r = y + z * (s * r + y);
	128	r += TH * s;
	129	w = x + r;
	130	if (ix >= 0x3ffe5942)
	131	{
	132	v = (__float128) iy;
	133	w = (v - 2.0Q * (x - (w * w / (w + v) - r)));
	134	if (sign < 0)
	135	w = -w;
	136	return w;
	137	}
	138	if (iy == 1)
	139	return w;
	140	else
	141	{ /* if allow error up to 2 ulp,
	142	simply return -1.0/(x+r) here */
	143	/* compute -1.0/(x+r) accurately */
	144	u1.value = w;
	145	u1.words32.w2 = 0;
	146	u1.words32.w3 = 0;
	147	v = r - (u1.value - x); /* u1+v = r+x */
	148	z = -1.0 / w;
	149	u.value = z;
	150	u.words32.w2 = 0;
	151	u.words32.w3 = 0;
152	s = 1.0 + u.value * u1.value;
153	return u.value + z * (s + u.value * v);
154	}
155	}
156
157
158
159	\f
160
161
162
163	/* s_tanl.c -- long double version of s_tan.c.
164	* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
165	*/
166
167	/* @(#)s_tan.c 5.1 93/09/24 */
168	/*
169	* ====================================================
170	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
171	*
172	* Developed at SunPro, a Sun Microsystems, Inc. business.
173	* Permission to use, copy, modify, and distribute this
174	* software is freely granted, provided that this notice
175	* is preserved.
176	* ====================================================
177	*/
178
179	/* tanl(x)
180	* Return tangent function of x.
181	*
182	* kernel function:
183	* __kernel_tanq ... tangent function on [-pi/4,pi/4]
184	* __ieee754_rem_pio2q ... argument reduction routine
185	*
186	* Method.
187	* Let S,C and T denote the sin, cos and tan respectively on
188	* [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
189	* in [-pi/4 , +pi/4], and let n = k mod 4.
190	* We have
191	*
192	* n sin(x) cos(x) tan(x)
193	* ----------------------------------------------------------
194	* 0 S C T
195	* 1 C -S -1/T
196	* 2 -S -C T
197	* 3 -C S -1/T
198	* ----------------------------------------------------------
199	*
200	* Special cases:
201	* Let trig be any of sin, cos, or tan.
202	* trig(+-INF) is NaN, with signals;
203	* trig(NaN) is that NaN;
204	*
205	* Accuracy:
206	* TRIG(x) returns trig(x) nearly rounded
207	*/
208
209
210	__float128
211	tanq (__float128 x)
212	{
213	__float128 y[2],z=0.0Q;
214	int64_t n, ix;
215
216	/* High word of x. */
217	GET_FLT128_MSW64(ix,x);
218
219	/* \|x\| ~< pi/4 */
220	ix &= 0x7fffffffffffffffLL;
0aa903b3	221	if(ix <= 0x3ffe921fb54442d1LL) return __quadmath_kernel_tanq(x,z,1);
87969c8c	222
	223	/* tanl(Inf or NaN) is NaN */
	224	else if (ix>=0x7fff000000000000LL) {
	225	if (ix == 0x7fff000000000000LL) {
	226	GET_FLT128_LSW64(n,x);
	227	}
	228	return x-x; /* NaN */
	229	}
	230
	231	/* argument reduction needed */
	232	else {
0aa903b3	233	n = __quadmath_rem_pio2q(x,y);
	234	/* 1 -- n even, -1 -- n odd */
	235	return __quadmath_kernel_tanq(y[0],y[1],1-((n&1)<<1));
87969c8c	236	}
87969c8c	237	}