1 /* s_log1pf.c -- float version of s_log1p.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
13 * ====================================================
17 #include <math_private.h>
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 */
31 static const float zero
= 0.0;
36 float hfsq
,f
,c
,s
,z
,R
,u
;
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 */
48 if(ax
<0x31000000) { /* |x| < 2**-29 */
49 math_force_eval(two25
+x
); /* raise inexact */
50 if (ax
<0x24800000) /* |x| < 2**-54 */
53 return x
- x
*x
*(float)0.5;
55 if(hx
>0||hx
<=((int32_t)0xbe95f61f)) {
56 k
=0;f
=x
;hu
=1;} /* -0.2929<x<0.41422 */
58 if (hx
>= 0x7f800000) return x
+x
;
65 c
= (k
>0)? (float)1.0-(u
-x
):x
-(u
-(float)1.0);
75 SET_FLOAT_WORD(u
,hu
|0x3f800000);/* normalize u */
78 SET_FLOAT_WORD(u
,hu
|0x3f000000); /* normalize u/2 */
79 hu
= (0x00800000-hu
)>>2;
84 if(hu
==0) { /* |f| < 2**-20 */
87 else {c
+= k
*ln2_lo
; return k
*ln2_hi
+c
;}
89 R
= hfsq
*((float)1.0-(float)0.66666666666666666*f
);
90 if(k
==0) return f
-R
; else
91 return k
*ln2_hi
-((R
-(k
*ln2_lo
+c
))-f
);
95 R
= z
*(Lp1
+z
*(Lp2
+z
*(Lp3
+z
*(Lp4
+z
*(Lp5
+z
*(Lp6
+z
*Lp7
))))));
96 if(k
==0) return f
-(hfsq
-s
*(hfsq
+R
)); else
97 return k
*ln2_hi
-((hfsq
-(s
*(hfsq
+R
)+(k
*ln2_lo
+c
)))-f
);
99 weak_alias (__log1pf
, log1pf
)