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

/*							s_atanl.c
 *
 *	Inverse circular tangent for 128-bit long double precision
 *      (arctangent)
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, atanq();
 *
 * y = atanq( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
 *
 * The function uses a rational approximation of the form
 * t + t^3 P(t^2)/Q(t^2), optimized for |t| < 0.09375.
 *
 * The argument is reduced using the identity
 *    arctan x - arctan u  =  arctan ((x-u)/(1 + ux))
 * and an 83-entry lookup table for arctan u, with u = 0, 1/8, ..., 10.25.
 * Use of the table improves the execution speed of the routine.
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -19, 19       4e5       1.7e-34     5.4e-35
 *
 *
 * WARNING:
 *
 * This program uses integer operations on bit fields of floating-point
 * numbers.  It does not work with data structures other than the
 * structure assumed.
 *
 */

/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>

    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, see
    <http://www.gnu.org/licenses/>.  */

#include "quadmath-imp.h"

/* arctan(k/8), k = 0, ..., 82 */
static const __float128 atantbl[84] = {
  0.0000000000000000000000000000000000000000E0Q,
  1.2435499454676143503135484916387102557317E-1Q, /* arctan(0.125)  */
  2.4497866312686415417208248121127581091414E-1Q,
  3.5877067027057222039592006392646049977698E-1Q,
  4.6364760900080611621425623146121440202854E-1Q,
  5.5859931534356243597150821640166127034645E-1Q,
  6.4350110879328438680280922871732263804151E-1Q,
  7.1882999962162450541701415152590465395142E-1Q,
  7.8539816339744830961566084581987572104929E-1Q,
  8.4415398611317100251784414827164750652594E-1Q,
  8.9605538457134395617480071802993782702458E-1Q,
  9.4200004037946366473793717053459358607166E-1Q,
  9.8279372324732906798571061101466601449688E-1Q,
  1.0191413442663497346383429170230636487744E0Q,
  1.0516502125483736674598673120862998296302E0Q,
  1.0808390005411683108871567292171998202703E0Q,
  1.1071487177940905030170654601785370400700E0Q,
  1.1309537439791604464709335155363278047493E0Q,
  1.1525719972156675180401498626127513797495E0Q,
  1.1722738811284763866005949441337046149712E0Q,
  1.1902899496825317329277337748293183376012E0Q,
  1.2068173702852525303955115800565576303133E0Q,
  1.2220253232109896370417417439225704908830E0Q,
  1.2360594894780819419094519711090786987027E0Q,
  1.2490457723982544258299170772810901230778E0Q,
  1.2610933822524404193139408812473357720101E0Q,
  1.2722973952087173412961937498224804940684E0Q,
  1.2827408797442707473628852511364955306249E0Q,
  1.2924966677897852679030914214070816845853E0Q,
  1.3016288340091961438047858503666855921414E0Q,
  1.3101939350475556342564376891719053122733E0Q,
  1.3182420510168370498593302023271362531155E0Q,
  1.3258176636680324650592392104284756311844E0Q,
  1.3329603993374458675538498697331558093700E0Q,
  1.3397056595989995393283037525895557411039E0Q,
  1.3460851583802539310489409282517796256512E0Q,
  1.3521273809209546571891479413898128509842E0Q,
  1.3578579772154994751124898859640585287459E0Q,
  1.3633001003596939542892985278250991189943E0Q,
  1.3684746984165928776366381936948529556191E0Q,
  1.3734007669450158608612719264449611486510E0Q,
  1.3780955681325110444536609641291551522494E0Q,
  1.3825748214901258580599674177685685125566E0Q,
  1.3868528702577214543289381097042486034883E0Q,
  1.3909428270024183486427686943836432060856E0Q,
  1.3948567013423687823948122092044222644895E0Q,
  1.3986055122719575950126700816114282335732E0Q,
  1.4021993871854670105330304794336492676944E0Q,
  1.4056476493802697809521934019958079881002E0Q,
  1.4089588955564736949699075250792569287156E0Q,
  1.4121410646084952153676136718584891599630E0Q,
  1.4152014988178669079462550975833894394929E0Q,
  1.4181469983996314594038603039700989523716E0Q,
  1.4209838702219992566633046424614466661176E0Q,
  1.4237179714064941189018190466107297503086E0Q,
  1.4263547484202526397918060597281265695725E0Q,
  1.4288992721907326964184700745371983590908E0Q,
  1.4313562697035588982240194668401779312122E0Q,
  1.4337301524847089866404719096698873648610E0Q,
  1.4360250423171655234964275337155008780675E0Q,
  1.4382447944982225979614042479354815855386E0Q,
  1.4403930189057632173997301031392126865694E0Q,
  1.4424730991091018200252920599377292525125E0Q,
  1.4444882097316563655148453598508037025938E0Q,
  1.4464413322481351841999668424758804165254E0Q,
  1.4483352693775551917970437843145232637695E0Q,
  1.4501726582147939000905940595923466567576E0Q,
  1.4519559822271314199339700039142990228105E0Q,
  1.4536875822280323362423034480994649820285E0Q,
  1.4553696664279718992423082296859928222270E0Q,
  1.4570043196511885530074841089245667532358E0Q,
  1.4585935117976422128825857356750737658039E0Q,
  1.4601391056210009726721818194296893361233E0Q,
  1.4616428638860188872060496086383008594310E0Q,
  1.4631064559620759326975975316301202111560E0Q,
  1.4645314639038178118428450961503371619177E0Q,
  1.4659193880646627234129855241049975398470E0Q,
  1.4672716522843522691530527207287398276197E0Q,
  1.4685896086876430842559640450619880951144E0Q,
  1.4698745421276027686510391411132998919794E0Q,
  1.4711276743037345918528755717617308518553E0Q,
  1.4723501675822635384916444186631899205983E0Q,
  1.4735431285433308455179928682541563973416E0Q, /* arctan(10.25) */
  1.5707963267948966192313216916397514420986E0Q  /* pi/2 */
};


