brine_fp/log.rs
1use super::consts::*;
2use super::unsigned::UnsignedNumeric;
3use super::signed::SignedNumeric;
4
5// Modified from the original to support precise numbers instead of floats
6// The original C code, the long comment, and the constants
7// below are from FreeBSD's /usr/src/lib/msun/src/e_log.c
8// and came with this notice. The go code is a simpler
9// version of the original C.
10//
11// ====================================================
12// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
13//
14// Developed at SunPro, a Sun Microsystems, Inc. business.
15// Permission to use, copy, modify, and distribute this
16// software is freely granted, provided that this notice
17// is preserved.
18// ====================================================
19//
20// __ieee754_log(x)
21// Return the logarithm of x
22//
23// Method :
24// 1. Argument Reduction: find k and f such that
25// x = 2**k * (1+f),
26// where sqrt(2)/2 < 1+f < sqrt(2) .
27//
28// 2. Approximation of log(1+f).
29// Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
30// = 2s + 2/3 s**3 + 2/5 s**5 + .....,
31// = 2s + s*R
32// We use a special Reme algorithm on [0,0.1716] to generate
33// a polynomial of degree 14 to approximate R. The maximum error
34// of this polynomial approximation is bounded by 2**-58.45. In
35// other words,
36// 2 4 6 8 10 12 14
37// R(z) ~ L1*s +L2*s +L3*s +L4*s +L5*s +L6*s +L7*s
38// (the values of L1 to L7 are listed in the program) and
39// | 2 14 | -58.45
40// | L1*s +...+L7*s - R(z) | <= 2
41// | |
42// Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
43// In order to guarantee error in log below 1ulp, we compute log by
44// log(1+f) = f - s*(f - R) (if f is not too large)
45// log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
46//
47// 3. Finally, log(x) = k*Ln2 + log(1+f).
48// = k*Ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*Ln2_lo)))
49// Here Ln2 is split into two floating point number:
50// Ln2_hi + Ln2_lo,
51// where n*Ln2_hi is always exact for |n| < 2000.
52//
53// Special cases:
54// log(x) is NaN with signal if x < 0 (including -INF) ;
55// log(+INF) is +INF; log(0) is -INF with signal;
56// log(NaN) is that NaN with no signal.
57//
58// Accuracy:
59// according to an error analysis, the error is always less than
60// 1 ulp (unit in the last place).
61//
62// Constants:
63// The hexadecimal values are the intended ones for the following
64// constants. The decimal values may be used, provided that the
65// compiler will convert from decimal to binary accurately enough
66// to produce the hexadecimal values shown.
67// Frexp breaks f into a normalized fraction
68// and an integral power of two.
69// It returns frac and exp satisfying f == frac × 2**exp,
70// with the absolute value of frac in the interval [½, 1).
71
72
73impl UnsignedNumeric {
74
75 /// Log returns the natural logarithm of x.
76 ///
77 /// Special cases are:
78 /// Log(+Inf) = +Inf
79 /// Log(0) = -Inf
80 /// Log(x < 0) = NaN
81 pub fn log(&self) -> Option<SignedNumeric> {
82 if self.eq(&ZERO_PREC) {
83 return None;
84 }
85
86 if self.eq(&ONE_PREC) {
87 return Some(SignedNumeric {
88 value: ZERO_PREC.clone(),
89 is_negative: false,
90 });
91 }
92
93 let (f1_init, ki_init) = self.frexp()?;
94
95 let (f1, ki) = if f1_init.less_than(&SQRT2OVERTWO) {
96 let new_f1 = f1_init.checked_mul(&TWO_PREC)?;
97 let new_k1 = ki_init.checked_sub(1)?;
98 (new_f1, new_k1)
99 } else {
100 (f1_init, ki_init)
101 };
102
103 let f = f1.signed().checked_sub(&UnsignedNumeric::one().signed())?;
104
105 let s_divisor = UnsignedNumeric { value: two() }.signed().checked_add(&f)?;
106 let s = &f.checked_div(&s_divisor)?;
107 let s2 = s.checked_mul(s)?.value;
108 let s4 = s2.checked_mul(&s2)?;
109 // s2 * (L1 + s4*(L3+s4*(L5+s4*L7)))
110 let t1 = s2.checked_mul(&L1.checked_add(&s4.checked_mul(
111 &L3.checked_add(&s4.checked_mul(&L5.checked_add(&s4.checked_mul(&L7)?)?)?)?
112 )?)?)?;
113
114 // s4 * (L2 + s4*(L4+s4*L6))
115 let t2 = s4.checked_mul(
116 &L2.checked_add(&s4.checked_mul(&L4.checked_add(&s4.checked_mul(&L6)?)?)?)?,
117 )?;
118
119 let r = t1.checked_add(&t2)?;
120 let hfsq = f
121 .checked_mul(&f)?
122 .checked_div(&UnsignedNumeric { value: two() }.signed())?;
123 let k = SignedNumeric {
124 value: UnsignedNumeric::new(u128::try_from(ki.abs()).ok()?),
125 is_negative: ki < 0,
126 };
127
128 // k*Ln2Hi - ((hfsq - (s*(hfsq+R) + k*Ln2Lo)) - f)
129 let kl2hi = k
130 .checked_mul(&LN2HI.signed())?
131 .checked_div(&LN2HI_SCALE.signed())?;
132 let shfsqr = s.checked_mul(&hfsq.checked_add(&r.signed())?)?;
133 let kl2lo = k
134 .checked_mul(&LN2LO.signed())?
135 .checked_div(&LN2LO_SCALE.signed())?;
136
137 let shfsqr_kl2lo = shfsqr.checked_add(&kl2lo)?;
138 let hfsq_shfsqr_kl2lo = hfsq.checked_sub(&shfsqr_kl2lo)?;
139 let f_hfsq_shfsqr_kl2lo = hfsq_shfsqr_kl2lo.checked_sub(&f)?;
140
141 kl2hi.checked_sub(&f_hfsq_shfsqr_kl2lo)
142 }
143
144}
145
146#[cfg(test)]
147mod tests {
148 use super::*;
149 use crate::InnerUint;
150
151 #[test]
152 fn test_log() {
153 let precision = InnerUint::from(5_000_000_000_u128); // correct to at least 9 decimal places
154
155 let test = UnsignedNumeric::new(9);
156 let log = test.log().unwrap().value;
157 let expected = UnsignedNumeric::new(21972245773362196)
158 .checked_div(&UnsignedNumeric::new(10000000000000000))
159 .unwrap();
160 assert!(log.almost_eq(&expected, precision));
161
162 let test2 = UnsignedNumeric::new(2);
163 assert!(test2.log().unwrap().value.almost_eq(
164 &UnsignedNumeric::new(6931471805599453)
165 .checked_div(&UnsignedNumeric::new(10000000000000000))
166 .unwrap(),
167 precision
168 ));
169
170 let test3 = &UnsignedNumeric::new(12)
171 .checked_div(&UnsignedNumeric::new(10))
172 .unwrap();
173 assert!(test3.log().unwrap().value.almost_eq(
174 &UnsignedNumeric::new(1823215567939546)
175 .checked_div(&UnsignedNumeric::new(10000000000000000))
176 .unwrap(),
177 precision
178 ));
179
180 let test5 = &UnsignedNumeric::new(15)
181 .checked_div(&UnsignedNumeric::new(10))
182 .unwrap();
183 assert!(test5.log().unwrap().value.almost_eq(
184 &UnsignedNumeric::new(4054651081081644)
185 .checked_div(&UnsignedNumeric::new(10000000000000000))
186 .unwrap(),
187 precision
188 ));
189
190 let test6 = UnsignedNumeric::new(4)
191 .checked_div(&UnsignedNumeric::new(1000000))
192 .unwrap();
193 assert!(test6.log().unwrap().value.almost_eq(
194 &UnsignedNumeric::new(12429216196844383)
195 .checked_div(&UnsignedNumeric::new(1000000000000000))
196 .unwrap(),
197 precision
198 ));
199 }
200}