gonzui

    1: /* s_log1pf.c -- float version of s_log1p.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
   20: ln2_hi =   6.9313812256e-01,    /* 0x3f317180 */
   21: ln2_lo =   9.0580006145e-06,    /* 0x3717f7d1 */
   22: two25 =    3.355443200e+07,     /* 0x4c000000 */
   23: Lp1 = 6.6666668653e-01, /* 3F2AAAAB */
   24: Lp2 = 4.0000000596e-01, /* 3ECCCCCD */
   25: Lp3 = 2.8571429849e-01, /* 3E924925 */
   26: Lp4 = 2.2222198546e-01, /* 3E638E29 */
   27: Lp5 = 1.8183572590e-01, /* 3E3A3325 */
   28: Lp6 = 1.5313838422e-01, /* 3E1CD04F */
   29: Lp7 = 1.4798198640e-01; /* 3E178897 */
   30: 
   31: static const float zero = 0.0;
   32: 
   33: float
   34: log1pf(float x)
   35: {
   36:         float hfsq,f,c,s,z,R,u;
   37:         int32_t k,hx,hu,ax;
   38: 
   39:         GET_FLOAT_WORD(hx,x);
   40:         ax = hx&0x7fffffff;
   41: 
   42:         k = 1;
   43:         if (hx < 0x3ed413d7) {                 /* x < 0.41422  */
   44:             if(ax>=0x3f800000) {               /* x <= -1.0 */
   45:                 if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=+inf */
   46:                 else return (x-x)/(x-x);      /* log1p(x<-1)=NaN */
   47:             }
   48:             if(ax<0x31000000) {                        /* |x| < 2**-29 */
   49:                 if(two25+x>zero                       /* raise inexact */
   50:                     &&ax<0x24800000)           /* |x| < 2**-54 */
   51:                     return x;
   52:                 else
   53:                     return x - x*x*(float)0.5;
   54:             }
   55:             if(hx>0||hx<=((int32_t)0xbe95f61f)) {
   56:                 k=0;f=x;hu=1;}        /* -0.2929<x<0.41422 */
   57:         } 
   58:         if (hx >= 0x7f800000) return x+x;
   59:         if(k!=0) {
   60:             if(hx<0x5a000000) {
   61:                 u  = (float)1.0+x; 
   62:                 GET_FLOAT_WORD(hu,u);
   63:                 k  = (hu>>23)-127;
   64:                 /* correction term */
   65:                 c  = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0);
   66:                 c /= u;
   67:             } else {
   68:                 u  = x;
   69:                 GET_FLOAT_WORD(hu,u);
   70:                 k  = (hu>>23)-127;
   71:                 c  = 0;
   72:             }
   73:             hu &= 0x007fffff;
   74:             if(hu<0x3504f7) {
   75:                 SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */
   76:             } else {
   77:                 k += 1; 
   78:                 SET_FLOAT_WORD(u,hu|0x3f000000);      /* normalize u/2 */
   79:                 hu = (0x00800000-hu)>>2;
   80:             }
   81:             f = u-(float)1.0;
   82:         }
   83:         hfsq=(float)0.5*f*f;
   84:         if(hu==0) {    /* |f| < 2**-20 */
   85:             if(f==zero) if(k==0) return zero;  
   86:                         else {c += k*ln2_lo; return k*ln2_hi+c;}
   87:             R = hfsq*((float)1.0-(float)0.66666666666666666*f);
   88:             if(k==0) return f-R; else
   89:                     return k*ln2_hi-((R-(k*ln2_lo+c))-f);
   90:         }
   91:         s = f/((float)2.0+f); 
   92:         z = s*s;
   93:         R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
   94:         if(k==0) return f-(hfsq-s*(hfsq+R)); else
   95:                  return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
   96: }