spirix 0.1.1

Two's complement floating-point arithmetic library
Documentation
use crate::core::integer::*;
use crate::core::undefined::*;
use crate::{Integer, Scalar, ScalarConstants};
use core::ops::*;
use i256::I256;
use num_traits::{AsPrimitive, WrappingAdd, WrappingMul, WrappingNeg, WrappingSub};
#[allow(private_bounds)]
impl<
        F: Integer
            + FullInt
            + Shl<isize, Output = F>
            + Shr<isize, Output = F>
            + Shl<F, Output = F>
            + Shr<F, Output = F>
            + Shl<E, Output = F>
            + Shr<E, Output = F>
            + WrappingNeg
            + WrappingAdd
            + WrappingMul
            + WrappingSub,
        E: Integer
            + FullInt
            + Shl<isize, Output = E>
            + Shr<isize, Output = E>
            + Shl<E, Output = E>
            + Shr<E, Output = E>
            + Shl<F, Output = E>
            + Shr<F, Output = E>
            + WrappingNeg
            + WrappingAdd
            + WrappingMul
            + WrappingSub,
    > Scalar<F, E>
where
    Scalar<F, E>: ScalarConstants,
    u8: AsPrimitive<F>,
    u16: AsPrimitive<F>,
    u32: AsPrimitive<F>,
    u64: AsPrimitive<F>,
    u128: AsPrimitive<F>,
    usize: AsPrimitive<F>,
    i8: AsPrimitive<F>,
    i16: AsPrimitive<F>,
    i32: AsPrimitive<F>,
    i64: AsPrimitive<F>,
    i128: AsPrimitive<F>,
    isize: AsPrimitive<F>,
    I256: From<F>,
    u8: AsPrimitive<E>,
    u16: AsPrimitive<E>,
    u32: AsPrimitive<E>,
    u64: AsPrimitive<E>,
    u128: AsPrimitive<E>,
    usize: AsPrimitive<E>,
    i8: AsPrimitive<E>,
    i16: AsPrimitive<E>,
    i32: AsPrimitive<E>,
    i64: AsPrimitive<E>,
    i128: AsPrimitive<E>,
    isize: AsPrimitive<E>,
    I256: From<E>,
{
    /// Subtracts another Scalar from this Scalar
    ///
    /// # Description
    ///
    /// Performs subtraction between two Scalars, handling special cases according to mathematical principles. Returns a finite Scalar unless the result exceeds representable range, in which case it may return an exploded or vanished Scalar.
    ///
    /// Subtraction process:
    /// 0. Checks for any escaped Scalars (vanished, exploded or undefined) and handles these cases
    /// 1. Checks for Zeros and returns the appropriate result
    /// 2. Aligns fractions by shifting the smaller value right based on exponent difference
    /// 3. Subtracts the aligned values
    /// 4. Normalizes the result and adjusts exponent accordingly
    /// 5. Escapes for underflow if necessary, overflow is naturally handled by escaped exponent alignment
    ///
    /// # Returns
    ///
    /// - `[℘ ]` ➔ `[℘ ]` First undefined state encountered
    /// - `[↑]` - `[↑]` ➔ `[℘ ↑-↑]` Undefined exploded minus exploded state
    /// - `[↑]` - `[#]` ➔ `[℘ ↑-]` Undefined exploded minus finite state
    /// - `[#]` - `[↑]` ➔ `[℘ -↑]` Undefined finite minus exploded state
    /// - `[↓]` - `[↓]` ➔ `[℘ ↓-↓]` Undefined vanished minus vanished state
    /// - `[↓]` - `[#]` ➔ `[-#]` Negative of the finite Scalar
    /// - `[#]` - `[↓]` ➔ `[#]` The finite Scalar
    /// - `[0]` - `[#]` ➔ `[-#]` Negative of the Scalar
    /// - `[#]` - `[0]` ➔ `[#]` The Scalar
    /// - `[0]` - `[0]` ➔ `[0]` Zero
    /// - `[#]` - `[#]` ➔ `[#]` or `[↑]` or `[↓]` A finite, exploded or vanished Scalar
    ///
    /// # Examples
    ///
    /// ```rust
    /// use spirix::{Scalar, ScalarF5E3};
    ///
    /// // Subtracting finite Scalars let a = Scalar::<i32, i8>::from(42_i32); let b = ScalarF5E3::from(12_i32); let diff = a - b; assert!(diff == 30_i32);
    ///
    /// // Subtracting with Zero assert!(a - 0_i32 == a); assert!(0_i32 - a == -a);
    ///
    /// // Subtraction that produces Zero let pos = ScalarF5E3::from(7_i32); assert!((pos - pos).is_zero());
    ///
    /// // Subtraction with vanished Scalars let tiny: ScalarF5E3 = ScalarF5E3::MIN_POS / 4_i32; assert!(tiny.vanished()); assert!(a - tiny == a); // Vanished value treated as Zero
    ///
    /// // Subtraction with exploded Scalars let huge: ScalarF5E3 = ScalarF5E3::MAX * 4_i32; assert!(huge.exploded()); assert!((a - huge).is_undefined());
    ///
    /// // Subtracting two exploded Scalars assert!((huge - huge).is_undefined());
    /// ```
    pub(crate) fn scalar_subtract_scalar(&self, scalar: &Self) -> Self {
        if !self.is_normal() || !scalar.is_normal() {
            if self.is_undefined() {
                return *self;
            }
            if scalar.is_undefined() {
                return *scalar;
            }
            // Infinity absorbs everything. [∞] is signless so −[∞] is a no-op; [∞]−X and X−[∞] both yield [∞].
            if self.is_infinite() || scalar.is_infinite() {
                return Self::INFINITY;
            }
            // Zero identity: X−[0] = X and [0]−X = −X.
            if scalar.is_zero() {
                return *self;
            }
            if self.is_zero() {
                return -scalar;
            }
            if self.exploded() && scalar.exploded() {
                return Self {
                    fraction: TRANSFINITE_MINUS_TRANSFINITE.prefix.sa(),
                    exponent: Self::ambiguous_exponent(),
                };
            }
            if self.vanished() && scalar.vanished() {
                return Self {
                    fraction: VANISHED_MINUS_VANISHED.prefix.sa(),
                    exponent: Self::ambiguous_exponent(),
                };
            }
            if self.exploded() {
                if scalar.vanished() {
                    return *self;
                }
                return Self {
                    fraction: TRANSFINITE_MINUS_FINITE.prefix.sa(),
                    exponent: Self::ambiguous_exponent(),
                };
            }
            if scalar.exploded() {
                if self.vanished() {
                    return -scalar;
                }
                return Self {
                    fraction: FINITE_MINUS_TRANSFINITE.prefix.sa(),
                    exponent: Self::ambiguous_exponent(),
                };
            }
            if self.vanished() {
                return -scalar;
            }
            if scalar.vanished() {
                return *self;
            }
            return *self;
        }

        let big_is_self = self.exponent.into_unsigned() > scalar.exponent.into_unsigned();
        let (big, small) = if big_is_self {
            (self, scalar)
        } else {
            (scalar, self)
        };
        let exp_diff = big.exponent.wrapping_sub(&small.exponent);
        // x86 implementation note: Spirix's signless design wants to drop small only when shift ≥ FRAC. But Wide = 2*FRAC bits, and the inflated form reaches 2^FRAC magnitude at the ±1.0 boundary, so big_f<<shift ± small_f needs 2*FRAC+2 bits — one more than Rust/x86 provides at any power-of-2 width. Tightening the threshold to FRAC-1 keeps the sub in signed 2*FRAC-bit arithmetic with no sign branching. Cost: ≤1 ULP when shift == FRAC-1. In Verilog this tightening is unnecessary.
        let small_lost = if big_is_self { *self } else { -scalar };
        if exp_diff.is_negative() {
            return small_lost;
        }
        let shift: isize = exp_diff.saturate();
        if shift >= Self::fraction_bits().wrapping_sub(1) {
            return small_lost;
        }

        let big_f = big.fraction.inflate(true).w_shl(shift);
        let small_f = small.fraction.inflate(true);
        // Pure two's complement subtract. With shift<=FRAC-2 both operands fit with room to spare, so signed wrap cannot occur regardless of sign.
        let result = if big_is_self {
            big_f.w_sub(small_f)
        } else {
            small_f.w_sub(big_f)
        };
        if result.w_is_zero() {
            return Self {
                fraction: F::zero(),
                exponent: Self::ambiguous_exponent(),
            };
        }
        let leading = result.leading_same();
        let fb = Self::fraction_bits();
        let delta: isize = fb.wrapping_sub(leading);
        // AMBIG=0 native: view small.exponent as an unsigned cycle position via `cycle_widen` (2x widening into <E as Inflate>::Wide), add the normalization delta (sign-extended since it can be negative), and bounds-check against the cycle's normal range [min_pos, max_pos].
        let small_pos = small.exponent.cycle_widen();
        let delta_e: E = delta.as_();
        let w_delta = delta_e.sign_extend();
        let offset_pos = small_pos.w_add(w_delta);
        let max_pos = Self::max_exponent().cycle_widen();
        let min_pos = Self::min_exponent().cycle_widen();
        if offset_pos > max_pos {
            return Self {
                fraction: result.w_shl(leading.wrapping_sub(1)).w_shr(fb).deflate(),
                exponent: Self::ambiguous_exponent(),
            };
        }
        if offset_pos < min_pos {
            return Self {
                fraction: result.w_shl(leading.wrapping_sub(2)).w_shr(fb).deflate(),
                exponent: Self::ambiguous_exponent(),
            };
        }
        let offset: E = offset_pos.deflate();
        // Main path extraction: `result << L >> FRAC` composed as a net shift of (L - FRAC). Writing it directly avoids the Rust shift-overflow semantics that bite when L == wide_bits (result = -1, full sign extension).
        let shl_amount = leading.wrapping_sub(fb);
        let canonical = if shl_amount >= 0 {
            result.w_shl(shl_amount)
        } else {
            result.w_shr(shl_amount.wrapping_neg())
        };
        Self {
            fraction: canonical.deflate(),
            exponent: offset,
        }
    }
}