[thirdparty/glibc.git] / sysdeps / ieee754 / ldbl-128ibm / e_asinl.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, see
    <http://www.gnu.org/licenses/>.  */

/* __ieee754_asin(x)
 * Method :
 *	Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
 *	we approximate asin(x) on [0,0.5] by
 *		asin(x) = x + x*x^2*R(x^2)
 *      Between .5 and .625 the approximation is
 *              asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
 *	For x in [0.625,1]
 *		asin(x) = pi/2-2*asin(sqrt((1-x)/2))
 *	Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
 *	then for x>0.98
 *		asin(x) = pi/2 - 2*(s+s*z*R(z))
 *			= pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
 *	For x<=0.98, let pio4_hi = pio2_hi/2, then
 *		f = hi part of s;
 *		c = sqrt(z) - f = (z-f*f)/(s+f) 	...f+c=sqrt(z)
 *	and
 *		asin(x) = pi/2 - 2*(s+s*z*R(z))
 *			= pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
 *			= pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
 *
 * Special cases:
 *	if x is NaN, return x itself;
 *	if |x|>1, return NaN with invalid signal.
 *
 */


#include <float.h>
#include <math.h>
#include <math_private.h>
long double sqrtl (long double);

static const long double
  one = 1.0L,
  huge = 1.0e+300L,
  pio2_hi = 1.5707963267948966192313216916397514420986L,
  pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
  pio4_hi = 7.8539816339744830961566084581987569936977E-1L,

	/* coefficient for R(x^2) */

  /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
     0 <= x <= 0.5
     peak relative error 1.9e-35  */
  pS0 = -8.358099012470680544198472400254596543711E2L,
  pS1 =  3.674973957689619490312782828051860366493E3L,
  pS2 = -6.730729094812979665807581609853656623219E3L,
  pS3 =  6.643843795209060298375552684423454077633E3L,
  pS4 = -3.817341990928606692235481812252049415993E3L,
  pS5 =  1.284635388402653715636722822195716476156E3L,
  pS6 = -2.410736125231549204856567737329112037867E2L,
  pS7 =  2.219191969382402856557594215833622156220E1L,
  pS8 = -7.249056260830627156600112195061001036533E-1L,
  pS9 =  1.055923570937755300061509030361395604448E-3L,

  qS0 = -5.014859407482408326519083440151745519205E3L,
  qS1 =  2.430653047950480068881028451580393430537E4L,
  qS2 = -4.997904737193653607449250593976069726962E4L,
  qS3 =  5.675712336110456923807959930107347511086E4L,
  qS4 = -3.881523118339661268482937768522572588022E4L,
  qS5 =  1.634202194895541569749717032234510811216E4L,
  qS6 = -4.151452662440709301601820849901296953752E3L,
  qS7 =  5.956050864057192019085175976175695342168E2L,
  qS8 = -4.175375777334867025769346564600396877176E1L,
  /* 1.000000000000000000000000000000000000000E0 */

  /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
     -0.0625 <= x <= 0.0625
     peak relative error 3.3e-35  */
  rS0 = -5.619049346208901520945464704848780243887E0L,
  rS1 =  4.460504162777731472539175700169871920352E1L,
  rS2 = -1.317669505315409261479577040530751477488E2L,
  rS3 =  1.626532582423661989632442410808596009227E2L,
  rS4 = -3.144806644195158614904369445440583873264E1L,
  rS5 = -9.806674443470740708765165604769099559553E1L,
  rS6 =  5.708468492052010816555762842394927806920E1L,
  rS7 =  1.396540499232262112248553357962639431922E1L,
  rS8 = -1.126243289311910363001762058295832610344E1L,
  rS9 = -4.956179821329901954211277873774472383512E-1L,
  rS10 =  3.313227657082367169241333738391762525780E-1L,

  sS0 = -4.645814742084009935700221277307007679325E0L,
  sS1 =  3.879074822457694323970438316317961918430E1L,
  sS2 = -1.221986588013474694623973554726201001066E2L,
  sS3 =  1.658821150347718105012079876756201905822E2L,
  sS4 = -4.804379630977558197953176474426239748977E1L,
  sS5 = -1.004296417397316948114344573811562952793E2L,
  sS6 =  7.530281592861320234941101403870010111138E1L,
  sS7 =  1.270735595411673647119592092304357226607E1L,
  sS8 = -1.815144839646376500705105967064792930282E1L,
  sS9 = -7.821597334910963922204235247786840828217E-2L,
  /*  1.000000000000000000000000000000000000000E0 */

 asinr5625 =  5.9740641664535021430381036628424864397707E-1L;


