kernel/utilities/
math.rs

1// Licensed under the Apache License, Version 2.0 or the MIT License.
2// SPDX-License-Identifier: Apache-2.0 OR MIT
3// Copyright Tock Contributors 2022.
4
5//! Helper functions for common mathematical operations.
6
7/// Represents an integral power-of-two as an exponent.
8#[derive(Copy, Clone, Debug, PartialEq, PartialOrd, Eq, Ord)]
9pub struct PowerOfTwo(u32);
10
11impl PowerOfTwo {
12    /// Returns the base-2 exponent as a numeric type
13    pub fn exp<R>(self) -> R
14    where
15        R: From<u32>,
16    {
17        From::from(self.0)
18    }
19
20    /// Converts a number two the nearest [`PowerOfTwo`] less-than-or-equal to it.
21    pub fn floor<F: Into<u32>>(f: F) -> PowerOfTwo {
22        PowerOfTwo(log_base_two(f.into()))
23    }
24
25    /// Converts a number two the nearest [`PowerOfTwo`] greater-than-or-equal to
26    /// it.
27    pub fn ceiling<F: Into<u32>>(f: F) -> PowerOfTwo {
28        PowerOfTwo(log_base_two(f.into().next_power_of_two()))
29    }
30
31    /// Creates a new [`PowerOfTwo`] representing the number zero.
32    pub fn zero() -> PowerOfTwo {
33        PowerOfTwo(0)
34    }
35
36    /// Converts a [`PowerOfTwo`] to a number.
37    pub fn as_num<F: From<u32>>(self) -> F {
38        (1 << self.0).into()
39    }
40}
41
42/// Get log base 2 of a number.
43///
44/// Note: this is the floor of the result. Also, an input of 0 results in an
45/// output of 0.
46pub fn log_base_two(num: u32) -> u32 {
47    num.checked_ilog2().unwrap_or(0)
48}
49
50/// Log base 2 of 64 bit unsigned integers.
51pub fn log_base_two_u64(num: u64) -> u32 {
52    num.checked_ilog2().unwrap_or(0)
53}
54
55// f32 log10 function adapted from [micromath](https://github.com/NeoBirth/micromath)
56const EXPONENT_MASK: u32 = 0b01111111_10000000_00000000_00000000;
57const EXPONENT_BIAS: u32 = 127;
58
59/// Return the absolute value of the floating point number.
60pub fn abs(n: f32) -> f32 {
61    f32::from_bits(n.to_bits() & 0x7FFF_FFFF)
62}
63
64fn extract_exponent_bits(x: f32) -> u32 {
65    (x.to_bits() & EXPONENT_MASK).overflowing_shr(23).0
66}
67
68fn extract_exponent_value(x: f32) -> i32 {
69    (extract_exponent_bits(x) as i32) - EXPONENT_BIAS as i32
70}
71
72fn ln_1to2_series_approximation(x: f32) -> f32 {
73    // Idea from https://stackoverflow.com/a/44232045/. Modified to not be
74    // restricted to int range and only values of x above 1.0 and got rid of
75    // most of the slow conversions, should work for all positive values of x.
76
77    // x may essentially be 1.0 but, as clippy notes, these kinds of floating
78    // point comparisons can fail when the bit pattern is not the same.
79    if abs(x - 1.0_f32) < f32::EPSILON {
80        return 0.0_f32;
81    }
82    let x_less_than_1: bool = x < 1.0;
83    // Note: we could use the fast inverse approximation here found in
84    // super::inv::inv_approx, but the precision of such an approximation is
85    // assumed not good enough.
86    let x_working: f32 = if x_less_than_1 { 1.0 / x } else { x };
87    // According to the SO post:
88    // ln(x) = ln((2^n)*y)= ln(2^n) + ln(y) = ln(2) * n + ln(y)
89    // Get exponent value.
90    let base2_exponent: u32 = extract_exponent_value(x_working) as u32;
91    let divisor: f32 = f32::from_bits(x_working.to_bits() & EXPONENT_MASK);
92    // Supposedly normalizing between 1.0 and 2.0.
93    let x_working: f32 = x_working / divisor;
94    // Approximate polynomial generated from maple in the post using Remez
95    // Algorithm: https://en.wikipedia.org/wiki/Remez_algorithm.
96    let ln_1to2_polynomial: f32 = -1.741_793_9_f32
97        + (2.821_202_6_f32
98            + (-1.469_956_8_f32 + (0.447_179_55_f32 - 0.056_570_85_f32 * x_working) * x_working)
99                * x_working)
100            * x_working;
101    // ln(2) * n + ln(y)
102    let result: f32 = (base2_exponent as f32) * core::f32::consts::LN_2 + ln_1to2_polynomial;
103    if x_less_than_1 { -result } else { result }
104}
105
106/// Compute the base 10 logarithm of `f`.
107pub fn log10(x: f32) -> f32 {
108    // Using change of base log10(x) = ln(x)/ln(10)
109    let ln10_recip = core::f32::consts::LOG10_E;
110    let fract_base_ln = ln10_recip;
111    let value_ln = ln_1to2_series_approximation(x);
112    value_ln * fract_base_ln
113}
kernel/utilities/math.rs

kernel/utilities/
math.rs
