/build/cargo-vendor-dir/libm-0.2.8/src/math/logf.rs
Line | Count | Source (jump to first uncovered line) |
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 | 0 | pub fn logf(mut x: f32) -> f32 { |
26 | 0 | let x1p25 = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25 |
27 | 0 |
|
28 | 0 | let mut ix = x.to_bits(); |
29 | 0 | let mut k = 0i32; |
30 | 0 |
|
31 | 0 | if (ix < 0x00800000) || ((ix >> 31) != 0) { |
32 | | /* x < 2**-126 */ |
33 | 0 | if ix << 1 == 0 { |
34 | 0 | return -1. / (x * x); /* log(+-0)=-inf */ |
35 | 0 | } |
36 | 0 | if (ix >> 31) != 0 { |
37 | 0 | return (x - x) / 0.; /* log(-#) = NaN */ |
38 | 0 | } |
39 | 0 | /* subnormal number, scale up x */ |
40 | 0 | k -= 25; |
41 | 0 | x *= x1p25; |
42 | 0 | ix = x.to_bits(); |
43 | 0 | } else if ix >= 0x7f800000 { |
44 | 0 | return x; |
45 | 0 | } else if ix == 0x3f800000 { |
46 | 0 | return 0.; |
47 | 0 | } |
48 | | |
49 | | /* reduce x into [sqrt(2)/2, sqrt(2)] */ |
50 | 0 | ix += 0x3f800000 - 0x3f3504f3; |
51 | 0 | k += ((ix >> 23) as i32) - 0x7f; |
52 | 0 | ix = (ix & 0x007fffff) + 0x3f3504f3; |
53 | 0 | x = f32::from_bits(ix); |
54 | 0 |
|
55 | 0 | let f = x - 1.; |
56 | 0 | let s = f / (2. + f); |
57 | 0 | let z = s * s; |
58 | 0 | let w = z * z; |
59 | 0 | let t1 = w * (LG2 + w * LG4); |
60 | 0 | let t2 = z * (LG1 + w * LG3); |
61 | 0 | let r = t2 + t1; |
62 | 0 | let hfsq = 0.5 * f * f; |
63 | 0 | let dk = k as f32; |
64 | 0 | s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI |
65 | 0 | } |