long double
__ieee754_asinl (long double x)
{
  long double a, t, w, p, q, c, r, s;
  int flag;

  if (__glibc_unlikely (isnan (x)))
    return x + x;
  flag = 0;
  a = __builtin_fabsl (x);
  if (a == 1.0L)	/* |x|>= 1 */
    return x * pio2_hi + x * pio2_lo;	/* asin(1)=+-pi/2 with inexact */
  else if (a >= 1.0L)
    return (x - x) / (x - x);	/* asin(|x|>1) is NaN */
  else if (a < 0.5L)
    {
      if (a < 6.938893903907228e-18L) /* |x| < 2**-57 */
	{
	  if (fabsl (x) < LDBL_MIN)
	    {
	      long double force_underflow = x * x;
	      math_force_eval (force_underflow);
	    }
	  long double force_inexact =  huge + x;
	  math_force_eval (force_inexact);
	  return x;		/* return x with inexact if x!=0 */
	}
      else
	{
	  t = x * x;
	  /* Mark to use pS, qS later on.  */
	  flag = 1;
	}
    }
  else if (a < 0.625L)
    {
      t = a - 0.5625;
      p = ((((((((((rS10 * t
		    + rS9) * t
		   + rS8) * t
		  + rS7) * t
		 + rS6) * t
		+ rS5) * t
	       + rS4) * t
	      + rS3) * t
	     + rS2) * t
	    + rS1) * t
	   + rS0) * t;

      q = ((((((((( t
		    + sS9) * t
		  + sS8) * t
		 + sS7) * t
		+ sS6) * t
	       + sS5) * t
	      + sS4) * t
	     + sS3) * t
	    + sS2) * t
	   + sS1) * t
	+ sS0;
      t = asinr5625 + p / q;
      if (x > 0.0L)
	return t;
      else
	return -t;
    }
  else
    {
      /* 1 > |x| >= 0.625 */
      w = one - a;
      t = w * 0.5;
    }

  p = (((((((((pS9 * t
	       + pS8) * t
	      + pS7) * t
	     + pS6) * t
	    + pS5) * t
	   + pS4) * t
	  + pS3) * t
	 + pS2) * t
	+ pS1) * t
       + pS0) * t;

  q = (((((((( t
	      + qS8) * t
	     + qS7) * t
	    + qS6) * t
	   + qS5) * t
	  + qS4) * t
	 + qS3) * t
	+ qS2) * t
       + qS1) * t
    + qS0;

  if (flag) /* 2^-57 < |x| < 0.5 */
    {
      w = p / q;
      return x + x * w;
    }

  s = __ieee754_sqrtl (t);
  if (a > 0.975L)
    {
      w = p / q;
      t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
    }
  else
    {
      w = ldbl_high (s);
      c = (t - w * w) / (s + w);
      r = p / q;
      p = 2.0 * s * r - (pio2_lo - 2.0 * c);
      q = pio4_hi - 2.0 * w;
      t = pio4_hi - (p - q);
    }

  if (x > 0.0L)
    return t;
  else
    return -t;
}
strong_alias (__ieee754_asinl, __asinl_finite)
Commit	Line	Data
f964490f RM	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>
0ac5ae23 UD	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 the
f964490f RM	18	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
59ba27a6 PE	31	License along with this library; if not, see
59ba27a6 PE	32	<http://www.gnu.org/licenses/>. */
f964490f RM	33
	34	/* __ieee754_asin(x)
	35	* Method :
	36	* Since asin(x) = x + x^3/6 + x^53/40 + x^715/336 + ...
	37	* we approximate asin(x) on [0,0.5] by
	38	* asin(x) = x + xx^2R(x^2)
	39	* Between .5 and .625 the approximation is
	40	* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
	41	* For x in [0.625,1]
	42	* asin(x) = pi/2-2*asin(sqrt((1-x)/2))
	43	* Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
	44	* then for x>0.98
	45	* asin(x) = pi/2 - 2(s+sz*R(z))
	46	* = pio2_hi - (2(s+sz*R(z)) - pio2_lo)
	47	* For x<=0.98, let pio4_hi = pio2_hi/2, then
	48	* f = hi part of s;
	49	* c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z)
	50	* and
	51	* asin(x) = pi/2 - 2(s+sz*R(z))
	52	* = pio4_hi+(pio4-2s)-(2szR(z)-pio2_lo)
	53	* = pio4_hi+(pio4-2f)-(2szR(z)-(pio2_lo+2c))
	54	*
	55	* Special cases:
	56	* if x is NaN, return x itself;
	57	* if \|x\|>1, return NaN with invalid signal.
	58	*
	59	*/
	60
	61