/* arctan t = t + t^3 p(t^2) / q(t^2)
   |t| <= 0.09375
   peak relative error 5.3e-37 */

static const __float128
  p0 = -4.283708356338736809269381409828726405572E1Q,
  p1 = -8.636132499244548540964557273544599863825E1Q,
  p2 = -5.713554848244551350855604111031839613216E1Q,
  p3 = -1.371405711877433266573835355036413750118E1Q,
  p4 = -8.638214309119210906997318946650189640184E-1Q,
  q0 = 1.285112506901621042780814422948906537959E2Q,
  q1 = 3.361907253914337187957855834229672347089E2Q,
  q2 = 3.180448303864130128268191635189365331680E2Q,
  q3 = 1.307244136980865800160844625025280344686E2Q,
  q4 = 2.173623741810414221251136181221172551416E1Q;
  /* q5 = 1.000000000000000000000000000000000000000E0 */

static const __float128 huge = 1.0e4930Q;

__float128
atanq (__float128 x)
{
  int k, sign;
  __float128 t, u, p, q;
  ieee854_float128 s;

  s.value = x;
  k = s.words32.w0;
  if (k & 0x80000000)
    sign = 1;
  else
    sign = 0;

  /* Check for IEEE special cases.  */
  k &= 0x7fffffff;
  if (k >= 0x7fff0000)
    {
      /* NaN. */
      if ((k & 0xffff) | s.words32.w1 | s.words32.w2 | s.words32.w3)
	return (x + x);

      /* Infinity. */
      if (sign)
	return -atantbl[83];
      else
	return atantbl[83];
    }

  if (k <= 0x3fc50000) /* |x| < 2**-58 */
    {
      math_check_force_underflow (x);
      /* Raise inexact. */
      if (huge + x > 0.0)
	return x;
    }

  if (k >= 0x40720000) /* |x| > 2**115 */
    {
      /* Saturate result to {-,+}pi/2 */
      if (sign)
	return -atantbl[83];
      else
	return atantbl[83];
    }

  if (sign)
      x = -x;

  if (k >= 0x40024800) /* 10.25 */
    {
      k = 83;
      t = -1.0/x;
    }
  else
    {
      /* Index of nearest table element.
	 Roundoff to integer is asymmetrical to avoid cancellation when t < 0
         (cf. fdlibm). */
      k = 8.0 * x + 0.25;
      u = 0.125Q * k;
      /* Small arctan argument.  */
      t = (x - u) / (1.0 + x * u);
    }

  /* Arctan of small argument t.  */
  u = t * t;
  p =     ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
  q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
  u = t * u * p / q  +  t;

  /* arctan x = arctan u  +  arctan t */
  u = atantbl[k] + u;
  if (sign)
    return (-u);
  else
    return u;
}
Commit	Line	Data
87969c8c	1	/* s_atanl.c
87969c8c	2	*
e409716d	3	* Inverse circular tangent for 128-bit long double precision
87969c8c	4	* (arctangent)
	5	*
	6	*
	7	*
	8	* SYNOPSIS:
	9	*
e409716d	10	* long double x, y, atanq();
87969c8c	11	*
e409716d	12	* y = atanq( x );
87969c8c	13	*
	14	*
	15	*
	16	* DESCRIPTION:
	17	*
	18	* Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
	19	*
	20	* The function uses a rational approximation of the form
	21	* t + t^3 P(t^2)/Q(t^2), optimized for \|t\| < 0.09375.
	22	*
	23	* The argument is reduced using the identity
	24	* arctan x - arctan u = arctan ((x-u)/(1 + ux))
	25	* and an 83-entry lookup table for arctan u, with u = 0, 1/8, ..., 10.25.
	26	* Use of the table improves the execution speed of the routine.
	27	*
	28	*
	29	*
	30	* ACCURACY:
	31	*
	32	* Relative error:
	33	* arithmetic domain # trials peak rms
	34	* IEEE -19, 19 4e5 1.7e-34 5.4e-35
	35	*
	36	*
	37	* WARNING:
	38	*
	39	* This program uses integer operations on bit fields of floating-point
	40	* numbers. It does not work with data structures other than the
	41	* structure assumed.
	42	*
	43	*/
	44
