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    if num == 0 {
50        0
51    } else {
52        31 - num.leading_zeros()
53    }
54}
55
56/// Log base 2 of 64 bit unsigned integers.
57pub fn log_base_two_u64(num: u64) -> u32 {
58    if num == 0 {
59        0
60    } else {
61        63 - num.leading_zeros()
62    }
63}
64
65// f32 log10 function adapted from [micromath](https://github.com/NeoBirth/micromath)
66const EXPONENT_MASK: u32 = 0b01111111_10000000_00000000_00000000;
67const EXPONENT_BIAS: u32 = 127;
68
69/// Return the absolute value of the floating point number.
70pub fn abs(n: f32) -> f32 {
71    f32::from_bits(n.to_bits() & 0x7FFF_FFFF)
72}
73
74fn extract_exponent_bits(x: f32) -> u32 {
75    (x.to_bits() & EXPONENT_MASK).overflowing_shr(23).0
76}
77
78fn extract_exponent_value(x: f32) -> i32 {
79    (extract_exponent_bits(x) as i32) - EXPONENT_BIAS as i32
80}
81
82fn ln_1to2_series_approximation(x: f32) -> f32 {
83    // Idea from https://stackoverflow.com/a/44232045/. Modified to not be
84    // restricted to int range and only values of x above 1.0 and got rid of
85    // most of the slow conversions, should work for all positive values of x.
86
87    // x may essentially be 1.0 but, as clippy notes, these kinds of floating
88    // point comparisons can fail when the bit pattern is not the same.
89    if abs(x - 1.0_f32) < f32::EPSILON {
90        return 0.0_f32;
91    }
92    let x_less_than_1: bool = x < 1.0;
93    // Note: we could use the fast inverse approximation here found in
94    // super::inv::inv_approx, but the precision of such an approximation is
95    // assumed not good enough.
96    let x_working: f32 = if x_less_than_1 { 1.0 / x } else { x };
97    // According to the SO post:
98    // ln(x) = ln((2^n)*y)= ln(2^n) + ln(y) = ln(2) * n + ln(y)
99    // Get exponent value.
100    let base2_exponent: u32 = extract_exponent_value(x_working) as u32;
101    let divisor: f32 = f32::from_bits(x_working.to_bits() & EXPONENT_MASK);
102    // Supposedly normalizing between 1.0 and 2.0.
103    let x_working: f32 = x_working / divisor;
104    // Approximate polynomial generated from maple in the post using Remez
105    // Algorithm: https://en.wikipedia.org/wiki/Remez_algorithm.
106    let ln_1to2_polynomial: f32 = -1.741_793_9_f32
107        + (2.821_202_6_f32
108            + (-1.469_956_8_f32 + (0.447_179_55_f32 - 0.056_570_85_f32 * x_working) * x_working)
109                * x_working)
110            * x_working;
111    // ln(2) * n + ln(y)
112    let result: f32 = (base2_exponent as f32) * f32::consts::LN_2 + ln_1to2_polynomial;
113    if x_less_than_1 {
114        -result
115    } else {
116        result
117    }
118}
119
120/// Compute the base 10 logarithm of `f`.
121pub fn log10(x: f32) -> f32 {
122    // Using change of base log10(x) = ln(x)/ln(10)
123    let ln10_recip = f32::consts::LOG10_E;
124    let fract_base_ln = ln10_recip;
125    let value_ln = ln_1to2_series_approximation(x);
126    value_ln * fract_base_ln
127}
kernel/utilities/math.rs

kernel/utilities/
math.rs
