Just.PreciseMath/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs

namespace Just.PreciseMath;

public readonly partial struct DoubleDouble :
    IAdditiveIdentity<DoubleDouble, DoubleDouble>,
    IMultiplicativeIdentity<DoubleDouble, DoubleDouble>,
    IUnaryPlusOperators<DoubleDouble, DoubleDouble>,
    IUnaryNegationOperators<DoubleDouble, DoubleDouble>,
    IAdditionOperators<DoubleDouble, DoubleDouble, DoubleDouble>,
    IAdditionOperators<DoubleDouble, double, DoubleDouble>,
    ISubtractionOperators<DoubleDouble, DoubleDouble, DoubleDouble>,
    ISubtractionOperators<DoubleDouble, double, DoubleDouble>,
    IMultiplyOperators<DoubleDouble, DoubleDouble, DoubleDouble>,
    IMultiplyOperators<DoubleDouble, double, DoubleDouble>,
    IDivisionOperators<DoubleDouble, DoubleDouble, DoubleDouble>,
    IDivisionOperators<DoubleDouble, double, DoubleDouble>
{
    /// <summary>Returns the operand unchanged.</summary>
    public static DoubleDouble operator +(DoubleDouble value)
    {
        return value;
    }

    /// <summary>Negates the value, including the high zero's sign, retaining canonical NaN and zero residuals.</summary>
    public static DoubleDouble operator -(DoubleDouble value)
    {
        // Negation preserves normalization; only NaN and zero residuals need canonicalization.
        return new DoubleDouble(double.IsNaN(value._high) ? double.NaN : -value._high,
            value._low == 0.0 ? 0.0 : -value._low);
    }

    /// <summary>Adds normalized expansions, retaining low-sum residuals under cancellation.</summary>
    public static DoubleDouble operator +(DoubleDouble left, DoubleDouble right)
    {
        if (!IsFinite(left) || !IsFinite(right) || (left._high == 0.0 && right._high == 0.0))
        {
            return new DoubleDouble(left._high + right._high);
        }
        if (!PreciseMathHelper.IsAdditionWithinFastRange(left._high) || !PreciseMathHelper.IsAdditionWithinFastRange(right._high))
        {
            return PreciseMathHelper.AddBoundary(left, right);
        }

        return PreciseMathHelper.AddFinite(left._high, left._low, right._high, right._low);
    }

    /// <summary>Subtracts normalized expansions.</summary>
    public static DoubleDouble operator -(DoubleDouble left, DoubleDouble right)
    {
        if (!IsFinite(left) || !IsFinite(right) || (left._high == 0.0 && right._high == 0.0))
        {
            return new DoubleDouble(left._high - right._high);
        }
        if (!PreciseMathHelper.IsAdditionWithinFastRange(left._high) || !PreciseMathHelper.IsAdditionWithinFastRange(right._high))
        {
            return PreciseMathHelper.AddBoundary(left, -right);
        }

        // With a canonical left low (never -0), the low TwoSum absorbs the
        // negated right zero low without changing either output component's bits.
        // Avoid unary negation's intermediate canonicalization on this finite path.
        return PreciseMathHelper.AddFinite(left._high, left._low, -right._high, -right._low);
    }

    /// <summary>Multiplies expansions using an FMA product residual and cross terms.</summary>
    /// <remarks>Results are approximate double-double values, not universally correctly rounded.</remarks>
    public static DoubleDouble operator *(DoubleDouble left, DoubleDouble right)
    {
        if (!IsFinite(left) || !IsFinite(right) || left._high == 0.0 || right._high == 0.0)
        {
            return new DoubleDouble(left._high * right._high);
        }
        if (!PreciseMathHelper.IsMultiplicationWithinFastRange(left._high, right._high))
        {
            return PreciseMathHelper.MultiplyBoundary(left, right);
        }

        (double product, double error) = PreciseMathHelper.TwoMultiply(left._high, right._high);
        error = Math.FusedMultiplyAdd(left._high, right._low, error);
        error = Math.FusedMultiplyAdd(left._low, right._high, error);
        error = Math.FusedMultiplyAdd(left._low, right._low, error);
        // With u = 2^-53, normalized inputs and the exponent-sum guard give
        // |error| < 4u*|product|, including rounding at the subnormal floor.
        // The product is normal and nonzero; its corrected sum remains finite.
        (double high, double low) = PreciseMathHelper.TwoQuickAdd(product, error);
        return new DoubleDouble(high, low == 0.0 ? 0.0 : low);
    }

    /// <summary>Divides expansions using a quotient estimate and two residual corrections.</summary>
    /// <remarks>Zero and nonfinite operands follow binary64 rules; precision decreases near underflow.</remarks>
    public static DoubleDouble operator /(DoubleDouble left, DoubleDouble right)
    {
        if (!IsFinite(left) || !IsFinite(right) || left._high == 0.0 || right._high == 0.0)
        {
            return new DoubleDouble(left._high / right._high);
        }
        if (!PreciseMathHelper.IsDivisionWithinFastRange(left._high) || !PreciseMathHelper.IsDivisionWithinFastRange(right._high))
        {
            return PreciseMathHelper.DivideBoundary(left, right);
        }

        double quotient = left._high / right._high;
        DoubleDouble remainder = left - (right * quotient);
        double correction = remainder._high / right._high;
        remainder -= right * correction;
        double finalCorrection = remainder._high / right._high;
        return FromComponents(quotient, correction) + finalCorrection;
    }

    /// <summary>Applies the expansion operation without discarding the low component.</summary>
    public static DoubleDouble operator +(DoubleDouble left, double right)
    {
        return PreciseMathHelper.AddScalar(left._high, left._low, right);
    }

    /// <summary>Applies the expansion operation without discarding the low component.</summary>
    public static DoubleDouble operator +(double left, DoubleDouble right)
    {
        return right + left;
    }

    /// <summary>Applies the expansion operation without discarding the low component.</summary>
    public static DoubleDouble operator -(DoubleDouble left, double right)
    {
        return PreciseMathHelper.AddScalar(left._high, left._low, -right);
    }

    /// <summary>Applies the expansion operation without discarding the low component.</summary>
    public static DoubleDouble operator -(double left, DoubleDouble right)
    {
        // Negate the components, not the result: exact cancellation must yield +0.
        return PreciseMathHelper.AddScalar(-right._high, -right._low, left);
    }

    /// <summary>Applies the expansion operation without discarding the low component.</summary>
    public static DoubleDouble operator *(DoubleDouble left, double right)
    {
        if (!IsFinite(left) || !double.IsFinite(right) || left._high == 0.0 || right == 0.0)
        {
            return new DoubleDouble(left._high * right);
        }
        if (!PreciseMathHelper.IsMultiplicationWithinFastRange(left._high, right))
        {
            return PreciseMathHelper.MultiplyBoundary(left, right);
        }

        (double product, double error) = PreciseMathHelper.TwoMultiply(left._high, right);
        error = Math.FusedMultiplyAdd(left._low, right, error);
        // Normalized input bounds the correction by O(u * product). The exponent
        // guard keeps the high product normal and its product residual representable.
        (double high, double low) = PreciseMathHelper.TwoQuickAdd(product, error);
        return new DoubleDouble(high, low == 0.0 ? 0.0 : low);
    }

    /// <summary>Applies the expansion operation without discarding the low component.</summary>
    public static DoubleDouble operator *(double left, DoubleDouble right)
    {
        return right * left;
    }

    /// <summary>Applies the expansion operation without discarding the low component.</summary>
    public static DoubleDouble operator /(DoubleDouble left, double right)
    {
        if (!IsFinite(left) || !double.IsFinite(right) || left._high == 0.0 || right == 0.0)
        {
            return new DoubleDouble(left._high / right);
        }
        if (!PreciseMathHelper.IsDivisionWithinFastRange(left._high) || !PreciseMathHelper.IsDivisionWithinFastRange(right))
        {
            return PreciseMathHelper.DivideBoundary(left, right);
        }

        double quotient = left._high / right;
        double remainder = Math.FusedMultiplyAdd(-quotient, right, left._high);
        double correction = (remainder + left._low) / right;
        // The exponent guard keeps the quotient normal. The correction is
        // O(u * quotient), so QuickTwoSum is ordered; one correction gives O(u^2) error.
        (double high, double low) = PreciseMathHelper.TwoQuickAdd(quotient, correction);
        return new DoubleDouble(high, low == 0.0 ? 0.0 : low);
    }

    /// <summary>Applies the expansion operation without discarding the low component.</summary>
    public static DoubleDouble operator /(double left, DoubleDouble right)
    {
        if (!double.IsFinite(left) || !IsFinite(right) || left == 0.0 || right._high == 0.0)
        {
            return new DoubleDouble(left / right._high);
        }
        if (!PreciseMathHelper.IsDivisionWithinFastRange(left) || !PreciseMathHelper.IsDivisionWithinFastRange(right._high))
        {
            return PreciseMathHelper.DivideBoundary(left, right);
        }

        double quotient = left / right._high;
        double remainder = Math.FusedMultiplyAdd(-quotient, right._high, left);
        remainder = Math.FusedMultiplyAdd(-quotient, right._low, remainder);
        double correction = remainder / right._high;
        // Using the high denominator in the correction adds only O(u^2) error.
        // As above, the guarded quotient dominates its correction in magnitude.
        (double high, double low) = PreciseMathHelper.TwoQuickAdd(quotient, correction);
        return new DoubleDouble(high, low == 0.0 ? 0.0 : low);
    }
}