kernel/utilities/
math.rs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
// Licensed under the Apache License, Version 2.0 or the MIT License.
// SPDX-License-Identifier: Apache-2.0 OR MIT
// Copyright Tock Contributors 2022.

//! Helper functions for common mathematical operations.

use core::f32;

/// Get closest power of two greater than the given number.
pub fn closest_power_of_two(mut num: u32) -> u32 {
    num -= 1;
    num |= num >> 1;
    num |= num >> 2;
    num |= num >> 4;
    num |= num >> 8;
    num |= num >> 16;
    num += 1;
    num
}

/// Represents an integral power-of-two as an exponent.
#[derive(Copy, Clone, Debug, PartialEq, PartialOrd, Eq, Ord)]
pub struct PowerOfTwo(u32);

impl PowerOfTwo {
    /// Returns the base-2 exponent as a numeric type
    pub fn exp<R>(self) -> R
    where
        R: From<u32>,
    {
        From::from(self.0)
    }

    /// Converts a number two the nearest [`PowerOfTwo`] less-than-or-equal to it.
    pub fn floor<F: Into<u32>>(f: F) -> PowerOfTwo {
        PowerOfTwo(log_base_two(f.into()))
    }

    /// Converts a number two the nearest [`PowerOfTwo`] greater-than-or-equal to
    /// it.
    pub fn ceiling<F: Into<u32>>(f: F) -> PowerOfTwo {
        PowerOfTwo(log_base_two(closest_power_of_two(f.into())))
    }

    /// Creates a new [`PowerOfTwo`] representing the number zero.
    pub fn zero() -> PowerOfTwo {
        PowerOfTwo(0)
    }

    /// Converts a [`PowerOfTwo`] to a number.
    pub fn as_num<F: From<u32>>(self) -> F {
        (1 << self.0).into()
    }
}

/// Get log base 2 of a number.
///
/// Note: this is the floor of the result. Also, an input of 0 results in an
/// output of 0.
pub fn log_base_two(num: u32) -> u32 {
    if num == 0 {
        0
    } else {
        31 - num.leading_zeros()
    }
}

/// Log base 2 of 64 bit unsigned integers.
pub fn log_base_two_u64(num: u64) -> u32 {
    if num == 0 {
        0
    } else {
        63 - num.leading_zeros()
    }
}

// f32 log10 function adapted from [micromath](https://github.com/NeoBirth/micromath)
const EXPONENT_MASK: u32 = 0b01111111_10000000_00000000_00000000;
const EXPONENT_BIAS: u32 = 127;

/// Return the absolute value of the floating point number.
pub fn abs(n: f32) -> f32 {
    f32::from_bits(n.to_bits() & 0x7FFF_FFFF)
}

fn extract_exponent_bits(x: f32) -> u32 {
    (x.to_bits() & EXPONENT_MASK).overflowing_shr(23).0
}

fn extract_exponent_value(x: f32) -> i32 {
    (extract_exponent_bits(x) as i32) - EXPONENT_BIAS as i32
}

fn ln_1to2_series_approximation(x: f32) -> f32 {
    // Idea from https://stackoverflow.com/a/44232045/. Modified to not be
    // restricted to int range and only values of x above 1.0 and got rid of
    // most of the slow conversions, should work for all positive values of x.

    // x may essentially be 1.0 but, as clippy notes, these kinds of floating
    // point comparisons can fail when the bit pattern is not the same.
    if abs(x - 1.0_f32) < f32::EPSILON {
        return 0.0_f32;
    }
    let x_less_than_1: bool = x < 1.0;
    // Note: we could use the fast inverse approximation here found in
    // super::inv::inv_approx, but the precision of such an approximation is
    // assumed not good enough.
    let x_working: f32 = if x_less_than_1 { 1.0 / x } else { x };
    // According to the SO post:
    // ln(x) = ln((2^n)*y)= ln(2^n) + ln(y) = ln(2) * n + ln(y)
    // Get exponent value.
    let base2_exponent: u32 = extract_exponent_value(x_working) as u32;
    let divisor: f32 = f32::from_bits(x_working.to_bits() & EXPONENT_MASK);
    // Supposedly normalizing between 1.0 and 2.0.
    let x_working: f32 = x_working / divisor;
    // Approximate polynomial generated from maple in the post using Remez
    // Algorithm: https://en.wikipedia.org/wiki/Remez_algorithm.
    let ln_1to2_polynomial: f32 = -1.741_793_9_f32
        + (2.821_202_6_f32
            + (-1.469_956_8_f32 + (0.447_179_55_f32 - 0.056_570_85_f32 * x_working) * x_working)
                * x_working)
            * x_working;
    // ln(2) * n + ln(y)
    let result: f32 = (base2_exponent as f32) * f32::consts::LN_2 + ln_1to2_polynomial;
    if x_less_than_1 {
        -result
    } else {
        result
    }
}

/// Compute the base 10 logarithm of `f`.
pub fn log10(x: f32) -> f32 {
    // Using change of base log10(x) = ln(x)/ln(10)
    let ln10_recip = f32::consts::LOG10_E;
    let fract_base_ln = ln10_recip;
    let value_ln = ln_1to2_series_approximation(x);
    value_ln * fract_base_ln
}