c98f0ea6	45	/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
87969c8c	46
	47	This library is free software; you can redistribute it and/or
	48	modify it under the terms of the GNU Lesser General Public
	49	License as published by the Free Software Foundation; either
	50	version 2.1 of the License, or (at your option) any later version.
	51
	52	This library is distributed in the hope that it will be useful,
	53	but WITHOUT ANY WARRANTY; without even the implied warranty of
	54	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	55	Lesser General Public License for more details.
	56
	57	You should have received a copy of the GNU Lesser General Public
e409716d	58	License along with this library; if not, see
e409716d	59	<http://www.gnu.org/licenses/>. */
87969c8c	60
	61	#include "quadmath-imp.h"
	62
	63	/* arctan(k/8), k = 0, ..., 82 */
	64	static const __float128 atantbl[84] = {
	65	0.0000000000000000000000000000000000000000E0Q,
	66	1.2435499454676143503135484916387102557317E-1Q, /* arctan(0.125) */
	67	2.4497866312686415417208248121127581091414E-1Q,
	68	3.5877067027057222039592006392646049977698E-1Q,
	69	4.6364760900080611621425623146121440202854E-1Q,
	70	5.5859931534356243597150821640166127034645E-1Q,
	71	6.4350110879328438680280922871732263804151E-1Q,
	72	7.1882999962162450541701415152590465395142E-1Q,
	73	7.8539816339744830961566084581987572104929E-1Q,
	74	8.4415398611317100251784414827164750652594E-1Q,
	75	8.9605538457134395617480071802993782702458E-1Q,
	76	9.4200004037946366473793717053459358607166E-1Q,
	77	9.8279372324732906798571061101466601449688E-1Q,
	78	1.0191413442663497346383429170230636487744E0Q,
	79	1.0516502125483736674598673120862998296302E0Q,
	80	1.0808390005411683108871567292171998202703E0Q,
	81	1.1071487177940905030170654601785370400700E0Q,
	82	1.1309537439791604464709335155363278047493E0Q,
	83	1.1525719972156675180401498626127513797495E0Q,
	84	1.1722738811284763866005949441337046149712E0Q,
	85	1.1902899496825317329277337748293183376012E0Q,
	86	1.2068173702852525303955115800565576303133E0Q,
	87	1.2220253232109896370417417439225704908830E0Q,
	88	1.2360594894780819419094519711090786987027E0Q,
	89	1.2490457723982544258299170772810901230778E0Q,
	90	1.2610933822524404193139408812473357720101E0Q,
	91	1.2722973952087173412961937498224804940684E0Q,
	92	1.2827408797442707473628852511364955306249E0Q,
	93	1.2924966677897852679030914214070816845853E0Q,
	94	1.3016288340091961438047858503666855921414E0Q,
	95	1.3101939350475556342564376891719053122733E0Q,
	96	1.3182420510168370498593302023271362531155E0Q,
	97	1.3258176636680324650592392104284756311844E0Q,
	98	1.3329603993374458675538498697331558093700E0Q,
	99	1.3397056595989995393283037525895557411039E0Q,
	100	1.3460851583802539310489409282517796256512E0Q,
	101	1.3521273809209546571891479413898128509842E0Q,
	102	1.3578579772154994751124898859640585287459E0Q,
	103	1.3633001003596939542892985278250991189943E0Q,
	104	1.3684746984165928776366381936948529556191E0Q,
	105	1.3734007669450158608612719264449611486510E0Q,
	106	1.3780955681325110444536609641291551522494E0Q,
	107	1.3825748214901258580599674177685685125566E0Q,
	108	1.3868528702577214543289381097042486034883E0Q,
	109	1.3909428270024183486427686943836432060856E0Q,
	110	1.3948567013423687823948122092044222644895E0Q,
	111	1.3986055122719575950126700816114282335732E0Q,
	112	1.4021993871854670105330304794336492676944E0Q,
	113	1.4056476493802697809521934019958079881002E0Q,
	114	1.4089588955564736949699075250792569287156E0Q,
	115	1.4121410646084952153676136718584891599630E0Q,
	116	1.4152014988178669079462550975833894394929E0Q,
	117	1.4181469983996314594038603039700989523716E0Q,
	118	1.4209838702219992566633046424614466661176E0Q,
	119	1.4237179714064941189018190466107297503086E0Q,
	120	1.4263547484202526397918060597281265695725E0Q,
	121	1.4288992721907326964184700745371983590908E0Q,
	122	1.4313562697035588982240194668401779312122E0Q,
	123	1.4337301524847089866404719096698873648610E0Q,
124	1.4360250423171655234964275337155008780675E0Q,
125	1.4382447944982225979614042479354815855386E0Q,
126	1.4403930189057632173997301031392126865694E0Q,
127	1.4424730991091018200252920599377292525125E0Q,
128	1.4444882097316563655148453598508037025938E0Q,
129	1.4464413322481351841999668424758804165254E0Q,
130	1.4483352693775551917970437843145232637695E0Q,
131	1.4501726582147939000905940595923466567576E0Q,
132	1.4519559822271314199339700039142990228105E0Q,
133	1.4536875822280323362423034480994649820285E0Q,
134	1.4553696664279718992423082296859928222270E0Q,
135	1.4570043196511885530074841089245667532358E0Q,
136	1.4585935117976422128825857356750737658039E0Q,
137	1.4601391056210009726721818194296893361233E0Q,
138	1.4616428638860188872060496086383008594310E0Q,
139	1.4631064559620759326975975316301202111560E0Q,
140	1.4645314639038178118428450961503371619177E0Q,
141	1.4659193880646627234129855241049975398470E0Q,
142	1.4672716522843522691530527207287398276197E0Q,
143	1.4685896086876430842559640450619880951144E0Q,
144	1.4698745421276027686510391411132998919794E0Q,
145	1.4711276743037345918528755717617308518553E0Q,
146	1.4723501675822635384916444186631899205983E0Q,
147	1.4735431285433308455179928682541563973416E0Q, /* arctan(10.25) */
148	1.5707963267948966192313216916397514420986E0Q /* pi/2 */
149	};
150
151
152	/* arctan t = t + t^3 p(t^2) / q(t^2)
153	\|t\| <= 0.09375
154	peak relative error 5.3e-37 */
155
156	static const __float128
157	p0 = -4.283708356338736809269381409828726405572E1Q,
158	p1 = -8.636132499244548540964557273544599863825E1Q,
159	p2 = -5.713554848244551350855604111031839613216E1Q,
160	p3 = -1.371405711877433266573835355036413750118E1Q,
161	p4 = -8.638214309119210906997318946650189640184E-1Q,
162	q0 = 1.285112506901621042780814422948906537959E2Q,
163	q1 = 3.361907253914337187957855834229672347089E2Q,
164	q2 = 3.180448303864130128268191635189365331680E2Q,
165	q3 = 1.307244136980865800160844625025280344686E2Q,
166	q4 = 2.173623741810414221251136181221172551416E1Q;
167	/* q5 = 1.000000000000000000000000000000000000000E0 */
168
c98f0ea6	169	static const __float128 huge = 1.0e4930Q;
87969c8c	170
	171	__float128
	172	atanq (__float128 x)
	173	{
	174	int k, sign;
	175	__float128 t, u, p, q;
	176	ieee854_float128 s;
	177
	178	s.value = x;
	179	k = s.words32.w0;
	180	if (k & 0x80000000)
	181	sign = 1;
	182	else
	183	sign = 0;
	184
	185	/* Check for IEEE special cases. */
	186	k &= 0x7fffffff;
	187	if (k >= 0x7fff0000)
	188	{
	189	/* NaN. */
	190	if ((k & 0xffff) \| s.words32.w1 \| s.words32.w2 \| s.words32.w3)
	191	return (x + x);
	192
	193	/* Infinity. */
	194	if (sign)
	195	return -atantbl[83];
	196	else
	197	return atantbl[83];
	198	}
	199
89a213c9	200	if (k <= 0x3fc50000) /* \|x\| < 2*-58 /
89a213c9	201	{
c98f0ea6	202	math_check_force_underflow (x);
89a213c9	203	/* Raise inexact. */
	204	if (huge + x > 0.0)
	205	return x;
	206	}
	207
	208	if (k >= 0x40720000) /* \|x\| > 2*115 /
	209	{
	210	/* Saturate result to {-,+}pi/2 */
	211	if (sign)
	212	return -atantbl[83];
	213	else
	214	return atantbl[83];
	215	}
	216
87969c8c	217	if (sign)
	218	x = -x;
	219
	220	if (k >= 0x40024800) /* 10.25 */
	221	{
	222	k = 83;
	223	t = -1.0/x;
	224	}
	225	else
	226	{
	227	/* Index of nearest table element.
	228	Roundoff to integer is asymmetrical to avoid cancellation when t < 0
	229	(cf. fdlibm). */
e409716d	230	k = 8.0 * x + 0.25;
87969c8c	231	u = 0.125Q * k;
	232	/* Small arctan argument. */
	233	t = (x - u) / (1.0 + x * u);
	234	}
	235
	236	/* Arctan of small argument t. */
	237	u = t * t;
	238	p = ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
	239	q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
	240	u = t * u * p / q + t;
	241
	242	/* arctan x = arctan u + arctan t */
	243	u = atantbl[k] + u;
	244	if (sign)
	245	return (-u);
	246	else
	247	return u;
	248	}