Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : ASINL.C

/* asinl.c
*
* Inverse circular sine, long double precision
*
*
*
* SYNOPSIS:
*
* double x, y, asinl();
*
* y = asinl( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between -pi/2 and +pi/2 whose sine is x.
*
* A rational function of the form x + x**3 P(x**2)/Q(x**2)
* is used for |x| in the interval [0, 0.5]. If |x| > 0.5 it is
* transformed by the identity
*
* asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1, 1 30000 2.7e-19 4.8e-20
*
*
* ERROR MESSAGES:
*
* message condition value returned
* asin domain |x| > 1 0.0
*
*/
/* acosl()
*
* Inverse circular cosine, long double precision
*
*
*
* SYNOPSIS:
*
* double x, y, acosl();
*
* y = acosl( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between -pi/2 and +pi/2 whose cosine
* is x.
*
* Analytically, acos(x) = pi/2 - asin(x). However if |x| is
* near 1, there is cancellation error in subtracting asin(x)
* from pi/2. Hence if x < -0.5,
*
* acos(x) = pi - 2.0 * asin( sqrt((1+x)/2) );
*
* or if x > +0.5,
*
* acos(x) = 2.0 * asin( sqrt((1-x)/2) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1, 1 30000 1.4e-19 3.5e-20
*
*
* ERROR MESSAGES:
*
* message condition value returned
* asin domain |x| > 1 0.0
*/

/* asin.c */

/*
Cephes Math Library Release 2.2: December, 1990
Copyright 1984, 1990 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

#include "mconf.h"

#ifdef UNK
static long double P[] = {
3.7769340062433674871612E-3L,
-6.1212919176969202969441E-1L,
5.9303993515791417710775E0L,
-1.8631697621590161441592E1L,
2.3314603132141795720634E1L,
-1.0087146579384916260197E1L,
};
static long double Q[] = {
/* 1.0000000000000000000000E0L,*/
-1.5684335624873146511217E1L,
7.8702951549021104258866E1L,
-1.7078401170625864261444E2L,
1.6712291455718995937376E2L,
-6.0522879476309497128868E1L,
};
#endif

#ifdef IBMPC
static short P[] = {
0x59d1,0x3509,0x7009,0xf786,0x3ff6,
0xbe97,0x93e6,0x7fab,0x9cb4,0xbffe,
0x8bf5,0x6810,0xd4dc,0xbdc5,0x4001,
0x9bd4,0x8d86,0xb77b,0x950d,0xc003,
0x3b0f,0x9e25,0x4ea5,0xba84,0x4003,
0xea38,0xc6a9,0xf3cf,0xa164,0xc002,
};
static short Q[] = {
/*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
0x1229,0x8516,0x09e9,0xfaf3,0xc002,
0xb5c3,0xf36f,0xe943,0x9d67,0x4005,
0xe11a,0xbe0f,0xb4fd,0xaac8,0xc006,
0x4c69,0x1355,0x7754,0xa71f,0x4006,
0xded7,0xa9fe,0x6db7,0xf217,0xc004,
};
#endif

#ifdef MIEEE
static long P[] = {
0x3ff60000,0xf7867009,0x350959d1,
0xbffe0000,0x9cb47fab,0x93e6be97,
0x40010000,0xbdc5d4dc,0x68108bf5,
0xc0030000,0x950db77b,0x8d869bd4,
0x40030000,0xba844ea5,0x9e253b0f,
0xc0020000,0xa164f3cf,0xc6a9ea38,
};
static long Q[] = {
/*0x3fff0000,0x80000000,0x00000000,*/
0xc0020000,0xfaf309e9,0x85161229,
0x40050000,0x9d67e943,0xf36fb5c3,
0xc0060000,0xaac8b4fd,0xbe0fe11a,
0x40060000,0xa71f7754,0x13554c69,
0xc0040000,0xf2176db7,0xa9feded7,
};
#endif

long double asinl(x)
long double x;
{
long double a, p, z, zz;
long double ldexpl(), sqrtl(), polevll(), p1evll();
short sign, flag;
extern long double PIO2L;

if( x > 0 )
{
sign = 1;
a = x;
}
else
{
sign = -1;
a = -x;
}

if( a > 1.0L )
{
mtherr( "asinl", DOMAIN );
return( 0.0L );
}

if( a < 1.0e-8L )
{
z = a;
goto done;
}

if( a > 0.5L )
{
zz = 0.5L -a;
zz = ldexpl( zz + 0.5L, -1 );
z = sqrtl( zz );
flag = 1;
}
else
{
z = a;
zz = z * z;
flag = 0;
}

p = zz * polevll( zz, P, 5)/p1evll( zz, Q, 5);
z = z * p + z;
if( flag != 0 )
{
z = z + z;
z = PIO2L - z;
}
done:
if( sign < 0 )
z = -z;
return(z);
}


extern long double PIO2L, PIL;

long double acosl(x)
long double x;
{
long double asinl(), sqrtl();

if( x < -1.0L )
goto domerr;

if( x < -0.5L)
return( PIL - 2.0L * asinl( sqrtl(0.5L*(1.0L+x)) ) );

if( x > 1.0L )
{
domerr: mtherr( "acosl", DOMAIN );
return( 0.0L );
}

if( x > 0.5L )
return( 2.0L * asinl( sqrtl(0.5L*(1.0L-x) ) ) );

return( PIO2L - asinl(x) );
}