Back to index

glibc  2.9
e_logf.c
Go to the documentation of this file.
00001 /* e_logf.c -- float version of e_log.c.
00002  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
00003  */
00004 
00005 /*
00006  * ====================================================
00007  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00008  *
00009  * Developed at SunPro, a Sun Microsystems, Inc. business.
00010  * Permission to use, copy, modify, and distribute this
00011  * software is freely granted, provided that this notice
00012  * is preserved.
00013  * ====================================================
00014  */
00015 
00016 #if defined(LIBM_SCCS) && !defined(lint)
00017 static char rcsid[] = "$NetBSD: e_logf.c,v 1.4 1995/05/10 20:45:54 jtc Exp $";
00018 #endif
00019 
00020 #include "math.h"
00021 #include "math_private.h"
00022 
00023 #ifdef __STDC__
00024 static const float
00025 #else
00026 static float
00027 #endif
00028 ln2_hi =   6.9313812256e-01,       /* 0x3f317180 */
00029 ln2_lo =   9.0580006145e-06,       /* 0x3717f7d1 */
00030 two25 =    3.355443200e+07, /* 0x4c000000 */
00031 Lg1 = 6.6666668653e-01,     /* 3F2AAAAB */
00032 Lg2 = 4.0000000596e-01,     /* 3ECCCCCD */
00033 Lg3 = 2.8571429849e-01, /* 3E924925 */
00034 Lg4 = 2.2222198546e-01, /* 3E638E29 */
00035 Lg5 = 1.8183572590e-01, /* 3E3A3325 */
00036 Lg6 = 1.5313838422e-01, /* 3E1CD04F */
00037 Lg7 = 1.4798198640e-01; /* 3E178897 */
00038 
00039 #ifdef __STDC__
00040 static const float zero   =  0.0;
00041 #else
00042 static float zero   =  0.0;
00043 #endif
00044 
00045 #ifdef __STDC__
00046        float __ieee754_logf(float x)
00047 #else
00048        float __ieee754_logf(x)
00049        float x;
00050 #endif
00051 {
00052        float hfsq,f,s,z,R,w,t1,t2,dk;
00053        int32_t k,ix,i,j;
00054 
00055        GET_FLOAT_WORD(ix,x);
00056 
00057        k=0;
00058        if (ix < 0x00800000) {                    /* x < 2**-126  */
00059            if ((ix&0x7fffffff)==0)
00060               return -two25/(x-x);        /* log(+-0)=-inf */
00061            if (ix<0) return (x-x)/(x-x);  /* log(-#) = NaN */
00062            k -= 25; x *= two25; /* subnormal number, scale up x */
00063            GET_FLOAT_WORD(ix,x);
00064        }
00065        if (ix >= 0x7f800000) return x+x;
00066        k += (ix>>23)-127;
00067        ix &= 0x007fffff;
00068        i = (ix+(0x95f64<<3))&0x800000;
00069        SET_FLOAT_WORD(x,ix|(i^0x3f800000));      /* normalize x or x/2 */
00070        k += (i>>23);
00071        f = x-(float)1.0;
00072        if((0x007fffff&(15+ix))<16) {      /* |f| < 2**-20 */
00073            if(f==zero) {
00074              if(k==0) return zero;  else {dk=(float)k;
00075                                       return dk*ln2_hi+dk*ln2_lo;}
00076            }
00077            R = f*f*((float)0.5-(float)0.33333333333333333*f);
00078            if(k==0) return f-R; else {dk=(float)k;
00079                    return dk*ln2_hi-((R-dk*ln2_lo)-f);}
00080        }
00081        s = f/((float)2.0+f);
00082        dk = (float)k;
00083        z = s*s;
00084        i = ix-(0x6147a<<3);
00085        w = z*z;
00086        j = (0x6b851<<3)-ix;
00087        t1= w*(Lg2+w*(Lg4+w*Lg6));
00088        t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
00089        i |= j;
00090        R = t2+t1;
00091        if(i>0) {
00092            hfsq=(float)0.5*f*f;
00093            if(k==0) return f-(hfsq-s*(hfsq+R)); else
00094                    return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
00095        } else {
00096            if(k==0) return f-s*(f-R); else
00097                    return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
00098        }
00099 }