reciprocal calc function

0-source/Just.PreciseMath/DDMath.cs

@@ -12,6 +12,42 @@ public static partial class DDMath
         return DoubleDouble.Abs(value);
     }
 
+    /// <summary>Returns the reciprocal with double-double precision.</summary>
+    /// <remarks>
+    /// 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 <c>1.0 / value</c>.
+    /// 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.
+    /// </remarks>
+    [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);
+    }
+
     /// <summary>Returns the nonnegative square root with double-double precision.</summary>
     /// <remarks>
     /// Uses power-of-two scaling and an FMA-based Newton correction. Results are

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;
+    }
+}

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,

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