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

kernel/utilities/
math.rs
