1 | /* sf_log1p.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 "fdlibm.h" |
---|
17 | |
---|
18 | #ifdef __STDC__ |
---|
19 | static const float |
---|
20 | #else |
---|
21 | static float |
---|
22 | #endif |
---|
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 */ |
---|
33 | |
---|
34 | #ifdef __STDC__ |
---|
35 | static const float zero = 0.0; |
---|
36 | #else |
---|
37 | static float zero = 0.0; |
---|
38 | #endif |
---|
39 | |
---|
40 | #ifdef __STDC__ |
---|
41 | float log1pf(float x) |
---|
42 | #else |
---|
43 | float log1pf(x) |
---|
44 | float x; |
---|
45 | #endif |
---|
46 | { |
---|
47 | float hfsq,f,c,s,z,R,u; |
---|
48 | __int32_t k,hx,hu,ax; |
---|
49 | |
---|
50 | GET_FLOAT_WORD(hx,x); |
---|
51 | ax = hx&0x7fffffff; |
---|
52 | |
---|
53 | k = 1; |
---|
54 | if (!FLT_UWORD_IS_FINITE(hx)) return x+x; |
---|
55 | if (hx < 0x3ed413d7) { /* x < 0.41422 */ |
---|
56 | if(ax>=0x3f800000) { /* x <= -1.0 */ |
---|
57 | if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=+inf */ |
---|
58 | else return (x-x)/(x-x); /* log1p(x<-1)=NaN */ |
---|
59 | } |
---|
60 | if(ax<0x31000000) { /* |x| < 2**-29 */ |
---|
61 | if(two25+x>zero /* raise inexact */ |
---|
62 | &&ax<0x24800000) /* |x| < 2**-54 */ |
---|
63 | return x; |
---|
64 | else |
---|
65 | return x - x*x*(float)0.5; |
---|
66 | } |
---|
67 | if(hx>0||hx<=((__int32_t)0xbe95f61f)) { |
---|
68 | k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */ |
---|
69 | } |
---|
70 | if(k!=0) { |
---|
71 | if(hx<0x5a000000) { |
---|
72 | u = (float)1.0+x; |
---|
73 | GET_FLOAT_WORD(hu,u); |
---|
74 | k = (hu>>23)-127; |
---|
75 | /* correction term */ |
---|
76 | c = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0); |
---|
77 | c /= u; |
---|
78 | } else { |
---|
79 | u = x; |
---|
80 | GET_FLOAT_WORD(hu,u); |
---|
81 | k = (hu>>23)-127; |
---|
82 | c = 0; |
---|
83 | } |
---|
84 | hu &= 0x007fffff; |
---|
85 | if(hu<0x3504f7) { |
---|
86 | SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */ |
---|
87 | } else { |
---|
88 | k += 1; |
---|
89 | SET_FLOAT_WORD(u,hu|0x3f000000); /* normalize u/2 */ |
---|
90 | hu = (0x00800000-hu)>>2; |
---|
91 | } |
---|
92 | f = u-(float)1.0; |
---|
93 | } |
---|
94 | hfsq=(float)0.5*f*f; |
---|
95 | if(hu==0) { /* |f| < 2**-20 */ |
---|
96 | if(f==zero) { if(k==0) return zero; |
---|
97 | else {c += k*ln2_lo; return k*ln2_hi+c;}} |
---|
98 | R = hfsq*((float)1.0-(float)0.66666666666666666*f); |
---|
99 | if(k==0) return f-R; else |
---|
100 | return k*ln2_hi-((R-(k*ln2_lo+c))-f); |
---|
101 | } |
---|
102 | s = f/((float)2.0+f); |
---|
103 | z = s*s; |
---|
104 | R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); |
---|
105 | if(k==0) return f-(hfsq-s*(hfsq+R)); else |
---|
106 | return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f); |
---|
107 | } |
---|
108 | |
---|
109 | #ifdef _DOUBLE_IS_32BITS |
---|
110 | |
---|
111 | #ifdef __STDC__ |
---|
112 | double log1p(double x) |
---|
113 | #else |
---|
114 | double log1p(x) |
---|
115 | double x; |
---|
116 | #endif |
---|
117 | { |
---|
118 | return (double) log1pf((float) x); |
---|
119 | } |
---|
120 | |
---|
121 | #endif /* defined(_DOUBLE_IS_32BITS) */ |
---|