1 | /* origin: FreeBSD /usr/src/lib/msun/src/e_logf.c */ |
2 | /* |
3 | * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
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 | const LN2_HI: f32 = 6.9313812256e-01; /* 0x3f317180 */ |
17 | const LN2_LO: f32 = 9.0580006145e-06; /* 0x3717f7d1 */ |
18 | /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ |
19 | const LG1: f32 = 0.66666662693; /* 0xaaaaaa.0p-24*/ |
20 | const LG2: f32 = 0.40000972152; /* 0xccce13.0p-25 */ |
21 | const LG3: f32 = 0.28498786688; /* 0x91e9ee.0p-25 */ |
22 | const LG4: f32 = 0.24279078841; /* 0xf89e26.0p-26 */ |
23 | |
24 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
25 | pub fn logf(mut x: f32) -> f32 { |
26 | let x1p25 = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25 |
27 | |
28 | let mut ix = x.to_bits(); |
29 | let mut k = 0i32; |
30 | |
31 | if (ix < 0x00800000) || ((ix >> 31) != 0) { |
32 | /* x < 2**-126 */ |
33 | if ix << 1 == 0 { |
34 | return -1. / (x * x); /* log(+-0)=-inf */ |
35 | } |
36 | if (ix >> 31) != 0 { |
37 | return (x - x) / 0.; /* log(-#) = NaN */ |
38 | } |
39 | /* subnormal number, scale up x */ |
40 | k -= 25; |
41 | x *= x1p25; |
42 | ix = x.to_bits(); |
43 | } else if ix >= 0x7f800000 { |
44 | return x; |
45 | } else if ix == 0x3f800000 { |
46 | return 0.; |
47 | } |
48 | |
49 | /* reduce x into [sqrt(2)/2, sqrt(2)] */ |
50 | ix += 0x3f800000 - 0x3f3504f3; |
51 | k += ((ix >> 23) as i32) - 0x7f; |
52 | ix = (ix & 0x007fffff) + 0x3f3504f3; |
53 | x = f32::from_bits(ix); |
54 | |
55 | let f = x - 1.; |
56 | let s = f / (2. + f); |
57 | let z = s * s; |
58 | let w = z * z; |
59 | let t1 = w * (LG2 + w * LG4); |
60 | let t2 = z * (LG1 + w * LG3); |
61 | let r = t2 + t1; |
62 | let hfsq = 0.5 * f * f; |
63 | let dk = k as f32; |
64 | s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI |
65 | } |
66 | |