ec0ce0d3	62	#include <float.h>
1ed0291c RH	63	#include <math.h>
1ed0291c RH	64	#include <math_private.h>
f964490f RM	65	long double sqrtl (long double);
f964490f RM	66
f964490f	67	static const long double
f964490f RM	68	one = 1.0L,
	69	huge = 1.0e+300L,
	70	pio2_hi = 1.5707963267948966192313216916397514420986L,
	71	pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
	72	pio4_hi = 7.8539816339744830961566084581987569936977E-1L,
	73
	74	/* coefficient for R(x^2) */
	75
	76	/* asin(x) = x + x^3 pS(x^2) / qS(x^2)
	77	0 <= x <= 0.5
	78	peak relative error 1.9e-35 */
	79	pS0 = -8.358099012470680544198472400254596543711E2L,
	80	pS1 = 3.674973957689619490312782828051860366493E3L,
	81	pS2 = -6.730729094812979665807581609853656623219E3L,
	82	pS3 = 6.643843795209060298375552684423454077633E3L,
	83	pS4 = -3.817341990928606692235481812252049415993E3L,
	84	pS5 = 1.284635388402653715636722822195716476156E3L,
	85	pS6 = -2.410736125231549204856567737329112037867E2L,
	86	pS7 = 2.219191969382402856557594215833622156220E1L,
	87	pS8 = -7.249056260830627156600112195061001036533E-1L,
	88	pS9 = 1.055923570937755300061509030361395604448E-3L,
	89
	90	qS0 = -5.014859407482408326519083440151745519205E3L,
	91	qS1 = 2.430653047950480068881028451580393430537E4L,
	92	qS2 = -4.997904737193653607449250593976069726962E4L,
	93	qS3 = 5.675712336110456923807959930107347511086E4L,
	94	qS4 = -3.881523118339661268482937768522572588022E4L,
	95	qS5 = 1.634202194895541569749717032234510811216E4L,
	96	qS6 = -4.151452662440709301601820849901296953752E3L,
	97	qS7 = 5.956050864057192019085175976175695342168E2L,
	98	qS8 = -4.175375777334867025769346564600396877176E1L,
	99	/* 1.000000000000000000000000000000000000000E0 */
	100
	101	/* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
	102	-0.0625 <= x <= 0.0625
	103	peak relative error 3.3e-35 */
	104	rS0 = -5.619049346208901520945464704848780243887E0L,
	105	rS1 = 4.460504162777731472539175700169871920352E1L,
	106	rS2 = -1.317669505315409261479577040530751477488E2L,
	107	rS3 = 1.626532582423661989632442410808596009227E2L,
	108	rS4 = -3.144806644195158614904369445440583873264E1L,
	109	rS5 = -9.806674443470740708765165604769099559553E1L,
	110	rS6 = 5.708468492052010816555762842394927806920E1L,
	111	rS7 = 1.396540499232262112248553357962639431922E1L,
	112	rS8 = -1.126243289311910363001762058295832610344E1L,
	113	rS9 = -4.956179821329901954211277873774472383512E-1L,
	114	rS10 = 3.313227657082367169241333738391762525780E-1L,
	115
	116	sS0 = -4.645814742084009935700221277307007679325E0L,
	117	sS1 = 3.879074822457694323970438316317961918430E1L,
	118	sS2 = -1.221986588013474694623973554726201001066E2L,
	119	sS3 = 1.658821150347718105012079876756201905822E2L,
	120	sS4 = -4.804379630977558197953176474426239748977E1L,
	121	sS5 = -1.004296417397316948114344573811562952793E2L,
	122	sS6 = 7.530281592861320234941101403870010111138E1L,
	123	sS7 = 1.270735595411673647119592092304357226607E1L,
	124	sS8 = -1.815144839646376500705105967064792930282E1L,
	125	sS9 = -7.821597334910963922204235247786840828217E-2L,
	126	/* 1.000000000000000000000000000000000000000E0 */
	127
	128	asinr5625 = 5.9740641664535021430381036628424864397707E-1L;
	129
	130
	131
