diff --git a/0-source/Just.PreciseMath/DDMath.cs b/0-source/Just.PreciseMath/DDMath.cs
index 303ec2a..dd936ae 100644
--- a/0-source/Just.PreciseMath/DDMath.cs
+++ b/0-source/Just.PreciseMath/DDMath.cs
@@ -55,6 +55,59 @@ public static partial class DDMath
return DoubleDouble.FromComponents(Math.ScaleB(root.High, rootExponent), Math.ScaleB(root.Low, rootExponent));
}
+ /// Returns the reciprocal square root with double-double precision.
+ ///
+ /// Uses power-of-two scaling and a compensated Newton step without double-double
+ /// division. Results are approximate, not guaranteed correctly rounded; positive
+ /// finite inputs are tested against a relative error bound of 2^-100 across the
+ /// binary64 range. Positive/negative zero maps to positive/negative infinity,
+ /// positive infinity to positive zero, and negative nonzero values or NaN to
+ /// canonical NaN.
+ ///
+ [Pure]
+ public static DoubleDouble InvSqrt(DoubleDouble value)
+ {
+ if (double.IsNaN(value.High) || value.High < 0.0)
+ {
+ return DoubleDouble.NaN;
+ }
+ if (value.High == 0.0)
+ {
+ return new DoubleDouble(Math.CopySign(double.PositiveInfinity, value.High));
+ }
+ if (double.IsPositiveInfinity(value.High))
+ {
+ return DoubleDouble.Zero;
+ }
+
+ // As in Sqrt, the even exponent is in [-1074, 1022], high is in [1, 4),
+ // and any underflow in a sparse scaled low is below the error bound.
+ int exponent = Math.ILogB(value.High) & ~1;
+ double high = Math.ScaleB(value.High, -exponent);
+ double low = Math.ScaleB(value.Low, -exponent);
+ double estimate = 1.0 / Math.Sqrt(high);
+
+ // With u=2^-53, the seed has O(u) relative error. Preserve the square's
+ // FMA residual before cancellation in 1 - (high + low)*estimate^2:
+ // using only the rounded square would leave O(u) error after refinement.
+ // All leading products are bounded and normal. Omitting low*squareError
+ // contributes only O(u^2), as does one Newton step's remaining error.
+ double square = estimate * estimate;
+ double squareError = Math.FusedMultiplyAdd(estimate, estimate, -square);
+ double residual = Math.FusedMultiplyAdd(-high, square, 1.0);
+ residual = Math.FusedMultiplyAdd(-high, squareError, residual);
+ residual = Math.FusedMultiplyAdd(-low, square, residual);
+ double correction = (0.5 * estimate) * residual;
+ DoubleDouble inverseRoot = DoubleDouble.FromComponents(estimate, correction);
+
+ // The result high is normal and finite over the entire positive input
+ // range. Rescale in the opposite direction to Sqrt and canonicalize any
+ // zero low, including underflow of an exceptionally sparse correction.
+ int inverseExponent = -(exponent / 2);
+ return DoubleDouble.FromComponents(Math.ScaleB(inverseRoot.High, inverseExponent),
+ Math.ScaleB(inverseRoot.Low, inverseExponent));
+ }
+
// For finite, nonzero normalized significands and exponents bounded by the
// integer-power domain (|exponent| < 2^42). Keep the exponent separate until
// the final result so an intermediate cannot overflow before reciprocation.
diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathInvSqrtTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathInvSqrtTests.cs
new file mode 100644
index 0000000..c029044
--- /dev/null
+++ b/1-tests/Just.PreciseMath.Tests/PreciseMathInvSqrtTests.cs
@@ -0,0 +1,186 @@
+using System.Numerics;
+using Shouldly;
+using Xunit;
+
+namespace Just.PreciseMath.Tests;
+
+public class PreciseMathInvSqrtTests
+{
+ [Theory]
+ [InlineData(2.0)]
+ [InlineData(3.0)]
+ [InlineData(5.0)]
+ [InlineData(1e-308)]
+ [InlineData(1e308)]
+ public void IrrationalInverseRootsRetainMoreThanBinary64Precision(double value)
+ {
+ DoubleDouble input = new(value);
+ DoubleDouble actual = DDMath.InvSqrt(input);
+ actual.Low.ShouldNotBe(0.0);
+ AssertInvSqrtBound(input, actual);
+ }
+
+ [Fact]
+ public void SpecialValuesMatchReciprocalSquareRootAndRemainCanonical()
+ {
+ double[] values = [0.0, -0.0, double.PositiveInfinity, double.NegativeInfinity,
+ double.NaN, -1.0, -double.Epsilon, double.MinValue];
+ foreach (double value in values)
+ {
+ DoubleDouble actual = DDMath.InvSqrt(new DoubleDouble(value));
+ // Binary64 is an independent oracle only for these special values.
+ double expected = 1.0 / Math.Sqrt(value);
+ BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(
+ BitConverter.DoubleToInt64Bits(double.IsNaN(expected) ? double.NaN : expected));
+ BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
+ DoubleDouble.IsCanonical(actual).ShouldBeTrue();
+ }
+
+ DoubleDouble negative = DoubleDouble.FromComponents(-1.0, Math.ScaleB(1.0, -54));
+ DoubleDouble.IsNaN(DDMath.InvSqrt(negative)).ShouldBeTrue();
+ }
+
+ [Fact]
+ public void PowersOfFourHaveExactInverseRootsAcrossTheFiniteRange()
+ {
+ // 1/sqrt(2^(2k)) = 2^-k exactly, including the minimum subnormal input.
+ for (int exponent = -1074; exponent <= 1022; exponent += 2)
+ {
+ DoubleDouble actual = DDMath.InvSqrt(new DoubleDouble(Math.ScaleB(1.0, exponent)));
+ actual.High.ShouldBe(Math.ScaleB(1.0, -(exponent / 2)));
+ BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
+ }
+ }
+
+ [Fact]
+ public void LowComponentsOnEitherSideOfOneAffectTheInverseRoot()
+ {
+ foreach (double low in new[] { Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) })
+ {
+ DoubleDouble input = DoubleDouble.FromComponents(1.0, low);
+ DoubleDouble actual = DDMath.InvSqrt(input);
+ actual.High.ShouldBe(1.0);
+ Math.Sign(actual.Low).ShouldBe(-Math.Sign(low));
+ AssertInvSqrtBound(input, actual);
+ }
+ }
+
+ [Fact]
+ public void PositiveFiniteInputsMeetExactRelativeErrorBound()
+ {
+ // Every input exponent, both parities, dense/sparse lows of either sign,
+ // and neighbors of binade transitions. The oracle uses the exact stored
+ // input sum, including any rounding during public input construction.
+ Random random = new(271828);
+ for (int exponent = -1074; exponent <= 1023; ++exponent)
+ {
+ double high = Math.ScaleB(1.0 + random.NextDouble(), exponent);
+ double low = Math.ScaleB(random.NextDouble(), exponent - 53);
+ double power = Math.ScaleB(1.0, exponent);
+ foreach (double boundary in new[] { Math.BitDecrement(power), power, Math.BitIncrement(power) })
+ {
+ if (boundary > 0.0)
+ {
+ DoubleDouble input = new(boundary);
+ AssertInvSqrtBound(input, DDMath.InvSqrt(input));
+ }
+ }
+ foreach (double residual in new[] { 0.0, low, -low, double.Epsilon, -double.Epsilon })
+ {
+ DoubleDouble input = DoubleDouble.FromComponents(high, residual);
+ if (input.High > 0.0)
+ {
+ AssertInvSqrtBound(input, DDMath.InvSqrt(input));
+ }
+ }
+ }
+
+ DoubleDouble[] boundaries = [new(1e-308), new(double.Epsilon),
+ new(Math.BitDecrement(Math.ScaleB(1.0, -1022))), new(Math.ScaleB(1.0, -1022)),
+ new(Math.BitIncrement(Math.ScaleB(1.0, -1022))), new(double.MaxValue),
+ DoubleDouble.FromComponents(double.MaxValue, Math.BitDecrement(Math.ScaleB(1.0, 970))),
+ DoubleDouble.FromComponents(double.MaxValue, -Math.ScaleB(1.0, 970)),
+ DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -53)),
+ DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -54)),
+ DoubleDouble.FromComponents(4.0, Math.ScaleB(1.0, -51)),
+ DoubleDouble.FromComponents(4.0, -Math.ScaleB(1.0, -52)),
+ DoubleDouble.FromComponents(0.0, double.Epsilon)];
+ foreach (DoubleDouble input in boundaries)
+ {
+ AssertInvSqrtBound(input, DDMath.InvSqrt(input));
+ }
+ }
+
+ [Fact]
+ public void DenseLowNormalizationBoundariesMeetTheBound()
+ {
+ // Exercise adjacent lows at half an ulp, both signs, for each exponent.
+ Random random = new(161803);
+ for (int exponent = -1074; exponent <= 1023; ++exponent)
+ {
+ double high = Math.ScaleB(1.0 + random.NextDouble(), exponent);
+ double halfUlp = Math.ScaleB(1.0, exponent - 53);
+ foreach (double magnitude in new[] { Math.BitDecrement(halfUlp), halfUlp, Math.BitIncrement(halfUlp) })
+ {
+ foreach (double low in new[] { magnitude, -magnitude })
+ {
+ DoubleDouble input = DoubleDouble.FromComponents(high, low);
+ if (DoubleDouble.IsFinite(input) && input.High > 0.0)
+ {
+ AssertInvSqrtBound(input, DDMath.InvSqrt(input));
+ }
+ }
+ }
+ }
+ }
+
+ [Fact]
+ public void RepresentableSparseCorrectionsAreNotDiscarded()
+ {
+ // (1+d)^(-1/2) = 1-d/2+O(d^2). For these dyadic d, the quadratic
+ // term is below half an ulp of d/2, even at the minimum subnormal.
+ foreach (int exponent in new[] { -100, -500, -1000, -1073 })
+ {
+ foreach (double sign in new[] { -1.0, 1.0 })
+ {
+ double low = Math.ScaleB(sign, exponent);
+ DoubleDouble actual = DDMath.InvSqrt(DoubleDouble.FromComponents(1.0, low));
+ actual.High.ShouldBe(1.0);
+ actual.Low.ShouldBe(Math.ScaleB(-sign, exponent - 1));
+ DoubleDouble.IsCanonical(actual).ShouldBeTrue();
+ }
+ }
+ }
+
+ private static void AssertInvSqrtBound(DoubleDouble input, DoubleDouble actual)
+ {
+ DoubleDouble.IsFinite(actual).ShouldBeTrue();
+ DoubleDouble.IsCanonical(actual).ShouldBeTrue();
+ (actual.High > 0.0).ShouldBeTrue();
+
+ // Exact dyadic oracle: x = X*2^-1074, y = Y*2^-1074. For positive
+ // x and y, |y/(1/sqrt(x)) - 1| <= t iff (1-t)^2 <= x*y^2 <= (1+t)^2.
+ // Cross-multiply with t=2^-100; no DD product, division or rounded root.
+ BigInteger x = Units(input.High) + Units(input.Low);
+ BigInteger y = Units(actual.High) + Units(actual.Low);
+ BigInteger scale = BigInteger.One << 100;
+ BigInteger product = (x * y * y) << 200;
+ BigInteger lower = ((scale - 1) * (scale - 1)) << 3222;
+ BigInteger upper = ((scale + 1) * (scale + 1)) << 3222;
+ (product >= lower && product <= upper).ShouldBeTrue(
+ $"InvSqrt 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/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs
index a9cf5c4..f48c1f2 100644
--- a/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs
+++ b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs
@@ -116,6 +116,8 @@ public class PreciseMathSqrtTests
new(Math.BitIncrement(Math.ScaleB(1.0, -1022))), new(double.MaxValue),
DoubleDouble.FromComponents(double.MaxValue, Math.ScaleB(1.0, 969)),
DoubleDouble.FromComponents(double.MaxValue, -Math.ScaleB(1.0, 969)),
+ DoubleDouble.FromComponents(double.MaxValue, Math.BitDecrement(Math.ScaleB(1.0, 970))),
+ DoubleDouble.FromComponents(double.MaxValue, -Math.ScaleB(1.0, 970)),
DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -53)),
DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -54)),
DoubleDouble.FromComponents(4.0, Math.ScaleB(1.0, -51)),
@@ -126,6 +128,84 @@ public class PreciseMathSqrtTests
}
}
+ [Fact]
+ public void RootRoundingMidpointsAndAdjacentLowsMeetTheBound()
+ {
+ // Exact squares of root midpoints 1 + 2^-53 and 2 - 2^-53:
+ // (1 + 2^-52) + 2^-106 and (4 - 2^-51) + 2^-106.
+ // Bracket each with adjacent low doubles, then exercise exponent parity
+ // and rescaling. At tiny exponents the input itself loses low precision;
+ // the oracle always checks the exact stored input, not the unscaled square.
+ double midpointLow = Math.ScaleB(1.0, -106);
+ for (int exponent = -1074; exponent <= 1022; ++exponent)
+ {
+ foreach (double high in new[] { Math.BitIncrement(1.0), Math.BitDecrement(4.0) })
+ {
+ foreach (double low in new[] { Math.BitDecrement(midpointLow), midpointLow, Math.BitIncrement(midpointLow) })
+ {
+ DoubleDouble input = DoubleDouble.FromComponents(Math.ScaleB(high, exponent), Math.ScaleB(low, exponent));
+ if (DoubleDouble.IsFinite(input) && input.High > 0.0)
+ {
+ AssertSqrtBound(input, DDMath.Sqrt(input));
+ }
+ }
+ }
+ }
+ }
+
+ [Fact]
+ public void PublicFactoryNormalizesLegacyZeroHighInputsBeforeTakingTheRoot()
+ {
+ // Legacy raw (0, low) examples are normalized at the public boundary.
+ foreach (double low in new[] { double.Epsilon, 1e-308, 2.0, double.MaxValue })
+ {
+ DoubleDouble input = DoubleDouble.FromComponents(0.0, low);
+ AssertSqrtBound(input, DDMath.Sqrt(input));
+ }
+ DoubleDouble.IsNaN(DDMath.Sqrt(DoubleDouble.FromComponents(0.0, -1.0))).ShouldBeTrue();
+ }
+
+ [Fact]
+ public void DenseLowNormalizationBoundariesMeetTheBound()
+ {
+ // Half-ulp lows and their neighbors test both sides of input normalization.
+ Random random = new(161803);
+ for (int exponent = -1074; exponent <= 1023; ++exponent)
+ {
+ double high = Math.ScaleB(1.0 + random.NextDouble(), exponent);
+ double halfUlp = Math.ScaleB(1.0, exponent - 53);
+ foreach (double magnitude in new[] { Math.BitDecrement(halfUlp), halfUlp, Math.BitIncrement(halfUlp) })
+ {
+ foreach (double low in new[] { magnitude, -magnitude })
+ {
+ DoubleDouble input = DoubleDouble.FromComponents(high, low);
+ if (DoubleDouble.IsFinite(input) && input.High > 0.0)
+ {
+ AssertSqrtBound(input, DDMath.Sqrt(input));
+ }
+ }
+ }
+ }
+ }
+
+ [Fact]
+ public void RepresentableSparseCorrectionsAreNotDiscarded()
+ {
+ // sqrt(1+d) = 1+d/2+O(d^2). For these dyadic d, the quadratic term
+ // is below half an ulp of d/2, including when d/2 is the minimum subnormal.
+ foreach (int exponent in new[] { -100, -500, -1000, -1073 })
+ {
+ foreach (double sign in new[] { -1.0, 1.0 })
+ {
+ double low = Math.ScaleB(sign, exponent);
+ DoubleDouble actual = DDMath.Sqrt(DoubleDouble.FromComponents(1.0, low));
+ actual.High.ShouldBe(1.0);
+ actual.Low.ShouldBe(Math.ScaleB(sign, exponent - 1));
+ DoubleDouble.IsCanonical(actual).ShouldBeTrue();
+ }
+ }
+ }
+
private static void AssertSqrtBound(DoubleDouble input, DoubleDouble actual)
{
DoubleDouble.IsFinite(actual).ShouldBeTrue();
diff --git a/README.md b/README.md
index 108b563..1c7cf88 100644
--- a/README.md
+++ b/README.md
@@ -88,6 +88,12 @@ The `DDMath` static class provides:
to retain extended precision, including for subnormal inputs, without squaring
an unscaled estimate near the exponent limits. Signed zero and positive infinity
are preserved; negative nonzero inputs and NaN return canonical NaN.
+- `InvSqrt(DoubleDouble)`: computes the reciprocal square root with power-of-two
+ scaling and a compensated Newton step, avoiding double-double division and its
+ allocating boundary paths. The estimate's squared-product residual is retained
+ with FMA. Signed zeros map to correspondingly signed infinities; positive infinity
+ maps to positive zero; negative nonzero inputs and NaN return canonical NaN.
+ This is a dedicated algorithm, not a claim of measured speedup over `1.0 / Sqrt(x)`.
- `Pow(DoubleDouble, int)`: exponentiation by squaring with a separately tracked
binary exponent. Supports the full `int` domain, including `int.MinValue`, and
reciprocates a bounded significand before final scaling for negative exponents.
@@ -115,6 +121,7 @@ The `DDMath` static class provides:
using Just.PreciseMath;
DoubleDouble root = DDMath.Sqrt(new DoubleDouble(2.0));
+DoubleDouble inverseRoot = DDMath.InvSqrt(new DoubleDouble(2.0));
DoubleDouble magnitude = DDMath.Abs(-root);
DoubleDouble smallPower = DDMath.Pow(new DoubleDouble(2.0), -1024);
DoubleDouble exponential = DDMath.Exp(new DoubleDouble(1.0));
@@ -124,10 +131,12 @@ DoubleDouble preciseExponent = DoubleDouble.FromComponents(0.5, 1e-30);
DoubleDouble precisePower = DDMath.Pow(new DoubleDouble(2.0), preciseExponent);
```
-Square-root tests compare the exact component sum against a `2^-100` relative
-error bound using integer inequalities. They include samples at every binary64
+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
-squares. This is a tested approximate-accuracy contract, not exhaustive coverage
+squares/powers of four. Square-root tests also bracket exact root-rounding midpoints;
+both suites exercise half-ulp low-component normalization boundaries.
+This is a tested approximate-accuracy contract, not exhaustive coverage
of all component pairs or a guarantee of correctly rounded results.
Logarithm tests compare exact component sums with independently generated