machina_softfloat/ops/
div.rs

1// SPDX-License-Identifier: MIT
2// IEEE 754 floating-point division.
3
4use crate::env::{ExcFlags, 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 division: a / b.
13pub fn div<F: FloatFormat>(a: F, b: F, env: &mut FloatEnv) -> F {
14    let pa = unpack::<F>(a);
15    let pb = unpack::<F>(b);
16    div_parts::<F>(&pa, &pb, env)
17}
18
19fn div_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 / Inf = NaN (INVALID)
34    if a.cls == FloatClass::Inf {
35        if b.cls == FloatClass::Inf {
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
47    // x / Inf = 0
48    if b.cls == FloatClass::Inf {
49        let mut r = FloatParts {
50            sign: result_sign,
51            exp: 0,
52            frac: 0,
53            cls: FloatClass::Zero,
54        };
55        return round_pack::<F>(&mut r, env);
56    }
57
58    // 0 / 0 = NaN (INVALID)
59    if a.cls == FloatClass::Zero {
60        if b.cls == FloatClass::Zero {
61            return return_nan::<F>(env);
62        }
63        let mut r = FloatParts {
64            sign: result_sign,
65            exp: 0,
66            frac: 0,
67            cls: FloatClass::Zero,
68        };
69        return round_pack::<F>(&mut r, env);
70    }
71
72    // x / 0 = Inf (DIVBYZERO)
73    if b.cls == FloatClass::Zero {
74        env.raise(ExcFlags::DIVBYZERO);
75        let mut r = FloatParts {
76            sign: result_sign,
77            exp: 0,
78            frac: 0,
79            cls: FloatClass::Inf,
80        };
81        return round_pack::<F>(&mut r, env);
82    }
83
84    // Both normal: divide fractions, subtract exponents.
85    let exp = a.exp - b.exp;
86
87    // We need a.frac / b.frac with enough precision.
88    // Both fracs have integer bit at position 126.
89    //
90    // Strategy: use iterative long division via shifting.
91    // We want ~128 bits of quotient. Since both operands
92    // are ~127-bit numbers (bit 126 set), the quotient
93    // is close to 1.xxx (at most 2.xxx if a.frac >= b.frac
94    // before any normalization).
95    //
96    // Shift the dividend left to get more quotient bits.
97    // dividend = a.frac, divisor = b.frac.
98    // If a.frac >= b.frac, quotient integer bit is 1 and
99    // exp stays. Otherwise, shift dividend left by 1 and
100    // decrement exp.
101
102    let (frac, exp) = div_frac(a.frac, b.frac, exp);
103
104    let mut result = FloatParts {
105        sign: result_sign,
106        exp,
107        frac,
108        cls: FloatClass::Normal,
109    };
110    round_pack::<F>(&mut result, env)
111}
112
113/// Divide two u128 fractions (integer bit at position 126).
114/// Returns (quotient_frac, adjusted_exp).
115///
116/// Algorithm: ensure dividend >= divisor (adjust exp),
117/// then iterative trial subtraction to produce ~128 bits
118/// of quotient with integer bit at position 126.
119fn div_frac(a_frac: u128, b_frac: u128, mut exp: i32) -> (u128, i32) {
120    if b_frac == 0 {
121        return (0, exp);
122    }
123
124    let mut rem = a_frac;
125
126    // If a_frac < b_frac, shift numerator left by 1
127    // and decrement exponent so quotient >= 1.
128    if rem < b_frac {
129        rem <<= 1;
130        exp -= 1;
131    }
132
133    // Now rem >= b_frac. The integer part of the
134    // quotient is 1. Produce 127 fractional bits
135    // via trial subtraction (total: 128 bits with
136    // integer bit at position 127).
137    let mut q: u128 = 0;
138    for _ in 0..128 {
139        q <<= 1;
140        if rem >= b_frac {
141            q |= 1;
142            rem -= b_frac;
143        }
144        rem <<= 1;
145    }
146
147    // Sticky bit from remaining remainder.
148    if rem != 0 {
149        q |= 1;
150    }
151
152    // q has integer bit at position 127. Shift right
153    // by 1 to place it at position 126 (our convention).
154    let sticky = q & 1;
155    let frac = (q >> 1) | sticky;
156
157    (frac, exp)
158}
159
160// ---------------------------------------------------------------
161// Convenience methods
162// ---------------------------------------------------------------
163
164macro_rules! impl_div {
165    ($ty:ty) => {
166        impl $ty {
167            pub fn div(self, other: Self, env: &mut FloatEnv) -> Self {
168                div::<Self>(self, other, env)
169            }
170        }
171    };
172}
173
174impl_div!(Float16);
175impl_div!(BFloat16);
176impl_div!(Float32);
177impl_div!(Float64);
178impl_div!(Float128);
179impl_div!(FloatX80);
machina_softfloat/ops/div.rs

machina_softfloat/ops/
div.rs
