lib/libm/asin.c

   1 /*                                                      asin.c
   2  *
   3  *      Inverse circular sine
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, asin();
  10  *
  11  * y = asin( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns radian angle between -pi/2 and +pi/2 whose sine is x.
  18  *
  19  * A rational function of the form x + x**3 P(x**2)/Q(x**2)
  20  * is used for |x| in the interval [0, 0.5].  If |x| > 0.5 it is
  21  * transformed by the identity
  22  *
  23  *    asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
  24  *
  25  *
  26  * ACCURACY:
  27  *
  28  *                      Relative error:
  29  * arithmetic   domain     # trials      peak         rms
  30  *    DEC      -1, 1        40000       2.6e-17     7.1e-18
  31  *    IEEE     -1, 1        10^6        1.9e-16     5.4e-17
  32  *
  33  *
  34  * ERROR MESSAGES:
  35  *
  36  *   message         condition      value returned
  37  * asin domain        |x| > 1           NAN
  38  *
  39  */
  40 \f/*                                                     acos()
  41  *
  42  *      Inverse circular cosine
  43  *
  44  *
  45  *
  46  * SYNOPSIS:
  47  *
  48  * double x, y, acos();
  49  *
  50  * y = acos( x );
  51  *
  52  *
  53  *
  54  * DESCRIPTION:
  55  *
  56  * Returns radian angle between 0 and pi whose cosine
  57  * is x.
  58  *
  59  * Analytically, acos(x) = pi/2 - asin(x).  However if |x| is
  60  * near 1, there is cancellation error in subtracting asin(x)
  61  * from pi/2.  Hence if x < -0.5,
  62  *
  63  *    acos(x) =  pi - 2.0 * asin( sqrt((1+x)/2) );
  64  *
  65  * or if x > +0.5,
  66  *
  67  *    acos(x) =  2.0 * asin(  sqrt((1-x)/2) ).
  68  *
  69  *
  70  * ACCURACY:
  71  *
  72  *                      Relative error:
  73  * arithmetic   domain     # trials      peak         rms
  74  *    DEC       -1, 1       50000       3.3e-17     8.2e-18
  75  *    IEEE      -1, 1       10^6        2.2e-16     6.5e-17
  76  *
  77  *
  78  * ERROR MESSAGES:
  79  *
  80  *   message         condition      value returned
  81  * asin domain        |x| > 1           NAN
  82  */
  83 \f
  84 /*                                                      asin.c  */
  85
  86 /*
  87 Cephes Math Library Release 2.8:  June, 2000
  88 Copyright 1984, 1995, 2000 by Stephen L. Moshier
  89 */
  90
  91 #include "mconf.h"
  92
  93 /* arcsin(x)  =  x + x^3 P(x^2)/Q(x^2)
  94    0 <= x <= 0.625
  95    Peak relative error = 1.2e-18  */
  96 #if UNK
  97 static double P[6] = {
  98  4.253011369004428248960E-3,
  99 -6.019598008014123785661E-1,
 100  5.444622390564711410273E0,
 101 -1.626247967210700244449E1,
 102  1.956261983317594739197E1,
 103 -8.198089802484824371615E0,
 104 };
 105 static double Q[5] = {
 106 /* 1.000000000000000000000E0, */
 107 -1.474091372988853791896E1,
 108  7.049610280856842141659E1,
 109 -1.471791292232726029859E2,
 110  1.395105614657485689735E2,
 111 -4.918853881490881290097E1,
 112 };
 113 #endif
 114 #if DEC
 115 static short P[24] = {
 116 0036213,0056330,0057244,0053234,
 117 0140032,0015011,0114762,0160255,
 118 0040656,0035130,0136121,0067313,
 119 0141202,0014616,0170474,0101731,
 120 0041234,0100076,0151674,0111310,
 121 0141003,0025540,0033165,0077246,
 122 };
 123 static short Q[20] = {
 124 /* 0040200,0000000,0000000,0000000, */
 125 0141153,0155310,0055360,0072530,
 126 0041614,0177001,0027764,0101237,
 127 0142023,0026733,0064653,0133266,
 128 0042013,0101264,0023775,0176351,
 129 0141504,0140420,0050660,0036543,
 130 };
 131 #endif
 132 #if IBMPC
 133 static short P[24] = {
 134 0x8ad3,0x0bd4,0x6b9b,0x3f71,
 135 0x5c16,0x333e,0x4341,0xbfe3,
 136 0x2dd9,0x178a,0xc74b,0x4015,
 137 0x907b,0xde27,0x4331,0xc030,
 138 0x9259,0xda77,0x9007,0x4033,
 139 0xafd5,0x06ce,0x656c,0xc020,
 140 };
 141 static short Q[20] = {
 142 /* 0x0000,0x0000,0x0000,0x3ff0, */
 143 0x0eab,0x0b5e,0x7b59,0xc02d,
 144 0x9054,0x25fe,0x9fc0,0x4051,
 145 0x76d7,0x6d35,0x65bb,0xc062,
 146 0xbf9d,0x84ff,0x7056,0x4061,
 147 0x07ac,0x0a36,0x9822,0xc048,
 148 };
 149 #endif
 150 #if MIEEE
 151 static short P[24] = {
 152 0x3f71,0x6b9b,0x0bd4,0x8ad3,
 153 0xbfe3,0x4341,0x333e,0x5c16,
 154 0x4015,0xc74b,0x178a,0x2dd9,
 155 0xc030,0x4331,0xde27,0x907b,
 156 0x4033,0x9007,0xda77,0x9259,
 157 0xc020,0x656c,0x06ce,0xafd5,
 158 };
 159 static short Q[20] = {
 160 /* 0x3ff0,0x0000,0x0000,0x0000, */
 161 0xc02d,0x7b59,0x0b5e,0x0eab,
 162 0x4051,0x9fc0,0x25fe,0x9054,
 163 0xc062,0x65bb,0x6d35,0x76d7,
 164 0x4061,0x7056,0x84ff,0xbf9d,
 165 0xc048,0x9822,0x0a36,0x07ac,
 166 };
 167 #endif
 168
 169 /* arcsin(1-x) = pi/2 - sqrt(2x)(1+R(x))
 170    0 <= x <= 0.5
 171    Peak relative error = 4.2e-18  */
 172 #if UNK
 173 static double R[5] = {
 174  2.967721961301243206100E-3,
 175 -5.634242780008963776856E-1,
 176  6.968710824104713396794E0,
 177 -2.556901049652824852289E1,
 178  2.853665548261061424989E1,
 179 };
 180 static double S[4] = {
 181 /* 1.000000000000000000000E0, */
 182 -2.194779531642920639778E1,
 183  1.470656354026814941758E2,
 184 -3.838770957603691357202E2,
 185  3.424398657913078477438E2,
 186 };
 187 #endif
 188 #if DEC
 189 static short R[20] = {
 190 0036102,0077034,0142164,0174103,
 191 0140020,0036222,0147711,0044173,
 192 0040736,0177655,0153631,0171523,
 193 0141314,0106525,0060015,0055474,
 194 0041344,0045422,0003630,0040344,
 195 };
 196 static short S[16] = {
 197 /* 0040200,0000000,0000000,0000000, */
 198 0141257,0112425,0132772,0166136,
 199 0042023,0010315,0075523,0175020,
 200 0142277,0170104,0126203,0017563,
 201 0042253,0034115,0102662,0022757,
 202 };
 203 #endif
 204 #if IBMPC
 205 static short R[20] = {
 206 0x9f08,0x988e,0x4fc3,0x3f68,
 207 0x290f,0x59f9,0x0792,0xbfe2,
 208 0x3e6a,0xbaf3,0xdff5,0x401b,
 209 0xab68,0xac01,0x91aa,0xc039,
 210 0x081d,0x40f3,0x8962,0x403c,
 211 };
 212 static short S[16] = {
 213 /* 0x0000,0x0000,0x0000,0x3ff0, */
 214 0x5d8c,0xb6bf,0xf2a2,0xc035,
 215 0x7f42,0xaf6a,0x6219,0x4062,
 216 0x63ee,0x9590,0xfe08,0xc077,
 217 0x44be,0xb0b6,0x6709,0x4075,
 218 };
 219 #endif
 220 #if MIEEE
 221 static short R[20] = {
 222 0x3f68,0x4fc3,0x988e,0x9f08,
 223 0xbfe2,0x0792,0x59f9,0x290f,
 224 0x401b,0xdff5,0xbaf3,0x3e6a,
 225 0xc039,0x91aa,0xac01,0xab68,
 226 0x403c,0x8962,0x40f3,0x081d,
 227 };
 228 static short S[16] = {
 229 /* 0x3ff0,0x0000,0x0000,0x0000, */
 230 0xc035,0xf2a2,0xb6bf,0x5d8c,
 231 0x4062,0x6219,0xaf6a,0x7f42,
 232 0xc077,0xfe08,0x9590,0x63ee,
 233 0x4075,0x6709,0xb0b6,0x44be,
 234 };
 235 #endif
 236
 237 /* pi/2 = PIO2 + MOREBITS.  */
 238 #ifdef DEC
 239 #define MOREBITS 5.721188726109831840122E-18
 240 #else
 241 #define MOREBITS 6.123233995736765886130E-17
 242 #endif
 243
 244 #ifdef ANSIPROT
 245 extern double polevl ( double, void *, int );
 246 extern double p1evl ( double, void *, int );
 247 extern double sqrt ( double );
 248 double asin ( double );
 249 #else
 250 double sqrt(), polevl(), p1evl();
 251 double asin();
 252 #endif
 253 extern double PIO2, PIO4, NAN;
 254
 255 double asin(x)
 256 double x;
 257 {
 258 double a, p, z, zz;
 259 short sign;
 260
 261 if( x > 0 )
 262         {
 263         sign = 1;
 264         a = x;
 265         }
 266 else
 267         {
 268         sign = -1;
 269         a = -x;
 270         }
 271
 272 if( a > 1.0 )
 273         {
 274         mtherr( "asin", DOMAIN );
 275         return( NAN );
 276         }
 277
 278 if( a > 0.625 )
 279         {
 280         /* arcsin(1-x) = pi/2 - sqrt(2x)(1+R(x))  */
 281         zz = 1.0 - a;
 282         p = zz * polevl( zz, R, 4)/p1evl( zz, S, 4);
 283         zz = sqrt(zz+zz);
 284         z = PIO4 - zz;
 285         zz = zz * p - MOREBITS;
 286         z = z - zz;
 287         z = z + PIO4;
 288         }
 289 else
 290         {
 291         if( a < 1.0e-8 )
 292                 {
 293                 return(x);
 294                 }
 295         zz = a * a;
 296         z = zz * polevl( zz, P, 5)/p1evl( zz, Q, 5);
 297         z = a * z + a;
 298         }
 299 if( sign < 0 )
 300         z = -z;
 301 return(z);
 302 }
 303
 304
 305
 306 double acos(x)
 307 double x;
 308 {
 309 double z;
 310
 311 if( (x < -1.0) || (x > 1.0) )
 312         {
 313         mtherr( "acos", DOMAIN );
 314         return( NAN );
 315         }
 316 if( x > 0.5 )
 317         {
 318         return( 2.0 * asin(  sqrt(0.5 - 0.5*x) ) );
 319         }
 320 z = PIO4 - asin(x);
 321 z = z + MOREBITS;
 322 z = z + PIO4;
 323 return( z );
 324 }