]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/flt-32/s_log1pf.c
1 /* s_log1pf.c -- float version of s_log1p.c.
5 * ====================================================
6 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 * Developed at SunPro, a Sun Microsystems, Inc. business.
9 * Permission to use, copy, modify, and distribute this
10 * software is freely granted, provided that this notice
12 * ====================================================
17 #include <math-barriers.h>
18 #include <math_private.h>
19 #include <math-underflow.h>
20 #include <libc-diag.h>
23 ln2_hi
= 6.9313812256e-01, /* 0x3f317180 */
24 ln2_lo
= 9.0580006145e-06, /* 0x3717f7d1 */
25 two25
= 3.355443200e+07, /* 0x4c000000 */
26 Lp1
= 6.6666668653e-01, /* 3F2AAAAB */
27 Lp2
= 4.0000000596e-01, /* 3ECCCCCD */
28 Lp3
= 2.8571429849e-01, /* 3E924925 */
29 Lp4
= 2.2222198546e-01, /* 3E638E29 */
30 Lp5
= 1.8183572590e-01, /* 3E3A3325 */
31 Lp6
= 1.5313838422e-01, /* 3E1CD04F */
32 Lp7
= 1.4798198640e-01; /* 3E178897 */
34 static const float zero
= 0.0;
39 float hfsq
,f
,c
,s
,z
,R
,u
;
46 if (hx
< 0x3ed413d7) { /* x < 0.41422 */
47 if(ax
>=0x3f800000) { /* x <= -1.0 */
48 if(x
==(float)-1.0) return -two25
/zero
; /* log1p(-1)=-inf */
49 else return (x
-x
)/(x
-x
); /* log1p(x<-1)=NaN */
51 if(ax
<0x31000000) { /* |x| < 2**-29 */
52 math_force_eval(two25
+x
); /* raise inexact */
53 if (ax
<0x24800000) /* |x| < 2**-54 */
55 math_check_force_underflow (x
);
59 return x
- x
*x
*(float)0.5;
61 if(hx
>0||hx
<=((int32_t)0xbe95f61f)) {
62 k
=0;f
=x
;hu
=1;} /* -0.2929<x<0.41422 */
64 if (hx
>= 0x7f800000) return x
+x
;
71 c
= (k
>0)? (float)1.0-(u
-x
):x
-(u
-(float)1.0);
81 SET_FLOAT_WORD(u
,hu
|0x3f800000);/* normalize u */
84 SET_FLOAT_WORD(u
,hu
|0x3f000000); /* normalize u/2 */
85 hu
= (0x00800000-hu
)>>2;
90 if(hu
==0) { /* |f| < 2**-20 */
93 else {c
+= k
*ln2_lo
; return k
*ln2_hi
+c
;}
95 R
= hfsq
*((float)1.0-(float)0.66666666666666666*f
);
96 if(k
==0) return f
-R
; else
97 return k
*ln2_hi
-((R
-(k
*ln2_lo
+c
))-f
);
101 R
= z
*(Lp1
+z
*(Lp2
+z
*(Lp3
+z
*(Lp4
+z
*(Lp5
+z
*(Lp6
+z
*Lp7
))))));
103 return f
- (hfsq
- s
* (hfsq
+ R
));
106 /* With GCC 7 when compiling with -Os the compiler warns
107 that c might be used uninitialized. This can't be true
108 because k must be 0 for c to be uninitialized and we
109 handled that computation earlier without using c. */
110 DIAG_PUSH_NEEDS_COMMENT
;
111 DIAG_IGNORE_Os_NEEDS_COMMENT (7, "-Wmaybe-uninitialized");
112 return k
* ln2_hi
- ((hfsq
- (s
* (hfsq
+ R
)
113 + (k
* ln2_lo
+ c
))) - f
);
114 DIAG_POP_NEEDS_COMMENT
;