f964490f RM	132	long double
f964490f RM	133	__ieee754_asinl (long double x)
f964490f	134	{
4ebd120c	135	long double a, t, w, p, q, c, r, s;
31dc8730	136	int flag;
f964490f	137
d81f90cc	138	if (__glibc_unlikely (isnan (x)))
6f10289e	139	return x + x;
f964490f	140	flag = 0;
4ebd120c AM	141	a = __builtin_fabsl (x);
4ebd120c AM	142	if (a == 1.0L) /* \|x\|>= 1 */
31dc8730	143	return x * pio2_hi + x * pio2_lo; /* asin(1)=+-pi/2 with inexact */
4ebd120c	144	else if (a >= 1.0L)
31dc8730	145	return (x - x) / (x - x); /* asin(\|x\|>1) is NaN */
4ebd120c	146	else if (a < 0.5L)
f964490f	147	{
4ebd120c	148	if (a < 6.938893903907228e-18L) /* \|x\| < 2*-57 /
f964490f	149	{
ec0ce0d3 JM	150	if (fabsl (x) < LDBL_MIN)
	151	{
	152	long double force_underflow = x * x;
	153	math_force_eval (force_underflow);
	154	}
fded7ed6 JM	155	long double force_inexact = huge + x;
	156	math_force_eval (force_inexact);
	157	return x; /* return x with inexact if x!=0 */
f964490f RM	158	}
	159	else
	160	{
	161	t = x * x;
	162	/* Mark to use pS, qS later on. */
	163	flag = 1;
	164	}
	165	}
4ebd120c	166	else if (a < 0.625L)
f964490f	167	{
4ebd120c	168	t = a - 0.5625;
f964490f RM	169	p = ((((((((((rS10 * t
	170	+ rS9) * t
	171	+ rS8) * t
	172	+ rS7) * t
	173	+ rS6) * t
	174	+ rS5) * t
	175	+ rS4) * t
	176	+ rS3) * t
	177	+ rS2) * t
	178	+ rS1) * t
	179	+ rS0) * t;
	180
	181	q = ((((((((( t
	182	+ sS9) * t
	183	+ sS8) * t
	184	+ sS7) * t
	185	+ sS6) * t
	186	+ sS5) * t
	187	+ sS4) * t
	188	+ sS3) * t
	189	+ sS2) * t
	190	+ sS1) * t
	191	+ sS0;
	192	t = asinr5625 + p / q;
31dc8730	193	if (x > 0.0L)
f964490f RM	194	return t;
	195	else
	196	return -t;
	197	}
	198	else
	199	{
	200	/* 1 > \|x\| >= 0.625 */
4ebd120c	201	w = one - a;
f964490f RM	202	t = w * 0.5;
	203	}
	204
	205	p = (((((((((pS9 * t
	206	+ pS8) * t
	207	+ pS7) * t
	208	+ pS6) * t
	209	+ pS5) * t
	210	+ pS4) * t
	211	+ pS3) * t
	212	+ pS2) * t
	213	+ pS1) * t
	214	+ pS0) * t;
	215
	216	q = (((((((( t
	217	+ qS8) * t
	218	+ qS7) * t
	219	+ qS6) * t
	220	+ qS5) * t
	221	+ qS4) * t
	222	+ qS3) * t
	223	+ qS2) * t
	224	+ qS1) * t
	225	+ qS0;
	226
	227	if (flag) /* 2^-57 < \|x\| < 0.5 */
	228	{
	229	w = p / q;
	230	return x + x * w;
	231	}
	232
	233	s = __ieee754_sqrtl (t);
4ebd120c	234	if (a > 0.975L)
f964490f RM	235	{
	236	w = p / q;
	237	t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
	238	}
	239	else
	240	{
4ebd120c	241	w = ldbl_high (s);
f964490f RM	242	c = (t - w * w) / (s + w);
	243	r = p / q;
	244	p = 2.0 * s * r - (pio2_lo - 2.0 * c);
	245	q = pio4_hi - 2.0 * w;
	246	t = pio4_hi - (p - q);
	247	}
	248
31dc8730	249	if (x > 0.0L)
f964490f RM	250	return t;
	251	else
	252	return -t;
	253	}
0ac5ae23	254	strong_alias (__ieee754_asinl, __asinl_finite)