text: /*							asin.c
   *
   *	Inverse circular sine
   *
   *
   *
   * SYNOPSIS:
   *
   * double x, y, cephes_asin();
   *
   * y = cephes_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
   *
   *    cephes_asin(x) = pi/2 - 2 cephes_asin( cephes_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
   *
   */
  /*							cephes_acos()
   *
   *	Inverse circular cosine
   *
   *
   *
   * SYNOPSIS:
   *
   * double x, y, cephes_acos();
   *
   * y = cephes_acos( x );
   *
   *
   *
   * DESCRIPTION:
   *
   * Returns radian angle between 0 and pi whose cosine
   * is x.
   *
   * Analytically, cephes_acos(x) = pi/2 - cephes_asin(x).  However if |x| is
   * near 1, there is cancellation error in subtracting cephes_asin(x)
   * from pi/2.  Hence if x < -0.5,
   *
   *    cephes_acos(x) =	 pi - 2.0 * cephes_asin( cephes_sqrt((1+x)/2) );
   *
   * or if x > +0.5,
   *
   *    cephes_acos(x) =	 2.0 * cephes_asin(  cephes_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 "mconf.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 - cephes_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 cephes_polevl ( double, void *, int );
  extern double cephes_p1evl ( double, void *, int );
  extern double cephes_sqrt ( double );
  double cephes_asin ( double );
  #else
  double cephes_sqrt(), cephes_polevl(), cephes_p1evl();
  double cephes_asin();
  #endif
  extern double PIO2, PIO4, NAN;
  
  double cephes_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 )
  	{
  	cephes_mtherr( "asin", DOMAIN );
  	return( NAN );
  	}
  
  if( a > 0.625 )
  	{
  	/* arcsin(1-x) = pi/2 - cephes_sqrt(2x)(1+R(x))  */
  	zz = 1.0 - a;
  	p = zz * cephes_polevl( zz, R, 4)/cephes_p1evl( zz, S, 4);
  	zz = cephes_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 * cephes_polevl( zz, P, 5)/cephes_p1evl( zz, Q, 5);
  	z = a * z + a;
  	}
  if( sign < 0 )
  	z = -z;
  return(z);
  }
  
  
  
  double cephes_acos(x)
  double x;
  {
  double z;
  
  if( (x < -1.0) || (x > 1.0) )
  	{
  	cephes_mtherr( "acos", DOMAIN );
  	return( NAN );
  	}
  if( x > 0.5 )
  	{
  	return( 2.0 * cephes_asin(  cephes_sqrt(0.5 - 0.5*x) ) );
  	}
  z = PIO4 - cephes_asin(x);
  z = z + MOREBITS;
  z = z + PIO4;
  return( z );
  }
  