Back to index

glibc  2.9
s_log1pf.c
Go to the documentation of this file.
00001 /* s_log1pf.c -- float version of s_log1p.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: s_log1pf.c,v 1.4 1995/05/10 20:47:48 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 Lp1 = 6.6666668653e-01,     /* 3F2AAAAB */
00032 Lp2 = 4.0000000596e-01,     /* 3ECCCCCD */
00033 Lp3 = 2.8571429849e-01, /* 3E924925 */
00034 Lp4 = 2.2222198546e-01, /* 3E638E29 */
00035 Lp5 = 1.8183572590e-01, /* 3E3A3325 */
00036 Lp6 = 1.5313838422e-01, /* 3E1CD04F */
00037 Lp7 = 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 __log1pf(float x)
00047 #else
00048        float __log1pf(x)
00049        float x;
00050 #endif
00051 {
00052        float hfsq,f,c,s,z,R,u;
00053        int32_t k,hx,hu,ax;
00054 
00055        GET_FLOAT_WORD(hx,x);
00056        ax = hx&0x7fffffff;
00057 
00058        k = 1;
00059        if (hx < 0x3ed413d7) {                    /* x < 0.41422  */
00060            if(ax>=0x3f800000) {           /* x <= -1.0 */
00061               if(x==(float)-1.0) return -two25/(x-x); /* log1p(-1)=+inf */
00062               else return (x-x)/(x-x);    /* log1p(x<-1)=NaN */
00063            }
00064            if(ax<0x31000000) {                   /* |x| < 2**-29 */
00065               if(two25+x>zero                    /* raise inexact */
00066                    &&ax<0x24800000)              /* |x| < 2**-54 */
00067                   return x;
00068               else
00069                   return x - x*x*(float)0.5;
00070            }
00071            if(hx>0||hx<=((int32_t)0xbe95f61f)) {
00072               k=0;f=x;hu=1;}       /* -0.2929<x<0.41422 */
00073        }
00074        if (hx >= 0x7f800000) return x+x;
00075        if(k!=0) {
00076            if(hx<0x5a000000) {
00077               u  = (float)1.0+x;
00078               GET_FLOAT_WORD(hu,u);
00079                k  = (hu>>23)-127;
00080               /* correction term */
00081                c  = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0);
00082               c /= u;
00083            } else {
00084               u  = x;
00085               GET_FLOAT_WORD(hu,u);
00086                k  = (hu>>23)-127;
00087               c  = 0;
00088            }
00089            hu &= 0x007fffff;
00090            if(hu<0x3504f7) {
00091                SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */
00092            } else {
00093                k += 1;
00094               SET_FLOAT_WORD(u,hu|0x3f000000);   /* normalize u/2 */
00095                hu = (0x00800000-hu)>>2;
00096            }
00097            f = u-(float)1.0;
00098        }
00099        hfsq=(float)0.5*f*f;
00100        if(hu==0) {   /* |f| < 2**-20 */
00101            if(f==zero) {
00102               if(k==0) return zero;
00103               else {c += k*ln2_lo; return k*ln2_hi+c;}
00104            }
00105            R = hfsq*((float)1.0-(float)0.66666666666666666*f);
00106            if(k==0) return f-R; else
00107                    return k*ln2_hi-((R-(k*ln2_lo+c))-f);
00108        }
00109        s = f/((float)2.0+f);
00110        z = s*s;
00111        R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
00112        if(k==0) return f-(hfsq-s*(hfsq+R)); else
00113                return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
00114 }
00115 weak_alias (__log1pf, log1pf)