diff --git a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs
index 0468bfd..4c9f2a7 100644
--- a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs
+++ b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs
@@ -158,36 +158,16 @@ public class MulMod64Benchmarks
}
///
- /// Uses the GPU-friendly helper; dense inputs now
- /// cost 24.4–61.9 ns, while NearFullRange and SparseOperands finish near 7.27 ns and ZeroOperands stay baseline-fast at
- /// 3.59 ns.
+ /// Uses the production helper shared by CPU and GPU callers.
///
///
- /// Observed means: CrossWordBlend 24.422 ns (6.86×), MixedBitPattern 61.850 ns, NearFullRange 7.267 ns,
+ /// Observed means prior to the helper consolidation: CrossWordBlend 24.422 ns, MixedBitPattern 61.850 ns, NearFullRange 7.267 ns,
/// PrimeSizedModulus 55.114 ns, SparseOperands 7.261 ns, ZeroOperands 3.587 ns.
///
[Benchmark]
- public ulong GpuCompatibleBaseline()
+ public ulong UnifiedMulMod()
{
- // TODO: Callers in production should migrate to ULongExtensions.MulMod64 where GPU parity
- // is not required; that baseline stayed 6–17× faster for dense 64-bit operands.
- return Input.Left.MulMod64GpuCompatible(Input.Right, Input.Modulus);
- }
-
- ///
- /// Uses the deferred native-modulo helper ;
- /// excels when operands are tiny (2.01 ns on ZeroOperands) yet remains the slowest choice on dense data at 146–295 ns.
- ///
- ///
- /// Observed means: CrossWordBlend 146.559 ns (40.63×), MixedBitPattern 255.555 ns, NearFullRange 21.390 ns,
- /// PrimeSizedModulus 294.486 ns, SparseOperands 18.182 ns, ZeroOperands 2.009 ns.
- ///
- [Benchmark]
- public ulong GpuCompatibleDeferred()
- {
- // TODO: Retire the deferred GPU shim from runtime paths and keep it only for benchmarks;
- // ULongExtensions.MulMod64 avoids the 6×–82× slowdown on real-world operand mixes.
- return Input.Left.MulMod64GpuCompatibleDeferred(Input.Right, Input.Modulus);
+ return Input.Left.MulMod64(Input.Right, Input.Modulus);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs
index 06a1bf1..b08f418 100644
--- a/EvenPerfectBitScanner/Program.cs
+++ b/EvenPerfectBitScanner/Program.cs
@@ -1030,7 +1030,7 @@ internal static class Program
Console.WriteLine(" --use-order test primality via q order");
Console.WriteLine(" --workaround-mod avoid '%' operator on the GPU");
// mod-automaton removed
- Console.WriteLine(" --filter-p= process only p from previous run results (required for --mersenne=bydivisor)");
+ Console.WriteLine(" --filter-p= process only p from previous run results (required for --mersenne=bydivisor; primes must be strictly increasing)");
Console.WriteLine(" --write-batch-size= overwrite frequency of disk writes (default 100 lines)");
Console.WriteLine(" --gcd-filter enable early sieve based on GCD");
Console.WriteLine(" --help, -help, --?, -?, /? show this help message");
diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs
index 43850f8..4b01a18 100644
--- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs
+++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs
@@ -1,5 +1,6 @@
using System;
using FluentAssertions;
+using PerfectNumbers.Core;
using PerfectNumbers.Core.Cpu;
using Xunit;
@@ -7,20 +8,23 @@ namespace PerfectNumbers.Core.Tests;
public class MersenneNumberDivisorCpuTesterTests
{
+ private const ulong LargeDivisor1 = 4_000_169UL;
+ private const ulong LargeDivisor2 = 4_000_561UL;
+
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_tester_tracks_divisors_across_primes()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
- tester.ConfigureFromMaxPrime(11UL);
+ tester.ConfigureFromMaxPrime(83UL);
- tester.IsPrime(5UL, out bool divisorsExhausted).Should().BeTrue();
+ tester.IsPrime(59UL, out bool divisorsExhausted).Should().BeTrue();
divisorsExhausted.Should().BeTrue();
- tester.IsPrime(7UL, out divisorsExhausted).Should().BeTrue();
+ tester.IsPrime(61UL, out divisorsExhausted).Should().BeTrue();
divisorsExhausted.Should().BeTrue();
- tester.IsPrime(11UL, out divisorsExhausted).Should().BeFalse();
+ tester.IsPrime(67UL, out divisorsExhausted).Should().BeTrue();
divisorsExhausted.Should().BeTrue();
}
@@ -29,66 +33,94 @@ public class MersenneNumberDivisorCpuTesterTests
public void ByDivisor_prime_test_limit_disables_divisor_exhaustion()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
- tester.ConfigureFromMaxPrime(11UL);
+ tester.ConfigureFromMaxPrime(83UL);
- tester.IsPrime(5UL, out bool divisorsExhausted, TimeSpan.Zero).Should().BeTrue();
+ tester.IsPrime(59UL, out bool divisorsExhausted, TimeSpan.Zero).Should().BeTrue();
divisorsExhausted.Should().BeFalse();
}
[Fact]
[Trait("Category", "Fast")]
- public void ByDivisor_session_checks_divisors_across_primes()
+ public void ByDivisor_session_checks_divisors_across_exponents()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
- tester.ConfigureFromMaxPrime(13UL);
+ tester.ConfigureFromMaxPrime(197UL);
using var session = tester.CreateDivisorSession();
- ulong[] primes = { 5UL, 7UL, 11UL, 13UL };
- byte[] hits = new byte[primes.Length];
+ ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL, 333_380UL };
+ byte[] hits = new byte[exponents.Length];
- ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL));
- session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits);
- hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
+ MontgomeryDivisorData divisor1 = MontgomeryDivisorData.FromModulus(LargeDivisor1);
+ ulong cycle1 = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor1, divisor1);
+ session.CheckDivisor(LargeDivisor1, divisor1, cycle1, exponents, hits);
+
+ byte[] expected1 = new byte[exponents.Length];
+ for (int i = 0; i < exponents.Length; i++)
+ {
+ expected1[i] = exponents[i].Pow2ModWindowedCpu(divisor1.Modulus) == 1UL ? (byte)1 : (byte)0;
+ }
+
+ hits.Should().Equal(expected1);
Array.Fill(hits, (byte)0);
- ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL));
- session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits);
- hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 });
+ MontgomeryDivisorData divisor2 = MontgomeryDivisorData.FromModulus(LargeDivisor2);
+ ulong cycle2 = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor2, divisor2);
+ session.CheckDivisor(LargeDivisor2, divisor2, cycle2, exponents, hits);
+
+ byte[] expected2 = new byte[exponents.Length];
+ for (int i = 0; i < exponents.Length; i++)
+ {
+ expected2[i] = exponents[i].Pow2ModWindowedCpu(divisor2.Modulus) == 1UL ? (byte)1 : (byte)0;
+ }
+
+ hits.Should().Equal(expected2);
}
[Fact]
[Trait("Category", "Fast")]
- public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite()
+ public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor1_as_composite()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
- tester.ConfigureFromMaxPrime(43UL);
+ tester.ConfigureFromMaxPrime(197UL);
using var session = tester.CreateDivisorSession();
- ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
+ MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(LargeDivisor1);
+ ulong cycle = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor1, divisorData);
+
+ ulong[] exponents = { 59UL, 67UL, cycle, cycle * 2UL, cycle + 1UL };
byte[] hits = new byte[exponents.Length];
+ session.CheckDivisor(LargeDivisor1, divisorData, cycle, exponents, hits);
- MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
- ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
- session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
+ byte[] expected = new byte[exponents.Length];
+ for (int i = 0; i < exponents.Length; i++)
+ {
+ expected[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0;
+ }
- hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
+ hits.Should().Equal(expected);
}
[Fact]
[Trait("Category", "Fast")]
- public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite()
+ public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor2_as_composite()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
- tester.ConfigureFromMaxPrime(61UL);
+ tester.ConfigureFromMaxPrime(197UL);
using var session = tester.CreateDivisorSession();
- ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
+ MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(LargeDivisor2);
+ ulong cycle = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor2, divisorData);
+
+ ulong[] exponents = { 61UL, 73UL, cycle, cycle * 2UL, cycle - 1UL };
byte[] hits = new byte[exponents.Length];
+ session.CheckDivisor(LargeDivisor2, divisorData, cycle, exponents, hits);
- MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
- ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
- session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
+ byte[] expected = new byte[exponents.Length];
+ for (int i = 0; i < exponents.Length; i++)
+ {
+ expected[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0;
+ }
- hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
+ hits.Should().Equal(expected);
}
}
diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs
index aceaea7..3ec73c6 100644
--- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs
+++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs
@@ -12,12 +12,12 @@ namespace PerfectNumbers.Core.Tests;
public class MersenneNumberDivisorGpuTesterTests
{
[Theory]
- [InlineData(33UL, 7UL, true)]
- [InlineData(35UL, 31UL, true)]
- [InlineData(37UL, 223UL, true)]
- [InlineData(33UL, 13UL, false)]
- [InlineData(35UL, 73UL, false)]
- [InlineData(37UL, 227UL, false)]
+ [InlineData(117_652UL, 4_000_169UL, true)]
+ [InlineData(333_380UL, 4_000_561UL, true)]
+ [InlineData(59UL, 4_000_169UL, false)]
+ [InlineData(61UL, 4_000_561UL, false)]
+ [InlineData(333_379UL, 4_000_561UL, false)]
+ [InlineData(117_653UL, 4_000_169UL, false)]
[Trait("Category", "Fast")]
public void IsDivisible_returns_expected(ulong exponent, ulong divisor, bool expected)
{
@@ -152,23 +152,39 @@ public class MersenneNumberDivisorGpuTesterTests
[Fact]
[Trait("Category", "Fast")]
- public void ByDivisor_session_checks_divisors_across_primes()
+ public void ByDivisor_session_checks_divisors_across_exponents()
{
var tester = new MersenneNumberDivisorByDivisorGpuTester();
- tester.ConfigureFromMaxPrime(43UL);
+ tester.ConfigureFromMaxPrime(197UL);
using var session = tester.CreateDivisorSession();
- ulong[] primes = { 31UL, 37UL, 41UL, 43UL };
- byte[] hits = new byte[primes.Length];
+ ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL, 333_380UL };
+ byte[] hits = new byte[exponents.Length];
- ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
- session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
- hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0 });
+ MontgomeryDivisorData divisor1 = MontgomeryDivisorData.FromModulus(4_000_169UL);
+ ulong cycle1 = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisor1);
+ session.CheckDivisor(4_000_169UL, divisor1, cycle1, exponents, hits);
+
+ byte[] expected1 = new byte[exponents.Length];
+ for (int i = 0; i < exponents.Length; i++)
+ {
+ expected1[i] = exponents[i].Pow2ModWindowedCpu(divisor1.Modulus) == 1UL ? (byte)1 : (byte)0;
+ }
+
+ hits.Should().Equal(expected1);
Array.Fill(hits, (byte)0);
- ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLength(13367UL, MontgomeryDivisorData.FromModulus(13367UL));
- session.CheckDivisor(13367UL, MontgomeryDivisorData.FromModulus(13367UL), cycle13367, primes, hits);
- hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
+ MontgomeryDivisorData divisor2 = MontgomeryDivisorData.FromModulus(4_000_561UL);
+ ulong cycle2 = MersenneDivisorCycles.CalculateCycleLength(4_000_561UL, divisor2);
+ session.CheckDivisor(4_000_561UL, divisor2, cycle2, exponents, hits);
+
+ byte[] expected2 = new byte[exponents.Length];
+ for (int i = 0; i < exponents.Length; i++)
+ {
+ expected2[i] = exponents[i].Pow2ModWindowedCpu(divisor2.Modulus) == 1UL ? (byte)1 : (byte)0;
+ }
+
+ hits.Should().Equal(expected2);
}
[Fact]
@@ -180,122 +196,126 @@ public class MersenneNumberDivisorGpuTesterTests
GpuBatchSize = 2,
};
- tester.ConfigureFromMaxPrime(47UL);
+ tester.ConfigureFromMaxPrime(197UL);
using var session = tester.CreateDivisorSession();
- ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
- byte[] hits = new byte[primes.Length];
+ ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL };
+ byte[] hits = new byte[exponents.Length];
- ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
- session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
+ MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(4_000_169UL);
+ ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisor);
+ session.CheckDivisor(4_000_169UL, divisor, cycle, exponents, hits);
- hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
+ byte[] expected = new byte[exponents.Length];
+ for (int i = 0; i < exponents.Length; i++)
+ {
+ expected[i] = exponents[i].Pow2ModWindowedCpu(divisor.Modulus) == 1UL ? (byte)1 : (byte)0;
+ }
+
+ hits.Should().Equal(expected);
}
[Fact]
[Trait("Category", "Fast")]
- public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite()
+ public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor1_as_composite()
{
var tester = new MersenneNumberDivisorByDivisorGpuTester
{
GpuBatchSize = 5,
};
- tester.ConfigureFromMaxPrime(43UL);
+ tester.ConfigureFromMaxPrime(197UL);
using var session = tester.CreateDivisorSession();
- ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
+ MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(4_000_169UL);
+ ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisorData);
+
+ ulong[] exponents = { 59UL, 67UL, cycle, cycle * 2UL, cycle + 1UL };
byte[] hits = new byte[exponents.Length];
- MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
- ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
-
- session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
+ session.CheckDivisor(4_000_169UL, divisorData, cycle, exponents, hits);
var hostBufferField = typeof(MersenneNumberDivisorByDivisorGpuTester.DivisorScanSession)
.GetField("_hostBuffer", BindingFlags.NonPublic | BindingFlags.Instance)!;
ulong[] hostBuffer = (ulong[])hostBufferField.GetValue(session)!;
Span residues = hostBuffer.AsSpan(0, exponents.Length);
- byte[] computedHits = new byte[exponents.Length];
+ byte[] expected = new byte[exponents.Length];
for (int i = 0; i < residues.Length; i++)
{
- computedHits[i] = residues[i] == divisorData.MontgomeryOne ? (byte)1 : (byte)0;
+ expected[i] = residues[i] == 1UL ? (byte)1 : (byte)0;
}
- computedHits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
+ byte[] direct = new byte[exponents.Length];
+ for (int i = 0; i < exponents.Length; i++)
+ {
+ direct[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0;
+ }
+
+ expected.Should().Equal(direct);
+ hits.Should().Equal(direct);
}
[Fact]
[Trait("Category", "Fast")]
- public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite()
+ public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor2_as_composite()
{
var tester = new MersenneNumberDivisorByDivisorGpuTester
{
GpuBatchSize = 5,
};
- tester.ConfigureFromMaxPrime(61UL);
+ tester.ConfigureFromMaxPrime(197UL);
using var session = tester.CreateDivisorSession();
- ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
+ MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(4_000_561UL);
+ ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_561UL, divisorData);
+
+ ulong[] exponents = { 61UL, 73UL, cycle, cycle * 2UL, cycle - 1UL };
byte[] hits = new byte[exponents.Length];
- MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
- ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
-
- session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
+ session.CheckDivisor(4_000_561UL, divisorData, cycle, exponents, hits);
var hostBufferField = typeof(MersenneNumberDivisorByDivisorGpuTester.DivisorScanSession)
.GetField("_hostBuffer", BindingFlags.NonPublic | BindingFlags.Instance)!;
ulong[] hostBuffer = (ulong[])hostBufferField.GetValue(session)!;
Span residues = hostBuffer.AsSpan(0, exponents.Length);
- byte[] computedHits = new byte[exponents.Length];
+ byte[] expected = new byte[exponents.Length];
for (int i = 0; i < residues.Length; i++)
{
- computedHits[i] = residues[i] == divisorData.MontgomeryOne ? (byte)1 : (byte)0;
+ expected[i] = residues[i] == 1UL ? (byte)1 : (byte)0;
}
- computedHits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
+ byte[] direct = new byte[exponents.Length];
+ for (int i = 0; i < exponents.Length; i++)
+ {
+ direct[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0;
+ }
+
+ expected.Should().Equal(direct);
+ hits.Should().Equal(direct);
}
[Fact]
[Trait("Category", "Fast")]
- public void ByDivisor_gpu_tester_skips_divisors_excluded_by_small_cycle_generation()
+ public void ByDivisor_gpu_tester_generates_cycles_for_large_divisors()
{
- var cycles = MersenneDivisorCycles.Shared;
- var tableField = typeof(MersenneDivisorCycles).GetField("_table", BindingFlags.NonPublic | BindingFlags.Instance)!;
- var smallCyclesField = typeof(MersenneDivisorCycles).GetField("_smallCycles", BindingFlags.NonPublic | BindingFlags.Instance)!;
+ var tester = new MersenneNumberDivisorByDivisorGpuTester();
+ tester.ConfigureFromMaxPrime(197UL);
- var originalTable = (List<(ulong Divisor, ulong Cycle)>)tableField.GetValue(cycles)!;
- var originalSmall = (ulong[]?)smallCyclesField.GetValue(cycles);
+ using var session = tester.CreateDivisorSession();
+ ulong[] exponents = { 59UL, 67UL, 117_652UL };
+ byte[] hits = new byte[exponents.Length];
- try
+ MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(4_000_169UL);
+ ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisorData);
+
+ session.CheckDivisor(4_000_169UL, divisorData, cycle, exponents, hits);
+
+ byte[] expected = new byte[exponents.Length];
+ for (int i = 0; i < exponents.Length; i++)
{
- var patchedTable = new List<(ulong Divisor, ulong Cycle)>();
- ulong[] patchedSmall = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1];
-
- MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(191UL);
- patchedSmall[191] = MersenneDivisorCycles.CalculateCycleLength(191UL, divisorData);
-
- tableField.SetValue(cycles, patchedTable);
- smallCyclesField.SetValue(cycles, patchedSmall);
- DivisorCycleCache.Shared.RefreshSnapshot();
-
- var tester = new MersenneNumberDivisorByDivisorGpuTester();
- tester.ConfigureFromMaxPrime(19UL);
-
- typeof(MersenneNumberDivisorByDivisorGpuTester)
- .GetField("_divisorLimit", BindingFlags.NonPublic | BindingFlags.Instance)!
- .SetValue(tester, 200UL);
-
- tester.IsPrime(19UL, out bool divisorsExhausted).Should().BeTrue();
- divisorsExhausted.Should().BeTrue();
- }
- finally
- {
- tableField.SetValue(cycles, originalTable);
- smallCyclesField.SetValue(cycles, originalSmall);
- DivisorCycleCache.Shared.RefreshSnapshot();
- GpuContextPool.DisposeAll();
+ expected[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0;
}
+
+ hits.Should().Equal(expected);
}
[Fact]
@@ -307,16 +327,23 @@ public class MersenneNumberDivisorGpuTesterTests
BatchSize = 2,
};
- tester.ConfigureFromMaxPrime(47UL);
+ tester.ConfigureFromMaxPrime(197UL);
using var session = tester.CreateDivisorSession();
- ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
- byte[] hits = new byte[primes.Length];
+ ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL };
+ byte[] hits = new byte[exponents.Length];
- ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
- session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
+ MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(4_000_169UL);
+ ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisor);
+ session.CheckDivisor(4_000_169UL, divisor, cycle, exponents, hits);
- hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
+ byte[] expected = new byte[exponents.Length];
+ for (int i = 0; i < exponents.Length; i++)
+ {
+ expected[i] = exponents[i].Pow2ModWindowedCpu(divisor.Modulus) == 1UL ? (byte)1 : (byte)0;
+ }
+
+ hits.Should().Equal(expected);
}
}
diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
index 7269b31..fed7e76 100644
--- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
+++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
@@ -633,17 +633,23 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit)
{
- if (prime <= 1UL)
- {
- return 0UL;
- }
+ // The by-divisor scanner never emits exponents below the first odd prime, so the defensive guard stays commented out to
+ // avoid branching on production workloads.
+ // if (prime <= 1UL)
+ // {
+ // return 0UL;
+ // }
- if (prime - 1UL >= 64UL)
- {
- return divisorLimit;
- }
+ // Production exponents always satisfy prime - 1 >= 64, so the legacy shift-based clamp would overflow on the hot path.
+ // Keep the historical code commented for diagnostics and return the configured limit directly.
+ // if (prime - 1UL >= 64UL)
+ // {
+ // return divisorLimit;
+ // }
+ //
+ // return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
- return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
+ return divisorLimit;
}
private sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -652,7 +658,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
private bool _disposed;
private ulong[]? _primeDeltas;
private ulong[]? _cycleRemainders;
- private ulong[]? _montgomeryDeltas;
+ private ulong[]? _residues;
private int _bufferCapacity;
private MersenneNumberDivisorResidueGpuEvaluator? _gpuResidueEvaluator;
@@ -674,10 +680,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
int length = primes.Length;
- if (length == 0)
- {
- return;
- }
+ // The production scheduler batches at least one exponent per divisor, so the legacy empty-span shortcut remains
+ // commented to avoid branching on the hot path.
+ // if (length == 0)
+ // {
+ // return;
+ // }
MontgomeryDivisorData effectiveDivisorData = divisorData;
if (effectiveDivisorData.Modulus != divisor)
@@ -765,13 +773,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan);
EnsureGpuEvaluator();
- Span residueSpan = _montgomeryDeltas!.AsSpan(0, length);
+ Span residueSpan = _residues!.AsSpan(0, length);
_gpuResidueEvaluator!.ComputeResidues(primes, divisorData, residueSpan);
- ulong montgomeryOne = divisorData.MontgomeryOne;
for (int i = 0; i < length; i++)
{
- hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == montgomeryOne ? (byte)1 : (byte)0;
+ hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == 1UL ? (byte)1 : (byte)0;
}
}
@@ -801,7 +808,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
return;
}
- if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _montgomeryDeltas is not null)
+ // Sessions release pooled buffers between scans; a single null check covers every array because buffers are rented and returned together.
+ if (requiredLength <= _bufferCapacity && _primeDeltas is not null)
{
return;
}
@@ -809,6 +817,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
int newCapacity = _bufferCapacity;
if (newCapacity == 0)
{
+ // Fresh sessions start without buffers; seed the first rent with a single-slot array.
newCapacity = 1;
}
@@ -821,51 +830,42 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
_primeDeltas = ArrayPool.Shared.Rent(newCapacity);
_cycleRemainders = ArrayPool.Shared.Rent(newCapacity);
- _montgomeryDeltas = ArrayPool.Shared.Rent(newCapacity);
+ _residues = ArrayPool.Shared.Rent(newCapacity);
_bufferCapacity = newCapacity;
}
private void ReturnBuffers()
{
- if (_primeDeltas is not null)
+ // Buffers are lazily allocated per session; the _primeDeltas guard applies to every pooled array because sessions rent and release them together.
+ if (_primeDeltas is null)
{
- ArrayPool.Shared.Return(_primeDeltas, clearArray: false);
- _primeDeltas = null;
- }
-
- if (_cycleRemainders is not null)
- {
- ArrayPool.Shared.Return(_cycleRemainders, clearArray: false);
- _cycleRemainders = null;
- }
-
- if (_montgomeryDeltas is not null)
- {
- ArrayPool.Shared.Return(_montgomeryDeltas, clearArray: false);
- _montgomeryDeltas = null;
+ return;
}
+ ArrayPool.Shared.Return(_primeDeltas, clearArray: false);
+ ArrayPool.Shared.Return(_cycleRemainders!, clearArray: false);
+ ArrayPool.Shared.Return(_residues!, clearArray: false);
+ _primeDeltas = null;
+ _cycleRemainders = null;
+ _residues = null;
_bufferCapacity = 0;
}
private static void ComputePrimeDeltas(ReadOnlySpan primes, Span deltas)
{
- if (primes.Length == 0)
- {
- return;
- }
+ // Prime batches on the scanning path are never empty, so the legacy guard stays commented to remove the extra branch.
+ // if (primes.Length == 0)
+ // {
+ // return;
+ // }
ulong previous = primes[0];
deltas[0] = previous;
+ // Production batches feed strictly increasing primes; callers outside the scanning path must uphold the contract.
for (int i = 1; i < primes.Length; i++)
{
ulong current = primes[i];
- if (current <= previous)
- {
- throw new ArgumentOutOfRangeException(nameof(primes), "Primes must be strictly increasing.");
- }
-
deltas[i] = current - previous;
previous = current;
}
diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs
index c6480ee..bcacd5c 100644
--- a/PerfectNumbers.Core/DivisorCycleCache.cs
+++ b/PerfectNumbers.Core/DivisorCycleCache.cs
@@ -86,11 +86,10 @@ public sealed class DivisorCycleCache
}
cycles[i] = cached;
+ continue;
}
- else
- {
- missingBuffer[missingCount++] = i;
- }
+
+ missingBuffer[missingCount++] = i;
}
if (missingCount != 0)
@@ -104,6 +103,7 @@ public sealed class DivisorCycleCache
{
ComputeCyclesCpu(divisors, cycles, missingIndices);
}
+
}
if (rentedMissing is not null)
diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
index 21769b7..9d10387 100644
--- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
@@ -35,7 +35,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
_kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(CheckKernel));
private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) =>
- _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputeMontgomeryExponentKernel));
+ _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputeExponentKernel));
private Action, ArrayView, byte> GetRemainderDeltaKernel(Accelerator accelerator) =>
_remainderDeltaKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, byte>(ComputeRemainderDeltasKernel));
@@ -530,46 +530,60 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
{
- if (maxPrime <= 1UL)
- {
- return 0UL;
- }
+ // The production configuration seeds the scan with exponents above the first odd prime, so the defensive guard remains
+ // commented to document the invariant without branching.
+ // if (maxPrime <= 1UL)
+ // {
+ // return 0UL;
+ // }
- if (maxPrime - 1UL >= 64UL)
- {
- return ulong.MaxValue;
- }
+ // Real workloads push maxPrime well past the 64-bit shift boundary, so keep the historical implementation commented out
+ // and return ulong.MaxValue directly.
+ // if (maxPrime - 1UL >= 64UL)
+ // {
+ // return ulong.MaxValue;
+ // }
+ //
+ // return (1UL << (int)(maxPrime - 1UL)) - 1UL;
- return (1UL << (int)(maxPrime - 1UL)) - 1UL;
+ return ulong.MaxValue;
}
private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit)
{
- if (prime <= 1UL)
- {
- return 0UL;
- }
+ // The by-divisor search emits candidate exponents starting at the first odd prime, so the defensive guard remains
+ // commented to document that invariant.
+ // if (prime <= 1UL)
+ // {
+ // return 0UL;
+ // }
- if (prime - 1UL >= 64UL)
- {
- return divisorLimit;
- }
+ // Production exponents are always large enough that the shift-based clamp would overflow, so skip the legacy branch and
+ // return the configured divisor limit directly.
+ // if (prime - 1UL >= 64UL)
+ // {
+ // return divisorLimit;
+ // }
+ //
+ // return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
- return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
+ return divisorLimit;
}
private static void CheckKernel(Index1D index, ArrayView divisors, ArrayView exponents, ArrayView hits)
{
MontgomeryDivisorData divisor = divisors[index];
ulong modulus = divisor.Modulus;
- if (modulus <= 1UL || (modulus & 1UL) == 0UL)
- {
- hits[index] = 0;
- return;
- }
+ // The divisor staging path filters out non-prime and even moduli before they reach the device, so this guard
+ // would never fire for EvenPerfectBitScanner workloads.
+ // if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+ // {
+ // hits[index] = 0;
+ // return;
+ // }
ulong exponent = exponents[index];
- hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
+ hits[index] = exponent.Pow2ModWindowedGpu(modulus) == 1UL ? (byte)1 : (byte)0;
}
public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -609,16 +623,20 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
}
int length = primes.Length;
- if (length == 0)
- {
- return;
- }
+ // Production batches always contain at least one exponent, so the empty-span guard remains commented out to remove
+ // the redundant branch.
+ // if (length == 0)
+ // {
+ // return;
+ // }
- if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL)
- {
- hits.Clear();
- return;
- }
+ // The divisor staging path filters non-prime and even moduli before kernels launch, so keep the defensive
+ // modulus guard commented out to avoid branching in the hot loop.
+ // if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL)
+ // {
+ // hits.Clear();
+ // return;
+ // }
ArrayView1D exponentsView = _exponentsBuffer.View;
ArrayView1D resultsView = _resultsBuffer.View;
@@ -777,33 +795,34 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
mask[globalIndex] = 1;
}
- private static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results)
+ private static void ComputeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results)
{
ulong modulus = divisor.Modulus;
- if (modulus <= 1UL || (modulus & 1UL) == 0UL)
- {
- results[index] = 0UL;
- return;
- }
+ // Kernel inputs are constrained to odd prime moduli, so this defensive branch would never execute on the configured
+ // scanning path.
+ // if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+ // {
+ // results[index] = 0UL;
+ // return;
+ // }
ulong exponent = exponents[index];
- results[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true);
+ results[index] = exponent.Pow2ModWindowedGpu(modulus);
}
public ulong DivisorLimit
{
get
{
- lock (_sync)
- {
- if (!_isConfigured)
- {
- throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
- }
+ // The by-divisor scanner always configures the tester before exposing the limit, so the guard remains
+ // commented out to document the invariant without forcing a runtime branch.
+ // if (!_isConfigured)
+ // {
+ // throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
+ // }
- return _divisorLimit;
- }
+ return _divisorLimit;
}
}
@@ -825,10 +844,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
internal BatchResources(MersenneNumberDivisorByDivisorGpuTester owner, Accelerator accelerator, int capacity)
{
int cycleCapacity = DivisorCycleCache.Shared.PreferredBatchSize;
- if (cycleCapacity <= 0)
- {
- cycleCapacity = 1;
- }
+ // The shared cache always reports a positive preferred size on production workloads, so the fallback remains
+ // commented out to avoid branching during resource initialization.
+ // if (cycleCapacity <= 0)
+ // {
+ // cycleCapacity = 1;
+ // }
int actualCapacity = Math.Max(capacity, cycleCapacity);
diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs
index 7c944cc..5c290a2 100644
--- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs
+++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs
@@ -104,21 +104,25 @@ internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable
private static void BuildKernel(Index1D index, ArrayView moduli, ArrayView destination)
{
int idx = index;
- if (idx >= moduli.Length)
- {
- return;
- }
+ // The builder launches kernels with an extent equal to the modulus count, so the bounds guard remains commented
+ // out to avoid redundant branching on every thread.
+ // if (idx >= moduli.Length)
+ // {
+ // return;
+ // }
ulong modulus = moduli[idx];
- if (modulus <= 1UL || (modulus & 1UL) == 0UL)
- {
- destination[idx] = default;
- return;
- }
+ // All staged divisors are odd primes on the scanner path, so keep the defensive modulus guard commented out to
+ // remove the branch from the kernel.
+ // if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+ // {
+ // destination[idx] = default;
+ // return;
+ // }
ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue64(modulus);
- ulong montgomeryTwo = montgomeryOne.MulMod64GpuCompatible(2UL, modulus);
+ ulong montgomeryTwo = montgomeryOne.MulMod64(2UL, modulus);
ulong montgomeryTwoSquared = MontgomeryMultiplyGpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
destination[idx] = new MontgomeryDivisorData(
diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs
index 8fdd163..f6c7d9f 100644
--- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs
+++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs
@@ -103,12 +103,23 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable
{
int idx = index;
int length = (int)exponents.Length;
- if (idx >= length)
- {
- return;
- }
+ // Batched launches match the buffer length on production workloads, so the bounds guard stays commented out to remove
+ // the redundant comparison.
+ // if (idx >= length)
+ // {
+ // return;
+ // }
ulong exponent = exponents[idx];
- results[idx] = exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: true);
+ ulong modulus = divisorData.Modulus;
+ // Residue batches only include valid odd moduli, so the defensive check stays commented to avoid re-validating
+ // the invariant on every iteration.
+ // if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+ // {
+ // results[idx] = 0UL;
+ // return;
+ // }
+
+ results[idx] = exponent.Pow2ModWindowedGpu(modulus);
}
}
diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs
index b403218..2e066c1 100644
--- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs
+++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs
@@ -805,13 +805,15 @@ internal static partial class PrimeOrderGpuHeuristics
for (int i = 0; i < factorCount; i++)
{
ulong factor = factors[i];
- if (factor <= 1UL)
- {
- continue;
- }
+ // The factorization routine emits strictly greater-than-one factors for phi on this path, so the
+ // defensive continue would never execute.
+ // if (factor <= 1UL)
+ // {
+ // continue;
+ // }
ulong reduced = phi / factor;
- if (reduced.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL)
+ if (reduced.Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
{
return false;
}
@@ -827,7 +829,7 @@ internal static partial class PrimeOrderGpuHeuristics
if (residue == 1UL || residue == 7UL)
{
ulong half = phi >> 1;
- if (Pow2EqualsOneKernel(half, divisor))
+ if ((half).Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
{
order = half;
}
@@ -848,10 +850,11 @@ internal static partial class PrimeOrderGpuHeuristics
{
ulong primeFactor = factors[i];
int exponent = exponents[i];
- if (primeFactor <= 1UL)
- {
- continue;
- }
+ // Factorization never yields unity on the scanning path, so skip the redundant guard.
+ // if (primeFactor <= 1UL)
+ // {
+ // continue;
+ // }
for (int iteration = 0; iteration < exponent; iteration++)
{
@@ -861,7 +864,7 @@ internal static partial class PrimeOrderGpuHeuristics
}
ulong reduced = order / primeFactor;
- if (Pow2EqualsOneKernel(reduced, divisor))
+ if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
{
order = reduced;
continue;
@@ -887,12 +890,14 @@ internal static partial class PrimeOrderGpuHeuristics
ArrayView1D compositeStack,
ArrayView1D statusOut)
{
- if (order == 0UL)
- {
- return false;
- }
+ // Order values handed to the kernel are strictly positive on the production path, so the zero guard stays commented
+ // out to remove the redundant branch.
+ // if (order == 0UL)
+ // {
+ // return false;
+ // }
- if (!Pow2EqualsOneKernel(order, divisor))
+ if ((order).Pow2ModWindowedGpu(divisor.Modulus) != 1UL)
{
return false;
}
@@ -922,10 +927,11 @@ internal static partial class PrimeOrderGpuHeuristics
{
ulong primeFactor = factors[i];
int exponent = exponents[i];
- if (primeFactor <= 1UL)
- {
- continue;
- }
+ // Factorization never yields one on the by-divisor path, so the defensive continue stays commented out.
+ // if (primeFactor <= 1UL)
+ // {
+ // continue;
+ // }
ulong reduced = order;
for (int iteration = 0; iteration < exponent; iteration++)
@@ -936,7 +942,7 @@ internal static partial class PrimeOrderGpuHeuristics
}
reduced /= primeFactor;
- if (Pow2EqualsOneKernel(reduced, divisor))
+ if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
{
return false;
}
@@ -998,28 +1004,32 @@ internal static partial class PrimeOrderGpuHeuristics
long candidateCapacity = candidates.Length;
int candidateLimit = candidateCapacity < HeuristicCandidateLimit ? (int)candidateCapacity : HeuristicCandidateLimit;
int candidateCount = BuildCandidatesKernel(order, workFactors, workExponents, factorCount, candidates, stackIndex, stackExponent, stackProduct, candidateLimit);
- if (candidateCount == 0)
- {
- return false;
- }
+ // Candidate enumeration always yields at least one order on the configured workloads, so the zero-count guard remains
+ // commented out to eliminate the redundant branch.
+ // if (candidateCount == 0)
+ // {
+ // return false;
+ // }
SortCandidatesKernel(prime, previousOrder, hasPreviousOrder != 0, candidates, candidateCount);
int powBudget = maxPowChecks <= 0 ? candidateCount : maxPowChecks;
- if (powBudget <= 0)
- {
- powBudget = candidateCount;
- }
+ // The heuristics always reserve at least one check, so the zero-budget fallback remains commented out.
+ // if (powBudget <= 0)
+ // {
+ // powBudget = candidateCount;
+ // }
int powUsed = 0;
for (int i = 0; i < candidateCount && powUsed < powBudget; i++)
{
ulong candidate = candidates[i];
- if (candidate <= 1UL)
- {
- continue;
- }
+ // Candidate generation skips 1 by construction, so the defensive continue stays commented out.
+ // if (candidate <= 1UL)
+ // {
+ // continue;
+ // }
if (powUsed >= powBudget)
{
@@ -1027,7 +1037,7 @@ internal static partial class PrimeOrderGpuHeuristics
}
powUsed++;
- if (!Pow2EqualsOneKernel(candidate, divisor))
+ if ((candidate).Pow2ModWindowedGpu(divisor.Modulus) != 1UL)
{
continue;
}
@@ -1307,10 +1317,11 @@ internal static partial class PrimeOrderGpuHeuristics
{
ulong primeFactor = factors[i];
int exponent = exponents[i];
- if (primeFactor <= 1UL)
- {
- continue;
- }
+ // Factorization never yields one for candidate orders, so keep the guard commented out to avoid re-checking inputs.
+ // if (primeFactor <= 1UL)
+ // {
+ // continue;
+ // }
ulong reduced = candidate;
for (int iteration = 0; iteration < exponent; iteration++)
@@ -1327,7 +1338,7 @@ internal static partial class PrimeOrderGpuHeuristics
}
powUsed++;
- if (Pow2EqualsOneKernel(reduced, divisor))
+ if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
{
return false;
}
@@ -1349,17 +1360,21 @@ internal static partial class PrimeOrderGpuHeuristics
ArrayView1D compositeStack,
ArrayView1D statusOut)
{
- if (initial <= 1UL)
- {
- return true;
- }
+ // Composite candidates sent here are always greater than one during production scans, so the trivial-return branch stays
+ // commented out.
+ // if (initial <= 1UL)
+ // {
+ // return true;
+ // }
int stackCapacity = (int)compositeStack.Length;
- if (stackCapacity <= 0)
- {
- statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow;
- return false;
- }
+ // Pollard-Rho always reserves a positive stack on production workloads, so the overflow guard remains commented out
+ // to avoid branching here.
+ // if (stackCapacity <= 0)
+ // {
+ // statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow;
+ // return false;
+ // }
int stackTop = 0;
compositeStack[stackTop] = initial;
@@ -1369,10 +1384,12 @@ internal static partial class PrimeOrderGpuHeuristics
{
stackTop--;
ulong composite = compositeStack[stackTop];
- if (composite <= 1UL)
- {
- continue;
- }
+ // Composite values enqueued for factoring are always greater than one on production workloads, so the guard stays
+ // commented out to avoid extra branching.
+ // if (composite <= 1UL)
+ // {
+ // continue;
+ // }
if (!PeelSmallPrimesKernel(
composite,
@@ -1580,9 +1597,7 @@ internal static partial class PrimeOrderGpuHeuristics
bLocal >>= BitOperations.TrailingZeroCount(bLocal);
if (aLocal > bLocal)
{
- ulong temp = aLocal;
- aLocal = bLocal;
- bLocal = temp;
+ (aLocal, bLocal) = (bLocal, aLocal);
}
bLocal -= aLocal;
@@ -1593,11 +1608,6 @@ internal static partial class PrimeOrderGpuHeuristics
}
}
- private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
- {
- return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL;
- }
-
private static ulong CalculateByDoublingKernel(ulong prime)
{
ulong order = 1UL;
@@ -1619,10 +1629,11 @@ internal static partial class PrimeOrderGpuHeuristics
private static GpuPow2ModStatus TryPow2ModBatchInternal(ReadOnlySpan exponents, UInt128 prime, Span remainders)
{
- if (exponents.Length == 0)
- {
- return GpuPow2ModStatus.Success;
- }
+ // Residue batches are never empty on production workloads, so skip the zero-length success shortcut.
+ // if (exponents.Length == 0)
+ // {
+ // return GpuPow2ModStatus.Success;
+ // }
if (remainders.Length < exponents.Length)
{
@@ -1698,19 +1709,22 @@ internal static partial class PrimeOrderGpuHeuristics
private static void ComputePow2ModCpu(ReadOnlySpan exponents, ulong prime, in MontgomeryDivisorData divisorData, Span results)
{
int length = exponents.Length;
+ ulong modulus = divisorData.Modulus;
+ _ = prime; // modulus and divisorData both describe the same q; retain the parameter to preserve the public signature.
for (int i = 0; i < length; i++)
{
- results[i] = Pow2ModCpu(exponents[i], prime, divisorData);
+ results[i] = exponents[i].Pow2ModWindowedCpu(modulus);
}
}
private static bool TryComputeOnGpuWide(ReadOnlySpan exponents, UInt128 prime, Span results)
{
int length = exponents.Length;
- if (length == 0)
- {
- return true;
- }
+ // Production launches always carry work, so keep the empty-batch shortcut commented out.
+ // if (length == 0)
+ // {
+ // return true;
+ // }
GpuUInt128[]? rentedExponents = null;
GpuUInt128[]? rentedResults = null;
@@ -1764,6 +1778,8 @@ internal static partial class PrimeOrderGpuHeuristics
}
finally
{
+ // Small batches use the stack-allocated fast path, leaving the rented arrays null. Retain the guards so pooled
+ // buffers are only returned when they were actually leased.
if (rentedExponents is not null)
{
ArrayPool.Shared.Return(rentedExponents, clearArray: false);
@@ -1788,16 +1804,6 @@ internal static partial class PrimeOrderGpuHeuristics
}
}
- private static ulong Pow2ModCpu(ulong exponent, ulong modulus, in MontgomeryDivisorData divisorData)
- {
- if (modulus <= 1UL)
- {
- return 0UL;
- }
-
- return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false);
- }
-
private static Action, MontgomeryDivisorData, ArrayView1D> GetPow2ModKernel(Accelerator accelerator)
{
return Pow2ModKernelCache.GetOrAdd(accelerator, static accel =>
@@ -1810,58 +1816,17 @@ internal static partial class PrimeOrderGpuHeuristics
private static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders)
{
- ulong exponent = exponents[index];
ulong modulus = divisor.Modulus;
- // ulong nPrime = divisor.NPrime;
- // ulong result = divisor.MontgomeryOne;
- // ulong baseValue = divisor.MontgomeryTwo;
- // ulong remainingExponent = exponent;
+ // The by-divisor heuristics call this kernel only with odd prime moduli, so the defensive branch stays commented to
+ // document the invariant without spending a branch in the kernel.
+ // if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+ // {
+ // remainders[index] = 0UL;
+ // return;
+ // }
- // if (remainingExponent == 0UL)
- // {
- // remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime);
- // return;
- // }
-
- // while (remainingExponent != 0UL)
- // {
- // if ((remainingExponent & 1UL) != 0UL)
- // {
- // result = result.MontgomeryMultiply(baseValue, modulus, nPrime);
- // }
-
- // remainingExponent >>= 1;
- // if (remainingExponent == 0UL)
- // {
- // break;
- // }
-
- // // Reusing baseValue to hold base^2 so we avoid allocating another register for the squared term.
- // baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime);
- // }
-
- // remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime);
-
- // Plain modular exponentiation without Montgomery reduction kept for reference.
- ulong plainResult = 1UL;// % modulus;
- ulong plainBase = 2UL;// % modulus;
- while (exponent != 0UL)
- {
- if ((exponent & 1UL) != 0UL)
- {
- plainResult = (plainResult * plainBase) % modulus;
- }
-
- exponent >>= 1;
- if (exponent == 0UL)
- {
- break;
- }
-
- plainBase = (plainBase * plainBase) % modulus;
- }
-
- remainders[index] = plainResult;
+ ulong exponent = exponents[index];
+ remainders[index] = exponent.Pow2ModWindowedGpu(modulus);
}
private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator)
diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs
index 55d6e81..819c8de 100644
--- a/PerfectNumbers.Core/MersenneDivisorCycles.cs
+++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs
@@ -74,23 +74,27 @@ public class MersenneDivisorCycles
return true;
}
- if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL)
- {
- return false;
- }
+ // Exponent and divisor validation happens during scheduling, so the defensive guard remains commented out to keep the
+ // hot path branch-free.
+ // if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL)
+ // {
+ // return false;
+ // }
- if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
- {
- var small = Shared._smallCycles;
- if (small is not null)
- {
- ulong cached = small[(int)divisor];
- if (cached != 0UL)
- {
- return cached == exponent;
- }
- }
- }
+ // Production scans never revisit divisors in the small-cycle snapshot, so skip the legacy lookup here. Factoring and
+ // diagnostics can query TryMatchSmallCycleForFactoring when they need those cached values.
+ // if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
+ // {
+ // var small = Shared._smallCycles;
+ // if (small is not null)
+ // {
+ // ulong cached = small[(int)divisor];
+ // if (cached != 0UL)
+ // {
+ // return cached == exponent;
+ // }
+ // }
+ // }
if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong heuristicCycle) && heuristicCycle != 0UL)
{
@@ -98,6 +102,7 @@ public class MersenneDivisorCycles
}
Dictionary? cache = factorCache;
+ // Each by-divisor scan session starts with a null cache and fills it lazily, so keep the allocation guard in place.
if (cache is null)
{
cache = new Dictionary(8);
@@ -109,7 +114,7 @@ public class MersenneDivisorCycles
return cycleLength == exponent;
}
- return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
+ return exponent.Pow2ModWindowedCpu(divisorData.Modulus) == 1UL;
}
public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
@@ -119,35 +124,70 @@ public class MersenneDivisorCycles
return false;
}
- if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
+ // Candidate divisors on the production path sit well above the precomputed snapshot, so skip the small-cycle lookup
+ // here and leave factoring scenarios to call TryMatchSmallCycleForFactoring when necessary.
+ // if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
+ // {
+ // MersenneDivisorCycles shared = Shared;
+ // ulong[]? small = shared._smallCycles;
+ // if (small is not null)
+ // {
+ // ulong cached = small[(int)divisor];
+ // if (cached != 0UL)
+ // {
+ // return cached == exponent;
+ // }
+ // }
+ // }
+
+ return exponent.Pow2ModWindowedCpu(divisorData.Modulus) == 1UL;
+ }
+
+ public static bool TryMatchSmallCycleForFactoring(ulong divisor, ulong exponent, out ulong cycleLength)
+ {
+ if (divisor > PerfectNumberConstants.MaxQForDivisorCycles)
{
- MersenneDivisorCycles shared = Shared;
- ulong[]? small = shared._smallCycles;
- if (small is not null)
- {
- ulong cached = small[(int)divisor];
- if (cached != 0UL)
- {
- return cached == exponent;
- }
- }
+ // Factoring helpers probe arbitrary divisors, so anything above the snapshot limit legitimately misses here.
+ cycleLength = 0UL;
+ return false;
}
- return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
+ ulong[]? small = Shared._smallCycles;
+ if (small is null)
+ {
+ // Snapshot recycling during diagnostics can temporarily clear the table; surface the miss so callers recompute.
+ cycleLength = 0UL;
+ return false;
+ }
+
+ ulong cached = small[(int)divisor];
+ if (cached == 0UL)
+ {
+ // Only a subset of small divisors carry cached cycles; signal the miss so factoring falls back to the GPU helper.
+ cycleLength = 0UL;
+ return false;
+ }
+
+ cycleLength = cached;
+ return cached == exponent;
}
private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent)
{
- if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL)
- {
- return false;
- }
+ // Production workloads pre-filter invalid inputs, so skip the defensive guards here. Factoring and diagnostics should
+ // validate inputs before calling into the candidate path.
+ // if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL)
+ // {
+ // return false;
+ // }
UInt128 step = (UInt128)exponent << 1;
- if (step == UInt128.Zero)
- {
- return false;
- }
+ // The overflow bail-out is handled by factoring-specific call sites; keep the historical branch commented to document the
+ // legacy behavior.
+ // if (step == UInt128.Zero)
+ // {
+ // return false;
+ // }
UInt128 adjusted = (UInt128)divisor - UInt128.One;
return adjusted % step == UInt128.Zero;
diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs
index fca8662..cf45857 100644
--- a/PerfectNumbers.Core/ULongExtensions.cs
+++ b/PerfectNumbers.Core/ULongExtensions.cs
@@ -12,13 +12,17 @@ public static class ULongExtensions
private const int Pow2WindowSize = 8;
private const ulong Pow2WindowFallbackThreshold = 32UL;
private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1);
+ private const int Pow2WindowSizeGpu = 5;
+ private const int Pow2WindowOddPowerCountGpu = 1 << (Pow2WindowSizeGpu - 1);
public static ulong CalculateOrder(this ulong q)
{
- if (q <= 2UL)
- {
- return 0UL;
- }
+ // Order computations only run on admissible divisors, so the small-prime guard remains commented out to keep the
+ // path branch-free.
+ // if (q <= 2UL)
+ // {
+ // return 0UL;
+ // }
ulong order = q - 1UL, prime, temp;
uint[] smallPrimes = PrimesGenerator.SmallPrimes;
@@ -77,7 +81,7 @@ public static class ULongExtensions
public static ulong MultiplyShiftRightShiftFirst(ulong value, ulong multiplier, int shift)
{
ulong high = value >> shift;
- ulong mask = (1UL << shift) - 1UL;
+ ulong mask = (1UL << shift) - 1UL;
ulong low = value & mask;
UInt128 highContribution = (UInt128)high * multiplier;
@@ -288,26 +292,220 @@ public static class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- // TODO: Replace this fallback with the UInt128 Montgomery helper measured fastest in
- // MulMod64Benchmarks so CPU callers stop paying for triple modulo operations.
- public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus);
-
- [MethodImpl(MethodImplOptions.AggressiveInlining)]
- public static ulong MulMod64GpuCompatible(this ulong a, ulong b, ulong modulus)
+ public static ulong MulMod64(this ulong a, ulong b, ulong modulus)
{
- // TODO: Remove this GPU-compatible shim from production once callers migrate to MulMod64,
- // which the benchmarks show is roughly 6-7× faster on dense 64-bit inputs.
- GpuUInt128 state = new(a % modulus);
- return state.MulMod(b, modulus);
+ // Scanner inputs always provide odd prime moduli, so the legacy modulus guard stays commented out to avoid extra
+ // branching.
+ // if (modulus <= 1UL)
+ // {
+ // return 0UL;
+ // }
+
+ ulong reducedA = a >= modulus ? a % modulus : a;
+ ulong reducedB = b >= modulus ? b % modulus : b;
+ GpuUInt128 state = new(reducedA);
+ return state.MulMod(reducedB, modulus);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public static ulong MulMod64GpuCompatibleDeferred(this ulong a, ulong b, ulong modulus)
+ public static ulong Pow2ModWindowedCpu(this ulong exponent, ulong modulus)
{
- // TODO: Move this deferred helper to the benchmark suite; the baseline MulMod64 avoids the
- // 5-40× slowdown seen across real-world operand distributions.
- GpuUInt128 state = new(a);
- return state.MulModWithNativeModulo(b, modulus);
+ // The by-divisor pipeline only schedules odd prime moduli, so the defensive check below
+ // would never trigger on that path. Leaving it commented documents the invariant without
+ // forcing a branch in the hot loop.
+ // if (modulus <= 1UL)
+ // {
+ // return 0UL;
+ // }
+
+ // Exponents flow from the strictly positive candidate list on the by-divisor scanner, so
+ // this guard would be redundant there.
+ // if (exponent == 0UL)
+ // {
+ // return 1UL % modulus;
+ // }
+
+ // Even though the by-divisor pipeline only feeds very large exponents, the legacy diagnostic callers
+ // still hit this small-exponent fast path, so keep the fallback in place.
+ if (exponent <= Pow2WindowFallbackThreshold)
+ {
+ return Pow2ModBinaryFallback(exponent, modulus);
+ }
+
+ int bitLength = GetPortableBitLength(exponent);
+ int windowSize = GetWindowSize(bitLength);
+ int oddPowerCount = 1 << (windowSize - 1);
+
+ Span oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount
+ ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount]
+ : new ulong[oddPowerCount];
+ oddPowerStorage = oddPowerStorage[..oddPowerCount];
+ InitializePlainOddPowers(modulus, oddPowerStorage);
+
+ // With modulus guaranteed to exceed one on the scanning path, the modulo operation would
+ // always yield one, so keep the literal to highlight the intended base value.
+ ulong result = 1UL;
+ int index = bitLength - 1;
+
+ while (index >= 0)
+ {
+ if (((exponent >> index) & 1UL) == 0UL)
+ {
+ result = result.MulMod64(result, modulus);
+ index--;
+ continue;
+ }
+
+ int windowStart = index - windowSize + 1;
+ if (windowStart < 0)
+ {
+ windowStart = 0;
+ }
+
+ while (((exponent >> windowStart) & 1UL) == 0UL)
+ {
+ windowStart++;
+ }
+
+ int windowLength = index - windowStart + 1;
+ for (int square = 0; square < windowLength; square++)
+ {
+ result = result.MulMod64(result, modulus);
+ }
+
+ ulong windowValue = (exponent >> windowStart) & ((1UL << windowLength) - 1UL);
+ int tableIndex = (int)((windowValue - 1UL) >> 1);
+ result = result.MulMod64(oddPowerStorage[tableIndex], modulus);
+
+ index = windowStart - 1;
+ }
+
+ return result;
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static ulong Pow2ModBinaryFallback(ulong exponent, ulong modulus)
+ {
+ // The scanning path keeps modulus strictly greater than one, but other utility callers still
+ // rely on this fallback for tiny exponents, so leave the modulo in place for safety.
+ ulong remainingExponent = exponent;
+ // Production moduli are always odd primes greater than two, so normalize the seed values directly.
+ ulong result = 1UL;
+ ulong baseValue = 2UL;
+
+ while (remainingExponent != 0UL)
+ {
+ if ((remainingExponent & 1UL) != 0UL)
+ {
+ result = result.MulMod64(baseValue, modulus);
+ }
+
+ remainingExponent >>= 1;
+ if (remainingExponent == 0UL)
+ {
+ break;
+ }
+
+ baseValue = baseValue.MulMod64(baseValue, modulus);
+ }
+
+ return result;
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus)
+ {
+ // The GPU pipeline feeds only odd prime moduli, so the removed guard would never fire for
+ // EvenPerfectBitScanner workloads.
+ // if (modulus <= 1UL)
+ // {
+ // return 0UL;
+ // }
+
+ // The GPU pipeline only emits positive exponents on production workloads, so skip the zero guard to avoid an extra branch.
+ // Diagnostic harnesses must normalize their inputs before invoking this helper.
+ // if (exponent == 0UL)
+ // {
+ // return 1UL;
+ // }
+
+ // The GPU kernels mostly process large exponents as well, but the binary fallback still serves
+ // diagnostic harnesses that probe edge cases around the threshold.
+ if (exponent <= Pow2WindowFallbackThreshold)
+ {
+ return Pow2ModBinaryGpuFallback(exponent, modulus);
+ }
+
+ int bitLength = GetPortableBitLength(exponent);
+ int windowSize = GetGpuWindowSize(bitLength);
+ int oddPowerCount = 1 << (windowSize - 1);
+
+ PlainOddPowerTableGpu oddPowers = new PlainOddPowerTableGpu(modulus, oddPowerCount);
+
+ ulong result = 1UL;
+ int index = bitLength - 1;
+
+ while (index >= 0)
+ {
+ if (((exponent >> index) & 1UL) == 0UL)
+ {
+ result = result.MulMod64(result, modulus);
+ index--;
+ continue;
+ }
+
+ int windowStart = index - windowSize + 1;
+ if (windowStart < 0)
+ {
+ windowStart = 0;
+ }
+
+ while (((exponent >> windowStart) & 1UL) == 0UL)
+ {
+ windowStart++;
+ }
+
+ int windowLength = index - windowStart + 1;
+ for (int square = 0; square < windowLength; square++)
+ {
+ result = result.MulMod64(result, modulus);
+ }
+
+ ulong windowValue = (exponent >> windowStart) & ((1UL << windowLength) - 1UL);
+ int tableIndex = (int)((windowValue - 1UL) >> 1);
+ result = result.MulMod64(oddPowers.Get(tableIndex), modulus);
+
+ index = windowStart - 1;
+ }
+
+ return result;
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static ulong Pow2ModBinaryGpuFallback(ulong exponent, ulong modulus)
+ {
+ ulong remainingExponent = exponent;
+ // Production moduli are always odd primes greater than two, so normalize the seed values directly.
+ ulong result = 1UL;
+ ulong baseValue = 2UL;
+
+ while (remainingExponent != 0UL)
+ {
+ if ((remainingExponent & 1UL) != 0UL)
+ {
+ result = result.MulMod64(baseValue, modulus);
+ }
+
+ remainingExponent >>= 1;
+ if (remainingExponent == 0UL)
+ {
+ break;
+ }
+
+ baseValue = baseValue.MulMod64(baseValue, modulus);
+ }
+
+ return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -347,10 +545,16 @@ public static class ULongExtensions
ulong result = divisor.MontgomeryOne;
ulong nPrime = divisor.NPrime;
- Span oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount] : default;
+ Span oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount
+ ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount]
+ : default;
ulong[]? pooledOddPowers = null;
- if (oddPowerCount < PerfectNumberConstants.PooledArrayThreshold)
+ if (oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount)
+ {
+ oddPowerStorage = oddPowerStorage[..oddPowerCount];
+ }
+ else if (oddPowerCount < PerfectNumberConstants.PooledArrayThreshold)
{
oddPowerStorage = new ulong[oddPowerCount];
}
@@ -548,6 +752,187 @@ public static class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static int GetGpuWindowSize(int bitLength)
+ {
+ int window = GetWindowSize(bitLength);
+ return window > Pow2WindowSizeGpu ? Pow2WindowSizeGpu : window;
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static void InitializePlainOddPowers(ulong modulus, Span destination)
+ {
+ int oddPowerCount = destination.Length;
+ // The by-divisor kernels always request at least one odd power, but other helpers reuse this initializer
+ // for smaller workloads, so the early exit stays in place.
+ if (oddPowerCount == 0)
+ {
+ return;
+ }
+
+ ulong baseValue = 2UL % modulus;
+ destination[0] = baseValue;
+ if (oddPowerCount == 1)
+ {
+ return;
+ }
+
+ ulong square = baseValue.MulMod64(baseValue, modulus);
+ for (int i = 1; i < oddPowerCount; i++)
+ {
+ // Reuse baseValue to accumulate successive odd powers and avoid allocating another local.
+ baseValue = baseValue.MulMod64(square, modulus);
+ destination[i] = baseValue;
+ }
+ }
+
+ private struct PlainOddPowerTableGpu
+ {
+ private readonly int _count;
+ private ulong _value0;
+ private ulong _value1;
+ private ulong _value2;
+ private ulong _value3;
+ private ulong _value4;
+ private ulong _value5;
+ private ulong _value6;
+ private ulong _value7;
+ private ulong _value8;
+ private ulong _value9;
+ private ulong _value10;
+ private ulong _value11;
+ private ulong _value12;
+ private ulong _value13;
+ private ulong _value14;
+ private ulong _value15;
+
+ public PlainOddPowerTableGpu(ulong modulus, int requestedCount)
+ {
+ _count = requestedCount > Pow2WindowOddPowerCountGpu ? Pow2WindowOddPowerCountGpu : requestedCount;
+ _value0 = 0UL;
+ _value1 = 0UL;
+ _value2 = 0UL;
+ _value3 = 0UL;
+ _value4 = 0UL;
+ _value5 = 0UL;
+ _value6 = 0UL;
+ _value7 = 0UL;
+ _value8 = 0UL;
+ _value9 = 0UL;
+ _value10 = 0UL;
+ _value11 = 0UL;
+ _value12 = 0UL;
+ _value13 = 0UL;
+ _value14 = 0UL;
+ _value15 = 0UL;
+
+ // Window planning always requests at least one odd power, so keep the zero-count check commented out to avoid
+ // branching here.
+ // if (_count == 0)
+ // {
+ // return;
+ // }
+
+ ulong baseValue = 2UL;
+ _value0 = baseValue;
+ if (_count == 1)
+ {
+ return;
+ }
+
+ ulong square = baseValue.MulMod64(baseValue, modulus);
+ ulong current = baseValue;
+ for (int i = 1; i < _count; i++)
+ {
+ current = current.MulMod64(square, modulus);
+ SetValue(i, current);
+ }
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public ulong Get(int index)
+ {
+ // Windowed lookups never exceed the populated table bounds on production workloads, so leave the check commented
+ // out to remove the branch.
+ // if ((uint)index >= (uint)_count)
+ // {
+ // return 0UL;
+ // }
+
+ return index switch
+ {
+ 0 => _value0,
+ 1 => _value1,
+ 2 => _value2,
+ 3 => _value3,
+ 4 => _value4,
+ 5 => _value5,
+ 6 => _value6,
+ 7 => _value7,
+ 8 => _value8,
+ 9 => _value9,
+ 10 => _value10,
+ 11 => _value11,
+ 12 => _value12,
+ 13 => _value13,
+ 14 => _value14,
+ _ => _value15,
+ };
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private void SetValue(int index, ulong value)
+ {
+ switch (index)
+ {
+ case 1:
+ _value1 = value;
+ break;
+ case 2:
+ _value2 = value;
+ break;
+ case 3:
+ _value3 = value;
+ break;
+ case 4:
+ _value4 = value;
+ break;
+ case 5:
+ _value5 = value;
+ break;
+ case 6:
+ _value6 = value;
+ break;
+ case 7:
+ _value7 = value;
+ break;
+ case 8:
+ _value8 = value;
+ break;
+ case 9:
+ _value9 = value;
+ break;
+ case 10:
+ _value10 = value;
+ break;
+ case 11:
+ _value11 = value;
+ break;
+ case 12:
+ _value12 = value;
+ break;
+ case 13:
+ _value13 = value;
+ break;
+ case 14:
+ _value14 = value;
+ break;
+ case 15:
+ _value15 = value;
+ break;
+ }
+ }
+ }
+
private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span destination)
{
int oddPowerCount = destination.Length;
diff --git a/todo_list.txt b/todo_list.txt
index 3ed48dc..9a7cdab 100644
--- a/todo_list.txt
+++ b/todo_list.txt
@@ -122,7 +122,6 @@ PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:203:
PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:207: // TODO: When divisorCycle is zero (cache miss), compute only that single cycle on the device
PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:214: // TODO: Remove this no-cycle branch and require on-demand cycle computation (without
PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:240: // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the
-PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:262: : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the
PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:283: // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:356: // TODO: Replace this guard with an on-demand CPU/GPU cycle computation path that skips
PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:396: // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here
@@ -266,7 +265,6 @@ todo_list.txt:121:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTeste
todo_list.txt:122:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:207: // TODO: When divisorCycle is zero (cache miss), compute only that single cycle on the device
todo_list.txt:123:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:214: // TODO: Remove this no-cycle branch and require on-demand cycle computation (without
todo_list.txt:124:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:240: // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the
-todo_list.txt:125:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:262: : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the
todo_list.txt:126:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:283: // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
todo_list.txt:127:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:356: // TODO: Replace this guard with an on-demand CPU/GPU cycle computation path that skips
todo_list.txt:128:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:396: // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here
@@ -407,17 +405,6 @@ PerfectNumbers.Core\ModResidueTracker.cs:257: // TODO: Us
PerfectNumbers.Core\ModResidueTracker.cs:263: // TODO: Consult MersenneDivisorCycles.Shared (or similar caches) here so we reuse precomputed cycle lengths
PerfectNumbers.Core\ModResidueTracker.cs:295: // TODO: Route this through the planned UInt128 windowed powmod so the inner MulMod adopts the
PerfectNumbers.Core\ModResidueTracker.cs:314: // TODO: Replace this double-and-add fallback with the UInt128 intrinsic-backed multiplier once the
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:49: set => _useDivisorCycles = value; // TODO: Remove this flag once GPU runs unconditionally fetch cycle lengths so every divisor check benefits from the measured speedups.
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:298: // TODO: Stage MontgomeryDivisorData in a GPU-ready cache keyed by divisor so we can copy
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:372: // TODO: Replace the block lease with a direct single-cycle lookup so the GPU path
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:395: // TODO: Drop this block hand-off once the cache exposes a direct
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:409: // TODO: Replace this `%` with the cycle-cache helper that subtracts precomputed
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:425: // TODO: Trigger an immediate GPU computation for the single missing cycle without storing
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:476: // TODO: Preload ProcessEightBitWindows-friendly divisor metadata for GPU batches so this
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:563: hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; // TODO: Replace Pow2MontgomeryMod with the ProcessEightBitWindows helper once available so GPU by-divisor scans reuse the benchmarked fast pow2 ladder.
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:632: // TODO: Replace this guard with an on-demand cycle computation that runs on the device selected
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:787: // TODO: Replace this `%` with the divisor-cycle remainder helper so we reuse cached
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:858: results[index] = exponent.Pow2MontgomeryModMontgomery(divisor); // TODO: Swap this Montgomery-domain pow2 with the windowed helper slated to replace the slower ladder in Pow2Minus1Mod benchmarks.
PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:43: // TODO: Swap this per-run enumerator creation for the shared PrimeIterator so by-divisor scans reuse the
PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:219: // TODO: Swap Task.Factory.StartNew for the shared low-overhead work queue once the scanner's
PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:240: bool useDivisorCycles = tester.UseDivisorCycles; // TODO: Remove the conditional path once all testers always enable divisor cycles to keep every divisor scan on the faster cached-length track.