From 8c2cb14d2b9e2bfffdd2ef3bf365659878a9eda9 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 13 Sep 2026 03:31:39 +0000 Subject: [PATCH] Divide exactly, or to a chosen precision Divide computed its fractional part through a double, so it capped at roughly sixteen significant digits regardless of what the operands carried. On a type whose whole purpose is arbitrary precision that was the wrong ceiling. A quotient with a terminating decimal expansion is now produced exactly, however many digits that takes. A reduced fraction terminates in base ten when its denominator is a product of twos and fives, which is cheap to test and usually cheap to reject: 1 / 8 is 0.125, and 1 / 2^64 keeps all sixty-four decimal places rather than the seventeen a double allows. A repeating quotient is taken to the precision of the wider operand, never fewer than MinimumDivisionPrecision significant digits, with the last digit rounded half away from zero. Whatever the scaled division discards is worth less than one unit of the dropped digits, so it can never carry the rounding decision across the halfway mark; at most it turns an exact tie into something above it, which rounds the same way. A new overload takes the precision explicitly. This changes results. Three tests pinned the old behaviour and are updated: two divisions that were correct only as far as a double reached, and one that compared Exp against 1/e when both happened to be equally imprecise. Exp still routes through a double, so that comparison is now made at the precision Exp actually delivers. It can also produce fewer significant digits than before in one case. The old implementation computed the integer part of a quotient exactly and appended a double-precision fraction, so a quotient with a long integer part came out with more digits than the new default: 924880000e19 divided by -4713100000e-19 gave fifty-four digits where it now gives fifty. That was not a contract so much as a side effect, and an erratic one - the same code gave seventeen digits for a quotient near one, and hundreds for a quotient with a large exponent. A uniform floor with an explicit override is the more predictable trade, and callers who want the extra digits can ask for them. Cost, from the benchmark suite: division is about twice as fast at eight and thirty significant digits, because the double round trip through formatting and parsing is gone, and about three and a half times slower at two hundred, because it now computes two hundred digits instead of discarding them. Verified by checking every division in the differential fuzz harness against the exact rational value rather than against the old implementation: 4,944 terminating quotients exact to the digit and 32,347 repeating ones within half a unit of their last significant digit, inside 871,454 checks with no failures. Also moves the benchmark project's coverage exemption out of the shared .NET workflow and into the project itself. That workflow is duplicated per repository rather than generated from one template - KtsuBuild carries a near-identical copy - so an exclusion added there would have to be repeated everywhere and re-applied on every resync. SonarQubeExclude in the benchmark project travels with the project and needs no workflow to know benchmarks exist, so the shared workflow goes back to its canonical form. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_018eTdSeGPQHGUKf9V3c2yXw --- .github/workflows/dotnet.yml | 2 +- CLAUDE.md | 3 +- .../PreciseNumber.Benchmarks.csproj | 6 + PreciseNumber.Benchmarks/README.md | 10 +- PreciseNumber.Test/PreciseNumberTests.cs | 130 +++++++++++++++- PreciseNumber/PreciseNumber.cs | 141 +++++++++++++++++- README.md | 10 +- 7 files changed, 282 insertions(+), 20 deletions(-) diff --git a/.github/workflows/dotnet.yml b/.github/workflows/dotnet.yml index fd26737..e6c4e35 100644 --- a/.github/workflows/dotnet.yml +++ b/.github/workflows/dotnet.yml @@ -352,7 +352,7 @@ jobs: '/d:sonar.host.url=https://sonarcloud.io' '/d:sonar.projectBaseDir=${{ github.workspace }}' '/d:sonar.cs.vscoveragexml.reportsPaths=coverage/**/coverage.xml' - '/d:sonar.coverage.exclusions=**/*Test*.cs,**/*.Tests.cs,**/*.Tests/**/*,**/*.Benchmarks/**/*,**/obj/**/*,**/*.dll,**/NativeExports.cs' + '/d:sonar.coverage.exclusions=**/*Test*.cs,**/*.Tests.cs,**/*.Tests/**/*,**/obj/**/*,**/*.dll,**/NativeExports.cs' '/d:sonar.cs.vstest.reportsPaths=coverage/**/*.trx' '/d:sonar.exclusions=**/NativeExports.cs' ) diff --git a/CLAUDE.md b/CLAUDE.md index e4872c2..141b003 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -33,7 +33,8 @@ dotnet run -c Release --project PreciseNumber.Benchmarks -- --filter '*' --job s ### Key Design Patterns - Factory methods `CreateFromInteger()` and `CreateFromFloatingPoint()` handle type-specific conversion logic -- Arithmetic operations use `MakeCommonized()` to align exponents before calculation +- Addition, subtraction and modulus align exponents before calculating; multiplication and division work on the significands directly +- `Divide` is exact when the quotient terminates, and otherwise rounds to a precision that never falls below the wider operand or `MinimumDivisionPrecision`. `Exp` and non-integer `Pow` still route through `double` - The `sanitize` constructor parameter controls whether trailing zeros are removed (default: true) - Constants (`Zero`, `One`, `Pi`, `E`, `Tau`) are pre-computed static instances diff --git a/PreciseNumber.Benchmarks/PreciseNumber.Benchmarks.csproj b/PreciseNumber.Benchmarks/PreciseNumber.Benchmarks.csproj index 54f9cc8..c6aa6ec 100644 --- a/PreciseNumber.Benchmarks/PreciseNumber.Benchmarks.csproj +++ b/PreciseNumber.Benchmarks/PreciseNumber.Benchmarks.csproj @@ -16,6 +16,12 @@ internal copy here through InternalsVisibleTo - referencing it again would make every polyfilled member ambiguous. --> $(NoWarn);KTSU0001 + + true diff --git a/PreciseNumber.Benchmarks/README.md b/PreciseNumber.Benchmarks/README.md index cce8188..f9dd876 100644 --- a/PreciseNumber.Benchmarks/README.md +++ b/PreciseNumber.Benchmarks/README.md @@ -57,11 +57,11 @@ operation returns a new instance, so a change that avoids an intermediate value `Allocated` column before it shows up in `Mean`. A comparison that allocates at all is a regression — none of them should. -One entry is easy to misread. `Divide` barely moves between 8 and 200 digits, and at 200 digits -it comes out faster than `Add`. That is not division being efficient — it computes its fractional -part through a `double`, so it discards everything past roughly the 17th significant digit and -the result it constructs is small no matter how large the operands were. Read it as a measure of -how much precision the operation throws away, not how fast it is. +`Divide` is the one row whose cost is not driven by the operands alone. It produces a terminating +quotient exactly and a repeating one to a chosen precision, so the `Digits` column moves it twice +over: wider operands are more work to divide, and they also raise the precision the quotient is +taken to. A number that looks expensive at 200 digits is doing proportionally more work, not doing +the same work badly. Benchmark operands come from a fixed digit pattern rather than a random source, so two runs on the same machine measure the same work. Numbers are still only comparable within a single run on diff --git a/PreciseNumber.Test/PreciseNumberTests.cs b/PreciseNumber.Test/PreciseNumberTests.cs index e35d0c8..4ef448c 100644 --- a/PreciseNumber.Test/PreciseNumberTests.cs +++ b/PreciseNumber.Test/PreciseNumberTests.cs @@ -1031,8 +1031,10 @@ public void TestOperatorDivide() PreciseNumber number1 = PreciseNumber.CreateFromComponents(-2, 12345); PreciseNumber number2 = PreciseNumber.CreateFromComponents(-3, 678); PreciseNumber result = number1 / number2; - Assert.AreEqual(BigInteger.Parse("18207964601769911504"), result.Significand); - Assert.AreEqual(-17, result.Exponent); + // 123.45 / 0.678 does not terminate, so it comes out at the default precision. + Assert.AreEqual(BigInteger.Parse("18207964601769911504424778761061946902654867256637"), result.Significand); + Assert.AreEqual(-47, result.Exponent); + Assert.AreEqual(PreciseNumber.MinimumDivisionPrecision, result.SignificantDigits); } [TestMethod] @@ -1603,7 +1605,11 @@ public void TestExpWithPositivePower() public void TestExpWithNegativePower() { PreciseNumber result = PreciseNumber.Exp(-1.ToPreciseNumber()); - PreciseNumber expected = PreciseNumber.One / PreciseNumber.E; // e^-1 = 1/e + + // Exp routes through a double, so 1/e is taken to the precision Exp actually delivers + // rather than the full precision Divide is now capable of. + PreciseNumber expected = PreciseNumber.Divide(PreciseNumber.One, PreciseNumber.E, result.SignificantDigits); + Assert.AreEqual(expected, result); } @@ -2000,6 +2006,122 @@ public void TestCountDigitsMatchesDecimalText() } } + [TestMethod] + public void TestDivideIsExactWhenTheQuotientTerminates() + { + // A quotient terminates exactly when the reduced denominator is a product of twos and + // fives, and those come out exact however many digits that takes. + (int Numerator, int Denominator, string Expected)[] cases = + [ + (1, 2, "0.5"), + (1, 4, "0.25"), + (1, 5, "0.2"), + (1, 8, "0.125"), + (1, 10, "0.1"), + (1, 16, "0.0625"), + (1, 20, "0.05"), + (1, 25, "0.04"), + (3, 8, "0.375"), + (-3, 8, "-0.375"), + (7, 1, "7"), + ]; + + foreach ((int numerator, int denominator, string expected) in cases) + { + PreciseNumber quotient = numerator.ToPreciseNumber() / denominator.ToPreciseNumber(); + Assert.AreEqual(expected, quotient.ToString(CultureInfo.InvariantCulture), $"{numerator}/{denominator}"); + } + } + + [TestMethod] + public void TestDivideKeepsEveryDigitOfALongTerminatingQuotient() + { + // 2^-64 terminates, but only after 64 decimal places - far more than the default + // precision would allow, and far more than a double could carry. + PreciseNumber quotient = PreciseNumber.One / PreciseNumber.CreateFromComponents(0, BigInteger.Pow(2, 64)); + + Assert.AreEqual( + "0.0000000000000000000542101086242752217003726400434970855712890625", + quotient.ToString(CultureInfo.InvariantCulture)); + + // Exact means exact: multiplying back reproduces the dividend. + Assert.AreEqual(PreciseNumber.One, quotient * PreciseNumber.CreateFromComponents(0, BigInteger.Pow(2, 64))); + } + + [TestMethod] + public void TestDivideProducesTheRequestedPrecisionWhenTheQuotientRepeats() + { + PreciseNumber one = PreciseNumber.One; + PreciseNumber three = 3.ToPreciseNumber(); + + Assert.AreEqual("0.3", PreciseNumber.Divide(one, three, 1).ToString(CultureInfo.InvariantCulture)); + Assert.AreEqual("0.333", PreciseNumber.Divide(one, three, 3).ToString(CultureInfo.InvariantCulture)); + Assert.AreEqual(20, PreciseNumber.Divide(one, three, 20).SignificantDigits); + Assert.AreEqual(200, PreciseNumber.Divide(one, three, 200).SignificantDigits); + } + + [TestMethod] + public void TestDivideRoundsHalfAwayFromZero() + { + PreciseNumber two = 2.ToPreciseNumber(); + PreciseNumber three = 3.ToPreciseNumber(); + PreciseNumber five = 5.ToPreciseNumber(); + PreciseNumber nine = 9.ToPreciseNumber(); + + // 0.666... rounds up, 0.333... rounds down, and the sign does not change which way. + Assert.AreEqual("0.667", PreciseNumber.Divide(two, three, 3).ToString(CultureInfo.InvariantCulture)); + Assert.AreEqual("0.333", PreciseNumber.Divide(PreciseNumber.One, three, 3).ToString(CultureInfo.InvariantCulture)); + Assert.AreEqual("-0.667", PreciseNumber.Divide(-two, three, 3).ToString(CultureInfo.InvariantCulture)); + Assert.AreEqual("-0.333", PreciseNumber.Divide(-PreciseNumber.One, three, 3).ToString(CultureInfo.InvariantCulture)); + + // 0.555... cut after one digit sits just above the halfway mark, so it rounds away. + Assert.AreEqual("0.6", PreciseNumber.Divide(five, nine, 1).ToString(CultureInfo.InvariantCulture)); + Assert.AreEqual("-0.6", PreciseNumber.Divide(-five, nine, 1).ToString(CultureInfo.InvariantCulture)); + } + + [TestMethod] + public void TestDivideNeverReducesTheOperandsPrecision() + { + // An operand carrying more digits than the floor pulls the quotient up to match it. + PreciseNumber wide = PreciseNumber.Parse(new string('7', 120), CultureInfo.InvariantCulture); + PreciseNumber three = 3.ToPreciseNumber(); + + Assert.AreEqual(120, (wide / three).SignificantDigits); + Assert.AreEqual(120, (three / wide).SignificantDigits); + Assert.AreEqual(PreciseNumber.MinimumDivisionPrecision, (three / 7.ToPreciseNumber()).SignificantDigits); + } + + [TestMethod] + public void TestDivideRejectsNonPositivePrecision() + { + PreciseNumber one = PreciseNumber.One; + PreciseNumber three = 3.ToPreciseNumber(); + + Assert.ThrowsExactly(() => PreciseNumber.Divide(one, three, 0)); + Assert.ThrowsExactly(() => PreciseNumber.Divide(one, three, -1)); + } + + [TestMethod] + public void TestDivideZeroDividend() + { + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Zero / 3.ToPreciseNumber()); + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Zero / PreciseNumber.NegativeOne); + } + + [TestMethod] + public void TestDivideRoundTripsThroughMultiplication() + { + // For a terminating quotient the round trip is exact, which is the strongest statement + // that can be made about a division. + foreach (int denominator in new[] { 2, 4, 5, 8, 10, 16, 20, 25, 32, 50, 64, 100, 125, 128 }) + { + PreciseNumber divisor = denominator.ToPreciseNumber(); + PreciseNumber quotient = 123456789.ToPreciseNumber() / divisor; + + Assert.AreEqual(123456789.ToPreciseNumber(), quotient * divisor, $"123456789/{denominator}"); + } + } + [TestMethod] public void TestPow10IsCorrectAcrossCacheBoundaries() { @@ -2057,7 +2179,7 @@ public void TestDivideWithExponentsBeyondDoubleRange() PreciseNumber result = left / right; - Assert.AreEqual("2.3333333333333333", result.ToString(CultureInfo.InvariantCulture)); + Assert.AreEqual("2.3333333333333333333333333333333333333333333333333", result.ToString(CultureInfo.InvariantCulture)); } [TestMethod] diff --git a/PreciseNumber/PreciseNumber.cs b/PreciseNumber/PreciseNumber.cs index 912389c..f7512fb 100644 --- a/PreciseNumber/PreciseNumber.cs +++ b/PreciseNumber/PreciseNumber.cs @@ -44,6 +44,12 @@ public record PreciseNumber /// private const string InvalidFormatMessage = "Input string was not in a correct format."; + /// + /// The fewest significant digits produces + /// when a quotient does not terminate. + /// + public const int MinimumDivisionPrecision = 50; + /// /// Pre-computed powers of ten, grown on demand. Declared before any other static state so /// that the static constants below can rely on it while they are being initialized. @@ -1291,30 +1297,149 @@ public static PreciseNumber Multiply(PreciseNumber left, PreciseNumber right) /// The number to divide. /// The number to divide by. /// The result of the division. + /// Thrown when is zero. + /// + /// A quotient whose decimal expansion terminates is produced exactly, however many digits that + /// takes. One that repeats is produced to the precision of the wider operand, and never fewer + /// than significant digits, with the last digit rounded + /// half away from zero. Use to choose + /// that precision. + /// public static PreciseNumber Divide(PreciseNumber left, PreciseNumber right) { Ensure.NotNull(left); Ensure.NotNull(right); + // Dividing must not silently discard precision the operands already carry. + int significantDigits = Math.Max( + Math.Max(left.SignificantDigits, right.SignificantDigits), + MinimumDivisionPrecision); + + return Divide(left, right, significantDigits); + } + + /// + /// Divides one number by another, to a chosen number of significant digits. + /// + /// The number to divide. + /// The number to divide by. + /// + /// The number of significant digits to produce when the quotient does not terminate. A quotient + /// that does terminate is exact regardless of this value. + /// + /// The result of the division. + /// Thrown when is zero. + /// Thrown when is less than one. + public static PreciseNumber Divide(PreciseNumber left, PreciseNumber right, int significantDigits) + { + Ensure.NotNull(left); + Ensure.NotNull(right); + + if (significantDigits < 1) + { + throw new ArgumentOutOfRangeException(nameof(significantDigits), significantDigits, "At least one significant digit is required."); + } + if (right.Significand.IsZero) { throw new DivideByZeroException(); } - if (Compare(left, right) == 0) + if (left.Significand.IsZero) { - return One; + return Zero; } - (BigInteger commonLeft, BigInteger commonRight, _) = CommonizeSignificands(left, right); + BigInteger numerator = left.Significand; + BigInteger denominator = right.Significand; + int exponent = left.Exponent - right.Exponent; - BigInteger integerComponent = BigInteger.DivRem(commonLeft, commonRight, out BigInteger remainder); + // Carry the sign on the numerator so the denominator can be factorized as a positive value. + if (denominator.Sign < 0) + { + numerator = -numerator; + denominator = -denominator; + } - // The common power of ten cancels between the remainder and the divisor, so it is left out - // entirely; including it would overflow to infinity for large exponents. - double fractionalComponent = double.CreateTruncating(remainder) / double.CreateTruncating(commonRight); + return TryDivideExactly(numerator, denominator, exponent, out PreciseNumber? exact) + ? exact + : DivideToPrecision(numerator, denominator, exponent, significantDigits); + } + + /// + /// Divides exactly, when the quotient has a terminating decimal expansion. + /// + /// The numerator, carrying the sign of the quotient. + /// The denominator, which must be positive. + /// The exponent the quotient's significand sits at. + /// The exact quotient, when there is one. + /// true if the quotient terminates and is exact; otherwise false. + private static bool TryDivideExactly(BigInteger numerator, BigInteger denominator, int exponent, [NotNullWhen(true)] out PreciseNumber? result) + { + // A fraction terminates in base ten exactly when its denominator is 2^twos * 5^fives. Most + // denominators are rejected by the first remainder test, which is why this is worth trying + // before falling back to a rounded quotient. + int twos = (int)BigInteger.TrailingZeroCount(denominator); + BigInteger remaining = denominator >> twos; + + int fives = 0; + while ((remaining % 5).IsZero) + { + remaining /= 5; + fives++; + } + + if (!remaining.IsOne) + { + result = null; + return false; + } + + // 1 / (2^p * 5^q) == (2^(k-p) * 5^(k-q)) / 10^k, where k is the larger of p and q. + int scale = Math.Max(twos, fives); + BigInteger significand = numerator * BigInteger.Pow(2, scale - twos) * BigInteger.Pow(5, scale - fives); + + result = new PreciseNumber(exponent - scale, significand); + return true; + } + + /// + /// Divides to a fixed number of significant digits, rounding the last of them half away from zero. + /// + /// The numerator, carrying the sign of the quotient. + /// The denominator, which must be positive. + /// The exponent the quotient's significand sits at. + /// The number of significant digits to produce. + /// The rounded quotient. + private static PreciseNumber DivideToPrecision(BigInteger numerator, BigInteger denominator, int exponent, int significantDigits) + { + // A quotient has either digits(numerator) - digits(denominator) digits or one more, so + // scaling by this much leaves at least one digit past the ones being kept: the digit the + // rounding decision is made on. + int scale = significantDigits + 1 - CountDigits(numerator) + CountDigits(denominator); + BigInteger scaled = scale > 0 ? numerator * Pow10(scale) : numerator; + int scaledExponent = exponent - Math.Max(scale, 0); + + BigInteger quotient = scaled / denominator; + int excess = CountDigits(quotient) - significantDigits; + + if (excess <= 0) + { + return new PreciseNumber(scaledExponent, quotient); + } + + BigInteger divisor = Pow10(excess); + BigInteger kept = BigInteger.DivRem(quotient, divisor, out BigInteger dropped); + + // Round half away from zero. Whatever the division above discarded is worth less than one + // unit of the dropped digits, so it can never carry the comparison across the halfway mark; + // at most it turns an exact tie into something above it, which rounds the same way. + if (BigInteger.Abs(dropped) * 2 >= divisor) + { + kept += quotient.Sign; + } - return new PreciseNumber(0, integerComponent) + fractionalComponent.ToPreciseNumber(); + return new PreciseNumber(scaledExponent + excess, kept); } /// diff --git a/README.md b/README.md index 8b9ac84..f2d1abf 100644 --- a/README.md +++ b/README.md @@ -351,9 +351,17 @@ You can control precision using: - **ReduceSignificance()**: Reduces to a specific number of significant digits +- **Divide(left, right, significantDigits)**: Chooses the precision of a quotient + +Division produces a terminating quotient exactly, however many digits that takes — `1 / 8` is +`0.125`, and `1 / 2^64` keeps all 64 decimal places. A repeating quotient is produced to the +precision of the wider operand, never fewer than `MinimumDivisionPrecision` (50) significant +digits, with the last digit rounded half away from zero. Pass an explicit precision to the +three-argument overload when you want something other than that. + ## Limitations -- Operations that inherently require approximation (like certain roots or logarithms) fall back to `double` precision for calculation +- `Exp()`, and `Pow()` with a non-integer power, are computed through `double` and are therefore limited to its precision. Addition, subtraction, multiplication and division are not - Conversion to standard types may throw `OverflowException` if the value is too large