diff --git a/0-source/Just.PreciseMath/DDMath.cs b/0-source/Just.PreciseMath/DDMath.cs index dd936ae..3d6b398 100644 --- a/0-source/Just.PreciseMath/DDMath.cs +++ b/0-source/Just.PreciseMath/DDMath.cs @@ -12,6 +12,42 @@ public static partial class DDMath return DoubleDouble.Abs(value); } + /// Returns the reciprocal with double-double precision. + /// + /// Specializes scalar/DD division for a numerator of one, omitting only its + /// redundant numerator checks. The arithmetic sequence and boundary handling + /// are unchanged, preserving the high and low component bits of 1.0 / value. + /// Results are approximate, not guaranteed correctly rounded; + /// finite results are tested against a relative error bound of 2^-100 plus one + /// minimum binary64 subnormal. Precision decreases near underflow. + /// Signed zeros map to signed infinities, signed infinities to signed zeros, + /// and NaN to canonical NaN. Finite overflow produces signed infinity. + /// + [Pure] + public static DoubleDouble Reciprocal(DoubleDouble value) + { + const double one = 1.0; + if (!DoubleDouble.IsFinite(value) || value._high == 0.0) + { + return new DoubleDouble(one / value._high); + } + if (!PreciseMathHelper.IsDivisionWithinFastRange(value._high)) + { + return PreciseMathHelper.DivideBoundary(one, value); + } + + double quotient = one / value._high; + double remainder = Math.FusedMultiplyAdd(-quotient, value._high, one); + remainder = Math.FusedMultiplyAdd(-quotient, value._low, remainder); + double correction = remainder / value._high; + // For normalized input, using the high denominator in the correction + // adds only O(u^2) error, with u=2^-53. The denominator exponent guard + // keeps the quotient normal and finite; it dominates the correction, + // so QuickTwoSum is ordered and its sum remains finite. + (double high, double low) = PreciseMathHelper.TwoQuickAdd(quotient, correction); + return new DoubleDouble(high, low == 0.0 ? 0.0 : low); + } + /// Returns the nonnegative square root with double-double precision. /// /// Uses power-of-two scaling and an FMA-based Newton correction. Results are diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathReciprocalTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathReciprocalTests.cs new file mode 100644 index 0000000..ab63d4a --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathReciprocalTests.cs @@ -0,0 +1,165 @@ +using System.Numerics; +using Shouldly; +using Xunit; + +namespace Just.PreciseMath.Tests; + +public class PreciseMathReciprocalTests +{ + [Theory] + [InlineData(3.0)] + [InlineData(7.0)] + [InlineData(-3.0)] + [InlineData(-7.0)] + public void OrdinaryReciprocalsRetainMoreThanBinary64Precision(double high) + { + foreach (double low in new[] { 0.0, Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) }) + { + DoubleDouble input = DoubleDouble.FromComponents(high, low); + DoubleDouble actual = DDMath.Reciprocal(input); + actual.Low.ShouldNotBe(0.0); + AssertReciprocalBound(input, actual); + } + } + + [Fact] + public void SpecialValuesFollowBinary64ReciprocalAndRemainCanonical() + { + double[] values = [0.0, -0.0, double.PositiveInfinity, double.NegativeInfinity, + double.NaN, BitConverter.Int64BitsToDouble(0x7ff0000000000001L)]; + foreach (double value in values) + { + DoubleDouble actual = DDMath.Reciprocal(new DoubleDouble(value)); + // Binary64 is an independent oracle for special values only. + double expected = 1.0 / value; + BitConverter.DoubleToInt64Bits(actual.High).ShouldBe( + BitConverter.DoubleToInt64Bits(double.IsNaN(expected) ? double.NaN : expected)); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); + DoubleDouble.IsCanonical(actual).ShouldBeTrue(); + } + } + + [Fact] + public void PowersOfTwoHaveExactReciprocalsOrSignedOverflowAcrossTheRange() + { + // 1/(s*2^k) = s*2^-k. All finite results here are exactly representable, + // including subnormals. Exponents below -1023 overflow the result. + for (int exponent = -1074; exponent <= 1023; ++exponent) + { + foreach (double sign in new[] { -1.0, 1.0 }) + { + DoubleDouble actual = DDMath.Reciprocal(new DoubleDouble(Math.ScaleB(sign, exponent))); + actual.High.ShouldBe(Math.ScaleB(sign, -exponent)); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); + DoubleDouble.IsCanonical(actual).ShouldBeTrue(); + } + } + } + + [Fact] + public void ReciprocalsMatchScalarDivisionBitsAcrossTheRange() + { + // Compatibility is bitwise, not just the same accuracy tolerance. Sample + // every exponent, both signs, dense/sparse lows and binade neighbors; + // this includes both ends of the scalar division fast-path guard. + Random random = new(314159); + for (int exponent = -1074; exponent <= 1023; ++exponent) + { + double high = Math.ScaleB(1.0 + random.NextDouble(), exponent); + double halfUlp = Math.ScaleB(1.0, exponent - 53); + double power = Math.ScaleB(1.0, exponent); + foreach (double sign in new[] { -1.0, 1.0 }) + { + foreach (double low in new[] { 0.0, halfUlp, -halfUlp, + Math.BitDecrement(halfUlp), -Math.BitDecrement(halfUlp), + Math.BitIncrement(halfUlp), -Math.BitIncrement(halfUlp), + double.Epsilon, -double.Epsilon }) + { + AssertMatchesDivision(DoubleDouble.FromComponents(sign * high, low)); + } + foreach (double boundary in new[] { Math.BitDecrement(power), power, Math.BitIncrement(power) }) + { + AssertMatchesDivision(new DoubleDouble(sign * boundary)); + } + } + } + + DoubleDouble[] specials = [DoubleDouble.Zero, new(-0.0), DoubleDouble.NaN, + new(double.PositiveInfinity), new(double.NegativeInfinity), + new(double.MaxValue), new(double.MinValue), + DoubleDouble.FromComponents(double.MaxValue, Math.BitDecrement(Math.ScaleB(1.0, 970))), + DoubleDouble.FromComponents(double.MinValue, -Math.BitDecrement(Math.ScaleB(1.0, 970)))]; + foreach (DoubleDouble input in specials) + { + AssertMatchesDivision(input); + } + } + + [Fact] + public void RepresentableSparseCorrectionsAreRetained() + { + // 1/(1+d) = 1-d+O(d^2). The omitted tail is less than half an ulp + // of d for these exact dyadics, including the minimum subnormal. + foreach (int exponent in new[] { -100, -500, -1000, -1074 }) + { + foreach (double sign in new[] { -1.0, 1.0 }) + { + foreach (double lowSign in new[] { -1.0, 1.0 }) + { + double low = Math.ScaleB(lowSign, exponent); + DoubleDouble input = DoubleDouble.FromComponents(sign, low); + DoubleDouble actual = DDMath.Reciprocal(input); + actual.High.ShouldBe(sign); + actual.Low.ShouldBe(-low); + AssertReciprocalBound(input, actual); + } + } + } + } + + private static void AssertMatchesDivision(DoubleDouble input) + { + DoubleDouble expected = 1.0 / input; + DoubleDouble actual = DDMath.Reciprocal(input); + string context = $"Input: ({input.High:R}, {input.Low:R})"; + BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High), context); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low), context); + DoubleDouble.IsCanonical(actual).ShouldBeTrue(); + if (DoubleDouble.IsFinite(input) && input.High != 0.0 && DoubleDouble.IsFinite(actual)) + { + // Division establishes compatibility, not accuracy: check the exact + // rational inequality independently for every finite matrix result. + AssertReciprocalBound(input, actual); + } + } + + private static void AssertReciprocalBound(DoubleDouble input, DoubleDouble actual) + { + DoubleDouble.IsFinite(actual).ShouldBeTrue(); + DoubleDouble.IsCanonical(actual).ShouldBeTrue(); + Math.Sign(actual.High).ShouldBe(Math.Sign(input.High)); + + // Exact dyadic oracle: x = X*2^-1074, y = Y*2^-1074. + // |y - 1/x| <= |1/x|*2^-100 + 2^-1074 is equivalent to + // |X*Y - 2^2148|*2^100 <= 2^2148 + |X|*2^100. + BigInteger x = Units(input.High) + Units(input.Low); + BigInteger y = Units(actual.High) + Units(actual.Low); + BigInteger scale = BigInteger.One << 2148; + BigInteger error = BigInteger.Abs((x * y) - scale) << 100; + (error <= scale + (BigInteger.Abs(x) << 100)).ShouldBeTrue( + $"Reciprocal bound failed for ({input.High:R}, {input.Low:R}): ({actual.High:R}, {actual.Low:R})"); + } + + 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; + } +} diff --git a/AGENTS.md b/AGENTS.md index a8bb0db..87a3362 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -52,8 +52,9 @@ Follow `.editorconfig`, not incidental style in unfinished code. - Use file-scoped namespaces, explicit types rather than `var`, and block-bodied methods. Preserve the configured expression-bodied property/accessor preferences. - Use `_camelCase` for non-public instance fields, including internal fields; - `s_camelCase` for non-public mutable static fields; PascalCase for constants - and static readonly fields. Do not rename internal fields to remove underscores. + `s_camelCase` for non-public mutable static fields; PascalCase for member constants + and static readonly fields. Method-local constants follow local-variable camelCase. + Do not rename internal fields to remove underscores. - Preserve parentheses that make mathematical grouping readable. `IDE0047` is intentionally disabled; do not re-enable it or remove grouping as style cleanup. - Document public APIs and non-obvious numerical preconditions. Explain algorithms, diff --git a/README.md b/README.md index 1c7cf88..4d8ba4b 100644 --- a/README.md +++ b/README.md @@ -84,6 +84,13 @@ The `DDMath` static class provides: - `Abs(DoubleDouble)`: preserves both components, maps either signed zero to positive zero and either infinity to positive infinity, and returns canonical NaN. It shares the existing `DoubleDouble.Abs` implementation. +- `Reciprocal(DoubleDouble)`: returns exactly the same high and low component bits + as `1.0 / value`. It specializes scalar/DD division for a numerator of one, + omitting only redundant numerator checks while preserving both divisions, + the FMA sequence, and normalization. Signed zeros map to signed infinities, + signed infinities to signed zeros, and NaN to canonical NaN. The allocating exact boundary + path handles extreme exponents; finite overflow produces signed infinity. + No speedup over scalar/DD division has been measured. - `Sqrt(DoubleDouble)`: uses power-of-two scaling and an FMA-based Newton correction to retain extended precision, including for subnormal inputs, without squaring an unscaled estimate near the exponent limits. Signed zero and positive infinity @@ -122,6 +129,7 @@ using Just.PreciseMath; DoubleDouble root = DDMath.Sqrt(new DoubleDouble(2.0)); DoubleDouble inverseRoot = DDMath.InvSqrt(new DoubleDouble(2.0)); +DoubleDouble reciprocal = DDMath.Reciprocal(new DoubleDouble(3.0)); DoubleDouble magnitude = DDMath.Abs(-root); DoubleDouble smallPower = DDMath.Pow(new DoubleDouble(2.0), -1024); DoubleDouble exponential = DDMath.Exp(new DoubleDouble(1.0)); @@ -131,6 +139,13 @@ DoubleDouble preciseExponent = DoubleDouble.FromComponents(0.5, 1e-30); DoubleDouble precisePower = DDMath.Pow(new DoubleDouble(2.0), preciseExponent); ``` +Reciprocal tests check bitwise equivalence with `1.0 / value` and independently +check exact rational error against `2^-100` relative plus one minimum binary64 +subnormal. They sample every binary64 exponent, both signs, dense and sparse lows, +binade neighbors, and special values; selected powers of two and sparse corrections +also have exact component checks. This preserves division's approximate-accuracy +contract, not a guarantee of correctly rounded results. + Square-root and inverse-square-root tests compare the exact component sum against a `2^-100` relative error bound using integer inequalities. They include samples at every binary64 exponent, boundary neighbors, both signs of the low component, and exact binary