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

kernel/utilities/
math.rs
