summaryrefslogtreecommitdiff
path: root/libm/double/asin.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/asin.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/double/asin.c')
-rw-r--r--libm/double/asin.c324
1 files changed, 0 insertions, 324 deletions
diff --git a/libm/double/asin.c b/libm/double/asin.c
deleted file mode 100644
index 1f83eccc8..000000000
--- a/libm/double/asin.c
+++ /dev/null
@@ -1,324 +0,0 @@
-/* asin.c
- *
- * Inverse circular sine
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, asin();
- *
- * y = asin( 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
- * DEC -1, 1 40000 2.6e-17 7.1e-18
- * IEEE -1, 1 10^6 1.9e-16 5.4e-17
- *
- *
- * ERROR MESSAGES:
- *
- * message condition value returned
- * asin domain |x| > 1 NAN
- *
- */
- /* acos()
- *
- * Inverse circular cosine
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, acos();
- *
- * y = acos( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns radian angle between 0 and pi 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
- * DEC -1, 1 50000 3.3e-17 8.2e-18
- * IEEE -1, 1 10^6 2.2e-16 6.5e-17
- *
- *
- * ERROR MESSAGES:
- *
- * message condition value returned
- * asin domain |x| > 1 NAN
- */
-
-/* asin.c */
-
-/*
-Cephes Math Library Release 2.8: June, 2000
-Copyright 1984, 1995, 2000 by Stephen L. Moshier
-*/
-
-#include <math.h>
-
-/* arcsin(x) = x + x^3 P(x^2)/Q(x^2)
- 0 <= x <= 0.625
- Peak relative error = 1.2e-18 */
-#if UNK
-static double P[6] = {
- 4.253011369004428248960E-3,
--6.019598008014123785661E-1,
- 5.444622390564711410273E0,
--1.626247967210700244449E1,
- 1.956261983317594739197E1,
--8.198089802484824371615E0,
-};
-static double Q[5] = {
-/* 1.000000000000000000000E0, */
--1.474091372988853791896E1,
- 7.049610280856842141659E1,
--1.471791292232726029859E2,
- 1.395105614657485689735E2,
--4.918853881490881290097E1,
-};
-#endif
-#if DEC
-static short P[24] = {
-0036213,0056330,0057244,0053234,
-0140032,0015011,0114762,0160255,
-0040656,0035130,0136121,0067313,
-0141202,0014616,0170474,0101731,
-0041234,0100076,0151674,0111310,
-0141003,0025540,0033165,0077246,
-};
-static short Q[20] = {
-/* 0040200,0000000,0000000,0000000, */
-0141153,0155310,0055360,0072530,
-0041614,0177001,0027764,0101237,
-0142023,0026733,0064653,0133266,
-0042013,0101264,0023775,0176351,
-0141504,0140420,0050660,0036543,
-};
-#endif
-#if IBMPC
-static short P[24] = {
-0x8ad3,0x0bd4,0x6b9b,0x3f71,
-0x5c16,0x333e,0x4341,0xbfe3,
-0x2dd9,0x178a,0xc74b,0x4015,
-0x907b,0xde27,0x4331,0xc030,
-0x9259,0xda77,0x9007,0x4033,
-0xafd5,0x06ce,0x656c,0xc020,
-};
-static short Q[20] = {
-/* 0x0000,0x0000,0x0000,0x3ff0, */
-0x0eab,0x0b5e,0x7b59,0xc02d,
-0x9054,0x25fe,0x9fc0,0x4051,
-0x76d7,0x6d35,0x65bb,0xc062,
-0xbf9d,0x84ff,0x7056,0x4061,
-0x07ac,0x0a36,0x9822,0xc048,
-};
-#endif
-#if MIEEE
-static short P[24] = {
-0x3f71,0x6b9b,0x0bd4,0x8ad3,
-0xbfe3,0x4341,0x333e,0x5c16,
-0x4015,0xc74b,0x178a,0x2dd9,
-0xc030,0x4331,0xde27,0x907b,
-0x4033,0x9007,0xda77,0x9259,
-0xc020,0x656c,0x06ce,0xafd5,
-};
-static short Q[20] = {
-/* 0x3ff0,0x0000,0x0000,0x0000, */
-0xc02d,0x7b59,0x0b5e,0x0eab,
-0x4051,0x9fc0,0x25fe,0x9054,
-0xc062,0x65bb,0x6d35,0x76d7,
-0x4061,0x7056,0x84ff,0xbf9d,
-0xc048,0x9822,0x0a36,0x07ac,
-};
-#endif
-
-/* arcsin(1-x) = pi/2 - sqrt(2x)(1+R(x))
- 0 <= x <= 0.5
- Peak relative error = 4.2e-18 */
-#if UNK
-static double R[5] = {
- 2.967721961301243206100E-3,
--5.634242780008963776856E-1,
- 6.968710824104713396794E0,
--2.556901049652824852289E1,
- 2.853665548261061424989E1,
-};
-static double S[4] = {
-/* 1.000000000000000000000E0, */
--2.194779531642920639778E1,
- 1.470656354026814941758E2,
--3.838770957603691357202E2,
- 3.424398657913078477438E2,
-};
-#endif
-#if DEC
-static short R[20] = {
-0036102,0077034,0142164,0174103,
-0140020,0036222,0147711,0044173,
-0040736,0177655,0153631,0171523,
-0141314,0106525,0060015,0055474,
-0041344,0045422,0003630,0040344,
-};
-static short S[16] = {
-/* 0040200,0000000,0000000,0000000, */
-0141257,0112425,0132772,0166136,
-0042023,0010315,0075523,0175020,
-0142277,0170104,0126203,0017563,
-0042253,0034115,0102662,0022757,
-};
-#endif
-#if IBMPC
-static short R[20] = {
-0x9f08,0x988e,0x4fc3,0x3f68,
-0x290f,0x59f9,0x0792,0xbfe2,
-0x3e6a,0xbaf3,0xdff5,0x401b,
-0xab68,0xac01,0x91aa,0xc039,
-0x081d,0x40f3,0x8962,0x403c,
-};
-static short S[16] = {
-/* 0x0000,0x0000,0x0000,0x3ff0, */
-0x5d8c,0xb6bf,0xf2a2,0xc035,
-0x7f42,0xaf6a,0x6219,0x4062,
-0x63ee,0x9590,0xfe08,0xc077,
-0x44be,0xb0b6,0x6709,0x4075,
-};
-#endif
-#if MIEEE
-static short R[20] = {
-0x3f68,0x4fc3,0x988e,0x9f08,
-0xbfe2,0x0792,0x59f9,0x290f,
-0x401b,0xdff5,0xbaf3,0x3e6a,
-0xc039,0x91aa,0xac01,0xab68,
-0x403c,0x8962,0x40f3,0x081d,
-};
-static short S[16] = {
-/* 0x3ff0,0x0000,0x0000,0x0000, */
-0xc035,0xf2a2,0xb6bf,0x5d8c,
-0x4062,0x6219,0xaf6a,0x7f42,
-0xc077,0xfe08,0x9590,0x63ee,
-0x4075,0x6709,0xb0b6,0x44be,
-};
-#endif
-
-/* pi/2 = PIO2 + MOREBITS. */
-#ifdef DEC
-#define MOREBITS 5.721188726109831840122E-18
-#else
-#define MOREBITS 6.123233995736765886130E-17
-#endif
-
-#ifdef ANSIPROT
-extern double polevl ( double, void *, int );
-extern double p1evl ( double, void *, int );
-extern double sqrt ( double );
-double asin ( double );
-#else
-double sqrt(), polevl(), p1evl();
-double asin();
-#endif
-extern double PIO2, PIO4, NAN;
-
-double asin(x)
-double x;
-{
-double a, p, z, zz;
-short sign;
-
-if( x > 0 )
- {
- sign = 1;
- a = x;
- }
-else
- {
- sign = -1;
- a = -x;
- }
-
-if( a > 1.0 )
- {
- mtherr( "asin", DOMAIN );
- return( NAN );
- }
-
-if( a > 0.625 )
- {
- /* arcsin(1-x) = pi/2 - sqrt(2x)(1+R(x)) */
- zz = 1.0 - a;
- p = zz * polevl( zz, R, 4)/p1evl( zz, S, 4);
- zz = sqrt(zz+zz);
- z = PIO4 - zz;
- zz = zz * p - MOREBITS;
- z = z - zz;
- z = z + PIO4;
- }
-else
- {
- if( a < 1.0e-8 )
- {
- return(x);
- }
- zz = a * a;
- z = zz * polevl( zz, P, 5)/p1evl( zz, Q, 5);
- z = a * z + a;
- }
-if( sign < 0 )
- z = -z;
-return(z);
-}
-
-
-
-double acos(x)
-double x;
-{
-double z;
-
-if( (x < -1.0) || (x > 1.0) )
- {
- mtherr( "acos", DOMAIN );
- return( NAN );
- }
-if( x > 0.5 )
- {
- return( 2.0 * asin( sqrt(0.5 - 0.5*x) ) );
- }
-z = PIO4 - asin(x);
-z = z + MOREBITS;
-z = z + PIO4;
-return( z );
-}