diff --git a/0-source/Just.PreciseMath/DDMath.cs b/0-source/Just.PreciseMath/DDMath.cs index 810c69d..ce2a445 100644 --- a/0-source/Just.PreciseMath/DDMath.cs +++ b/0-source/Just.PreciseMath/DDMath.cs @@ -12,31 +12,18 @@ 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 Square(DoubleDouble value) + { + return DoubleDouble.Square(value); + } + + /// [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); - } - - return PreciseMathHelper.DivideScalarFinite(one, value); + return DoubleDouble.Reciprocal(value); } /// diff --git a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs index 37085ed..e66ae4e 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs @@ -6,6 +6,43 @@ public readonly partial struct DoubleDouble : IMultiplyOperators, IDivisionOperators { + /// Returns the square using a specialized double-double product. + /// + /// Combines the equal cross terms before FMA accumulation, retaining low squared. + /// Results are approximate and may differ in low bits from value * value; + /// finite results are tested against 2^-100 relative error plus one minimum + /// binary64 subnormal, not a guarantee of correctly rounded components. + /// Either zero maps to positive zero, either infinity to positive infinity, + /// and NaN to canonical NaN. Uses the allocating multiplication fallback at + /// exponent boundaries. No performance improvement is guaranteed. + /// + [Pure] + public static DoubleDouble Square(DoubleDouble value) + { + if (!IsFinite(value) || value._high == 0.0) + { + return new DoubleDouble(value._high * value._high); + } + if (!PreciseMathHelper.IsMultiplicationWithinFastRange(value._high, value._high)) + { + return PreciseMathHelper.MultiplyBoundary(value, value); + } + + (double product, double error) = PreciseMathHelper.TwoSquare(value._high); + // The guard bounds the high exponent to [-450, 450], so doubling high + // is exact and finite. FMA avoids rounding the cross product separately. + error = Math.FusedMultiplyAdd(value._high + value._high, value._low, error); + // As for multiplication, normalized input makes |error| < 4u*product; + // the positive high product is normal and the corrected sum is finite. + (double high, double low) = PreciseMathHelper.TwoQuickAdd(product, error); + // Incorporate low squared AFTER the cross term has been absorbed into + // high. Otherwise an exactly representable tail can round away before + // cancellation, e.g. (1 - 2^-54)^2 loses its 2^-108 residual. + low = Math.FusedMultiplyAdd(value._low, value._low, low); + (high, low) = PreciseMathHelper.TwoQuickAdd(high, low); + return new DoubleDouble(high, low == 0.0 ? 0.0 : low); + } + /// Returns the operand unchanged. public static DoubleDouble operator +(DoubleDouble value) { diff --git a/0-source/Just.PreciseMath/DoubleDouble.Formatting.cs b/0-source/Just.PreciseMath/DoubleDouble.Formatting.cs index 5ca76b0..ba5145f 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.Formatting.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Formatting.cs @@ -22,6 +22,37 @@ public readonly partial struct DoubleDouble : IFormattable return ToString(null, provider); } + /// Returns the exact decimal expansion of the stored value using invariant culture. + /// + /// Uses allocating integer arithmetic, without rounding or scientific notation. + /// Omits unnecessary fractional zeros and preserves negative zero; nonfinite + /// values use NaN, Infinity, and -Infinity. Sparse pairs can require over a + /// thousand characters. Parse with invariant culture to recover the components; + /// this is an exact-value format, not a shortest-round-trip format. + /// + public string ToStringExact() + { + if (!double.IsFinite(_high) || _high == 0.0) + { + return _high.ToString(CultureInfo.InvariantCulture); + } + + (BigInteger numerator, BigInteger denominator) = ConversionFraction(); + // A sum of binary64 values has a power-of-two denominator. Multiplying + // by the matching power of five makes an integer coefficient over 10^n. + int decimalPlaces = (int)(denominator.GetBitLength() - 1); + BigInteger coefficient = BigInteger.Abs(numerator) * BigInteger.Pow(5, decimalPlaces); + string digits = coefficient.ToString(CultureInfo.InvariantCulture); + int end = digits.Length; + while (decimalPlaces > 0 && digits[end - 1] == '0') + { + end--; + decimalPlaces--; + } + string sign = numerator.Sign < 0 ? "-" : string.Empty; + return sign + FormattingFixed(digits[..end], decimalPlaces, "."); + } + /// /// 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 diff --git a/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs b/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs index 0c0894d..58feb7b 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs @@ -5,6 +5,33 @@ public readonly partial struct DoubleDouble /// Gets the binary radix of the components. public static int Radix => 2; + /// 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); + } + + return PreciseMathHelper.DivideScalarFinite(one, value); + } + /// Returns the absolute value, preserving both components and canonicalizing NaN. public static DoubleDouble Abs(DoubleDouble value) { diff --git a/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs b/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs index cb07fa4..067436f 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs @@ -244,7 +244,8 @@ public readonly partial struct DoubleDouble : IRootFunctions /// Returns the nonnegative square root with double-double precision. /// - /// Uses power-of-two scaling and an FMA-based Newton correction. Results are + /// Uses power-of-two scaling and an FMA-based Newton correction. An exponent- + /// tracked correction retains sparse lows endangered by downscaling. Results are /// approximate, not guaranteed correctly rounded; finite positive inputs are /// tested against a relative error bound of 2^-100 across the binary64 range. /// Signed zero and positive infinity are preserved. Negative nonzero values @@ -264,11 +265,19 @@ public readonly partial struct DoubleDouble : IRootFunctions // Choose an even exponent, rounding negative odd exponents down too. // For positive finite normalized inputs, exponent is in [-1074, 1022] - // and the scaled high is in [1, 4). Scaling the high is exact; any - // underflow in a very sparse low is far below the relative error bound. + // and the scaled high is in [1, 4). Scaling the high is exact. int exponent = Math.ILogB(value.High) & ~1; double high = Math.ScaleB(value.High, -exponent); double low = Math.ScaleB(value.Low, -exponent); + // Dividing the scaled low by 2*estimate (in [2, 4]) can round at the + // subnormal floor before output rescaling restores its magnitude. In + // that sparse domain, refine high alone and apply the original low with + // separate exponents afterward. Never add both the scaled and original low. + bool separateLow = exponent > 0 && value.Low != 0.0 && Math.Abs(low) < Math.ScaleB(1.0, -1020); + if (separateLow) + { + low = 0.0; + } double estimate = Math.Sqrt(high); // Estimate is in [1, 2]. FMA avoids rounding its square before @@ -282,7 +291,10 @@ public readonly partial struct DoubleDouble : IRootFunctions // Normalize again to canonicalize a zero low, including any underflow // of an exceptionally sparse correction during rescaling. int rootExponent = exponent / 2; - return DoubleDouble.FromComponents(Math.ScaleB(root.High, rootExponent), Math.ScaleB(root.Low, rootExponent)); + root = DoubleDouble.FromComponents(Math.ScaleB(root.High, rootExponent), Math.ScaleB(root.Low, rootExponent)); + // In this branch |originalLow/originalHigh| < 2^-1020. The omitted + // quadratic term is negligible even relative to the low correction. + return separateLow ? root + RootLowCorrection(value, root, 2) : root; } /// Returns the reciprocal square root with double-double precision. diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs index ce09c43..4208975 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs @@ -978,6 +978,114 @@ public class DoubleDoubleArithmeticTests } } + [Theory] + [InlineData(-450)] + [InlineData(0)] + [InlineData(450)] + public void SquareRetainsExactCrossAndLowSquaredTerms(int exponent) + { + // (h - h*2^-54)^2 = h^2 - h^2*2^-53 + h^2*2^-108. + // (h + h*2^-53)^2 = h^2 + h^2*2^-52 + h^2*2^-106. + // Both sums fit exactly in two components, including fast-range endpoints. + double high = Math.ScaleB(1.0, exponent); + double highSquared = Math.ScaleB(1.0, 2 * exponent); + DoubleDouble below = DoubleDouble.FromComponents(high, -Math.ScaleB(1.0, exponent - 54)); + DoubleDouble above = DoubleDouble.FromComponents(high, Math.ScaleB(1.0, exponent - 53)); + foreach (double sign in new[] { -1.0, 1.0 }) + { + CheckBoundary(DoubleDouble.Square(sign * below), Math.BitDecrement(highSquared), Math.ScaleB(1.0, (2 * exponent) - 108)); + CheckBoundary(DDMath.Square(sign * above), Math.BitIncrement(highSquared), Math.ScaleB(1.0, (2 * exponent) - 106)); + } + } + + [Fact] + public void SquareHandlesSpecialValuesAndExactBinaryPowers() + { + foreach (double value in new[] { 0.0, -0.0, double.NaN, double.PositiveInfinity, double.NegativeInfinity }) + { + CheckBits(DoubleDouble.Square(new DoubleDouble(value)), value * value); + CheckBits(DDMath.Square(new DoubleDouble(value)), value * value); + } + for (int exponent = -1074; exponent <= 1023; exponent++) + { + // A binary power squared needs only exponent arithmetic, not a DD oracle. + double high = Math.ScaleB(1.0, exponent); + double expected = Math.ScaleB(1.0, 2 * exponent); + CheckBits(DoubleDouble.Square(new DoubleDouble(high)), expected); + CheckBits(DoubleDouble.Square(new DoubleDouble(-high)), expected); + } + } + + [Fact] + public void SquareMeetsExactRationalBoundsAcrossRangesAndDispatchTransitions() + { + Random random = new(57721); + for (int exponent = -1074; exponent <= 1023; exponent++) + { + double high = Math.ScaleB(1.0 + random.NextDouble(), exponent); + double low = Math.ScaleB(random.NextDouble(), exponent - 53); + foreach (double residual in new[] { 0.0, low, -low, double.Epsilon, -double.Epsilon }) + { + AssertSquare(DoubleDouble.FromComponents(high, residual)); + } + double power = Math.ScaleB(1.0, exponent); + AssertSquare(new DoubleDouble(Math.BitDecrement(power))); + AssertSquare(new DoubleDouble(power)); + AssertSquare(new DoubleDouble(Math.BitIncrement(power))); + } + } + + [Fact] + public void SquareUsesTheCompleteInputForOverflowAndUnderflowClassification() + { + // Neighborhoods of sqrt(overflow midpoint) and sqrt(epsilon/2). + // The exact integer oracle, not Math.Sqrt or another DD operation, + // decides which side of each threshold the actual stored input lies on. + double overflowHigh = Math.ScaleB(1.0, 512); + // At low=-2^457, the square exceeds the midpoint by low^2; + // the next more-negative low puts it strictly below the midpoint. + double overflowLow = -Math.ScaleB(1.0, 457); + double underflowHigh = Math.ScaleB(Math.Sqrt(0.5), -537); + double underflowLow = Math.ScaleB(1.0, -591); + foreach (double high in new[] { Math.BitDecrement(overflowHigh), overflowHigh, Math.BitIncrement(overflowHigh) }) + { + foreach (double low in new[] { Math.BitDecrement(overflowLow), overflowLow, Math.BitIncrement(overflowLow), 0.0 }) + { + AssertSquare(DoubleDouble.FromComponents(high, low)); + } + } + foreach (double high in new[] { Math.BitDecrement(underflowHigh), underflowHigh, Math.BitIncrement(underflowHigh) }) + { + foreach (double low in new[] { -underflowLow, 0.0, underflowLow }) + { + AssertSquare(DoubleDouble.FromComponents(high, low)); + } + } + } + + private static void AssertSquare(DoubleDouble value) + { + BigInteger input = Units(value); + BigInteger numerator = input * input; + BigInteger denominator = BigInteger.One << 1074; + BigInteger overflow = Units(double.MaxValue) + Units(Math.ScaleB(1.0, 970)); + DoubleDouble actual = DoubleDouble.Square(value); + DoubleDouble.IsCanonical(actual).ShouldBeTrue(); + DoubleDouble.IsNegative(actual).ShouldBeFalse(); + if (numerator >= overflow * denominator) + { + CheckBits(actual, double.PositiveInfinity); + } + else + { + AssertRelative(actual, numerator, denominator); + // The absolute error floor must not hide a wrong zero/nonzero result. + DoubleDouble.IsZero(actual).ShouldBe((numerator << 1) <= denominator); + } + CheckBoundary(DoubleDouble.Square(-value), actual.High, actual.Low); + CheckBoundary(DDMath.Square(value), actual.High, actual.Low); + } + private static (double[] Scalars, DoubleDouble[] Values) ScalarArithmeticCompatibilityCases() { List scalars = [0.0, -0.0, 1.0, -1.0, 2.0, -2.0, 3.0, -3.0, 100.0, -100.0, diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleFormattingTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleFormattingTests.cs index fab29ed..f7a0860 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleFormattingTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleFormattingTests.cs @@ -4,6 +4,79 @@ namespace Just.PreciseMath.Tests; public class DoubleDoubleFormattingTests { + [Theory] + [InlineData(1.25, 0.0, "1.25")] + [InlineData(-1.25, 0.0, "-1.25")] + [InlineData(1000.0, 0.0, "1000")] + [InlineData(1.0, 5.551115123125783e-17, "1.000000000000000055511151231257827021181583404541015625")] + [InlineData(9007199254740992.0, 1.0, "9007199254740993")] + public void ExactFormattingPreservesTheCompleteDyadicValue(double high, double low, string expected) + { + // Exact binary fractions: the nonzero fractional low is 2^-54. + DoubleDouble value = DoubleDouble.FromComponents(high, low); + value.ToStringExact().ShouldBe(expected); + } + + [Fact] + public void ExactFormattingUsesInvariantSpecialSymbolsAndPreservesZeroSigns() + { + CultureInfo original = CultureInfo.CurrentCulture; + try + { + CultureInfo.CurrentCulture = CultureInfo.GetCultureInfo("fr-FR"); + new DoubleDouble(-1.25).ToStringExact().ShouldBe("-1.25"); + DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, + DoubleDouble.NaN, DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity]; + string[] expected = ["0", "-0", "NaN", "Infinity", "-Infinity"]; + for (int i = 0; i < values.Length; i++) + { + string text = values[i].ToStringExact(); + text.ShouldBe(expected[i]); + DoubleDouble parsed = DoubleDouble.Parse(text, CultureInfo.InvariantCulture); + BitConverter.DoubleToInt64Bits(parsed.High).ShouldBe(BitConverter.DoubleToInt64Bits(values[i].High)); + BitConverter.DoubleToInt64Bits(parsed.Low).ShouldBe(0L); + } + } + finally + { + CultureInfo.CurrentCulture = original; + } + } + + [Fact] + public void ExactFormattingPreservesSparsePairsBeyondTheStandardPrecisionLimit() + { + DoubleDouble[] values = [new(double.Epsilon), new(double.MaxValue), DoubleDouble.Pi, + DoubleDouble.FromComponents(1.0, double.Epsilon), + DoubleDouble.FromComponents(1.0, -double.Epsilon), + DoubleDouble.FromComponents(double.MaxValue, double.Epsilon), + DoubleDouble.FromComponents(double.MaxValue, -double.Epsilon), + DoubleDouble.FromComponents(double.MaxValue, Math.BitDecrement(Math.ScaleB(1.0, 970)))]; + foreach (DoubleDouble value in values) + { + AssertExactText(value); + AssertExactText(-value); + } + DoubleDouble.FromComponents(1.0, double.Epsilon).ToStringExact().Length.ShouldBeGreaterThan(999); + } + + [Fact] + public void ExactFormattingMatchesIndependentIntegerValuesAcrossEveryExponent() + { + Random random = new(271828); + for (int exponent = -1074; exponent <= 1023; exponent++) + { + double high = Math.ScaleB(1.0 + random.NextDouble(), exponent); + double denseLow = Math.ScaleB(random.NextDouble(), exponent - 53); + foreach (double low in new[] { 0.0, denseLow, -denseLow, double.Epsilon, -double.Epsilon }) + { + DoubleDouble value = DoubleDouble.FromComponents(high, low); + AssertExactText(value); + AssertExactText(-value); + } + } + } + [Theory] [InlineData("G0")] [InlineData("g0")] @@ -138,4 +211,40 @@ public class DoubleDoubleFormattingTests { Should.Throw(() => DoubleDouble.One.ToString(format, CultureInfo.InvariantCulture)); } + + private static void AssertExactText(DoubleDouble value) + { + string text = value.ToStringExact(); + text.Length.ShouldBeLessThanOrEqualTo(2048); + text.ShouldNotContain("E"); + text.ShouldNotContain("e"); + int point = text.IndexOf('.', StringComparison.Ordinal); + int places = point < 0 ? 0 : text.Length - point - 1; + if (point >= 0) + { + text.ShouldNotEndWith("0"); + text.ShouldNotEndWith("."); + } + string digits = point < 0 ? text : text.Remove(point, 1); + BigInteger coefficient = BigInteger.Parse(digits, CultureInfo.InvariantCulture); + // Independently decode the input as integer multiples of 2^-1074 and + // compare with the printed integer coefficient / 10^places exactly. + BigInteger units = ExactFormattingUnits(value.High) + ExactFormattingUnits(value.Low); + (coefficient << 1074).ShouldBe(units * BigInteger.Pow(10, places)); + DoubleDouble parsed = DoubleDouble.Parse(text, CultureInfo.InvariantCulture); + BitConverter.DoubleToInt64Bits(parsed.High).ShouldBe(BitConverter.DoubleToInt64Bits(value.High)); + BitConverter.DoubleToInt64Bits(parsed.Low).ShouldBe(BitConverter.DoubleToInt64Bits(value.Low)); + } + + private static BigInteger ExactFormattingUnits(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; + } } diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs index 94fbcb8..327bc39 100644 --- a/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs @@ -202,6 +202,83 @@ public class PreciseMathSqrtTests } } + [Theory] + [InlineData(-1.0)] + [InlineData(1.0)] + public void ScalingPreservesRepresentableSparseRootCorrections(double sign) + { + // sqrt(2^1000 + d) = 2^500 + d/2^501 + O(d^2/2^1500). + // Scaling d=+/-2^-174 down by 2^1000 erases it, but the root's + // +/-2^-675 low is representable. Prove its rounding interval exactly. + DoubleDouble input = DoubleDouble.FromComponents(Math.ScaleB(1.0, 1000), Math.ScaleB(sign, -174)); + AssertSparseRoot(input, Math.ScaleB(1.0, 500), Math.ScaleB(sign, -675)); + } + + [Fact] + public void SparseRootCorrectionsSurviveScalingAndCorrectionDivisionBoundaries() + { + // Exact squares with both even and odd high exponents. Include scaled + // lows on both sides of the separate-correction dispatch, subnormal + // scaled lows, fully erased scaled lows, and subnormal output lows. + foreach (int exponent in new[] { 0, 2, 100, 1000, 1022 }) + { + foreach (double rootMantissa in new[] { 1.0, 1.5 }) + { + double high = Math.ScaleB(rootMantissa * rootMantissa, exponent); + double rootHigh = Math.ScaleB(rootMantissa, exponent / 2); + foreach (int gap in new[] { -1019, -1020, -1021, -1074, -1075, -1174 }) + { + double low = Math.ScaleB(1.0, exponent + gap); + double rootLow = Math.ScaleB(low, -(exponent / 2)) / (2.0 * rootMantissa); + if (rootLow == 0.0) + { + continue; // This matrix targets representable nonzero corrections. + } + foreach (double sign in new[] { -1.0, 1.0 }) + { + // AssertSparseRoot independently proves this candidate's + // rounding interval; the first-order formula is not an oracle. + AssertSparseRoot(DoubleDouble.FromComponents(high, sign * low), rootHigh, sign * rootLow); + } + } + } + } + } + + [Fact] + public void SparseRootDispatchNeighborsRetainBothSignsWithoutDoubleCounting() + { + // With high=2^1000, lows around 2^-20 scale to the dispatch at 2^-1020. + // Keep the predicted output normal so candidate scaling is exact. + double threshold = Math.ScaleB(1.0, -20); + foreach (double low in new[] { Math.BitDecrement(threshold), threshold, Math.BitIncrement(threshold) }) + { + foreach (double sign in new[] { -1.0, 1.0 }) + { + DoubleDouble input = DoubleDouble.FromComponents(Math.ScaleB(1.0, 1000), sign * low); + AssertSparseRoot(input, Math.ScaleB(1.0, 500), Math.ScaleB(sign * low, -501)); + } + } + } + + private static void AssertSparseRoot(DoubleDouble input, double high, double low) + { + BigInteger x = (Units(input.High) + Units(input.Low)) << 1074; + BigInteger center = Units(high) + Units(low); + BigInteger previous = Units(high) + Units(Math.BitDecrement(low)); + BigInteger following = Units(high) + Units(Math.BitIncrement(low)); + // Strict bounds avoid relying on a tie rule or a rounded sqrt oracle. + ((previous + center) * (previous + center) < (x << 2)).ShouldBeTrue(); + ((center + following) * (center + following) > (x << 2)).ShouldBeTrue(); + DoubleDouble actual = DoubleDouble.Sqrt(input); + actual.High.ShouldBe(high); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(low)); + DoubleDouble facade = DDMath.Sqrt(input); + BitConverter.DoubleToInt64Bits(facade.High).ShouldBe(BitConverter.DoubleToInt64Bits(actual.High)); + BitConverter.DoubleToInt64Bits(facade.Low).ShouldBe(BitConverter.DoubleToInt64Bits(actual.Low)); + AssertSqrtBound(input, actual); + } + private static void AssertSqrtBound(DoubleDouble input, DoubleDouble actual) { DoubleDouble.IsFinite(actual).ShouldBeTrue(); diff --git a/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticBenchmarks.cs b/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticBenchmarks.cs index 8da7594..536177a 100644 --- a/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticBenchmarks.cs +++ b/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticBenchmarks.cs @@ -87,6 +87,18 @@ public class ArithmeticBenchmarks return _doubleDoubleLeft * _doubleDoubleRight; } + [Benchmark(Baseline = true), BenchmarkCategory("Squaring")] + public DoubleDouble DoubleDoubleSquareViaMultiply() + { + return _doubleDoubleLeft * _doubleDoubleLeft; + } + + [Benchmark, BenchmarkCategory("Squaring")] + public DoubleDouble DoubleDoubleSquare() + { + return DoubleDouble.Square(_doubleDoubleLeft); + } + [Benchmark(Baseline = true), BenchmarkCategory("Division")] public double DoubleDivide() { diff --git a/README.md b/README.md index 036e13e..fd2293d 100644 --- a/README.md +++ b/README.md @@ -109,6 +109,15 @@ 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. +- `Square(DoubleDouble)`: forwards to the specialized `DoubleDouble.Square` kernel. + It combines the equal cross terms with FMA, normalizes, then incorporates the + low-component square without losing selected exactly representable tails. + Either signed zero maps to positive zero, either infinity to positive infinity, + and NaN to canonical NaN. It shares multiplication's allocating boundary fallback + and tested `2^-100` relative error plus one minimum binary64 subnormal bound. + Results can differ in low bits from `value * value`; neither route is universally + correctly rounded. No speedup has been measured, and existing kernels have not + been rewritten to use it. - `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 redundant numerator checks and sharing the finite scalar-numerator kernel, @@ -118,7 +127,11 @@ The `DDMath` static class provides: 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 + an unscaled estimate near the exponent limits. Sparse lows endangered by input + scaling or correction division are applied separately with tracked exponents: + for example, `sqrt((2^1000, 2^-174))` retains the low component `2^-675`. + This exceptional correction can reach the allocating boundary machinery; + ordinary inputs retain the original scaled correction path. 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 @@ -195,6 +208,7 @@ The `DDMath` static class provides: using Just.PreciseMath; DoubleDouble root = DoubleDouble.Sqrt(new DoubleDouble(2.0)); +DoubleDouble square = DoubleDouble.Square(root); DoubleDouble inverseRoot = DoubleDouble.InvSqrt(new DoubleDouble(2.0)); DoubleDouble cubeRoot = DoubleDouble.Cbrt(new DoubleDouble(-8.0)); DoubleDouble distance = DoubleDouble.Hypot(new DoubleDouble(3.0), new DoubleDouble(4.0)); @@ -227,6 +241,9 @@ squares/powers of four. Square-root tests also bracket exact root-rounding midpo 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. +Square-root tests separately verify selected sparse output residuals with exact +squared midpoint inequalities, including either low sign, dispatch neighbors, +and corrections that become representable only after rescaling. Cube-root and nth-root tests check `2^-100` relative error for finite nonzero inputs with degree magnitude at least two. Small degrees use exact integer-power @@ -306,6 +323,14 @@ sampled error bounds does not guarantee correct range decisions for every input. - `TryFormat(Span, ...)` implements `ISpanFormattable` with the same formats. It currently allocates via `ToString`; insufficient space returns `false`, writes zero characters, and leaves the destination unchanged. +- `ToStringExact()` emits the complete decimal expansion of the stored component + sum, without rounding, scientific notation, or unnecessary fractional zeros. + Unlike ordinary `ToString`, it always uses **invariant culture**. It preserves + negative zero and emits `NaN`, `Infinity`, or `-Infinity` for special values. + It allocates through exact integer arithmetic; sparse pairs can require over + a thousand characters. Use `DoubleDouble.Parse(text, CultureInfo.InvariantCulture)` + to recover the components. This is not a shortest-round-trip format and does not + change the standard formats' 999-digit precision limit or the parser's length limit. `DoubleDouble` implements `ISignedNumber`, including the inherited `INumberBase` contracts: binary radix, classification, absolute value, magnitude @@ -374,7 +399,7 @@ bool success = DoubleDouble.TryParse("1.25e-2".AsSpan(), CultureInfo.InvariantCu Natural `DDMath.Log`, the complete exponential family on `DoubleDouble`, and all three `DDMath.Pow` overloads are implemented. Logarithms in other bases, generic-math interfaces beyond `ISignedNumber`, -`IFloatingPointConstants`, `IRootFunctions`, and `IExponentialFunctions`, additional text formats/general +`IFloatingPointConstants`, `IRootFunctions`, and `IExponentialFunctions`, additional text formats/shortest- round-trip formatting, and non-arithmetic 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 @@ -400,6 +425,8 @@ on CI workflow runs. BenchmarkDotNet measures arithmetic throughput, dependent-chain latency, and allocations, including comparisons of `DoubleDouble`, `decimal`, and `double`, mixed scalar operations, and exponent-boundary paths. These are performance measurements, not accuracy tests. +Squaring cases compare `DoubleDouble.Square(value)` with `value * value`; their +presence is not evidence of a speedup. After the Release build above, run from the repository root: