Skip to main content

machina_softfloat/ops/
mul.rs

1// SPDX-License-Identifier: MIT
2// IEEE 754 floating-point multiplication.
3
4use crate::env::FloatEnv;
5use crate::parts::{
6    nan_propagate, return_nan, round_pack, unpack, FloatClass, FloatParts,
7};
8use crate::types::{
9    BFloat16, Float128, Float16, Float32, Float64, FloatFormat, FloatX80,
10};
11
12/// Floating-point multiplication: a * b.
13pub fn mul<F: FloatFormat>(a: F, b: F, env: &mut FloatEnv) -> F {
14    let pa = unpack::<F>(a);
15    let pb = unpack::<F>(b);
16    mul_parts::<F>(&pa, &pb, env)
17}
18
19fn mul_parts<F: FloatFormat>(
20    a: &FloatParts,
21    b: &FloatParts,
22    env: &mut FloatEnv,
23) -> F {
24    let result_sign = a.sign ^ b.sign;
25
26    // NaN propagation
27    if a.is_nan() || b.is_nan() {
28        let mut r = nan_propagate(a, b, env);
29        r.sign = result_sign;
30        return round_pack::<F>(&mut r, env);
31    }
32
33    // Inf * 0 = NaN (INVALID)
34    if a.cls == FloatClass::Inf {
35        if b.cls == FloatClass::Zero {
36            return return_nan::<F>(env);
37        }
38        let mut r = FloatParts {
39            sign: result_sign,
40            exp: 0,
41            frac: 0,
42            cls: FloatClass::Inf,
43        };
44        return round_pack::<F>(&mut r, env);
45    }
46    if b.cls == FloatClass::Inf {
47        if a.cls == FloatClass::Zero {
48            return return_nan::<F>(env);
49        }
50        let mut r = FloatParts {
51            sign: result_sign,
52            exp: 0,
53            frac: 0,
54            cls: FloatClass::Inf,
55        };
56        return round_pack::<F>(&mut r, env);
57    }
58
59    // Zero handling
60    if a.cls == FloatClass::Zero || b.cls == FloatClass::Zero {
61        let mut r = FloatParts {
62            sign: result_sign,
63            exp: 0,
64            frac: 0,
65            cls: FloatClass::Zero,
66        };
67        return round_pack::<F>(&mut r, env);
68    }
69
70    // Both normal: multiply fractions and add exponents.
71    let exp = a.exp + b.exp;
72
73    // u128 * u128 -> upper 128 bits.
74    // Our fracs have the integer bit at position 126.
75    // Product's integer bits would be at position 252
76    // in a 256-bit result. We need the upper 128 bits.
77    let (hi, lo) = mul_u128(a.frac, b.frac);
78
79    // The product has the integer-bit pair at bit
80    // 252 of the 256-bit result = bit 124 of `hi`.
81    // We need to left-align so the integer bit is at 126.
82    // Shift left by 2, but carry in bits from lo.
83    let frac = (hi << 2) | (lo >> 126);
84    let sticky = if lo & ((1u128 << 126) - 1) != 0 { 1 } else { 0 };
85    let frac = frac | sticky;
86
87    let mut result = FloatParts {
88        sign: result_sign,
89        exp,
90        frac,
91        cls: FloatClass::Normal,
92    };
93    round_pack::<F>(&mut result, env)
94}
95
96/// Multiply two u128 values, return (hi, lo) of the 256-bit
97/// result. Uses 64-bit half decomposition.
98fn mul_u128(a: u128, b: u128) -> (u128, u128) {
99    let a_lo = a as u64 as u128;
100    let a_hi = (a >> 64) as u64 as u128;
101    let b_lo = b as u64 as u128;
102    let b_hi = (b >> 64) as u64 as u128;
103
104    let ll = a_lo * b_lo;
105    let lh = a_lo * b_hi;
106    let hl = a_hi * b_lo;
107    let hh = a_hi * b_hi;
108
109    // Accumulate the middle terms.
110    let mid = (ll >> 64)
111        + (lh & 0xFFFF_FFFF_FFFF_FFFF)
112        + (hl & 0xFFFF_FFFF_FFFF_FFFF);
113
114    let lo =
115        (ll & 0xFFFF_FFFF_FFFF_FFFF) | ((mid & 0xFFFF_FFFF_FFFF_FFFF) << 64);
116    let hi = hh + (lh >> 64) + (hl >> 64) + (mid >> 64);
117
118    (hi, lo)
119}
120
121// ---------------------------------------------------------------
122// Convenience methods
123// ---------------------------------------------------------------
124
125macro_rules! impl_mul {
126    ($ty:ty) => {
127        impl $ty {
128            pub fn mul(self, other: Self, env: &mut FloatEnv) -> Self {
129                mul::<Self>(self, other, env)
130            }
131        }
132    };
133}
134
135impl_mul!(Float16);
136impl_mul!(BFloat16);
137impl_mul!(Float32);
138impl_mul!(Float64);
139impl_mul!(Float128);
140impl_mul!(FloatX80);