la-stack 0.4.6

Fast, stack-allocated linear algebra for fixed dimensions
Documentation 
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

#![forbid(unsafe_code)]

//! Shared binary64 rounding primitives for certified arithmetic.
//!
//! Representation and rounding use the IEEE 754 model in `REFERENCES.md`
//! \[9-10\]. The interval and reduction certificates share these primitives.

/// Return the exact error in a rounded binary64 sum.
///
/// This is `FastTwoSum` with operands ordered by magnitude. With IEEE-754
/// round-to-nearest and gradual underflow, `rounded + error` equals the
/// exact-real sum whenever the rounded sum is finite. Ordering prevents an
/// intermediate overflow even at the finite range boundary; see
/// `REFERENCES.md` \[17\], Theorem 5.1. Callers supply finite operands and their
/// rounded sum.
#[inline]
pub(crate) const fn two_sum_error(left: f64, right: f64, rounded: f64) -> f64 {
    let (large, small) = if left.abs() >= right.abs() {
        (left, right)
    } else {
        (right, left)
    };
    let virtual_small = rounded - large;
    small - virtual_small
}

/// Decompose a nonzero finite binary64 magnitude as `significand × 2^exponent`.
#[inline]
pub(crate) const fn decompose_magnitude(value: f64) -> (u128, i64) {
    let magnitude_bits = value.to_bits() & 0x7fff_ffff_ffff_ffff;
    let biased_exponent = ((magnitude_bits >> 52) & 0x7ff).cast_signed();
    let fraction = magnitude_bits & 0x000f_ffff_ffff_ffff;

    if biased_exponent == 0 {
        (fraction as u128, -1074)
    } else {
        (
            (fraction | (1_u64 << 52)) as u128,
            biased_exponent - 1023 - 52,
        )
    }
}

/// Compare two positive values represented as `significand × 2^exponent`.
#[inline]
const fn compare_binary_magnitudes(
    left_significand: u128,
    left_exponent: i64,
    right_significand: u128,
    right_exponent: i64,
) -> i8 {
    let left_zeros = left_significand.trailing_zeros() as i64;
    let right_zeros = right_significand.trailing_zeros() as i64;
    let normalized_left = left_significand >> left_zeros.cast_unsigned();
    let normalized_right = right_significand >> right_zeros.cast_unsigned();
    let normalized_left_exponent = left_exponent + left_zeros;
    let normalized_right_exponent = right_exponent + right_zeros;

    let left_top =
        normalized_left_exponent + (u128::BITS - normalized_left.leading_zeros() - 1) as i64;
    let right_top =
        normalized_right_exponent + (u128::BITS - normalized_right.leading_zeros() - 1) as i64;
    if left_top < right_top {
        return -1;
    }
    if left_top > right_top {
        return 1;
    }

    let common_exponent = if normalized_left_exponent < normalized_right_exponent {
        normalized_left_exponent
    } else {
        normalized_right_exponent
    };
    let aligned_left =
        normalized_left << (normalized_left_exponent - common_exponent).cast_unsigned();
    let aligned_right =
        normalized_right << (normalized_right_exponent - common_exponent).cast_unsigned();
    if aligned_left < aligned_right {
        -1
    } else if aligned_left > aligned_right {
        1
    } else {
        0
    }
}

/// Compare the exact-real product `left × right` with its rounded result.
///
/// Nonzero finite operands contribute at most 53 significand bits each, so the
/// exact product fits in `u128`. Integer comparison selects the outward
/// endpoint even when the rounded product is zero; see `REFERENCES.md` \[9-10\].
#[inline]
pub(crate) const fn compare_product_with_rounded(left: f64, right: f64, rounded: f64) -> i8 {
    let negative = left.is_sign_negative() != right.is_sign_negative();
    if rounded == 0.0 {
        return if negative { -1 } else { 1 };
    }

    let (left_significand, left_exponent) = decompose_magnitude(left);
    let (right_significand, right_exponent) = decompose_magnitude(right);
    let exact_significand = left_significand * right_significand;
    let exact_exponent = left_exponent + right_exponent;
    let (rounded_significand, rounded_exponent) = decompose_magnitude(rounded);
    let magnitude_relation = compare_binary_magnitudes(
        exact_significand,
        exact_exponent,
        rounded_significand,
        rounded_exponent,
    );

    if negative {
        -magnitude_relation
    } else {
        magnitude_relation
    }
}

#[cfg(test)]
mod tests {
    use super::compare_binary_magnitudes;

    #[test]
    fn binary_magnitude_comparison_orders_distinct_top_exponents() {
        assert_eq!(compare_binary_magnitudes(1, 1, 1, 0), 1);
        assert_eq!(compare_binary_magnitudes(1, 0, 1, 1), -1);
    }
}