diff --git a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs new file mode 100644 index 0000000..221f2f0 --- /dev/null +++ b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs @@ -0,0 +1,198 @@ +namespace Just.PreciseMath; + +public readonly partial struct DoubleDouble : + IAdditiveIdentity, + IMultiplicativeIdentity, + IUnaryPlusOperators, + IUnaryNegationOperators, + IAdditionOperators, + IAdditionOperators, + ISubtractionOperators, + ISubtractionOperators, + IMultiplyOperators, + IMultiplyOperators, + IDivisionOperators, + IDivisionOperators +{ + /// Returns the operand unchanged. + public static DoubleDouble operator +(DoubleDouble value) + { + return value; + } + + /// Negates the value, including the high zero's sign, retaining canonical NaN and zero residuals. + 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); + } + + /// Adds normalized expansions, retaining low-sum residuals under cancellation. + 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 (Math.Max(Math.ILogB(left._high), Math.ILogB(right._high)) > 1020) + { + return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left) + PreciseMathHelper.ArithmeticUnits(right), BigInteger.One << 1074); + } + + (double high, double highError) = PreciseMathHelper.TwoAdd(left._high, right._high); + (double low, double lowError) = PreciseMathHelper.TwoAdd(left._low, right._low); + (double middle, double middleError) = PreciseMathHelper.TwoAdd(highError, low); + (double sum, double sumError) = PreciseMathHelper.TwoAdd(high, middle); + return FromComponents(sum, sumError + (middleError + lowError)); + } + + /// Subtracts normalized expansions. + public static DoubleDouble operator -(DoubleDouble left, DoubleDouble right) + { + return left + (-right); + } + + /// Multiplies expansions using an FMA product residual and cross terms. + /// Results are approximate double-double values, not universally correctly rounded. + 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); + } + int exponent = Math.ILogB(left._high) + Math.ILogB(right._high); + if (exponent < -900 || exponent > 900) + { + return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left) * PreciseMathHelper.ArithmeticUnits(right), BigInteger.One << 2148); + } + + (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); + return FromComponents(product, error); + } + + /// Divides expansions using a quotient estimate and two residual corrections. + /// Zero and nonfinite operands follow binary64 rules; precision decreases near underflow. + 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); + } + int leftExponent = Math.ILogB(left._high); + int rightExponent = Math.ILogB(right._high); + if (Math.Abs(leftExponent) > 450 || Math.Abs(rightExponent) > 450) + { + return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left), PreciseMathHelper.ArithmeticUnits(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; + } + + /// Applies the expansion operation without discarding the low component. + public static DoubleDouble operator +(DoubleDouble left, double right) + { + return PreciseMathHelper.AddScalar(left._high, left._low, right); + } + + /// Applies the expansion operation without discarding the low component. + public static DoubleDouble operator +(double left, DoubleDouble right) + { + return right + left; + } + + /// Applies the expansion operation without discarding the low component. + public static DoubleDouble operator -(DoubleDouble left, double right) + { + return PreciseMathHelper.AddScalar(left._high, left._low, -right); + } + + /// Applies the expansion operation without discarding the low component. + 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); + } + + /// Applies the expansion operation without discarding the low component. + 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); + } + int exponent = Math.ILogB(left._high) + Math.ILogB(right); + if (exponent < -900 || exponent > 900) + { + return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left) * PreciseMathHelper.ArithmeticUnits(right), BigInteger.One << 2148); + } + + (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); + } + + /// Applies the expansion operation without discarding the low component. + public static DoubleDouble operator *(double left, DoubleDouble right) + { + return right * left; + } + + /// Applies the expansion operation without discarding the low component. + 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); + } + int leftExponent = Math.ILogB(left._high); + int rightExponent = Math.ILogB(right); + if (Math.Abs(leftExponent) > 450 || Math.Abs(rightExponent) > 450) + { + return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left), PreciseMathHelper.ArithmeticUnits(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); + } + + /// Applies the expansion operation without discarding the low component. + 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); + } + int leftExponent = Math.ILogB(left); + int rightExponent = Math.ILogB(right._high); + if (Math.Abs(leftExponent) > 450 || Math.Abs(rightExponent) > 450) + { + return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left), PreciseMathHelper.ArithmeticUnits(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); + } +} diff --git a/0-source/Just.PreciseMath/DoubleDouble.Comparison.cs b/0-source/Just.PreciseMath/DoubleDouble.Comparison.cs new file mode 100644 index 0000000..c4f439c --- /dev/null +++ b/0-source/Just.PreciseMath/DoubleDouble.Comparison.cs @@ -0,0 +1,93 @@ +namespace Just.PreciseMath; + +public readonly partial struct DoubleDouble : IComparable, IComparable, + IComparisonOperators +{ + /// Returns the normalized components. + public void Decompose(out double high, out double low) + { + high = _high; + low = _low; + } + + /// Orders NaN before other values, and compares finite values using both components. + public int CompareTo(DoubleDouble other) + { + int order = _high.CompareTo(other._high); + return order != 0 ? order : _low.CompareTo(other._low); + } + + /// + public int CompareTo(object? obj) + { + if (obj is null) + { + return 1; + } + if (obj is DoubleDouble other) + { + return CompareTo(other); + } + throw new ArgumentException($"Object must be of type {nameof(DoubleDouble)}.", nameof(obj)); + } + + /// Applies binary64 IsNaN classification to the canonical high component. + public static bool IsNaN(DoubleDouble value) + { + return double.IsNaN(value._high); + } + + /// Applies binary64 IsFinite classification to the canonical high component. + public static bool IsFinite(DoubleDouble value) + { + return double.IsFinite(value._high); + } + + /// Applies binary64 IsInfinity classification to the canonical high component. + public static bool IsInfinity(DoubleDouble value) + { + return double.IsInfinity(value._high); + } + + /// Applies binary64 IsPositiveInfinity classification to the canonical high component. + public static bool IsPositiveInfinity(DoubleDouble value) + { + return double.IsPositiveInfinity(value._high); + } + + /// Applies binary64 IsNegativeInfinity classification to the canonical high component. + public static bool IsNegativeInfinity(DoubleDouble value) + { + return double.IsNegativeInfinity(value._high); + } + + /// Applies binary64 IsNegative classification to the canonical high component. + public static bool IsNegative(DoubleDouble value) + { + return double.IsNegative(value._high); + } + + /// Compares both components; NaN operands are unordered. + public static bool operator <(DoubleDouble left, DoubleDouble right) + { + return left._high < right._high || (left._high == right._high && left._low < right._low); + } + + /// Compares both components; NaN operands are unordered. + public static bool operator >(DoubleDouble left, DoubleDouble right) + { + return left._high > right._high || (left._high == right._high && left._low > right._low); + } + + /// Compares both components; NaN operands are unordered. + public static bool operator <=(DoubleDouble left, DoubleDouble right) + { + return left._high < right._high || (left._high == right._high && left._low <= right._low); + } + + /// Compares both components; NaN operands are unordered. + public static bool operator >=(DoubleDouble left, DoubleDouble right) + { + return left._high > right._high || (left._high == right._high && left._low >= right._low); + } +} diff --git a/0-source/Just.PreciseMath/DoubleDouble.Conversion.cs b/0-source/Just.PreciseMath/DoubleDouble.Conversion.cs new file mode 100644 index 0000000..deb6b12 --- /dev/null +++ b/0-source/Just.PreciseMath/DoubleDouble.Conversion.cs @@ -0,0 +1,312 @@ +namespace Just.PreciseMath; + +/// +/// Explicit integer conversions truncate toward zero; IConvertible integer conversions round +/// to nearest with ties to even. Both check the resulting integer's range and reject nonfinite +/// values with OverflowException. IConvertible reports TypeCode.Object; Boolean conversion is +/// false only for zero. Char, DateTime, and enum conversions throw InvalidCastException. +/// Numeric conversions ignore their format provider; string conversion uses it. +/// +public readonly partial struct DoubleDouble : IConvertible +{ + TypeCode IConvertible.GetTypeCode() + { + return TypeCode.Object; + } + + bool IConvertible.ToBoolean(IFormatProvider? provider) + { + return _high != 0.0; + } + + char IConvertible.ToChar(IFormatProvider? provider) + { + throw new InvalidCastException("DoubleDouble cannot be converted to Char."); + } + + DateTime IConvertible.ToDateTime(IFormatProvider? provider) + { + throw new InvalidCastException("DoubleDouble cannot be converted to DateTime."); + } + + byte IConvertible.ToByte(IFormatProvider? provider) + { + return (byte)ConversionRoundedInteger(); + } + + sbyte IConvertible.ToSByte(IFormatProvider? provider) + { + return (sbyte)ConversionRoundedInteger(); + } + + short IConvertible.ToInt16(IFormatProvider? provider) + { + return (short)ConversionRoundedInteger(); + } + + ushort IConvertible.ToUInt16(IFormatProvider? provider) + { + return (ushort)ConversionRoundedInteger(); + } + + int IConvertible.ToInt32(IFormatProvider? provider) + { + return (int)ConversionRoundedInteger(); + } + + uint IConvertible.ToUInt32(IFormatProvider? provider) + { + return (uint)ConversionRoundedInteger(); + } + + long IConvertible.ToInt64(IFormatProvider? provider) + { + return (long)ConversionRoundedInteger(); + } + + ulong IConvertible.ToUInt64(IFormatProvider? provider) + { + return (ulong)ConversionRoundedInteger(); + } + + decimal IConvertible.ToDecimal(IFormatProvider? provider) + { + return (decimal)this; + } + + double IConvertible.ToDouble(IFormatProvider? provider) + { + return (double)this; + } + + float IConvertible.ToSingle(IFormatProvider? provider) + { + return (float)this; + } + + object IConvertible.ToType(Type conversionType, IFormatProvider? provider) + { + ArgumentNullException.ThrowIfNull(conversionType); + if (conversionType == typeof(DoubleDouble) || conversionType == typeof(object)) + { + return this; + } + if (!conversionType.IsEnum && Type.GetTypeCode(conversionType) is TypeCode code && code is not (TypeCode.Object or TypeCode.Empty or TypeCode.DBNull)) + { + return Convert.ChangeType(this, code, provider); + } + throw new InvalidCastException($"DoubleDouble cannot be converted to {conversionType.Name}."); + } + + private BigInteger ConversionRoundedInteger() + { + (BigInteger numerator, BigInteger denominator) = ConversionFraction(); + return ConversionRoundQuotient(numerator, denominator); + } + /// Constructs an exact representation of a 32-bit integer. + public DoubleDouble(int value) : this((double)value) + { + } + + /// Constructs an exact representation of a 64-bit integer. + public DoubleDouble(long value) + { + double high = value; + // The rounded high is integral but may be +2^63 for long.MaxValue. + // Int128 holds it and the exact residual without allocating. + double low = (double)((Int128)value - (Int128)high); + this = new DoubleDouble(high, low); + } + + /// Converts a 32-bit integer exactly. + public static explicit operator DoubleDouble(int value) + { + return new DoubleDouble(value); + } + + /// Converts a 64-bit integer exactly. + public static explicit operator DoubleDouble(long value) + { + return new DoubleDouble(value); + } + + /// Truncates toward zero; throws OverflowException when the truncated value is out of range or nonfinite. + public static explicit operator int(DoubleDouble value) + { + return (int)value.ConversionInteger(); + } + + /// Truncates toward zero; throws OverflowException when the truncated value is out of range or nonfinite. + public static explicit operator long(DoubleDouble value) + { + return (long)value.ConversionInteger(); + } + + /// Converts a binary64 value without losing information, preserving signed zero. + public static explicit operator DoubleDouble(double value) + { + return new DoubleDouble(value); + } + + /// Converts a binary32 value exactly, preserving signed zero. + public static explicit operator DoubleDouble(float value) + { + return new DoubleDouble((double)value); + } + + /// Rounds to binary64, nearest with ties to even; preserves nonfinite values and signed zero. + public static explicit operator double(DoubleDouble value) + { + // Normalization already rounds the complete sum to the high component. + return value._high; + } + + /// Rounds the exact component sum directly to binary32, nearest with ties to even. + public static explicit operator float(DoubleDouble value) + { + // Canonical nonfinite values and zeros also have a zero low component. + // With no residual, the binary64-to-binary32 cast already rounds once. + if (value._low == 0.0) + { + return (float)value._high; + } + (BigInteger numerator, BigInteger denominator) = value.ConversionFraction(); + return (float)ConversionRoundBinary(numerator, denominator, 24, -149); + } + + /// + /// Constructs from the exact decimal coefficient and scale, rounding the high component + /// and then its exact residual to binary64, each with ties to even. + /// + public DoubleDouble(decimal value) + { + Span bits = stackalloc int[4]; + decimal.GetBits(value, bits); + BigInteger numerator = (uint)bits[0] + ((BigInteger)(uint)bits[1] << 32) + ((BigInteger)(uint)bits[2] << 64); + BigInteger denominator = BigInteger.Pow(10, (bits[3] >> 16) & 0xff); + bool negative = bits[3] < 0; + if (negative) + { + numerator = -numerator; + } + double high = ConversionRoundBinary(numerator, denominator, 53, -1074); + (BigInteger highNumerator, BigInteger highDenominator) = ConversionDoubleFraction(high); + double low = ConversionRoundBinary((numerator * highDenominator) - (highNumerator * denominator), denominator * highDenominator, 53, -1074); + // Rounding the residual can reach a midpoint: normalize the resulting pair. + this = FromComponents(numerator.IsZero && negative ? -0.0 : high, low); + } + + /// Converts decimal using the exact coefficient and scale, not a decimal round trip. + public static explicit operator DoubleDouble(decimal value) + { + return new DoubleDouble(value); + } + + /// + /// Rounds the exact sum to the greatest decimal scale (up to 28) whose coefficient fits + /// 96 bits, with ties to even. Nonfinite values or magnitudes above decimal.MaxValue throw OverflowException. + /// + public static explicit operator decimal(DoubleDouble value) + { + (BigInteger numerator, BigInteger denominator) = value.ConversionFraction(); + bool negative = double.IsNegative(value._high); + numerator = BigInteger.Abs(numerator); + BigInteger maximum = (BigInteger.One << 96) - 1; + if (numerator > maximum * denominator) + { + throw new OverflowException("The value is outside the decimal range."); + } + for (int scale = 28; scale >= 0; scale--) + { + BigInteger coefficient = ConversionRoundQuotient(numerator * BigInteger.Pow(10, scale), denominator); + if (coefficient <= maximum) + { + return new decimal(unchecked((int)(uint)(coefficient & uint.MaxValue)), + unchecked((int)(uint)((coefficient >> 32) & uint.MaxValue)), + unchecked((int)(uint)(coefficient >> 64)), negative, (byte)scale); + } + } + throw new OverflowException("The value is outside the decimal range."); + } + + // Signed numerator and positive denominator; symmetric nearest-even integer rounding. + private static BigInteger ConversionRoundQuotient(BigInteger numerator, BigInteger denominator) + { + BigInteger quotient = BigInteger.DivRem(BigInteger.Abs(numerator), denominator, out BigInteger remainder); + int comparison = (remainder << 1).CompareTo(denominator); + if (comparison > 0 || (comparison == 0 && !quotient.IsEven)) + { + quotient++; + } + return numerator.Sign < 0 ? -quotient : quotient; + } + + // Rounds a rational directly to a binary precision, with a minimum subnormal quantum. + // Used for binary64 decimal decomposition and binary32 output (returned exactly in binary64). + private static double ConversionRoundBinary(BigInteger numerator, BigInteger denominator, int precision, int minimumShift) + { + if (numerator.IsZero) + { + return 0.0; + } + bool negative = numerator.Sign < 0; + numerator = BigInteger.Abs(numerator); + int exponent = checked((int)(numerator.GetBitLength() - denominator.GetBitLength())); + bool below = exponent >= 0 ? numerator < (denominator << exponent) : (numerator << -exponent) < denominator; + if (below) + { + exponent--; + } + int shift = Math.Max(exponent - precision + 1, minimumShift); + BigInteger rounded = shift >= 0 + ? ConversionRoundQuotient(numerator, denominator << shift) + : ConversionRoundQuotient(numerator << -shift, denominator); + double result = Math.ScaleB((double)rounded, shift); + return negative ? -result : result; + } + + private BigInteger ConversionInteger() + { + (BigInteger numerator, BigInteger denominator) = ConversionFraction(); + return numerator / denominator; + } + + // Exact signed rational decomposition, with positive denominator; rejects nonfinite values. + private (BigInteger Numerator, BigInteger Denominator) ConversionFraction() + { + if (!double.IsFinite(_high)) + { + throw new OverflowException("A nonfinite value cannot be converted to a finite number."); + } + (BigInteger highNumerator, BigInteger highDenominator) = ConversionDoubleFraction(_high); + (BigInteger lowNumerator, BigInteger lowDenominator) = ConversionDoubleFraction(_low); + // Both denominators are powers of two: align to the larger one rather + // than multiplying them and inflating every subsequent exact operation. + int shift = (int)(highDenominator.GetBitLength() - lowDenominator.GetBitLength()); + return shift >= 0 + ? (highNumerator + (lowNumerator << shift), highDenominator) + : ((highNumerator << -shift) + lowNumerator, lowDenominator); + } + + // The input must be finite. Binary64 is an integer significand times a power of two. + private static (BigInteger Numerator, BigInteger Denominator) ConversionDoubleFraction(double value) + { + if (value == 0.0) + { + return (BigInteger.Zero, BigInteger.One); + } + ulong bits = BitConverter.DoubleToUInt64Bits(value); + int exponent = (int)((bits >> 52) & 0x7ff); + BigInteger significand = bits & 0x000fffffffffffffUL; + if (exponent != 0) + { + significand += BigInteger.One << 52; + } + int shift = exponent == 0 ? -1074 : exponent - 1075; + if ((bits >> 63) != 0) + { + significand = -significand; + } + return shift >= 0 ? (significand << shift, BigInteger.One) : (significand, BigInteger.One << -shift); + } +} diff --git a/0-source/Just.PreciseMath/DoubleDouble.Formatting.cs b/0-source/Just.PreciseMath/DoubleDouble.Formatting.cs new file mode 100644 index 0000000..5ca76b0 --- /dev/null +++ b/0-source/Just.PreciseMath/DoubleDouble.Formatting.cs @@ -0,0 +1,128 @@ +using System.Globalization; + +namespace Just.PreciseMath; + +public readonly partial struct DoubleDouble : IFormattable +{ + /// Formats with G32 and the current culture. + public override string ToString() + { + return ToString(null, null); + } + + /// Formats with the specified G, E, or F format and the current culture. + public string ToString(string? format) + { + return ToString(format, null); + } + + /// Formats with G32 and the supplied culture (or current culture when null). + public string ToString(IFormatProvider? provider) + { + return ToString(null, provider); + } + + /// + /// Formats the exact component sum, rounding to nearest with ties to even. + /// Supports G/g (significant digits, default 32; G0 also means 32), E/e + /// (fractional digits, default 6), and F/f (fractional digits, culture default). + /// Precision is limited to 0 through 999. Other standard or custom formats throw FormatException. + /// G uses scientific notation for exponents below -4 or at least the precision; + /// trailing fractional zeros are removed. G is not a shortest-round-trip format. + /// Signed zero and culture-specific signs, separators, and nonfinite symbols are preserved. + /// + public string ToString(string? format, IFormatProvider? formatProvider) + { + NumberFormatInfo info = NumberFormatInfo.GetInstance(formatProvider); + char specifier = string.IsNullOrEmpty(format) ? 'G' : format[0]; + char kind = char.ToUpperInvariant(specifier); + if (kind is not ('G' or 'E' or 'F')) + { + throw new FormatException("Only G, E, and F numeric formats are supported."); + } + int precision = kind == 'G' ? 32 : kind == 'E' ? 6 : info.NumberDecimalDigits; + if (format is { Length: > 1 }) + { + precision = 0; + foreach (char digit in format.AsSpan(1)) + { + if (digit is < '0' or > '9' || precision > 99) + { + throw new FormatException("Numeric precision must be between 0 and 999."); + } + precision = (precision * 10) + digit - '0'; + } + } + if (kind == 'G' && precision == 0) + { + precision = 32; + } + if (double.IsNaN(_high)) + { + return info.NaNSymbol; + } + if (double.IsInfinity(_high)) + { + return double.IsNegative(_high) ? info.NegativeInfinitySymbol : info.PositiveInfinitySymbol; + } + (BigInteger numerator, BigInteger denominator) = ConversionFraction(); + numerator = BigInteger.Abs(numerator); + string sign = double.IsNegative(_high) ? info.NegativeSign : string.Empty; + if (kind == 'F') + { + BigInteger rounded = ConversionRoundQuotient(numerator * BigInteger.Pow(10, precision), denominator); + return sign + FormattingFixed(rounded.ToString(CultureInfo.InvariantCulture), precision, info.NumberDecimalSeparator); + } + + int exponent = 0; + if (!numerator.IsZero) + { + // Decimal digit lengths give an estimate no more than one above floor(log10(n/d)). + exponent = numerator.ToString(CultureInfo.InvariantCulture).Length - denominator.ToString(CultureInfo.InvariantCulture).Length; + bool below = exponent >= 0 ? numerator < denominator * BigInteger.Pow(10, exponent) : numerator * BigInteger.Pow(10, -exponent) < denominator; + if (below) + { + exponent--; + } + } + int significantDigits = kind == 'E' ? precision + 1 : precision; + int shift = significantDigits - 1 - exponent; + BigInteger coefficient = shift >= 0 + ? ConversionRoundQuotient(numerator * BigInteger.Pow(10, shift), denominator) + : ConversionRoundQuotient(numerator, denominator * BigInteger.Pow(10, -shift)); + string digits = coefficient.ToString(CultureInfo.InvariantCulture); + if (digits.Length > significantDigits) + { + digits = digits[..^1]; + exponent++; + } + digits = digits.PadLeft(significantDigits, '0'); + if (kind == 'E' || exponent < -4 || exponent >= precision) + { + string fractional = kind == 'G' ? digits[1..].TrimEnd('0') : digits[1..]; + string mantissa = digits[..1] + (fractional.Length == 0 ? string.Empty : info.NumberDecimalSeparator + fractional); + string exponentSign = exponent < 0 ? info.NegativeSign : info.PositiveSign; + string exponentDigits = Math.Abs(exponent).ToString(kind == 'E' ? "D3" : "D2", CultureInfo.InvariantCulture); + return sign + mantissa + (char.IsLower(specifier) ? "e" : "E") + exponentSign + exponentDigits; + } + int decimalPlaces = significantDigits - 1 - exponent; + // Trim only fractional zeros, before inserting a potentially multi-character separator. + int end = digits.Length; + while (decimalPlaces > 0 && digits[end - 1] == '0') + { + end--; + decimalPlaces--; + } + return sign + FormattingFixed(digits[..end], decimalPlaces, info.NumberDecimalSeparator); + } + + private static string FormattingFixed(string digits, int decimalPlaces, string separator) + { + if (decimalPlaces <= 0) + { + return digits + new string('0', -decimalPlaces); + } + digits = digits.PadLeft(decimalPlaces + 1, '0'); + return digits[..^decimalPlaces] + separator + digits[^decimalPlaces..]; + } +} diff --git a/0-source/Just.PreciseMath/DoubleDouble.Parsing.cs b/0-source/Just.PreciseMath/DoubleDouble.Parsing.cs new file mode 100644 index 0000000..fd32c93 --- /dev/null +++ b/0-source/Just.PreciseMath/DoubleDouble.Parsing.cs @@ -0,0 +1,207 @@ +using System.Globalization; + +namespace Just.PreciseMath; + +/// +/// Parsing supports decimal/scientific notation with ASCII digits, surrounding whitespace, +/// culture-specific signs and decimal separator, and NaN/infinity symbols (case-insensitive). +/// Group separators, currency, hexadecimal notation, and NumberStyles options are not supported. +/// Inputs are limited to 4096 characters, including surrounding whitespace, to bound work. +/// The exact decimal coefficient and exponent are converted to a normalized high/low pair, +/// not through double or decimal. Components are rounded nearest, ties to even, then normalized; +/// a second rounding at the overflow midpoint is kept finite when the exact input is below it. +/// Overflow succeeds with signed infinity; underflow succeeds with signed zero. +/// This is not a shortest-round-trip parser/formatter contract. +/// +public readonly partial struct DoubleDouble : ISpanParsable +{ + private const int MaximumParsingLength = 4096; + + /// Parses decimal/scientific text, using the current culture when provider is null. + /// The input is null. + /// The input is malformed, unsupported, or longer than 4096 characters. + public static DoubleDouble Parse(string s, IFormatProvider? provider = null) + { + ArgumentNullException.ThrowIfNull(s); + return Parse(s.AsSpan(), provider); + } + + /// Parses decimal/scientific text, using the current culture when provider is null. + /// The input is malformed, unsupported, or longer than 4096 characters. + public static DoubleDouble Parse(ReadOnlySpan s, IFormatProvider? provider = null) + { + if (!TryParse(s, provider, out DoubleDouble result)) + { + throw new FormatException("Invalid or unsupported DoubleDouble text (maximum 4096 characters)."); + } + return result; + } + + /// Parses using the current culture. Returns false and Zero for null, invalid, or oversized input. + public static bool TryParse(string? s, out DoubleDouble result) + { + return TryParse(s, null, out result); + } + + /// Parses using the supplied culture (current when null). Returns false and Zero for null, invalid, or oversized input. + public static bool TryParse(string? s, IFormatProvider? provider, out DoubleDouble result) + { + return TryParse(s.AsSpan(), provider, out result); + } + + /// Parses using the current culture. Returns false and Zero for invalid or oversized input. + public static bool TryParse(ReadOnlySpan s, out DoubleDouble result) + { + return TryParse(s, null, out result); + } + + /// Parses using the supplied culture (current when null). Returns false and Zero for invalid or oversized input. + public static bool TryParse(ReadOnlySpan s, IFormatProvider? provider, out DoubleDouble result) + { + result = Zero; + if (s.Length > MaximumParsingLength) + { + return false; + } + s = s.Trim(); + if (s.IsEmpty) + { + return false; + } + NumberFormatInfo info = NumberFormatInfo.GetInstance(provider); + // Custom symbols can themselves start with a numeric sign. + if (s.Equals(info.NaNSymbol, StringComparison.OrdinalIgnoreCase)) + { + result = NaN; + return true; + } + if (s.Equals(info.PositiveInfinitySymbol, StringComparison.OrdinalIgnoreCase)) + { + result = new DoubleDouble(double.PositiveInfinity); + return true; + } + if (s.Equals(info.NegativeInfinitySymbol, StringComparison.OrdinalIgnoreCase)) + { + result = new DoubleDouble(double.NegativeInfinity); + return true; + } + bool negative = ParsingConsumeSign(ref s, info); + if (s.Equals(info.NaNSymbol, StringComparison.OrdinalIgnoreCase)) + { + result = NaN; + return true; + } + if (s.Equals(info.PositiveInfinitySymbol, StringComparison.OrdinalIgnoreCase)) + { + result = new DoubleDouble(negative ? double.NegativeInfinity : double.PositiveInfinity); + return true; + } + + BigInteger coefficient = BigInteger.Zero; + int significantDigits = 0; + int fractionalDigits = 0; + bool hasDigit = false; + bool hasSeparator = false; + int position = 0; + while (position < s.Length) + { + char digit = s[position]; + if (digit is >= '0' and <= '9') + { + coefficient = (coefficient * 10) + digit - '0'; + hasDigit = true; + if (!coefficient.IsZero) + { + significantDigits++; + } + if (hasSeparator) + { + fractionalDigits++; + } + position++; + } + else if (!hasSeparator && s[position..].StartsWith(info.NumberDecimalSeparator, StringComparison.Ordinal)) + { + hasSeparator = true; + position += info.NumberDecimalSeparator.Length; + } + else + { + break; + } + } + if (!hasDigit) + { + return false; + } + + int exponent = 0; + if (position < s.Length && s[position] is 'e' or 'E') + { + s = s[(position + 1)..]; + bool negativeExponent = ParsingConsumeSign(ref s, info); + if (s.IsEmpty) + { + return false; + } + foreach (char digit in s) + { + if (digit is < '0' or > '9') + { + return false; + } + // Saturate beyond any offset the bounded mantissa can cancel. + // Huge exponents still require validating every remaining digit. + exponent = Math.Min((exponent * 10) + digit - '0', MaximumParsingLength * 2); + } + exponent = negativeExponent ? -exponent : exponent; + } + else if (position != s.Length) + { + return false; + } + + int decimalExponent = exponent - fractionalDigits; + int magnitude = significantDigits + decimalExponent - 1; + if (coefficient.IsZero || magnitude < -324) + { + result = new DoubleDouble(negative ? -0.0 : 0.0); + } + else if (magnitude > 308) + { + result = new DoubleDouble(negative ? double.NegativeInfinity : double.PositiveInfinity); + } + else + { + // Boundary decades need exact rounding. The magnitude checks also + // keep enormous exponents from requesting unbounded powers of ten. + BigInteger numerator = negative ? -coefficient : coefficient; + BigInteger denominator = BigInteger.One; + if (decimalExponent >= 0) + { + numerator *= BigInteger.Pow(10, decimalExponent); + } + else + { + denominator = BigInteger.Pow(10, -decimalExponent); + } + result = PreciseMathHelper.ArithmeticFromRatio(numerator, denominator); + } + return true; + } + + // A sign is optional; do not consume whitespace between it and the number. + private static bool ParsingConsumeSign(ref ReadOnlySpan text, NumberFormatInfo info) + { + if (info.PositiveSign.Length != 0 && text.StartsWith(info.PositiveSign, StringComparison.Ordinal)) + { + text = text[info.PositiveSign.Length..]; + } + else if (info.NegativeSign.Length != 0 && text.StartsWith(info.NegativeSign, StringComparison.Ordinal)) + { + text = text[info.NegativeSign.Length..]; + return true; + } + return false; + } +} diff --git a/0-source/Just.PreciseMath/DoubleDouble.cs b/0-source/Just.PreciseMath/DoubleDouble.cs index 17f03ed..534bcaa 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.cs @@ -1,40 +1,99 @@ namespace Just.PreciseMath; /// -/// Represents higher precision floating point type +/// Represents a normalized, fixed-size sum of two binary64 values. /// -public readonly struct DoubleDouble : +/// +/// Finite components are nonoverlapping; the high component is the rounded sum and +/// the low component retains its residual. NaN and infinities have a positive-zero +/// low component. A zero low input preserves the high zero's sign; exact nonzero +/// cancellation produces positive zero. Exponent range is limited by binary64, +/// and precision decreases near underflow. Arithmetic is not guaranteed correctly +/// rounded. Equals treats NaNs as equal for collections, while operators do not. +/// +public readonly partial struct DoubleDouble : IEquatable, IEqualityOperators { internal readonly double _high; internal readonly double _low; - // No normalization, for internal use only + /// + /// Stores trusted components without normalization or validation. + /// + /// + /// Callers must supply a normalized finite pair, or a canonical NaN/infinity + /// with a positive-zero low component. Zero residuals must be positive zero. + /// Use FromComponents for arbitrary pairs or rounded arithmetic intermediates. + /// internal DoubleDouble(double high, double low) { _high = high; _low = low; } + /// + /// Creates the normalized sum of two arbitrary components. Nonfinite sums have + /// a positive-zero low component; NaNs are canonicalized to double.NaN. + /// + /// + /// Uses round-to-nearest, ties-to-even binary64 arithmetic. A zero low input + /// preserves the high zero's sign; exact nonzero cancellation yields positive + /// zero. Overflow produces infinity and precision decreases near underflow. + /// + public static DoubleDouble FromComponents(double high, double low) + { + if (low == 0.0) + { + return new DoubleDouble(high); + } + + double sum = high + low; + if (!double.IsFinite(sum)) + { + return new DoubleDouble(sum); + } + + // Magnitude-ordered QuickTwoSum, reusing the checked sum. Subtracting + // the larger input avoids TwoSum's possible intermediate overflow when + // a small opposite-sign input precedes a value near the binary64 limit. + double error = Math.Abs(high) >= Math.Abs(low) + ? low - (sum - high) + : high - (sum - low); + return new DoubleDouble(sum, error == 0.0 ? 0.0 : error); + } + /// /// Constructs new DoubleDouble from a given double. /// /// Initial high component - public DoubleDouble(double high) : this(high, 0.0) + public DoubleDouble(double high) : this(double.IsNaN(high) ? double.NaN : high, 0.0) { } #region Static constants + /// + /// Represents an additive identity value. + /// + public static DoubleDouble AdditiveIdentity => Zero; + /// + /// Represents a multiplicative identity value. + /// + public static DoubleDouble MultiplicativeIdentity => One; + /// /// Represents a value that is not a number (NaN). /// - public static DoubleDouble NaN => new(double.NaN, double.NaN); + public static DoubleDouble NaN => new(double.NaN); /// /// Represents a unit value. /// public static DoubleDouble One => new(1.0, 0); /// + /// Represents a negative unit value. + /// + public static DoubleDouble NegativeOne => new(-1.0, 0); + /// /// Represents a zero value. /// public static DoubleDouble Zero => new(); @@ -64,22 +123,37 @@ public readonly struct DoubleDouble : /// [Pure] - public bool Equals(DoubleDouble other) => _high == other._high && _low == other._low; + public bool Equals(DoubleDouble other) + { + return _high.Equals(other._high) && _low.Equals(other._low); + } /// [Pure] - public override bool Equals(object? obj) => obj is DoubleDouble other && this.Equals(other); + public override bool Equals(object? obj) + { + return obj is DoubleDouble other && Equals(other); + } /// [Pure] - public override int GetHashCode() => HashCode.Combine(_high, _low); + public override int GetHashCode() + { + return HashCode.Combine(_high, _low); + } /// - /// TODO: fill + /// Tests numerical equality; NaN operands are never equal. /// [Pure, MethodImpl(MethodImplOptions.AggressiveInlining)] - public static bool operator ==(DoubleDouble left, DoubleDouble right) => left.Equals(right); + public static bool operator ==(DoubleDouble left, DoubleDouble right) + { + return left._high == right._high && left._low == right._low; + } /// - /// TODO: fill + /// Tests numerical inequality; NaN operands are always unequal. /// [Pure, MethodImpl(MethodImplOptions.AggressiveInlining)] - public static bool operator !=(DoubleDouble left, DoubleDouble right) => !left.Equals(right); + public static bool operator !=(DoubleDouble left, DoubleDouble right) + { + return !(left == right); + } } diff --git a/0-source/Just.PreciseMath/PreciseMathHelper.cs b/0-source/Just.PreciseMath/PreciseMathHelper.cs index 7265b1d..5c909d3 100644 --- a/0-source/Just.PreciseMath/PreciseMathHelper.cs +++ b/0-source/Just.PreciseMath/PreciseMathHelper.cs @@ -2,6 +2,9 @@ namespace Just.PreciseMath; internal static class PreciseMathHelper { + // General TwoSum: no magnitude ordering required, but inputs, sum, and + // intermediate subtractions must stay finite. Arithmetic callers bound the + // exponents; arbitrary-component normalization uses magnitude ordering instead. [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static (double Res, double Err) TwoAdd(double a, double b) { @@ -10,6 +13,8 @@ internal static class PreciseMathHelper return (r, (a - (r - t)) + (b - t)); } + // QuickTwoSum requires |a| >= |b| and finite inputs/sum. It returns the + // rounded sum and its exact residual; callers canonicalize a zero residual. [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static (double Res, double Err) TwoQuickAdd(double a, double b) { @@ -18,13 +23,15 @@ internal static class PreciseMathHelper return (r, b - (r - a)); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static (double Res, double Err) TwoSubstract(double a, double b) + internal static (double Res, double Err) TwoSubtract(double a, double b) { double r = a - b; double t = r - a; return (r, (a - (r - t)) - (b + t)); } + // FMA gives the exact product residual when it is representable. Near + // underflow it rounds, and an overflowing product cannot use this transform. [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static (double Res, double Err) TwoMultiply(double a, double b) { @@ -33,10 +40,108 @@ internal static class PreciseMathHelper return (r, Math.FusedMultiplyAdd(a, b, -r)); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static (double Res, double Err) TwoSuare(double a) + internal static (double Res, double Err) TwoSquare(double a) { double r = a * a; return (r, Math.FusedMultiplyAdd(a, a, -r)); } + + // The first two arguments are normalized components (a negated zero low is + // also allowed). Sharing this path preserves both subtraction orders without + // constructing a temporary expansion for the scalar or the negated operand. + internal static DoubleDouble AddScalar(double high, double low, double value) + { + if (!double.IsFinite(high) || !double.IsFinite(value) || (high == 0.0 && value == 0.0)) + { + return new DoubleDouble(high + value); + } + if (Math.Max(Math.ILogB(high), Math.ILogB(value)) > 1020) + { + return ArithmeticFromRatio(ArithmeticUnits(high) + ArithmeticUnits(low) + ArithmeticUnits(value), BigInteger.One << 1074); + } + + (double sum, double error) = TwoAdd(high, value); + // Near high-component cancellation, Sterbenz makes the first sum exact, + // so error is zero and this retains low exactly. Otherwise its rounding + // contributes only O(u^2) relative error. The final TwoSum normalizes. + (double result, double residual) = TwoAdd(sum, error + low); + return new DoubleDouble(result, residual == 0.0 ? 0.0 : residual); + } + + // The boundary path uses bounded binary integers (at most about 4200 bits), not + // arbitrary-precision storage. It avoids overflow and double rounding in EFTs + // at the binary64 exponent limits. The common path remains allocation-free. + internal static BigInteger ArithmeticUnits(DoubleDouble value) + { + return ArithmeticUnits(value._high) + ArithmeticUnits(value._low); + } + + internal static BigInteger ArithmeticUnits(double value) + { + long bits = BitConverter.DoubleToInt64Bits(value); + int exponent = (int)((bits >> 52) & 0x7ff); + BigInteger significand = bits & 0xfffffffffffffL; + if (exponent != 0) + { + significand = (significand + (BigInteger.One << 52)) << (exponent - 1); + } + return bits < 0 ? -significand : significand; + } + + // Shared with decimal-text parsing. Denominator must be nonzero; round the + // high and exact residual, then normalize without a spurious second overflow. + internal static DoubleDouble ArithmeticFromRatio(BigInteger numerator, BigInteger denominator) + { + if (denominator.Sign < 0) + { + numerator = -numerator; + denominator = -denominator; + } + double high = ArithmeticRound(numerator, denominator); + if (!double.IsFinite(high) || high == 0.0) + { + return new DoubleDouble(high); + } + BigInteger residual = (numerator << 1074) - (ArithmeticUnits(high) * denominator); + double low = ArithmeticRound(residual, denominator << 1074); + if (double.IsInfinity(high + low)) + { + // The exact value rounded to a finite high, but rounding its residual + // can land on the overflow midpoint. Select the adjacent finite pair + // rather than overflow on this second rounding. True overflow already + // returned above; this loses at most one low-component ulp. + low = Math.BitDecrement(Math.Abs(low)) * Math.Sign(low); + } + return DoubleDouble.FromComponents(high, low); + } + + // Round an exact rational to binary64, ties-to-even, including subnormal and + // overflow boundaries. Denominator is positive; sign is retained on underflow. + private static double ArithmeticRound(BigInteger numerator, BigInteger denominator) + { + bool negative = numerator.Sign < 0; + numerator = BigInteger.Abs(numerator); + if (numerator.IsZero) + { + return 0.0; + } + int exponent = (int)(numerator.GetBitLength() - denominator.GetBitLength()); + bool below = exponent >= 0 ? numerator < (denominator << exponent) : (numerator << -exponent) < denominator; + if (below) + { + --exponent; + } + int shift = Math.Max(exponent - 52, -1074); + BigInteger dividend = shift < 0 ? numerator << -shift : numerator; + BigInteger divisor = shift > 0 ? denominator << shift : denominator; + BigInteger rounded = BigInteger.DivRem(dividend, divisor, out BigInteger remainder); + int comparison = (remainder << 1).CompareTo(divisor); + if (comparison > 0 || (comparison == 0 && !rounded.IsEven)) + { + ++rounded; + } + double result = Math.ScaleB((double)rounded, shift); + return negative ? -result : result; + } } diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs new file mode 100644 index 0000000..1fd6b98 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs @@ -0,0 +1,279 @@ +using System.Numerics; +using Shouldly; +using Xunit; + +namespace Just.PreciseMath.Tests; + +public class DoubleDoubleArithmeticTests +{ + [Fact] + public void CancellationRetainsBothLowSumTerms() + { + double small = Math.ScaleB(1.0, -54); + double tiny = Math.ScaleB(1.0, -108); + DoubleDouble left = DoubleDouble.FromComponents(1.0, small); + DoubleDouble right = DoubleDouble.FromComponents(-1.0, tiny); + Check(left + right, small, tiny); + Check(right + left, small, tiny); + Check(left - (-right), small, tiny); + Check(+left, 1.0, small); + Check(-left, -1.0, -small); + Check(left - left, 0.0, 0.0); + } + + [Fact] + public void ScalarOverloadsPreserveOperandOrderAndResiduals() + { + DoubleDouble value = DoubleDouble.FromComponents(2.0, Math.ScaleB(1.0, -80)); + Check(3.0 - value, 1.0, -value.Low); + Check(value - 3.0, -1.0, value.Low); + Check(value + 3.0, 5.0, value.Low); + Check(3.0 + value, 5.0, value.Low); + Check(value * 2.0, 4.0, 2.0 * value.Low); + Check(2.0 * value, 4.0, 2.0 * value.Low); + Check(value / 2.0, 1.0, value.Low / 2.0); + Check(6.0 / new DoubleDouble(2.0), 3.0, 0.0); + } + + [Fact] + public void ScalarLeftSubtractionAppliesTheRequestedOperandOrder() + { + // Review-1 §1: the former operator -(double, DoubleDouble) returned arg - lvalue, + // so 3.0 - DD(2.0) produced -1 instead of 1. + Check(3.0 - new DoubleDouble(2.0), 1.0, 0.0); + Check(new DoubleDouble(2.0) - 3.0, -1.0, 0.0); + Check(3.0 - DoubleDouble.FromComponents(2.0, Math.ScaleB(1.0, -80)), 1.0, -Math.ScaleB(1.0, -80)); + Check(DoubleDouble.FromComponents(2.0, Math.ScaleB(1.0, -80)) - 3.0, -1.0, Math.ScaleB(1.0, -80)); + } + + [Fact] + public void DivisionByAValueCarriedOnlyInTheLowComponentIsFinite() + { + // Review-2 §5: with high == 0 and low != 0 the former division returned + // Infinity because it divided by the zero high component. The public factory + // now folds such a pair into its high component, so the quotient is finite. + DoubleDouble denominator = DoubleDouble.FromComponents(0.0, 1.0); + Check(denominator, 1.0, 0.0); + Check(new DoubleDouble(4.0) / denominator, 4.0, 0.0); + Check(4.0 / denominator, 4.0, 0.0); + Check(new DoubleDouble(4.0) / DoubleDouble.FromComponents(0.0, -2.0), -2.0, 0.0); + } + + [Fact] + public void ScalarCancellationPreservesTheRemainingExpansionInBothOrders() + { + // At a normal binade boundary, 2^e - BitDecrement(2^e) = 2^(e-53). + // The low input becomes the representable residual of that exact difference. + foreach (int exponent in new[] { -967, -900, -450, 0, 450, 969, 1020, 1023 }) + { + foreach (double sign in new[] { -1.0, 1.0 }) + { + double high = Math.ScaleB(1.0, exponent); + double low = sign * Math.ScaleB(1.0, exponent - 107); + double scalar = sign * Math.BitDecrement(high); + double difference = sign * Math.ScaleB(1.0, exponent - 53); + DoubleDouble value = DoubleDouble.FromComponents(sign * high, low); + Check(value - scalar, difference, low); + Check(scalar - value, -difference, -low); + Check(value + (-scalar), difference, low); + Check((-scalar) + value, difference, low); + } + } + } + + [Fact] + public void MixedSpecialValuesIgnoreFiniteResidualsButPreserveResultSigns() + { + foreach (double sign in new[] { -1.0, 1.0 }) + { + DoubleDouble value = DoubleDouble.FromComponents(sign, sign * Math.ScaleB(1.0, -54)); + foreach (double scalar in new[] { double.NaN, double.NegativeInfinity, double.PositiveInfinity }) + { + CheckBits(value + scalar, sign + scalar); + CheckBits(scalar + value, scalar + sign); + CheckBits(value - scalar, sign - scalar); + CheckBits(scalar - value, scalar - sign); + CheckBits(value * scalar, sign * scalar); + CheckBits(scalar * value, scalar * sign); + CheckBits(value / scalar, sign / scalar); + CheckBits(scalar / value, scalar / sign); + } + foreach (double zero in new[] { 0.0, -0.0 }) + { + CheckBits(value * zero, sign * zero); + CheckBits(zero * value, zero * sign); + CheckBits(value / zero, sign / zero); + CheckBits(zero / value, zero / sign); + Check(value + zero, value.High, value.Low); + Check(zero + value, value.High, value.Low); + Check(value - zero, value.High, value.Low); + Check(zero - value, -value.High, -value.Low); + } + } + } + + [Fact] + public void ProductAndQuotientRetainExtraPrecision() + { + // (1 + 2^-52)(1 - 2^-52) = 1 - 2^-104 exactly. + Check(new DoubleDouble(1.0 + Math.ScaleB(1.0, -52)) * new DoubleDouble(1.0 - Math.ScaleB(1.0, -52)), + 1.0, -Math.ScaleB(1.0, -104)); + // Binary expansion of 1/3, rounding high then residual ties-to-even. + Check(DoubleDouble.One / 3.0, 0.3333333333333333, 1.850371707708594e-17); + } + + [Fact] + public void ScalarProductNormalizesACorrectionBeyondTheHighMidpoint() + { + // (1 + 2^-53)(1 + 2^-52) = 1 + 3*2^-53 + 2^-105, exactly. + // The rounded high advances twice above 1; its residual is still exact. + DoubleDouble value = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -53)); + double scalar = Math.BitIncrement(1.0); + double high = 1.0 + Math.ScaleB(1.0, -51); + double low = -Math.ScaleB(1.0, -53) + Math.ScaleB(1.0, -105); + Check(value * scalar, high, low); + Check(scalar * value, high, low); + Check(value * (-scalar), -high, -low); + Check((-scalar) * value, -high, -low); + } + + [Fact] + public void ScalarDivisionRetainsNumeratorAndDenominatorResiduals() + { + double low = Math.ScaleB(1.0, -80); + DoubleDouble value = DoubleDouble.FromComponents(1.0, low); + Check(value / 1.0, 1.0, low); + Check(value / (-1.0), -1.0, -low); + // Compare the reciprocal to its exact rational, not another DD operator. + BigInteger numerator = BigInteger.One << 2148; + AssertRelative(1.0 / value, numerator, Units(value)); + AssertRelative(-1.0 / value, -numerator, Units(value)); + } + + [Fact] + public void ExtremeFiniteOperationsDoNotOverflowIntermediates() + { + DoubleDouble maximum = new(double.MaxValue); + DoubleDouble third = maximum / 3.0; + DoubleDouble thirdPair = maximum / new DoubleDouble(3.0); + DoubleDouble thirdScalar = double.MaxValue / new DoubleDouble(3.0); + AssertRelative(third, Units(maximum), 3); + thirdPair.ShouldBe(third); + thirdScalar.ShouldBe(third); + Check(new DoubleDouble(double.Epsilon) / new DoubleDouble(double.Epsilon), 1.0, 0.0); + Check(new DoubleDouble(double.Epsilon) * new DoubleDouble(Math.ScaleB(1.0, 1023)), Math.ScaleB(1.0, -51), 0.0); + Check(new DoubleDouble(Math.ScaleB(1.0, -1022)) / 2.0, Math.ScaleB(1.0, -1023), 0.0); + Check(maximum * 2.0, double.PositiveInfinity, 0.0); + Check(maximum + maximum, double.PositiveInfinity, 0.0); + // High-only addition overflows, but the complete sum is exactly MaxValue. + DoubleDouble below = DoubleDouble.FromComponents(double.MaxValue, -Math.ScaleB(1.0, 969)); + Check(below + Math.ScaleB(1.0, 969), double.MaxValue, 0.0); + } + + [Theory] + [InlineData(1.0)] + [InlineData(-1.0)] + public void ExactOverflowMidpointStillOverflows(double sign) + { + DoubleDouble maximum = new(sign * double.MaxValue); + double halfUlp = sign * Math.ScaleB(1.0, 970); + Check(maximum + halfUlp, sign * double.PositiveInfinity, 0.0); + Check(maximum - (-halfUlp), sign * double.PositiveInfinity, 0.0); + Check(DoubleDouble.FromComponents(sign * double.MaxValue, halfUlp), sign * double.PositiveInfinity, 0.0); + } + + [Fact] + public void SpecialValueMatrixMatchesBinary64IncludingZeroSigns() + { + double[] values = [0.0, -0.0, 1.0, -1.0, double.PositiveInfinity, double.NegativeInfinity, double.NaN]; + foreach (double left in values) + { + foreach (double right in values) + { + DoubleDouble a = new(left); + DoubleDouble b = new(right); + CheckBits(a + b, left + right); + CheckBits(a - b, left - right); + CheckBits(a * b, left * right); + CheckBits(a / b, left / right); + CheckBits(a + right, left + right); + CheckBits(left + b, left + right); + CheckBits(a - right, left - right); + CheckBits(left - b, left - right); + CheckBits(a * right, left * right); + CheckBits(left * b, left * right); + CheckBits(a / right, left / right); + CheckBits(left / b, left / right); + } + } + } + + [Fact] + public void DeterministicArithmeticMeetsConservativeErrorBound() + { + Random random = new(1729); + for (int i = 0; i < 250; ++i) + { + DoubleDouble a = DoubleDouble.FromComponents(Math.ScaleB((random.NextDouble() * 2.0) - 1.0, random.Next(-400, 401)), + Math.ScaleB(random.NextDouble(), random.Next(-500, -450))); + DoubleDouble b = DoubleDouble.FromComponents(Math.ScaleB((random.NextDouble() * 2.0) - 1.0, random.Next(-400, 401)), + Math.ScaleB(random.NextDouble(), random.Next(-500, -450))); + BigInteger x = Units(a); + BigInteger y = Units(b); + AssertRelative(a + b, x + y, BigInteger.One); + AssertRelative(a - b, x - y, BigInteger.One); + AssertRelative(a * b, x * y, BigInteger.One << 1074); + AssertRelative(a / b, x << 1074, y); + } + } + + private static void Check(DoubleDouble value, double high, double low) + { + value.High.ShouldBe(high); + value.Low.ShouldBe(low); + } + + private static void CheckBits(DoubleDouble value, double expected) + { + if (double.IsNaN(expected)) + { + DoubleDouble.IsNaN(value).ShouldBeTrue(); + } + else + { + BitConverter.DoubleToInt64Bits(value.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected)); + } + BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(0L); + } + + // Independent oracle: every finite binary64 is an integer multiple of 2^-1074. + private static BigInteger Units(DoubleDouble value) + { + return Units(value.High) + Units(value.Low); + } + + private static BigInteger Units(double value) + { + long bits = BitConverter.DoubleToInt64Bits(value); + int exponent = (int)((bits >> 52) & 0x7ff); + BigInteger significand = bits & 0xfffffffffffffL; + if (exponent != 0) + { + significand += BigInteger.One << 52; + significand <<= exponent - 1; + } + return bits < 0 ? -significand : significand; + } + + private static void AssertRelative(DoubleDouble actual, BigInteger numerator, BigInteger denominator) + { + DoubleDouble.IsFinite(actual).ShouldBeTrue(); + BigInteger error = BigInteger.Abs((Units(actual) * denominator) - numerator); + // <= 2^-100 relative error plus one minimum subnormal (rounding floor). + (error <= (BigInteger.Abs(numerator) >> 100) + BigInteger.Abs(denominator)).ShouldBeTrue(); + if (actual.High != 0.0) + { + (Math.Abs(actual.Low) <= Math.ScaleB(1.0, Math.ILogB(actual.High) - 53)).ShouldBeTrue(); + } + } +} diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleBoundaryTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleBoundaryTests.cs new file mode 100644 index 0000000..1d0965a --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleBoundaryTests.cs @@ -0,0 +1,434 @@ +using System.Numerics; +using Shouldly; +using Xunit; + +namespace Just.PreciseMath.Tests; + +public class DoubleDoubleBoundaryTests +{ + [Theory] + [InlineData(1.0, "+")] + [InlineData(-1.0, "+")] + [InlineData(1.0, "-")] + [InlineData(-1.0, "-")] + public void AdditionBelowOverflowMidpointRemainsFinite(double sign, string operation) + { + // Exact magnitude = MaxValue + 2^970 - 2^916, strictly below + // the binary64 overflow midpoint. The right factory call reduces to + // (BitDecrement(2^969), 0), so scalar overloads share this case. + DoubleDouble left = DoubleDouble.FromComponents(sign * double.MaxValue, sign * Math.ScaleB(1.0, 969)); + DoubleDouble right = DoubleDouble.FromComponents(sign * Math.ScaleB(1.0, 969), -sign * Math.ScaleB(1.0, 916)); + if (operation == "-") + { + right = -right; + } + Rational expected = Expected(Exact(left), Exact(right), operation); + BelowOverflowMidpoint(expected).ShouldBeTrue(); + // A canonical finite pair meets the requested accuracy; infinity is + // not forced by the representational limit or the error contract. + DoubleDouble finiteWitness = DoubleDouble.FromComponents(sign * double.MaxValue, + sign * Math.BitDecrement(Math.ScaleB(1.0, 970))); + AssertAccurate(finiteWitness, expected, "finite witness"); + AssertOperation(left, right, operation); + } + + [Theory] + [InlineData(1.0, "+", false)] + [InlineData(-1.0, "+", false)] + [InlineData(1.0, "+", true)] + [InlineData(-1.0, "+", true)] + [InlineData(1.0, "-", false)] + [InlineData(-1.0, "-", false)] + public void ScalarAdditionBelowOverflowMidpointRemainsFinite(double sign, string operation, bool scalarLeft) + { + DoubleDouble pair = DoubleDouble.FromComponents(sign * double.MaxValue, sign * Math.ScaleB(1.0, 969)); + double scalar = sign * Math.BitDecrement(Math.ScaleB(1.0, 969)); + if (operation == "-") + { + scalar = -scalar; + } + Rational expected = Expected(Exact(pair), Exact(scalar), operation); + DoubleDouble actual = operation == "-" ? pair - scalar : scalarLeft ? scalar + pair : pair + scalar; + AssertAccurate(actual, expected, Describe(pair, new DoubleDouble(scalar), operation)); + } + + [Theory] + [InlineData(1.0)] + [InlineData(-1.0)] + public void MultiplicationBelowOverflowMidpointRemainsFinite(double sign) + { + // Exact magnitude = MaxValue + 2^970 - 3*2^914. + DoubleDouble left = DoubleDouble.FromComponents(sign * double.MaxValue, sign * Math.ScaleB(1.0, 969)); + DoubleDouble right = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -55)); + BelowOverflowMidpoint(Exact(left) * Exact(right)).ShouldBeTrue(); + AssertOperation(left, right, "*"); + } + + [Theory] + [InlineData(1.0)] + [InlineData(-1.0)] + public void DivisionWithLowNumeratorBelowOverflowMidpointRemainsFinite(double sign) + { + DoubleDouble left = DoubleDouble.FromComponents(sign * double.MaxValue, sign * Math.ScaleB(1.0, 969)); + DoubleDouble right = DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -55)); + Rational expected = Exact(left) / Exact(right); + BelowOverflowMidpoint(expected).ShouldBeTrue(); + AssertAccurate(left / right, expected, Describe(left, right, "/")); + } + + [Theory] + [InlineData(1.0, false)] + [InlineData(-1.0, false)] + [InlineData(1.0, true)] + [InlineData(-1.0, true)] + public void DivisionBelowOverflowMidpointRemainsFinite(double sign, bool scalarLeft) + { + DoubleDouble left = new(sign * double.MaxValue); + DoubleDouble right = DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -54)); + Rational expected = Exact(left) / Exact(right); + BelowOverflowMidpoint(expected).ShouldBeTrue(); + DoubleDouble actual = scalarLeft ? (sign * double.MaxValue) / right : left / right; + AssertAccurate(actual, expected, Describe(left, right, "/")); + } + + [Fact] + public void FactoryAvoidsIntermediateOverflowForFiniteOppositeSignSums() + { + // With U = 2^971 and M = MaxValue, M - 1.5U rounds to M - U. + // Smaller-first TwoSum computes (M - U) - (-1.5U) = M + 0.5U, + // which overflows even though the original exact sum is finite. + foreach (double sign in new[] { -1.0, 1.0 }) + { + double small = -sign * Math.ScaleB(3.0, 970); + double large = sign * double.MaxValue; + foreach (DoubleDouble actual in new[] { DoubleDouble.FromComponents(small, large), + DoubleDouble.FromComponents(large, small) }) + { + actual.High.ShouldBe(sign * Math.BitDecrement(double.MaxValue)); + actual.Low.ShouldBe(-sign * Math.ScaleB(1.0, 970)); + Exact(actual).CompareTo(Exact(small) + Exact(large)).ShouldBe(0); + AssertNormalized(actual); + } + } + } + + [Fact] + public void FactoryPreservesExactFiniteSumsAcrossExponentBoundaries() + { + double[] components = + [ + 0.0, -0.0, double.Epsilon, -double.Epsilon, + Math.BitDecrement(Math.ScaleB(1.0, -1022)), Math.ScaleB(1.0, -1022), + Math.ScaleB(1.0, -969), Math.ScaleB(1.0, -53), 1.0, + Math.BitIncrement(1.0), Math.ScaleB(1.0, 970), double.MaxValue, + Math.ScaleB(3.0, 970), -Math.ScaleB(3.0, 970), + -Math.ScaleB(1.0, -1022), -1.0, -double.MaxValue + ]; + foreach (double high in components) + { + foreach (double low in components) + { + Rational expected = Exact(high) + Exact(low); + if (!BelowOverflowMidpoint(expected)) + { + continue; + } + DoubleDouble actual = DoubleDouble.FromComponents(high, low); + Exact(actual).CompareTo(expected).ShouldBe(0, $"factory ({high:R}, {low:R})"); + AssertNormalized(actual); + DoubleDouble repeated = DoubleDouble.FromComponents(actual.High, actual.Low); + repeated.Equals(actual).ShouldBeTrue(); + repeated.GetHashCode().ShouldBe(actual.GetHashCode()); + } + } + } + + [Theory] + [InlineData("+")] + [InlineData("-")] + [InlineData("*")] + [InlineData("/")] + public void ArithmeticAcrossFastPathTransitionsMeetsExactRationalBound(string operation) + { + int[] exponents = [-1074, -1022, -970, -901, -900, -899, -451, -450, -449, + -54, -1, 0, 1, 54, 449, 450, 451, 899, 900, 901, 969, 1020, 1021, 1023]; + Random random = new(0x5eed); + foreach (int leftExponent in exponents) + { + foreach (int rightExponent in exponents) + { + for (int sample = 0; sample < 4; ++sample) + { + DoubleDouble left = Sample(random, leftExponent); + DoubleDouble right = Sample(random, rightExponent); + AssertOperation(left, right, operation); + // Exercise scalar overloads independently, not via equality + // with the corresponding potentially faulty DD operation. + AssertScalarOperations(left, right.High, operation); + } + } + } + } + + [Fact] + public void CancellationAcrossBinadesRetainsSmallResiduals() + { + int[] exponents = [-1022, -969, -450, 0, 450, 969, 1020, 1023]; + foreach (int exponent in exponents) + { + double high = Math.ScaleB(1.0, exponent); + foreach (int gap in new[] { 53, 54, 105, 106, 107, 200, 1000 }) + { + double low = Math.ScaleB(1.0, exponent - gap); + DoubleDouble left = DoubleDouble.FromComponents(high, low); + DoubleDouble right = DoubleDouble.FromComponents(-high, Math.ScaleB(1.0, exponent - gap - 54)); + Rational expected = Exact(left) + Exact(right); + AssertAccurate(left + right, expected, Describe(left, right, "+")); + AssertAccurate(right + left, expected, Describe(right, left, "+")); + AssertAccurate(left - (-right), expected, Describe(left, -right, "-")); + } + } + } + + [Fact] + public void FiniteComparisonsAgreeWithExactValuesAtAdjacentHighMidpoints() + { + List values = [new(0.0), new(-0.0)]; + foreach (int exponent in new[] { -1022, -970, -450, 0, 450, 970, 1023 }) + { + double high = Math.ScaleB(1.0, exponent); + double adjacent = Math.BitIncrement(high); + double midpointLow = Math.ScaleB(1.0, exponent - 53); + foreach (double sign in new[] { -1.0, 1.0 }) + { + values.Add(new DoubleDouble(sign * high)); + values.Add(DoubleDouble.FromComponents(sign * high, sign * midpointLow)); + values.Add(DoubleDouble.FromComponents(sign * adjacent, -sign * midpointLow)); + values.Add(DoubleDouble.FromComponents(sign * high, sign * Math.BitDecrement(midpointLow))); + values.Add(DoubleDouble.FromComponents(sign * high, sign * Math.BitIncrement(midpointLow))); + } + } + foreach (DoubleDouble left in values) + { + foreach (DoubleDouble right in values) + { + int order = Exact(left).CompareTo(Exact(right)); + string context = Describe(left, right, "compare"); + Math.Sign(left.CompareTo(right)).ShouldBe(Math.Sign(order), context); + (left < right).ShouldBe(order < 0, context); + (left > right).ShouldBe(order > 0, context); + (left <= right).ShouldBe(order <= 0, context); + (left >= right).ShouldBe(order >= 0, context); + (left == right).ShouldBe(order == 0, context); + (left != right).ShouldBe(order != 0, context); + left.Equals(right).ShouldBe(order == 0, context); + if (order == 0) + { + left.GetHashCode().ShouldBe(right.GetHashCode(), context); + } + } + } + } + + [Fact] + public void NonzeroUnderflowRetainsResultSign() + { + foreach (double sign in new[] { -1.0, 1.0 }) + { + DoubleDouble tiny = new(sign * double.Epsilon); + DoubleDouble[] zeros = [tiny * 0.25, 0.25 * tiny, tiny / 4.0, + tiny * new DoubleDouble(0.25), tiny / new DoubleDouble(4.0), + (sign * double.Epsilon) / new DoubleDouble(4.0)]; + foreach (DoubleDouble zero in zeros) + { + BitConverter.DoubleToInt64Bits(zero.High).ShouldBe(sign < 0.0 ? long.MinValue : 0L); + BitConverter.DoubleToInt64Bits(zero.Low).ShouldBe(0L); + } + } + } + + private static DoubleDouble Sample(Random random, int exponent) + { + double sign = random.Next(2) == 0 ? -1.0 : 1.0; + double high = Math.ScaleB(sign * (1.0 + random.NextDouble()), exponent); + double low = Math.ScaleB((random.NextDouble() * 2.0) - 1.0, exponent - random.Next(53, 121)); + return DoubleDouble.FromComponents(high, low); + } + + private static void AssertOperation(DoubleDouble left, DoubleDouble right, string operation) + { + Rational expected = Expected(Exact(left), Exact(right), operation); + if (!BelowOverflowMidpoint(expected)) + { + return; + } + DoubleDouble actual = operation switch + { + "+" => left + right, + "-" => left - right, + "*" => left * right, + "/" => left / right, + _ => throw new ArgumentOutOfRangeException(nameof(operation)) + }; + AssertAccurate(actual, expected, Describe(left, right, operation)); + AssertNormalized(actual); + } + + private static void AssertScalarOperations(DoubleDouble left, double right, string operation) + { + Rational forward = Expected(Exact(left), Exact(right), operation); + Rational reverse = Expected(Exact(right), Exact(left), operation); + if (BelowOverflowMidpoint(forward)) + { + DoubleDouble actual = operation switch + { + "+" => left + right, + "-" => left - right, + "*" => left * right, + "/" => left / right, + _ => throw new ArgumentOutOfRangeException(nameof(operation)) + }; + AssertAccurate(actual, forward, Describe(left, new DoubleDouble(right), operation)); + AssertNormalized(actual); + } + if (BelowOverflowMidpoint(reverse)) + { + DoubleDouble actual = operation switch + { + "+" => right + left, + "-" => right - left, + "*" => right * left, + "/" => right / left, + _ => throw new ArgumentOutOfRangeException(nameof(operation)) + }; + AssertAccurate(actual, reverse, Describe(new DoubleDouble(right), left, operation)); + AssertNormalized(actual); + } + } + + private static Rational Expected(Rational left, Rational right, string operation) + { + return operation switch + { + "+" => left + right, + "-" => left - right, + "*" => left * right, + "/" => left / right, + _ => throw new ArgumentOutOfRangeException(nameof(operation)) + }; + } + + private static bool BelowOverflowMidpoint(Rational value) + { + return value.Abs().CompareTo(Exact(double.MaxValue) + Exact(Math.ScaleB(1.0, 970))) < 0; + } + + private static void AssertNormalized(DoubleDouble value) + { + double.IsFinite(value.High).ShouldBeTrue(); + double.IsFinite(value.Low).ShouldBeTrue(); + (value.High + value.Low).ShouldBe(value.High); + if (value.Low == 0.0) + { + BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(0L); + } + } + + private static string Describe(DoubleDouble left, DoubleDouble right, string operation) + { + return $"({left.High:R}, {left.Low:R}) {operation} ({right.High:R}, {right.Low:R})"; + } + + private static void AssertAccurate(DoubleDouble actual, Rational expected, string context) + { + string diagnostic = $"{context}: actual ({actual.High:R}, {actual.Low:R})"; + double.IsFinite(actual.High).ShouldBeTrue(diagnostic); + double.IsFinite(actual.Low).ShouldBeTrue(diagnostic); + Rational error = (Exact(actual) - expected).Abs(); + // Conservative contract, not a correct-rounding assertion. All + // comparisons, including the subnormal floor, use exact rationals. + Rational tolerance = (expected.Abs() * new Rational(1, BigInteger.One << 100)) + + Exact(double.Epsilon); + error.CompareTo(tolerance).ShouldBeLessThanOrEqualTo(0, diagnostic); + } + + private static Rational Exact(DoubleDouble value) + { + return Exact(value.High) + Exact(value.Low); + } + + private static Rational Exact(double value) + { + // Decode IEEE-754 directly; no production helpers or conversions. + double.IsFinite(value).ShouldBeTrue(); + ulong bits = BitConverter.DoubleToUInt64Bits(value); + int biasedExponent = (int)((bits >> 52) & 0x7ff); + BigInteger significand = bits & 0x000f_ffff_ffff_ffffUL; + int exponent = -1074; + if (biasedExponent != 0) + { + significand += BigInteger.One << 52; + exponent = biasedExponent - 1075; + } + if ((bits >> 63) != 0) + { + significand = -significand; + } + return exponent >= 0 + ? new Rational(significand << exponent, BigInteger.One) + : new Rational(significand, BigInteger.One << -exponent); + } + + private readonly struct Rational + { + private readonly BigInteger _numerator; + private readonly BigInteger _denominator; + + public Rational(BigInteger numerator, BigInteger denominator) + { + if (denominator.IsZero) + { + throw new DivideByZeroException(); + } + BigInteger divisor = BigInteger.GreatestCommonDivisor(numerator, denominator); + _numerator = numerator / divisor * denominator.Sign; + _denominator = BigInteger.Abs(denominator / divisor); + } + + public Rational Abs() + { + return new Rational(BigInteger.Abs(_numerator), _denominator); + } + + public int CompareTo(Rational other) + { + return (_numerator * other._denominator).CompareTo(other._numerator * _denominator); + } + + public static Rational operator +(Rational left, Rational right) + { + return new Rational((left._numerator * right._denominator) + (right._numerator * left._denominator), + left._denominator * right._denominator); + } + + public static Rational operator -(Rational value) + { + return new Rational(-value._numerator, value._denominator); + } + + public static Rational operator -(Rational left, Rational right) + { + return left + (-right); + } + + public static Rational operator *(Rational left, Rational right) + { + return new Rational(left._numerator * right._numerator, left._denominator * right._denominator); + } + + public static Rational operator /(Rational left, Rational right) + { + return new Rational(left._numerator * right._denominator, left._denominator * right._numerator); + } + } +} diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleComparisonTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleComparisonTests.cs new file mode 100644 index 0000000..022696a --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleComparisonTests.cs @@ -0,0 +1,107 @@ +using Shouldly; +using Xunit; + +namespace Just.PreciseMath.Tests; + +public class DoubleDoubleComparisonTests +{ + [Fact] + public void RelationalOperatorsCoverEqualHighComponentsAndSpecialValues() + { + // Explicit numerical order, including opposite low-component signs and + // equivalent signed zeros. NaNs are tested separately as unordered. + DoubleDouble[] ordered = + [ + new(double.NegativeInfinity), new(-double.MaxValue), + DoubleDouble.FromComponents(-1.0, -double.Epsilon), new(-1.0), DoubleDouble.FromComponents(-1.0, double.Epsilon), + new(-double.Epsilon), new(-0.0), new(0.0), new(double.Epsilon), + DoubleDouble.FromComponents(1.0, -double.Epsilon), new(1.0), DoubleDouble.FromComponents(1.0, double.Epsilon), + new(double.MaxValue), new(double.PositiveInfinity) + ]; + for (int i = 0; i < ordered.Length; ++i) + { + for (int j = 0; j < ordered.Length; ++j) + { + bool bothZero = (i is 6 or 7) && (j is 6 or 7); + (ordered[i] < ordered[j]).ShouldBe(i < j && !bothZero); + (ordered[i] > ordered[j]).ShouldBe(i > j && !bothZero); + (ordered[i] <= ordered[j]).ShouldBe(i <= j || bothZero); + (ordered[i] >= ordered[j]).ShouldBe(i >= j || bothZero); + } + } + } + + [Fact] + public void OrderingUsesLowComponentAndSupportsCollections() + { + DoubleDouble one = DoubleDouble.One; + DoubleDouble above = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -100)); + (above > one).ShouldBeTrue(); + (one < above).ShouldBeTrue(); + (one <= above).ShouldBeTrue(); + (above >= one).ShouldBeTrue(); + above.CompareTo(one).ShouldBeGreaterThan(0); + new SortedSet { above, one }.Count.ShouldBe(2); + ((IComparable)one).CompareTo(null).ShouldBe(1); + Should.Throw(() => ((IComparable)one).CompareTo("1")); + } + + [Fact] + public void DistinctLowComponentsAreNotCollapsedByOrderingOrCollections() + { + // Review-1 §3: a = (1, 0) and b = (1, 1e-30) are distinct values. The former + // high-only CompareTo reported them equal, so a SortedSet retained only one. + DoubleDouble a = DoubleDouble.One; + DoubleDouble b = DoubleDouble.FromComponents(1.0, 1e-30); + b.High.ShouldBe(1.0); + b.Low.ShouldBe(1e-30); + (b > a).ShouldBeTrue(); + (a < b).ShouldBeTrue(); + b.CompareTo(a).ShouldBeGreaterThan(0); + b.Equals(a).ShouldBeFalse(); + new SortedSet { b, a }.Count.ShouldBe(2); + } + + [Fact] + public void NaNOperatorsAreUnorderedButCompareToProvidesTotalOrder() + { + DoubleDouble nan = DoubleDouble.NaN; + foreach (DoubleDouble other in new[] { nan, DoubleDouble.Zero, new DoubleDouble(-0.0), + new DoubleDouble(double.NegativeInfinity), new DoubleDouble(double.PositiveInfinity), + DoubleDouble.FromComponents(1.0, double.Epsilon), DoubleDouble.FromComponents(-1.0, -double.Epsilon) }) + { + (nan < other).ShouldBeFalse(); + (nan > other).ShouldBeFalse(); + (nan <= other).ShouldBeFalse(); + (nan >= other).ShouldBeFalse(); + (other < nan).ShouldBeFalse(); + (other > nan).ShouldBeFalse(); + (other <= nan).ShouldBeFalse(); + (other >= nan).ShouldBeFalse(); + } + nan.CompareTo(nan).ShouldBe(0); + nan.CompareTo(DoubleDouble.Zero).ShouldBeLessThan(0); + } + + [Theory] + [InlineData(1.0, 0.0)] + [InlineData(0.0, -1.0)] + [InlineData(double.NaN, 0.0)] + [InlineData(0.0, double.NaN)] + [InlineData(double.PositiveInfinity, 0.0)] + [InlineData(0.0, double.NegativeInfinity)] + [InlineData(double.PositiveInfinity, double.NegativeInfinity)] + public void ClassificationFollowsCanonicalHigh(double high, double low) + { + DoubleDouble value = DoubleDouble.FromComponents(high, low); + DoubleDouble.IsNaN(value).ShouldBe(double.IsNaN(value.High)); + DoubleDouble.IsFinite(value).ShouldBe(double.IsFinite(value.High)); + DoubleDouble.IsInfinity(value).ShouldBe(double.IsInfinity(value.High)); + DoubleDouble.IsPositiveInfinity(value).ShouldBe(double.IsPositiveInfinity(value.High)); + DoubleDouble.IsNegativeInfinity(value).ShouldBe(double.IsNegativeInfinity(value.High)); + DoubleDouble.IsNegative(value).ShouldBe(double.IsNegative(value.High)); + value.Decompose(out double actualHigh, out double actualLow); + actualHigh.ShouldBe(value.High); + actualLow.ShouldBe(value.Low); + } +} diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleConversionTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleConversionTests.cs new file mode 100644 index 0000000..5b00d46 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleConversionTests.cs @@ -0,0 +1,266 @@ +using System.Numerics; +using Shouldly; +using Xunit; + +namespace Just.PreciseMath.Tests; + +public class DoubleDoubleConversionTests +{ + [Theory] + [InlineData(0.0, 0.0, false)] + [InlineData(1.0, -1.0, false)] + [InlineData(0.0, double.Epsilon, true)] + [InlineData(0.0, -double.Epsilon, true)] + [InlineData(double.NaN, 0.0, true)] + [InlineData(double.PositiveInfinity, 0.0, true)] + [InlineData(double.NegativeInfinity, 0.0, true)] + public void BooleanConversionUsesNormalizedZero(double high, double low, bool expected) + { + IConvertible value = DoubleDouble.FromComponents(high, low); + value.ToBoolean(null).ShouldBe(expected); + value.ToType(typeof(bool), null).ShouldBe(expected); + } + + [Theory] + [InlineData(0.0)] + [InlineData(double.Epsilon)] + [InlineData(-double.Epsilon)] + [InlineData(double.MaxValue)] + [InlineData(-double.MaxValue)] + [InlineData(double.PositiveInfinity)] + [InlineData(double.NegativeInfinity)] + [InlineData(double.NaN)] + public void SingleComponentFloatConversionsMatchBinary64Casts(double high) + { + DoubleDouble value = new(high); + int expectedBits = BitConverter.SingleToInt32Bits((float)high); + BitConverter.SingleToInt32Bits((float)value).ShouldBe(expectedBits); + BitConverter.SingleToInt32Bits(((IConvertible)value).ToSingle(null)).ShouldBe(expectedBits); + } + + [Fact] + public void SingleComponentFloatConversionDoesNotAllocate() + { + DoubleDouble value = new(1.25); + float result = 0.0f; + for (int i = 0; i < 100; ++i) + { + result = (float)value; + } + long before = GC.GetAllocatedBytesForCurrentThread(); + for (int i = 0; i < 100; ++i) + { + result = (float)value; + } + long allocated = GC.GetAllocatedBytesForCurrentThread() - before; + result.ShouldBe(1.25f); + allocated.ShouldBe(0L); + } + + [Fact] + public void DecimalConstructionPreservesSignedScaledZero() + { + decimal zero = new(0, 0, 0, true, 28); + DoubleDouble value = new(zero); + BitConverter.DoubleToInt64Bits(value.High).ShouldBe(long.MinValue); + BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(0L); + ((IConvertible)value).ToBoolean(null).ShouldBeFalse(); + } + + [Fact] + public void IntegerInputsRemainExactAcrossBinary64RoundingBoundaries() + { + foreach (long input in new[] { 0L, 1L, -1L, long.MinValue, long.MinValue + 1, long.MaxValue - 1, long.MaxValue }) + { + CheckIntegerInput(input); + } + for (int exponent = 53; exponent < 63; ++exponent) + { + long center = 1L << exponent; + long halfUlp = 1L << (exponent - 53); + foreach (long offset in new[] { -halfUlp - 1, -halfUlp, -halfUlp + 1, halfUlp - 1, halfUlp, halfUlp + 1 }) + { + CheckIntegerInput(center + offset); + CheckIntegerInput(-center - offset); + } + } + Random random = new(1729); + for (int i = 0; i < 250; ++i) + { + CheckIntegerInput(random.NextInt64(long.MinValue, long.MaxValue)); + } + } + + [Fact] + public void IntegerConstructionDoesNotAllocate() + { + DoubleDouble value = default; + for (int i = 0; i < 100; ++i) + { + value = new DoubleDouble(long.MaxValue); + } + long before = GC.GetAllocatedBytesForCurrentThread(); + for (int i = 0; i < 100; ++i) + { + value = new DoubleDouble(long.MaxValue); + } + long allocated = GC.GetAllocatedBytesForCurrentThread() - before; + value.High.ShouldBe(Math.ScaleB(1.0, 63)); + value.Low.ShouldBe(-1.0); + allocated.ShouldBe(0L); + } + + private static void CheckIntegerInput(long input) + { + foreach (DoubleDouble value in new[] { new DoubleDouble(input), (DoubleDouble)input }) + { + // Both components of an integer input are integers. BigInteger + // recombines them exactly without rounding the sum to binary64. + (new BigInteger(value.High) + new BigInteger(value.Low)).ShouldBe(new BigInteger(input)); + value.High.ShouldBe((double)input); + ((long)value).ShouldBe(input); + } + } + + [Fact] + public void IntegerInputsPreserveEveryBit() + { + DoubleDouble value = (DoubleDouble)9007199254740993L; + value.High.ShouldBe(9007199254740992.0); + value.Low.ShouldBe(1.0); + new DoubleDouble(long.MaxValue).Low.ShouldBe(-1.0); + new DoubleDouble(long.MinValue).Low.ShouldBe(0.0); + ((DoubleDouble)int.MaxValue).High.ShouldBe(2147483647.0); + new DoubleDouble(1).High.ShouldBe(1.0); + } + + [Fact] + public void BinaryConversionsRoundOnceUsingBothComponents() + { + double midpoint = 1.0 + Math.ScaleB(1.0, -24); + ((float)DoubleDouble.FromComponents(midpoint, Math.ScaleB(1.0, -80))).ShouldBe(MathF.BitIncrement(1.0f)); + ((float)DoubleDouble.FromComponents(midpoint, -Math.ScaleB(1.0, -80))).ShouldBe(1.0f); + ((float)new DoubleDouble(midpoint)).ShouldBe(1.0f); + ((float)DoubleDouble.FromComponents(Math.ScaleB(1.0, -150), double.Epsilon)).ShouldBe(float.Epsilon); + ((float)new DoubleDouble(Math.ScaleB(1.0, -150))).ShouldBe(0.0f); + ((float)new DoubleDouble(double.MaxValue)).ShouldBe(float.PositiveInfinity); + float.IsNaN((float)DoubleDouble.NaN).ShouldBeTrue(); + ((double)new DoubleDouble(double.NegativeInfinity)).ShouldBe(double.NegativeInfinity); + ((double)DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54))).ShouldBe(1.0); + ((DoubleDouble)double.Epsilon).High.ShouldBe(double.Epsilon); + ((DoubleDouble)float.Epsilon).High.ShouldBe(Math.ScaleB(1.0, -149)); + BitConverter.DoubleToInt64Bits((double)(DoubleDouble)(-0.0)).ShouldBe(long.MinValue); + BitConverter.SingleToInt32Bits((float)(DoubleDouble)(-0.0f)).ShouldBe(int.MinValue); + } + + [Fact] + public void DecimalInputComputesTheExactBinaryResidual() + { + DoubleDouble tenth = new(0.1m); + tenth.High.ShouldBe(0.1); + // 1/10 - binary64(0.1) = -1/(5 * 2^55), rounded to binary64. + tenth.Low.ShouldBe(-5.551115123125783e-18); + DoubleDouble maximum = (DoubleDouble)decimal.MaxValue; + maximum.High.ShouldBe(Math.ScaleB(1.0, 96)); + maximum.Low.ShouldBe(-1.0); + ((decimal)maximum).ShouldBe(decimal.MaxValue); + ((decimal)new DoubleDouble(decimal.MinValue)).ShouldBe(decimal.MinValue); + ((decimal)new DoubleDouble(0.0000000000000000000000000001m)).ShouldBe(0.0000000000000000000000000001m); + ((decimal)tenth).ShouldBe(0.1m); + } + + [Fact] + public void DecimalOutputRoundsExactSumToNearestEven() + { + ((decimal)DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54))).ShouldBe(1.0000000000000000555111512313m); + ((decimal)new DoubleDouble(Math.ScaleB(1.0, -29))).ShouldBe(0.0000000018626451492309570312m); + ((decimal)DoubleDouble.FromComponents(Math.ScaleB(1.0, -29), double.Epsilon)).ShouldBe(0.0000000018626451492309570313m); + ((decimal)new DoubleDouble(double.Epsilon)).ShouldBe(0m); + Should.Throw(() => (decimal)new DoubleDouble(Math.ScaleB(1.0, 96))); + Should.Throw(() => (decimal)DoubleDouble.NaN); + Should.Throw(() => (decimal)new DoubleDouble(double.NegativeInfinity)); + } + + [Fact] + public void ConvertibleIntegersUseNearestEvenAndCheckAllTargetRanges() + { + IConvertible value = DoubleDouble.FromComponents(2.5, double.Epsilon); + value.ToByte(null).ShouldBe((byte)3); + value.ToSByte(null).ShouldBe((sbyte)3); + value.ToInt16(null).ShouldBe((short)3); + value.ToUInt16(null).ShouldBe((ushort)3); + value.ToInt32(null).ShouldBe(3); + value.ToUInt32(null).ShouldBe(3U); + value.ToInt64(null).ShouldBe(3L); + value.ToUInt64(null).ShouldBe(3UL); + ((IConvertible)new DoubleDouble(2.5)).ToInt32(null).ShouldBe(2); + ((IConvertible)new DoubleDouble(-2.5)).ToInt32(null).ShouldBe(-2); + ((IConvertible)DoubleDouble.FromComponents(-2.5, -double.Epsilon)).ToInt32(null).ShouldBe(-3); + ((IConvertible)DoubleDouble.FromComponents(Math.ScaleB(1.0, 64), -1.0)).ToUInt64(null).ShouldBe(ulong.MaxValue); + Should.Throw(() => ((IConvertible)new DoubleDouble(255.5)).ToByte(null)); + Should.Throw(() => ((IConvertible)new DoubleDouble(127.5)).ToSByte(null)); + Should.Throw(() => ((IConvertible)new DoubleDouble(32767.5)).ToInt16(null)); + Should.Throw(() => ((IConvertible)new DoubleDouble(65535.5)).ToUInt16(null)); + Should.Throw(() => ((IConvertible)new DoubleDouble(2147483647.5)).ToInt32(null)); + Should.Throw(() => ((IConvertible)new DoubleDouble(4294967295.5)).ToUInt32(null)); + Should.Throw(() => ((IConvertible)new DoubleDouble(Math.ScaleB(1.0, 63))).ToInt64(null)); + Should.Throw(() => ((IConvertible)new DoubleDouble(Math.ScaleB(1.0, 64))).ToUInt64(null)); + Should.Throw(() => ((IConvertible)new DoubleDouble(-1.0)).ToUInt64(null)); + Should.Throw(() => ((IConvertible)DoubleDouble.NaN).ToInt32(null)); + } + + [Fact] + public void ConvertibleDispatchPreservesTypeAndUsesNumericPolicies() + { + IConvertible value = new DoubleDouble(1.25); + value.GetTypeCode().ShouldBe(TypeCode.Object); + value.ToBoolean(null).ShouldBeTrue(); + ((IConvertible)DoubleDouble.Zero).ToBoolean(null).ShouldBeFalse(); + ((IConvertible)DoubleDouble.NaN).ToBoolean(null).ShouldBeTrue(); + value.ToDecimal(null).ShouldBe(1.25m); + value.ToDouble(null).ShouldBe(1.25); + value.ToSingle(null).ShouldBe(1.25f); + value.ToString(System.Globalization.CultureInfo.InvariantCulture).ShouldBe("1.25"); + value.ToType(typeof(DoubleDouble), null).ShouldBe(new DoubleDouble(1.25)); + value.ToType(typeof(object), null).ShouldBe(new DoubleDouble(1.25)); + value.ToType(typeof(int), null).ShouldBe(1); + value.ToType(typeof(string), System.Globalization.CultureInfo.InvariantCulture).ShouldBe("1.25"); + value.ToType(typeof(decimal), null).ShouldBe(1.25m); + value.ToType(typeof(double), null).ShouldBe(1.25); + value.ToType(typeof(float), null).ShouldBe(1.25f); + value.ToType(typeof(bool), null).ShouldBe(true); + Should.Throw(() => value.ToChar(null)); + Should.Throw(() => value.ToDateTime(null)); + Should.Throw(() => value.ToType(typeof(Guid), null)); + Should.Throw(() => value.ToType(typeof(DayOfWeek), null)); + Should.Throw(() => value.ToType(null!, null)); + } + + [Fact] + public void ChangeTypeRecognizesDoubleDoubleAndDelegatesNumericTargets() + { + // Review-1 §10: Convert.ChangeType(One, typeof(DoubleDouble)) threw + // InvalidCastException because ToType did not recognize its own type. + System.Globalization.CultureInfo invariant = System.Globalization.CultureInfo.InvariantCulture; + Convert.ChangeType(DoubleDouble.One, typeof(DoubleDouble), invariant).ShouldBe(DoubleDouble.One); + Convert.ChangeType(new DoubleDouble(1.25), typeof(int), invariant).ShouldBe(1); + Convert.ChangeType(new DoubleDouble(1.25), typeof(double), invariant).ShouldBe(1.25); + Convert.ChangeType(new DoubleDouble(1.25), typeof(string), invariant).ShouldBe("1.25"); + Should.Throw(() => Convert.ChangeType(DoubleDouble.One, typeof(Guid), invariant)); + } + + [Fact] + public void ExplicitIntegersTruncateTheCompleteExpansion() + { + ((int)DoubleDouble.FromComponents(1.0, -1e-30)).ShouldBe(0); + ((int)DoubleDouble.FromComponents(-1.0, 1e-30)).ShouldBe(0); + ((long)DoubleDouble.FromComponents(9007199254740992.0, 1.0)).ShouldBe(9007199254740993L); + ((long)DoubleDouble.FromComponents(9223372036854775808.0, -1.0)).ShouldBe(long.MaxValue); + ((long)new DoubleDouble(-9223372036854775808.0)).ShouldBe(long.MinValue); + ((int)DoubleDouble.FromComponents(2147483648.0, -0.25)).ShouldBe(int.MaxValue); + Should.Throw(() => (int)new DoubleDouble(2147483648.0)); + Should.Throw(() => (long)new DoubleDouble(9223372036854775808.0)); + Should.Throw(() => (int)DoubleDouble.NaN); + Should.Throw(() => (long)new DoubleDouble(double.PositiveInfinity)); + } +} diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleFormattingTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleFormattingTests.cs new file mode 100644 index 0000000..d801b56 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleFormattingTests.cs @@ -0,0 +1,119 @@ +using System.Globalization; +using Shouldly; +using Xunit; + +namespace Just.PreciseMath.Tests; + +public class DoubleDoubleFormattingTests +{ + [Fact] + public void GeneralFormattingTrimsOnlyFractionalZerosAtMaximumPrecision() + { + NumberFormatInfo provider = new() { NumberDecimalSeparator = "::" }; + new DoubleDouble(1.25).ToString("G999", provider).ShouldBe("1::25"); + new DoubleDouble(1000.0).ToString("G999", provider).ShouldBe("1000"); + new DoubleDouble(-0.0).ToString("G999", provider).ShouldBe("-0"); + new DoubleDouble(999.5).ToString("G3", provider).ShouldBe("1E+03"); + } + + [Fact] + public void GeneralFormattingRetainsLowComponentDigits() + { + DoubleDouble value = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54)); + value.ToString("G32", CultureInfo.InvariantCulture).ShouldBe("1.0000000000000000555111512312578"); + value.ToString("G", CultureInfo.InvariantCulture).ShouldBe("1.0000000000000000555111512312578"); + value.ToString("", CultureInfo.InvariantCulture).ShouldBe("1.0000000000000000555111512312578"); + ((IFormattable)value).ToString(null, CultureInfo.InvariantCulture).ShouldBe("1.0000000000000000555111512312578"); + new DoubleDouble(100000.0).ToString("G5", CultureInfo.InvariantCulture).ShouldBe("1E+05"); + new DoubleDouble(0.00001m).ToString("g3", CultureInfo.InvariantCulture).ShouldBe("1e-05"); + new DoubleDouble(0.0001m).ToString("G3", CultureInfo.InvariantCulture).ShouldBe("0.0001"); + new DoubleDouble(999.5).ToString("G3", CultureInfo.InvariantCulture).ShouldBe("1E+03"); + } + + [Fact] + public void FormattingDoesNotRouteThroughDecimalOrRestrictExponentRange() + { + // Review-1 §6 verified the former decimal-based formatter threw for 1e100, + // printed "0" for 1e-100, threw for NaN/infinity, and emitted for PI digits + // already wrong at binary64 precision. Expected digits are exact references: + // the PI pair rounded to 32 significant digits, ties to even + // (Python Fraction/Decimal at precision 120), and the exact binary64 values of + // the powers of ten. + DoubleDouble.PI.ToString("G32", CultureInfo.InvariantCulture).ShouldBe("3.1415926535897932384626433832795"); + new DoubleDouble(1e100).ToString("G", CultureInfo.InvariantCulture).ShouldBe("1.0000000000000000159028911097599E+100"); + new DoubleDouble(1e-100).ToString("G", CultureInfo.InvariantCulture).ShouldBe("1.0000000000000000199918998026029E-100"); + // The binary64 values are exactly representable as sums with a zero residual, + // so the G32 text must agree with the BCL formatter for the same value. + new DoubleDouble(1e100).ToString("G32", CultureInfo.InvariantCulture) + .ShouldBe((1e100).ToString("G32", CultureInfo.InvariantCulture)); + new DoubleDouble(1e-100).ToString("G32", CultureInfo.InvariantCulture) + .ShouldBe((1e-100).ToString("G32", CultureInfo.InvariantCulture)); + DoubleDouble.NaN.ToString(CultureInfo.InvariantCulture).ShouldBe("NaN"); + new DoubleDouble(double.PositiveInfinity).ToString(CultureInfo.InvariantCulture).ShouldBe("Infinity"); + new DoubleDouble(double.NegativeInfinity).ToString(CultureInfo.InvariantCulture).ShouldBe("-Infinity"); + } + + [Theory] + [InlineData(2.5, "F0", "2")] + [InlineData(3.5, "F0", "4")] + [InlineData(-2.5, "F0", "-2")] + [InlineData(1.25, "F1", "1.2")] + [InlineData(9.5, "E0", "1E+001")] + [InlineData(0.125, "e2", "1.25e-001")] + [InlineData(0.0, "E2", "0.00E+000")] + [InlineData(-0.0, "F2", "-0.00")] + public void FixedAndExponentialRoundToEven(double value, string format, string expected) + { + new DoubleDouble(value).ToString(format, CultureInfo.InvariantCulture).ShouldBe(expected); + } + + [Fact] + public void FormattingCoversFullBinary64Range() + { + new DoubleDouble(double.MaxValue).ToString("E5", CultureInfo.InvariantCulture).ShouldBe("1.79769E+308"); + new DoubleDouble(double.Epsilon).ToString("G6", CultureInfo.InvariantCulture).ShouldBe("4.94066E-324"); + DoubleDouble.FromComponents(1.0, double.Epsilon).ToString("F324", CultureInfo.InvariantCulture) + .ShouldBe("1." + new string('0', 323) + "5"); + DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54)).ToString("F30", CultureInfo.InvariantCulture) + .ShouldBe("1.000000000000000055511151231258"); + DoubleDouble.FromComponents(2.5, double.Epsilon).ToString("F0", CultureInfo.InvariantCulture).ShouldBe("3"); + } + + [Fact] + [System.Diagnostics.CodeAnalysis.SuppressMessage("Globalization", "CA1305:Specify IFormatProvider", Justification = "Tests intentionally exercise the current-culture overloads and derive their separator from CurrentCulture.")] + public void FormattingUsesProviderForSeparatorsSignsAndSpecialValues() + { + NumberFormatInfo provider = new() + { + NumberDecimalSeparator = ",", + NegativeSign = "minus", + PositiveSign = "plus", + NumberDecimalDigits = 3, + NaNSymbol = "not-number", + PositiveInfinitySymbol = "infinite", + NegativeInfinitySymbol = "minus-infinite" + }; + new DoubleDouble(-1.25).ToString("F", provider).ShouldBe("minus1,250"); + new DoubleDouble(125.0).ToString("E2", provider).ShouldBe("1,25Eplus002"); + new DoubleDouble(-1.25).ToString(provider).ShouldBe("minus1,25"); + DoubleDouble.NaN.ToString("G", provider).ShouldBe("not-number"); + new DoubleDouble(double.PositiveInfinity).ToString("F2", provider).ShouldBe("infinite"); + new DoubleDouble(double.NegativeInfinity).ToString("E", provider).ShouldBe("minus-infinite"); + new DoubleDouble(1.25).ToString("F2", CultureInfo.GetCultureInfo("fr-FR")).ShouldBe("1,25"); + new DoubleDouble(1.25).ToString().ShouldBe("1" + CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator + "25"); + new DoubleDouble(1.25).ToString("F1").ShouldBe("1" + CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator + "2"); + } + + [Theory] + [InlineData("R")] + [InlineData("N2")] + [InlineData("0.00")] + [InlineData("G1000")] + [InlineData("F-1")] + [InlineData("F 2")] + [InlineData("E999999999999999999999")] + public void UnsupportedOrUnboundedFormatsThrow(string format) + { + Should.Throw(() => DoubleDouble.One.ToString(format, CultureInfo.InvariantCulture)); + } +} diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleParsingTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleParsingTests.cs new file mode 100644 index 0000000..8458755 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleParsingTests.cs @@ -0,0 +1,273 @@ +using System.Globalization; +using System.Numerics; +using Shouldly; +using Xunit; + +namespace Just.PreciseMath.Tests; + +public class DoubleDoubleParsingTests +{ + [Theory] + [InlineData("9007199254740993", 9007199254740992.0, 1.0)] + [InlineData("-9007199254740993", -9007199254740992.0, -1.0)] + [InlineData("9.007199254740993e15", 9007199254740992.0, 1.0)] + [InlineData(" +900719925474099300E-2\t", 9007199254740992.0, 1.0)] + [InlineData("1.000000000000000055511151231257827021181583404541015625", 1.0, 5.551115123125783e-17)] + public void ParsingRetainsDecimalInformationBeyondBinary64(string text, double high, double low) + { + // Exact binary cases: 2^53 + 1 and 1 + 2^-54. Neither may be rounded + // through double or decimal before extracting the low component. + AssertParsed(text, CultureInfo.InvariantCulture, high, low); + } + + [Theory] + [InlineData("0.1", 0.1, -5.551115123125783e-18)] + [InlineData("-0.1", -0.1, 5.551115123125783e-18)] + [InlineData("1.23456789012345678901234567890123456789", 1.2345678901234567, 9.858021020478981e-17)] + [InlineData("1e-300", 1e-300, -2.5059094e-317)] + [InlineData("2.4703282292062328e-324", double.Epsilon, 0.0)] + public void NonDyadicDecimalInputsRetainTheRoundedResidual(string text, double high, double low) + { + // Reproduce independently with Python fractions: v = Fraction(text), + // h = float(v), l = float(v - Fraction(h)); canonicalize a zero low to +0. + AssertParsed(text, CultureInfo.InvariantCulture, high, low); + } + + [Theory] + [InlineData("G99")] + [InlineData("E99")] + [InlineData("F99")] + public void SupportedFormatterOutputCanPreserveAnExactlyRepresentablePair(string format) + { + double low = Math.ScaleB(1.0, -80); + DoubleDouble value = DoubleDouble.FromComponents(1.0, low); + CultureInfo culture = CultureInfo.GetCultureInfo("fr-FR"); + // This dyadic has a terminating decimal expansion within the chosen + // precision. No general G32 round-trip guarantee follows from this test. + AssertParsed(value.ToString(format, culture), culture, 1.0, low); + } + + [Theory] + [InlineData(".125", 0.125)] + [InlineData("125.e-3", 0.125)] + [InlineData("00000125E-3", 0.125)] + [InlineData("\r\n 1.25 \t", 1.25)] + [InlineData("0", 0.0)] + [InlineData("-0.000e999999999999999999999", -0.0)] + [InlineData("1e999999999999999999999", double.PositiveInfinity)] + [InlineData("-1e999999999999999999999", double.NegativeInfinity)] + [InlineData("1e-999999999999999999999", 0.0)] + [InlineData("-1e-999999999999999999999", -0.0)] + [InlineData("NaN", double.NaN)] + [InlineData("-nan", double.NaN)] + [InlineData("Infinity", double.PositiveInfinity)] + [InlineData("+infinity", double.PositiveInfinity)] + [InlineData("-Infinity", double.NegativeInfinity)] + public void ParsingHandlesGrammarAndSpecialValues(string text, double high) + { + AssertParsed(text, CultureInfo.InvariantCulture, high, 0.0); + } + + [Theory] + [InlineData("")] + [InlineData(" \t\r\n")] + [InlineData("+")] + [InlineData(".")] + [InlineData("e1")] + [InlineData("1e")] + [InlineData("1e+")] + [InlineData("1e--1")] + [InlineData("1e999999999999999999x")] + [InlineData("0e999999999999999999x")] + [InlineData("--1")] + [InlineData("+ 1")] + [InlineData("1 2")] + [InlineData("1.2.3")] + [InlineData("1e2e3")] + [InlineData("1,000")] + [InlineData("$1")] + [InlineData("(1)")] + [InlineData("0x10")] + [InlineData("1_000")] + [InlineData("１２３")] + [InlineData("1\0")] + [InlineData("NaNx")] + public void InvalidInputFailsWithoutLeavingAPartialResult(string text) + { + DoubleDouble.TryParse(text, CultureInfo.InvariantCulture, out DoubleDouble fromString).ShouldBeFalse(); + DoubleDouble.TryParse(text.AsSpan(), CultureInfo.InvariantCulture, out DoubleDouble fromSpan).ShouldBeFalse(); + AssertPositiveZero(fromString); + AssertPositiveZero(fromSpan); + Should.Throw(() => DoubleDouble.Parse(text, CultureInfo.InvariantCulture)); + Should.Throw(() => DoubleDouble.Parse(text.AsSpan(), CultureInfo.InvariantCulture)); + } + + [Fact] + public void NullStringFailsTryParseAndThrowsArgumentNullFromParse() + { + DoubleDouble.TryParse((string?)null, out DoubleDouble result).ShouldBeFalse(); + AssertPositiveZero(result); + DoubleDouble.TryParse((string?)null, CultureInfo.InvariantCulture, out result).ShouldBeFalse(); + AssertPositiveZero(result); + Should.Throw(() => DoubleDouble.Parse((string)null!, CultureInfo.InvariantCulture)); + } + + [Fact] + public void CultureSuppliesDecimalSeparatorAndBothExponentSigns() + { + NumberFormatInfo info = (NumberFormatInfo)CultureInfo.InvariantCulture.NumberFormat.Clone(); + info.NumberDecimalSeparator = "::"; + info.PositiveSign = "plus"; + info.NegativeSign = "minus"; + AssertParsed("minus1::25eplus2", info, -125.0, 0.0); + AssertParsed("plus125eminus2", info, 1.25, 0.0); + AssertParsed("12,5", CultureInfo.GetCultureInfo("fr-FR"), 12.5, 0.0); + DoubleDouble.TryParse("1.25", info, out DoubleDouble result).ShouldBeFalse(); + AssertPositiveZero(result); + } + + [Fact] + public void CultureSpecialSymbolsAreRecognizedBeforeConsumingTheirSigns() + { + NumberFormatInfo info = (NumberFormatInfo)CultureInfo.InvariantCulture.NumberFormat.Clone(); + info.NaNSymbol = "+missing"; + info.PositiveInfinitySymbol = "-unbounded"; + info.NegativeInfinitySymbol = "negative-limit"; + AssertParsed("+MISSING", info, double.NaN, 0.0); + AssertParsed("-UNBOUNDED", info, double.PositiveInfinity, 0.0); + AssertParsed("NEGATIVE-LIMIT", info, double.NegativeInfinity, 0.0); + } + + [Fact] + public void DefaultTryParseAndNullProvidersUseCurrentCulture() + { + CultureInfo previous = CultureInfo.CurrentCulture; + try + { + CultureInfo.CurrentCulture = CultureInfo.GetCultureInfo("fr-FR"); + const string Text = "1,25"; + DoubleDouble.Parse(Text, null).High.ShouldBe(1.25); + DoubleDouble.Parse(Text.AsSpan(), null).High.ShouldBe(1.25); + DoubleDouble.TryParse(Text, out DoubleDouble fromString).ShouldBeTrue(); + DoubleDouble.TryParse(Text.AsSpan(), out DoubleDouble fromSpan).ShouldBeTrue(); + fromString.High.ShouldBe(1.25); + fromSpan.High.ShouldBe(1.25); + AssertParsed(Text, null, 1.25, 0.0); + } + finally + { + CultureInfo.CurrentCulture = previous; + } + } + + [Fact] + public void GenericStringAndSpanParsingInterfacesAreImplemented() + { + ParseString("9007199254740993").Low.ShouldBe(1.0); + ParseSpan("9007199254740993").Low.ShouldBe(1.0); + } + + [Fact] + public void InputLengthIsBoundedBeforeIgnoringWhitespaceOrLeadingZeros() + { + AssertParsed(new string('0', 4095) + "1", CultureInfo.InvariantCulture, 1.0, 0.0); + AssertParsed("1" + new string('0', 4095), CultureInfo.InvariantCulture, double.PositiveInfinity, 0.0); + InvalidInputFailsWithoutLeavingAPartialResult(new string('0', 4096) + "1"); + InvalidInputFailsWithoutLeavingAPartialResult(new string(' ', 4096) + "1"); + } + + [Fact] + public void LongMantissasCanCancelLargeExponentsWithoutOverflowOrUnderflow() + { + AssertParsed("1" + new string('0', 4000) + "e-4000", CultureInfo.InvariantCulture, 1.0, 0.0); + AssertParsed("0." + new string('0', 4000) + "1e4001", CultureInfo.InvariantCulture, 1.0, 0.0); + } + + [Fact] + public void ExactDyadicsCoverSubnormalAndOverflowRoundingBoundaries() + { + foreach (int sign in new[] { -1, 1 }) + { + // Integer coefficients and powers of two define independent exact inputs. + AssertParsed(DyadicText(sign, -1074), CultureInfo.InvariantCulture, sign * double.Epsilon, 0.0); + AssertParsed(DyadicText(sign, -1075), CultureInfo.InvariantCulture, sign < 0 ? -0.0 : 0.0, 0.0); + AssertParsed(DyadicText(3 * sign, -1076), CultureInfo.InvariantCulture, sign * double.Epsilon, 0.0); + BigInteger maximum = (BigInteger.One << 1024) - (BigInteger.One << 971); + AssertParsed(DyadicText(sign * maximum, 0), CultureInfo.InvariantCulture, sign * double.MaxValue, 0.0); + BigInteger midpoint = maximum + (BigInteger.One << 970); + AssertParsed(DyadicText(sign * midpoint, 0), CultureInfo.InvariantCulture, + sign < 0 ? double.NegativeInfinity : double.PositiveInfinity, 0.0); + // Residual rounding reaches the overflow midpoint although the exact + // input is below it. The adjacent finite pair must be selected. + AssertParsed(DyadicText(sign * (midpoint - (BigInteger.One << 916)), 0), CultureInfo.InvariantCulture, + sign * double.MaxValue, sign * Math.BitDecrement(Math.ScaleB(1.0, 970))); + } + } + + [Theory] + [InlineData(2, false)] + [InlineData(3, true)] + public void LowComponentRoundingUsesExactResidualAndTiesToEven(int tail, bool roundUp) + { + // 1 + 2^-54 + tail*2^-108. At tail=2 the low is at its midpoint; + // at tail=3 it lies above the midpoint. The high remains exactly 1. + BigInteger coefficient = (BigInteger.One << 108) + (BigInteger.One << 54) + tail; + double low = Math.ScaleB(1.0, -54); + AssertParsed(DyadicText(coefficient, -108), CultureInfo.InvariantCulture, + 1.0, roundUp ? Math.BitIncrement(low) : low); + } + + private static T ParseString(string text) where T : IParsable + { + T.TryParse(text, CultureInfo.InvariantCulture, out T? result).ShouldBeTrue(); + result.ShouldBe(T.Parse(text, CultureInfo.InvariantCulture)); + return T.Parse(text, CultureInfo.InvariantCulture); + } + + private static T ParseSpan(ReadOnlySpan text) where T : ISpanParsable + { + T.TryParse(text, CultureInfo.InvariantCulture, out T? result).ShouldBeTrue(); + result.ShouldBe(T.Parse(text, CultureInfo.InvariantCulture)); + return T.Parse(text, CultureInfo.InvariantCulture); + } + + private static string DyadicText(BigInteger coefficient, int exponent) + { + // c*2^-k = (c*5^k)*10^-k. No production conversion or formatting helpers. + return exponent >= 0 + ? (coefficient << exponent).ToString(CultureInfo.InvariantCulture) + : (coefficient * BigInteger.Pow(5, -exponent)).ToString(CultureInfo.InvariantCulture) + + "e" + exponent.ToString(CultureInfo.InvariantCulture); + } + + private static void AssertPositiveZero(DoubleDouble value) + { + BitConverter.DoubleToInt64Bits(value.High).ShouldBe(0L); + BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(0L); + } + + private static void AssertParsed(string text, IFormatProvider? provider, double high, double low) + { + DoubleDouble.TryParse(text, provider, out DoubleDouble fromString).ShouldBeTrue(); + DoubleDouble.TryParse(text.AsSpan(), provider, out DoubleDouble fromSpan).ShouldBeTrue(); + DoubleDouble[] results = [DoubleDouble.Parse(text, provider), DoubleDouble.Parse(text.AsSpan(), provider), + fromString, fromSpan]; + foreach (DoubleDouble result in results) + { + result.High.ShouldBe(high); + result.Low.ShouldBe(low); + if (double.IsFinite(high)) + { + (result.High + result.Low).ShouldBe(result.High); + } + if (high == 0.0) + { + BitConverter.DoubleToInt64Bits(result.High).ShouldBe(BitConverter.DoubleToInt64Bits(high)); + } + if (low == 0.0) + { + BitConverter.DoubleToInt64Bits(result.Low).ShouldBe(0L); + } + } + } +} diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleRepresentationTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleRepresentationTests.cs new file mode 100644 index 0000000..945fb62 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleRepresentationTests.cs @@ -0,0 +1,100 @@ +using Shouldly; +using Xunit; + +namespace Just.PreciseMath.Tests; + +public class DoubleDoubleRepresentationTests +{ + [Fact] + public void ArbitraryComponentsUseThePublicFactoryNotAPublicPairConstructor() + { + typeof(DoubleDouble).GetConstructor([typeof(double), typeof(double)]).ShouldBeNull(); + DoubleDouble value = DoubleDouble.FromComponents(1.0, 1.0); + value.High.ShouldBe(2.0); + value.Low.ShouldBe(0.0); + } + + [Fact] + public void TrustedConstructorPreservesAlreadyNormalizedComponents() + { + DoubleDouble value = new(1.0, Math.ScaleB(1.0, -54)); + value.High.ShouldBe(1.0); + value.Low.ShouldBe(Math.ScaleB(1.0, -54)); + DoubleDouble negativeZero = new(-0.0, 0.0); + BitConverter.DoubleToInt64Bits(negativeZero.High).ShouldBe(long.MinValue); + BitConverter.DoubleToInt64Bits(negativeZero.Low).ShouldBe(0L); + } + + [Fact] + public void NegationPreservesNormalizationAndCanonicalComponents() + { + DoubleDouble[] values = [DoubleDouble.NaN, new(double.PositiveInfinity), new(double.NegativeInfinity), + new(0.0), new(-0.0), new(double.Epsilon), new(-double.Epsilon), + DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54))]; + foreach (DoubleDouble value in values) + { + DoubleDouble negated = -value; + double expectedHigh = DoubleDouble.IsNaN(value) ? double.NaN : -value.High; + double expectedLow = value.Low == 0.0 ? 0.0 : -value.Low; + BitConverter.DoubleToInt64Bits(negated.High).ShouldBe(BitConverter.DoubleToInt64Bits(expectedHigh)); + BitConverter.DoubleToInt64Bits(negated.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expectedLow)); + (-negated).Equals(value).ShouldBeTrue(); + } + } + + [Fact] + public void FactoryCanonicalizesNaNPayloadsAndZeroResidualSigns() + { + double nan = BitConverter.Int64BitsToDouble(0x7ff8000000000001L); + foreach (DoubleDouble value in new[] { new DoubleDouble(nan), DoubleDouble.FromComponents(nan, 0.0), + DoubleDouble.FromComponents(1.0, nan), DoubleDouble.FromComponents(double.PositiveInfinity, double.NegativeInfinity) }) + { + BitConverter.DoubleToInt64Bits(value.High).ShouldBe(BitConverter.DoubleToInt64Bits(double.NaN)); + BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(0L); + } + DoubleDouble negativeZero = DoubleDouble.FromComponents(-0.0, -0.0); + BitConverter.DoubleToInt64Bits(negativeZero.High).ShouldBe(long.MinValue); + BitConverter.DoubleToInt64Bits(negativeZero.Low).ShouldBe(0L); + DoubleDouble cancelled = DoubleDouble.FromComponents(1.0, -1.0); + BitConverter.DoubleToInt64Bits(cancelled.High).ShouldBe(0L); + BitConverter.DoubleToInt64Bits(cancelled.Low).ShouldBe(0L); + } + + [Theory] + [InlineData(0.0, 1.0, 1.0, 0.0)] + [InlineData(1.0, 1.0, 2.0, 0.0)] + [InlineData(1.0, -1.0, 0.0, 0.0)] + [InlineData(1e300, -1e300, 0.0, 0.0)] + [InlineData(double.Epsilon, double.Epsilon, 2 * double.Epsilon, 0.0)] + [InlineData(double.MaxValue, double.MaxValue, double.PositiveInfinity, 0.0)] + [InlineData(1.0, double.PositiveInfinity, double.PositiveInfinity, 0.0)] + [InlineData(double.NegativeInfinity, 1.0, double.NegativeInfinity, 0.0)] + public void FactoryNormalizesComponents(double high, double low, double expectedHigh, double expectedLow) + { + DoubleDouble value = DoubleDouble.FromComponents(high, low); + value.High.ShouldBe(expectedHigh); + value.Low.ShouldBe(expectedLow); + } + + [Fact] + public void NaNHasCanonicalComponentsAndCollectionEquality() + { + DoubleDouble value = DoubleDouble.FromComponents(double.PositiveInfinity, double.NegativeInfinity); + double.IsNaN(value.High).ShouldBeTrue(); + value.Low.ShouldBe(0.0); + value.Equals(DoubleDouble.NaN).ShouldBeTrue(); + value.GetHashCode().ShouldBe(DoubleDouble.NaN.GetHashCode()); + new HashSet { value }.Contains(DoubleDouble.NaN).ShouldBeTrue(); + (value == DoubleDouble.NaN).ShouldBeFalse(); + (value != DoubleDouble.NaN).ShouldBeTrue(); + } + + [Fact] + public void ZeroSignsArePreservedButEqual() + { + DoubleDouble negative = new(-0.0); + BitConverter.DoubleToInt64Bits(negative.High).ShouldBe(long.MinValue); + negative.Equals(DoubleDouble.Zero).ShouldBeTrue(); + negative.GetHashCode().ShouldBe(DoubleDouble.Zero.GetHashCode()); + } +} diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs index ac4fa90..004ea1e 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs @@ -5,6 +5,28 @@ namespace Just.PreciseMath.Tests; public class DoubleDoubleTests { + [Theory] + [InlineData("PI", 3.141592653589793, 1.2246467991473532e-16)] + [InlineData("E", 2.718281828459045, 1.4456468917292502e-16)] + [InlineData("LN2", 0.6931471805599453, 2.3190468138462996e-17)] + public void ConstantsHaveNearestBinary64Residuals(string name, double high, double low) + { + // Each residual is round_binary64(constant - exact_binary64(high)). + // Reproduced with Python decimal at precision 90: e = Decimal(1).exp(), + // ln(2) = Decimal(2).ln(), pi = 16*atan(1/5) - 4*atan(1/239), using + // atan(x) = sum((-1)^k*x^(2k+1)/(2k+1)) until |term| < 1e-95. + // low = float(reference - Decimal.from_float(float(reference))). + DoubleDouble value = name switch + { + "PI" => DoubleDouble.PI, + "E" => DoubleDouble.E, + "LN2" => DoubleDouble.LN2, + _ => throw new ArgumentOutOfRangeException(nameof(name)), + }; + value.High.ShouldBe(high); + value.Low.ShouldBe(low); + } + [Fact] public void OneHasExpectedComponents() { @@ -13,4 +35,52 @@ public class DoubleDoubleTests value.High.ShouldBe(1.0); value.Low.ShouldBe(0.0); } + + [Theory] + [InlineData(1.0, 0.0)] + [InlineData(10.0, 0.0)] + [InlineData(100.0, 0.0)] + [InlineData(-1.0, 0.0)] + [InlineData(-10.0, 0.0)] + [InlineData(-100.0, 0.0)] + [InlineData(3.141592653589793, 1.2246467991473532e-16)] + [InlineData(2.718281828459045, 1.4456468917292502e-16)] + [InlineData(0.6931471805599453, 2.3190468138462996e-17)] + public void AdditiveIdentity(double high, double low) + { + DoubleDouble value = new(high, low); + + DoubleDouble result = value + DoubleDouble.AdditiveIdentity; + DoubleDouble resultInversedOrder = DoubleDouble.AdditiveIdentity + value; + + result.High.ShouldBe(high); + result.Low.ShouldBe(low); + + resultInversedOrder.High.ShouldBe(high); + resultInversedOrder.Low.ShouldBe(low); + } + + [Theory] + [InlineData(1.0, 0.0)] + [InlineData(10.0, 0.0)] + [InlineData(100.0, 0.0)] + [InlineData(-1.0, 0.0)] + [InlineData(-10.0, 0.0)] + [InlineData(-100.0, 0.0)] + [InlineData(3.141592653589793, 1.2246467991473532e-16)] + [InlineData(2.718281828459045, 1.4456468917292502e-16)] + [InlineData(0.6931471805599453, 2.3190468138462996e-17)] + public void MultiplicativeIdentity(double high, double low) + { + DoubleDouble value = new(high, low); + + DoubleDouble result = value * DoubleDouble.MultiplicativeIdentity; + DoubleDouble resultInversedOrder = DoubleDouble.MultiplicativeIdentity * value; + + result.High.ShouldBe(high); + result.Low.ShouldBe(low); + + resultInversedOrder.High.ShouldBe(high); + resultInversedOrder.Low.ShouldBe(low); + } } diff --git a/AGENTS.md b/AGENTS.md index e775bd6..ff75047 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -25,8 +25,10 @@ Do not claim API completeness or accuracy beyond tested contracts. - Use the ignored, repository-local `.hermes/` directory for plans, checklists, investigation notes, and handoff context. Create or update notes as useful; keep them concise and revalidate them against current files. Optional - `.hermes/project-context.md` holds setup and review background. Do not store - secrets, force-add this directory, or make builds or tests depend on it. + `.hermes/README.md` indexes active notes and `.hermes/project-context.md` holds + condensed implementation context. Treat `.hermes/archive/` as historical, not + current instructions or status. Do not store secrets, force-add this directory, + or make builds or tests depend on it. ## Working rules @@ -67,10 +69,11 @@ Follow `.editorconfig`, not incidental style in unfinished code. Mathematically equivalent formulas can round differently. - State and verify algorithm preconditions, especially magnitude ordering for quick-sum transforms, normalization assumptions, and overflow/underflow limits. -- The internal two-component `DoubleDouble` constructor currently does not normalize. - Do not assume arbitrary pairs are canonical. Establish the intended normalization, - NaN, infinity, and signed-zero contracts before changing construction, equality, - hashing, ordering, or classification; keep those operations consistent. +- The internal two-component `DoubleDouble` constructor does not normalize or + validate. Use `FromComponents` for arbitrary pairs and reserve raw construction + for proven normalized/canonical results. Preserve the documented normalization, + NaN, infinity, and signed-zero contracts across construction, equality, hashing, + ordering, and classification; do not change one in isolation. - For affected operations, cover cancellation, widely separated magnitudes, zero and signed zero, subnormals, extreme finite values, infinities, and NaNs. Check intermediate overflow/underflow even when the final result is representable. diff --git a/README.md b/README.md index cb7292b..755e416 100644 --- a/README.md +++ b/README.md @@ -6,15 +6,112 @@ precision than a single `double` while using a fixed-size representation, rather than arbitrary-precision arithmetic. > **Work in progress.** The public API is incomplete and may change. Numerical -> accuracy has not been validated, and the library is not ready for production use. +> contracts are covered by regression tests, not an exhaustive accuracy +> certification. The library is not ready for production use. -## Planned scope +## DoubleDouble core -- Double-double arithmetic, constants, comparisons, conversions, and formatting. -- Common functions including `Abs`, `Sqrt`, `Pow`, `Exp`, and `Log`. -- Correctness tests against higher-precision references and performance benchmarks. +`DoubleDouble` stores a normalized high/low pair. Use `new DoubleDouble(value)` +for a single `double`, or `DoubleDouble.FromComponents(high, low)` for arbitrary +components. The factory normalizes finite sums and canonicalizes NaN/infinity +with a positive-zero low component. The two-component constructor is internal +and performs no normalization or validation; it is reserved for trusted, +already-normalized results. The constants `PI`, `E`, and `LN2` include binary64 +residuals checked against independently computed high-precision values. -These are development goals, not a list of currently supported features. +- Arithmetic: unary `+`/`-`, binary `+`, `-`, `*`, `/`, and both operand orders with + a `double`. Addition retains residuals under cancellation; multiplication uses + fused multiply-add; division uses residual corrections. Mixed `double` operators + use specialized scalar paths rather than promoting the scalar to `DoubleDouble`. + Their finite fast paths normalize once with a final sum transform; scalar + division uses one compensated quotient correction within the error contract below. +- Exponent boundaries: bounded `BigInteger` calculations avoid intermediate + overflow and underflow on the exceptional finite path. Ordinary arithmetic uses + floating-point transforms without allocations. The stored value remains two + doubles; this is not an arbitrary-precision API. +- Comparisons use both components. `Equals` treats NaNs as equal and signed zeros + as equal for collections. `CompareTo` orders NaN before other values. Numerical + equality and relational operators treat NaN as unordered, like `double`. +- Signed zero is preserved by single-value construction and unary negation. + `FromComponents` with a zero low input preserves the high zero's + sign. Exact cancellation of nonzero values yields positive zero. Arithmetic + special values follow binary64 rules. + +Arithmetic is approximate double-double arithmetic, **not a promise of correctly +rounded 106-bit results**. The deterministic rational-oracle tests check a +conservative error bound of `2^-100` relative plus one minimum binary64 subnormal, +with exact component checks for selected representable cases. Near underflow, +extended precision necessarily decreases; overflow produces infinity. Performance +has not been benchmarked, including the allocating exponent-boundary path. + +## Conversions and formatting + +- Explicit conversions support `double`, `float`, `int`, `long`, and `decimal` + in both directions. Integer inputs are exact. Decimal inputs use their exact + coefficient and scale to compute the high component and its residual. +- Integer casts truncate the complete expansion toward zero and throw + `OverflowException` for nonfinite or out-of-range results. `IConvertible` + integer conversions instead round to nearest, ties to even, with range checks. +- Binary32 output rounds the complete expansion directly, including low-component + decisions at midpoints. Decimal output rounds to the greatest fitting scale up + to 28; nonfinite values and magnitudes above `decimal.MaxValue` throw. +- `IConvertible` reports `TypeCode.Object`, supports conversion to itself, and + treats only numerical zero as false. Char, DateTime, and enum conversions are + unsupported and throw `InvalidCastException`. +- `ToString` formats the exact component sum, supports culture-sensitive + `G`/`g`, `E`/`e`, and `F`/`f`, and rounds ties to even. Precision is bounded to + 0–999; other standard and custom formats throw `FormatException`. Default + `G32` is **not** shortest-round-trip formatting. NaN, infinities, and signed + zero are supported without converting through decimal. + +Conversions, parsing, and formatting use allocating `BigInteger` intermediates +where needed to preserve precision; no additional dependency is required. + +## Parsing + +`Parse` and `TryParse` accept strings and `ReadOnlySpan` and implement +`IParsable` / `ISpanParsable`. This initial parser +preserves high/low precision rather than parsing through `double` or `decimal`. +It converts an exact decimal coefficient/exponent into rounded high and residual +components, then normalizes the pair. It does not promise universally correctly +rounded 106-bit results or a general `ToString` round trip. + +```csharp +using System.Globalization; +using Just.PreciseMath; + +DoubleDouble value = DoubleDouble.Parse("9007199254740993", CultureInfo.InvariantCulture); +// value.High == 9007199254740992.0; value.Low == 1.0 + +bool success = DoubleDouble.TryParse("1.25e-2".AsSpan(), CultureInfo.InvariantCulture, + out DoubleDouble parsed); +``` + +- Supported grammar: optional sign, ASCII decimal digits with an optional decimal + separator, and optional `e`/`E` exponent with sign and digits. At least one + mantissa digit is required; `.5` and `1.` are accepted with invariant culture. + Surrounding whitespace is allowed; internal whitespace is not. +- Signs and the decimal separator come from the supplied culture; a null or + omitted provider uses the current culture. Culture-specific NaN and infinity + symbols are recognized case-insensitively. Signed zero is preserved. +- Group separators, currency, parentheses, hexadecimal notation, digit separators, + and `NumberStyles` overloads are not supported. +- Input is limited to **4096 characters**, including surrounding whitespace. + Huge exponents are bounded before constructing powers of ten. Well-formed + overflow succeeds with signed infinity; underflow rounds to a subnormal or + signed zero. A second rounding just below the overflow midpoint stays finite. +- `Parse` throws `ArgumentNullException` for a null string and `FormatException` + for invalid, unsupported, or oversized input. `TryParse` returns `false` and + positive `Zero` for those inputs. + +## Deferred scope + +The planned `PreciseMath` static class and its `Abs`, `Sqrt`, `Pow`, `Exp`, and +`Log` functions are not implemented. Broader generic-math interfaces, expanded +parsing/round-trip formatting, and performance benchmarks remain deferred. +Replacing allocating arithmetic boundary fallbacks is also deferred; the current +`BigInteger` paths remain in place. That optimization does not require removing +`BigInteger` from conversions, parsing, formatting, or independent test oracles. ## Build and test