gonzui

    1: /* s_atanf.c -- float version of s_atan.c.
    2:  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
    3:  */
    4: 
    5: /*
    6:  * ====================================================
    7:  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    8:  *
    9:  * Developed at SunPro, a Sun Microsystems, Inc. business.
   10:  * Permission to use, copy, modify, and distribute this
   11:  * software is freely granted, provided that this notice 
   12:  * is preserved.
   13:  * ====================================================
   14:  */
   15: 
   16: #include "math.h"
   17: #include "math_private.h"
   18: 
   19: static const float atanhi[] = {
   20:   4.6364760399e-01, /* atan(0.5)hi 0x3eed6338 */
   21:   7.8539812565e-01, /* atan(1.0)hi 0x3f490fda */
   22:   9.8279368877e-01, /* atan(1.5)hi 0x3f7b985e */
   23:   1.5707962513e+00, /* atan(inf)hi 0x3fc90fda */
   24: };
   25: 
   26: static const float atanlo[] = {
   27:   5.0121582440e-09, /* atan(0.5)lo 0x31ac3769 */
   28:   3.7748947079e-08, /* atan(1.0)lo 0x33222168 */
   29:   3.4473217170e-08, /* atan(1.5)lo 0x33140fb4 */
   30:   7.5497894159e-08, /* atan(inf)lo 0x33a22168 */
   31: };
   32: 
   33: static const float aT[] = {
   34:   3.3333334327e-01, /* 0x3eaaaaaa */
   35:  -2.0000000298e-01, /* 0xbe4ccccd */
   36:   1.4285714924e-01, /* 0x3e124925 */
   37:  -1.1111110449e-01, /* 0xbde38e38 */
   38:   9.0908870101e-02, /* 0x3dba2e6e */
   39:  -7.6918758452e-02, /* 0xbd9d8795 */
   40:   6.6610731184e-02, /* 0x3d886b35 */
   41:  -5.8335702866e-02, /* 0xbd6ef16b */
   42:   4.9768779427e-02, /* 0x3d4bda59 */
   43:  -3.6531571299e-02, /* 0xbd15a221 */
   44:   1.6285819933e-02, /* 0x3c8569d7 */
   45: };
   46: 
   47: static const float 
   48: one   = 1.0,
   49: huge   = 1.0e30;
   50: 
   51: float
   52: atanf(float x)
   53: {
   54:         float w,s1,s2,z;
   55:         int32_t ix,hx,id;
   56: 
   57:         GET_FLOAT_WORD(hx,x);
   58:         ix = hx&0x7fffffff;
   59:         if(ix>=0x50800000) {   /* if |x| >= 2^34 */
   60:             if(ix>0x7f800000)
   61:                 return x+x;           /* NaN */
   62:             if(hx>0) return  atanhi[3]+atanlo[3];
   63:             else     return -atanhi[3]-atanlo[3];
   64:         } if (ix < 0x3ee00000) {       /* |x| < 0.4375 */
   65:             if (ix < 0x31000000) {     /* |x| < 2^-29 */
   66:                 if(huge+x>one) return x;      /* raise inexact */
   67:             }
   68:             id = -1;
   69:         } else {
   70:         x = fabsf(x);
   71:         if (ix < 0x3f980000) {         /* |x| < 1.1875 */
   72:             if (ix < 0x3f300000) {     /* 7/16 <=|x|<11/16 */
   73:                 id = 0; x = ((float)2.0*x-one)/((float)2.0+x); 
   74:             } else {                   /* 11/16<=|x|< 19/16 */
   75:                 id = 1; x  = (x-one)/(x+one); 
   76:             }
   77:         } else {
   78:             if (ix < 0x401c0000) {     /* |x| < 2.4375 */
   79:                 id = 2; x  = (x-(float)1.5)/(one+(float)1.5*x);
   80:             } else {                   /* 2.4375 <= |x| < 2^66 */
   81:                 id = 3; x  = -(float)1.0/x;
   82:             }
   83:         }}
   84:     /* end of argument reduction */
   85:         z = x*x;
   86:         w = z*z;
   87:     /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
   88:         s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
   89:         s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
   90:         if (id<0) return x - x*(s1+s2);
   91:         else {
   92:             z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
   93:             return (hx<0)? -z:z;
   94:         }
   95: }