Print this page




  88 #define one     xxx[0]
  89 #define huge    xxx[1]
  90 #define pio2_hi xxx[2]
  91 #define pio2_lo xxx[3]
  92 #define pio4_hi xxx[4]
  93 #define pS0     xxx[5]
  94 #define pS1     xxx[6]
  95 #define pS2     xxx[7]
  96 #define pS3     xxx[8]
  97 #define pS4     xxx[9]
  98 #define pS5     xxx[10]
  99 #define qS1     xxx[11]
 100 #define qS2     xxx[12]
 101 #define qS3     xxx[13]
 102 #define qS4     xxx[14]
 103 /* INDENT ON */
 104 
 105 double
 106 asin(double x) {
 107         double t, w, p, q, c, r, s;
 108         int hx, ix;
 109 
 110         hx = ((int *) &x)[HIWORD];
 111         ix = hx & 0x7fffffff;
 112         if (ix >= 0x3ff00000) {      /* |x| >= 1 */
 113                 if (((ix - 0x3ff00000) | ((int *) &x)[LOWORD]) == 0)
 114                         /* asin(1)=+-pi/2 with inexact */
 115                         return x * pio2_hi + x * pio2_lo;
 116                 else if (isnan(x))
 117 #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
 118                         return ix >= 0x7ff80000 ? x : (x - x) / (x - x);
 119                         /* assumes sparc-like QNaN */
 120 #else
 121                         return (x - x) / (x - x);       /* asin(|x|>1) is NaN */
 122 #endif
 123                 else
 124                         return _SVID_libm_err(x, x, 2);
 125         }
 126         else if (ix < 0x3fe00000) {  /* |x| < 0.5 */
 127                 if (ix < 0x3e400000) {       /* if |x| < 2**-27 */
 128                         if (huge + x > one)
 129                                 return x;       /* return x with inexact if
 130                                                  * x != 0 */
 131                 }
 132                 else
 133                         t = x * x;
 134                 p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 +
 135                         t * (pS4 + t * pS5)))));
 136                 q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4)));
 137                 w = p / q;
 138                 return x + x * w;
 139         }
 140         /* 1 > |x| >= 0.5 */
 141         w = one - fabs(x);
 142         t = w * 0.5;
 143         p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
 144         q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4)));
 145         s = sqrt(t);
 146         if (ix >= 0x3FEF3333) {      /* if |x| > 0.975 */
 147                 w = p / q;
 148                 t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
 149         }
 150         else {
 151                 w = s;
 152                 ((int *) &w)[LOWORD] = 0;


  88 #define one     xxx[0]
  89 #define huge    xxx[1]
  90 #define pio2_hi xxx[2]
  91 #define pio2_lo xxx[3]
  92 #define pio4_hi xxx[4]
  93 #define pS0     xxx[5]
  94 #define pS1     xxx[6]
  95 #define pS2     xxx[7]
  96 #define pS3     xxx[8]
  97 #define pS4     xxx[9]
  98 #define pS5     xxx[10]
  99 #define qS1     xxx[11]
 100 #define qS2     xxx[12]
 101 #define qS3     xxx[13]
 102 #define qS4     xxx[14]
 103 /* INDENT ON */
 104 
 105 double
 106 asin(double x) {
 107         double t, w, p, q, c, r, s;
 108         int hx, ix, i;
 109 
 110         hx = ((int *) &x)[HIWORD];
 111         ix = hx & 0x7fffffff;
 112         if (ix >= 0x3ff00000) {      /* |x| >= 1 */
 113                 if (((ix - 0x3ff00000) | ((int *) &x)[LOWORD]) == 0)
 114                         /* asin(1)=+-pi/2 with inexact */
 115                         return x * pio2_hi + x * pio2_lo;
 116                 else if (isnan(x))
 117 #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
 118                         return ix >= 0x7ff80000 ? x : (x - x) / (x - x);
 119                         /* assumes sparc-like QNaN */
 120 #else
 121                         return (x - x) / (x - x);       /* asin(|x|>1) is NaN */
 122 #endif
 123                 else
 124                         return _SVID_libm_err(x, x, 2);
 125         }
 126         else if (ix < 0x3fe00000) {  /* |x| < 0.5 */
 127                 if (ix < 0x3e400000) {       /* if |x| < 2**-27 */
 128                         if ((i = (int) x) == 0)
 129                                 return x;       /* return x with inexact if
 130                                                  * x != 0 */
 131                 }

 132                 t = x * x;
 133                 p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 +
 134                         t * (pS4 + t * pS5)))));
 135                 q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4)));
 136                 w = p / q;
 137                 return x + x * w;
 138         }
 139         /* 1 > |x| >= 0.5 */
 140         w = one - fabs(x);
 141         t = w * 0.5;
 142         p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
 143         q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4)));
 144         s = sqrt(t);
 145         if (ix >= 0x3FEF3333) {      /* if |x| > 0.975 */
 146                 w = p / q;
 147                 t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
 148         }
 149         else {
 150                 w = s;
 151                 ((int *) &w)[LOWORD] = 0;