Replace heuristic sqrt with integer algorithm

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

PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs

@@ -306,25 +306,36 @@ Cleanup:
 
     internal static ulong ComputeHeuristicSqrt(ulong n)
     {
-        ulong sqrt = (ulong)Math.Sqrt(n);
+        if (n < 2UL)
-        UInt128 square = (UInt128)sqrt * sqrt;
-
-        while (square > n)
         {
-            sqrt--;
+            return n;
-            square = (UInt128)sqrt * sqrt;
         }
 
-        ulong next = sqrt + 1UL;
+        ulong bit = 1UL << 62;
-        UInt128 nextSquare = (UInt128)next * next;
+        while (bit > n)
-        while (nextSquare <= n)
         {
-            sqrt = next;
+            bit >>= 2;
-            next++;
-            nextSquare = (UInt128)next * next;
         }
 
-        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<uint> GetCombinedDivisors(byte nMod10)

PerfectNumbers.Core/HeuristicPrimeTester.cs

@@ -431,25 +431,36 @@ Cleanup:
 
     internal static ulong ComputeHeuristicSqrt(ulong n)
     {
-        ulong sqrt = (ulong)Math.Sqrt(n);
+        if (n < 2UL)
-        UInt128 square = (UInt128)sqrt * sqrt;
-
-        while (square > n)
         {
-            sqrt--;
+            return n;
-            square = (UInt128)sqrt * sqrt;
         }
 
-        ulong next = sqrt + 1UL;
+        ulong bit = 1UL << 62;
-        UInt128 nextSquare = (UInt128)next * next;
+        while (bit > n)
-        while (nextSquare <= n)
         {
-            sqrt = next;
+            bit >>= 2;
-            next++;
-            nextSquare = (UInt128)next * next;
         }
 
-        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<byte> GetGroupBEndingOrder(byte nMod10) => nMod10 switch