From eeeff42be6d319559abefe7a00dc8694002673a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 29 Oct 2025 01:02:31 +0100 Subject: [PATCH] Replace heuristic sqrt with integer algorithm --- .../HeuristicCombinedPrimeTesterTests.cs | 47 +++++++++++++++++++ .../HeuristicCombinedPrimeTester.cs | 37 ++++++++++----- PerfectNumbers.Core/HeuristicPrimeTester.cs | 37 ++++++++++----- 3 files changed, 95 insertions(+), 26 deletions(-) diff --git a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs index 307fdc7..2d76983 100644 --- a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs +++ b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs @@ -1,4 +1,5 @@ using System; +using System.Numerics; using System.Collections.Generic; using PerfectNumbers.Core; using FluentAssertions; @@ -193,4 +194,50 @@ public class HeuristicCombinedPrimeTesterTests (groupBValue % 10U).Should().NotBe(3U); } } + + [Theory] + [Trait("Category", "Fast")] + [InlineData(0UL)] + [InlineData(1UL)] + [InlineData(2UL)] + [InlineData(3UL)] + [InlineData(4UL)] + [InlineData(5UL)] + [InlineData(8UL)] + [InlineData(9UL)] + [InlineData(15UL)] + [InlineData(16UL)] + [InlineData(17UL)] + [InlineData(24UL)] + [InlineData(25UL)] + [InlineData(26UL)] + [InlineData(255UL)] + [InlineData(256UL)] + [InlineData(257UL)] + [InlineData(1_000UL)] + [InlineData(1_048_575UL)] + [InlineData(1_048_576UL)] + [InlineData(1_048_577UL)] + [InlineData(4_000_000UL)] + [InlineData(4_000_001UL)] + [InlineData(4_000_004UL)] + [InlineData(65_535_999UL)] + [InlineData(65_536_000UL)] + [InlineData(65_536_001UL)] + [InlineData(4_294_967_295UL)] + [InlineData(4_294_967_296UL)] + [InlineData(4_294_967_297UL)] + [InlineData(18_446_744_073_709_551_615UL)] + [InlineData(18_446_744_073_709_551_614UL)] + [InlineData(18_446_744_073_709_551_613UL)] + public void ComputeHeuristicSqrt_returns_floor_square_root(ulong value) + { + ulong sqrt = HeuristicCombinedPrimeTester.ComputeHeuristicSqrt(value); + + bool squareBelowOrEqual = (UInt128)sqrt * sqrt <= (UInt128)value; + bool nextSquareAbove = (UInt128)(sqrt + 1UL) * (UInt128)(sqrt + 1UL) > (UInt128)value; + + squareBelowOrEqual.Should().BeTrue(); + nextSquareAbove.Should().BeTrue(); + } } diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 1ae2880..8d29efc 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -306,25 +306,36 @@ Cleanup: internal static ulong ComputeHeuristicSqrt(ulong n) { - ulong sqrt = (ulong)Math.Sqrt(n); - UInt128 square = (UInt128)sqrt * sqrt; - - while (square > n) + if (n < 2UL) { - sqrt--; - square = (UInt128)sqrt * sqrt; + return n; } - ulong next = sqrt + 1UL; - UInt128 nextSquare = (UInt128)next * next; - while (nextSquare <= n) + ulong bit = 1UL << 62; + while (bit > n) { - sqrt = next; - next++; - nextSquare = (UInt128)next * next; + bit >>= 2; } - return sqrt; + ulong result = 0UL; + ulong remainder = n; + while (bit != 0) + { + ulong candidate = result + bit; + if (remainder >= candidate) + { + remainder -= candidate; + result = (result >> 1) + bit; + } + else + { + result >>= 1; + } + + bit >>= 2; + } + + return result; } private static ReadOnlySpan GetCombinedDivisors(byte nMod10) diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index c3983a1..d8fb2cd 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -431,25 +431,36 @@ Cleanup: internal static ulong ComputeHeuristicSqrt(ulong n) { - ulong sqrt = (ulong)Math.Sqrt(n); - UInt128 square = (UInt128)sqrt * sqrt; - - while (square > n) + if (n < 2UL) { - sqrt--; - square = (UInt128)sqrt * sqrt; + return n; } - ulong next = sqrt + 1UL; - UInt128 nextSquare = (UInt128)next * next; - while (nextSquare <= n) + ulong bit = 1UL << 62; + while (bit > n) { - sqrt = next; - next++; - nextSquare = (UInt128)next * next; + bit >>= 2; } - return sqrt; + ulong result = 0UL; + ulong remainder = n; + while (bit != 0) + { + ulong candidate = result + bit; + if (remainder >= candidate) + { + remainder -= candidate; + result = (result >> 1) + bit; + } + else + { + result >>= 1; + } + + bit >>= 2; + } + + return result; } private static ReadOnlySpan GetGroupBEndingOrder(byte nMod10) => nMod10 switch