Compare commits

...
Author SHA1 Message Date
Marek Lemańczyk 7ea7cf527e Merge pull request #168 from mlemanczyk/codex/remove-redundant-wrapper-methods-in-mersenne-path
Eliminate IDisposable dependencies
2025-10-13 22:16:12 +02:00
Marek Lemańczyk c6008a65b4 Cache GPU remainder configs and fold constant modulos 2025-10-13 22:15:58 +02:00
Marek Lemańczyk 0b43880959 Remove IDisposable usage and using declarations 2025-10-13 21:55:16 +02:00
Marek Lemańczyk 305a489d6c Merge pull request #167 from mlemanczyk/codex/analyze-if-conditions-in-evenperfectbitscanner
Document unreachable guards in GPU by-divisor pipeline
2025-10-13 20:42:48 +02:00
Marek Lemańczyk f2eb6580c4 Comment unreachable configuration guards in by-divisor testers 2025-10-13 20:41:14 +02:00
Marek Lemańczyk b0ba6401ca Remove configuration locks from by-divisor testers 2025-10-13 20:41:07 +02:00
Marek Lemańczyk 6e89700f80 Clarify invariants across GPU divisor helpers 2025-10-13 20:41:00 +02:00
Marek Lemańczyk bccb2f43c9 Document unreachable guards in GPU by-divisor pipeline 2025-10-13 19:56:47 +02:00
Marek Lemańczyk 65e010ded5 Merge pull request #166 from mlemanczyk/codex/locate-montgomery-arithmetic-usage-in-gpu-kernels-ir9oyy
Align prime-order GPU pow2 kernel with plain modular helper
2025-10-13 18:21:54 +02:00
Marek Lemańczyk a255dc718d Restore snapshot missing-cycle guard 2025-10-13 18:21:22 +02:00
Marek Lemańczyk c81579865c Compute missing small divisor cycles on demand 2025-10-13 18:18:57 +02:00
Marek Lemańczyk c8d6143ab7 Simplify by-divisor buffer guards and document prime ordering 2025-10-13 18:12:06 +02:00
Marek Lemańczyk 594a38b9d7 Move by-divisor tests to large-cycle divisors 2025-10-13 18:05:58 +02:00
Marek Lemańczyk c0b6b29837 Document divisor invariants and refresh pow2 tests 2025-10-13 17:42:56 +02:00
Marek Lemańczyk af909e9e2d Simplify prime guards and adopt unified mulmod 2025-10-13 16:56:50 +02:00
Marek Lemańczyk 3f5904c340 Clarify windowed pow2 invariants and tidy Montgomery pooling 2025-10-13 14:30:06 +02:00
Marek Lemańczyk 79695b0c56 Switch GPU pow2mod helpers to plain windowed arithmetic 2025-10-13 14:21:35 +02:00
Marek Lemańczyk 1916196e89 Introduce windowed plain pow2mod helper for CPU checks 2025-10-13 14:21:28 +02:00
Marek Lemańczyk 73e27326a2 Use plain modular math for single-shot CPU checks 2025-10-13 14:21:22 +02:00
Marek Lemańczyk 092947e216 Remove GPU Montgomery usage tracker note 2025-10-13 12:19:01 +02:00
Marek Lemańczyk 716ffcb243 Route prime-order pow2 kernel through plain modular helper 2025-10-13 05:08:34 +02:00
Marek Lemańczyk e10394afa3 Merge pull request #165 from mlemanczyk/revert-164-codex/locate-montgomery-arithmetic-usage-in-gpu-kernels
Revert "Use standard modular math for single-use GPU exponent checks"
2025-10-13 05:00:42 +02:00
Marek Lemańczyk 8b02993bac Revert "Use standard modular math for single-use GPU exponent checks" 2025-10-13 04:59:42 +02:00
Marek Lemańczyk ad6a5c6e61 Merge pull request #164 from mlemanczyk/codex/locate-montgomery-arithmetic-usage-in-gpu-kernels
Use standard modular math for single-use GPU exponent checks
2025-10-13 04:59:36 +02:00
34 changed files with 1332 additions and 856 deletions

No files matched your search

@@ -158,36 +158,16 @@ public class MulMod64Benchmarks
}
/// <summary>
/// Uses the GPU-friendly <see cref="ULongExtensions.MulMod64GpuCompatible(ulong, ulong, ulong)"/> 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 <see cref="ULongExtensions.MulMod64(ulong, ulong, ulong)"/> helper shared by CPU and GPU callers.
/// </summary>
/// <remarks>
/// 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.
/// </remarks>
[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);
}
/// <summary>
/// Uses the deferred native-modulo helper <see cref="ULongExtensions.MulMod64GpuCompatibleDeferred(ulong, ulong, ulong)"/>;
/// excels when operands are tiny (2.01 ns on ZeroOperands) yet remains the slowest choice on dense data at 146–295 ns.
/// </summary>
/// <remarks>
/// 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.
/// </remarks>
[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)]
+2 -1
View File
@@ -540,7 +540,7 @@ public class ProgramTests
public void Main_displays_help_when_requested()
{
var main = typeof(Program).GetMethod("Main", BindingFlags.NonPublic | BindingFlags.Static)!;
using var writer = new StringWriter();
var writer = new StringWriter();
TextWriter original = Console.Out;
Console.SetOut(writer);
@@ -554,6 +554,7 @@ public class ProgramTests
}
string output = writer.ToString();
writer.Dispose();
output.Should().Contain("Usage:");
output.Should().Contain("--mersenne-device=cpu|gpu");
output.Should().Contain("--primes-device=cpu|gpu");
+16 -23
View File
@@ -720,7 +720,21 @@ internal static class Program
if (_useByDivisorMode)
{
RunByDivisorMode(byDivisorCandidates, threadCount, _primeTestLimit);
if (_byDivisorTester is null)
{
throw new InvalidOperationException("By-divisor tester has not been initialized.");
}
MersenneNumberDivisorByDivisorTester.Run(
byDivisorCandidates,
_byDivisorTester,
_byDivisorPreviousResults,
_byDivisorStartPrime,
static () => _lastCompositeP = true,
static () => _lastCompositeP = false,
PrintResult,
threadCount,
_primeTestLimit);
FlushBuffer();
StringBuilderPool.Return(_outputBuilder!);
return;
@@ -861,27 +875,6 @@ internal static class Program
return candidates;
}
private static void RunByDivisorMode(List<ulong> candidates, int threadCount, TimeSpan? primeTestLimit)
{
if (_byDivisorTester is null)
{
throw new InvalidOperationException("By-divisor tester has not been initialized.");
}
// TODO: Inline this helper by invoking MersenneNumberDivisorByDivisorTester.Run directly at call sites so
// the CLI avoids this extra frame during the hot composite-reporting path.
MersenneNumberDivisorByDivisorTester.Run(
candidates,
_byDivisorTester,
_byDivisorPreviousResults,
_byDivisorStartPrime,
static () => _lastCompositeP = true,
static () => _lastCompositeP = false,
PrintResult,
threadCount,
primeTestLimit);
}
private static unsafe void LoadResultsFile(string resultsFileName, Action<ulong, bool, bool> lineProcessorAction)
{
using FileStream readStream = new(resultsFileName, FileMode.Open, FileAccess.Read, FileShare.ReadWrite);
@@ -1030,7 +1023,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=<path> process only p from previous run results (required for --mersenne=bydivisor)");
Console.WriteLine(" --filter-p=<path> process only p from previous run results (required for --mersenne=bydivisor; primes must be strictly increasing)");
Console.WriteLine(" --write-batch-size=<value> 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");
@@ -18,13 +18,14 @@ public class GpuPrimeLimiterTests
var startedSecond = new TaskCompletionSource<bool>();
var sw = new Stopwatch();
using var first = GpuPrimeWorkLimiter.Acquire();
var first = GpuPrimeWorkLimiter.Acquire();
sw.Start();
var secondTask = Task.Run(() =>
{
using var second = GpuPrimeWorkLimiter.Acquire();
var second = GpuPrimeWorkLimiter.Acquire();
startedSecond.TrySetResult(true);
second.Dispose();
});
// Give secondTask a moment to attempt acquire (it should block).
@@ -34,23 +34,19 @@ public class NttMontgomeryTests
// Device execution
GpuContextPool.GpuContextLease context = GpuContextPool.Rent();
try
{
var accelerator = context.Accelerator;
using var buffer = accelerator.Allocate1D<GpuUInt128>(gpu.Length);
buffer.View.CopyFromCPU(ref gpu[0], gpu.Length);
NttGpuMath.SquareDevice(accelerator, buffer.View, modulus, primitiveRoot);
buffer.View.CopyToCPU(ref gpu[0], gpu.Length);
var accelerator = context.Accelerator;
var buffer = accelerator.Allocate1D<GpuUInt128>(gpu.Length);
buffer.View.CopyFromCPU(ref gpu[0], gpu.Length);
NttGpuMath.SquareDevice(accelerator, buffer.View, modulus, primitiveRoot);
buffer.View.CopyToCPU(ref gpu[0], gpu.Length);
for (int i = 0; i < cpu.Length; i++)
{
gpu[i].Should().Be(cpu[i]);
}
}
finally
for (int i = 0; i < cpu.Length; i++)
{
context.Dispose();
gpu[i].Should().Be(cpu[i]);
}
buffer.Dispose();
context.Dispose();
}
[Fact]
@@ -17,13 +17,14 @@ public class GpuWorkLimiterTests
var startedSecond = new TaskCompletionSource<bool>();
var sw = new Stopwatch();
using var first = GpuWorkLimiter.Acquire();
var first = GpuWorkLimiter.Acquire();
sw.Start();
var secondTask = Task.Run(() =>
{
using var second = GpuWorkLimiter.Acquire();
var second = GpuWorkLimiter.Acquire();
startedSecond.TrySetResult(true);
second.Dispose();
});
await Task.Delay(50);
@@ -14,8 +14,8 @@ public class MersenneDivisorCyclesGenerateGpuTests
{
MersenneDivisorCycles.GenerateGpu(path, maxDivisor: 50UL, batchSize: 16, skipCount: 0, nextPosition: 0);
var pairs = new List<(ulong divisor, ulong cycle)>();
using var stream = new FileStream(path, FileMode.Open, FileAccess.Read);
using var reader = new BinaryReader(stream);
var stream = new FileStream(path, FileMode.Open, FileAccess.Read);
var reader = new BinaryReader(stream);
while (stream.Position < stream.Length)
{
ulong d = reader.ReadUInt64();
@@ -41,6 +41,9 @@ public class MersenneDivisorCyclesGenerateGpuTests
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(pair.divisor);
pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLength(pair.divisor, divisorData));
}
reader.Dispose();
stream.Dispose();
}
finally
{
@@ -61,8 +64,8 @@ public class MersenneDivisorCyclesGenerateGpuTests
MersenneDivisorCycles.Generate(path, maxDivisor: 50UL, threads: 4);
var pairs = new List<(ulong divisor, ulong cycle)>();
using var stream = new FileStream(path, FileMode.Open, FileAccess.Read);
using var reader = new BinaryReader(stream);
var stream = new FileStream(path, FileMode.Open, FileAccess.Read);
var reader = new BinaryReader(stream);
while (stream.Position < stream.Length)
{
ulong divisor = reader.ReadUInt64();
@@ -95,6 +98,9 @@ public class MersenneDivisorCyclesGenerateGpuTests
var (_, completeCount) = MersenneDivisorCycles.FindLast(path);
completeCount.Should().Be(expected.Count);
reader.Dispose();
stream.Dispose();
}
finally
{
@@ -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,100 @@ 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];
var session = tester.CreateDivisorSession();
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);
session.Dispose();
}
[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 };
var session = tester.CreateDivisorSession();
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);
session.Dispose();
}
[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 };
var session = tester.CreateDivisorSession();
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);
session.Dispose();
}
}
@@ -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,41 @@ 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];
var session = tester.CreateDivisorSession();
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);
session.Dispose();
}
[Fact]
@@ -180,122 +198,134 @@ 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];
var session = tester.CreateDivisorSession();
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);
session.Dispose();
}
[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 };
var session = tester.CreateDivisorSession();
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<ulong> 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);
session.Dispose();
}
[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 };
var session = tester.CreateDivisorSession();
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<ulong> 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);
session.Dispose();
}
[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);
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);
session.Dispose();
}
[Fact]
@@ -307,16 +337,25 @@ 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];
var session = tester.CreateDivisorSession();
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);
session.Dispose();
}
}
@@ -3,6 +3,7 @@ using System.Buffers;
using System.Collections.Concurrent;
using System.Collections.Generic;
using System.Numerics;
using System.Threading;
using PerfectNumbers.Core;
using PerfectNumbers.Core.Gpu;
@@ -10,14 +11,11 @@ namespace PerfectNumbers.Core.Cpu;
public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester
{
private readonly object _sync = new();
private readonly ConcurrentBag<DivisorScanSession> _sessionPool = new();
private ulong _divisorLimit;
private ulong _lastStatusDivisor;
private int _lastStatusDivisor;
private bool _isConfigured;
private int _batchSize = 1_024;
private ulong _configuredMaxPrime;
private ByDivisorDeltasDevice _deltasDevice = ByDivisorDeltasDevice.Cpu;
private ByDivisorMontgomeryDevice _montgomeryDevice = ByDivisorMontgomeryDevice.Cpu;
@@ -34,79 +32,60 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
public ByDivisorDeltasDevice DeltasDevice
{
get => _deltasDevice;
set
{
lock (_sync)
{
_deltasDevice = value;
}
}
set => _deltasDevice = value;
}
public ByDivisorMontgomeryDevice MontgomeryDevice
{
get => _montgomeryDevice;
set
{
lock (_sync)
{
_montgomeryDevice = value;
}
}
set => _montgomeryDevice = value;
}
public void ConfigureFromMaxPrime(ulong maxPrime)
{
lock (_sync)
{
_configuredMaxPrime = maxPrime;
_divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
_lastStatusDivisor = 0UL;
_isConfigured = true;
}
_divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
_lastStatusDivisor = 0;
_isConfigured = true;
}
public ulong DivisorLimit
{
get
{
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
// The by-divisor driver configures the tester during startup, so exposing the divisor limit never
// re-enters this method before initialization completes.
// if (!_isConfigured)
// {
// throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
// }
return _divisorLimit;
}
return _divisorLimit;
}
}
public ulong GetAllowedMaxDivisor(ulong prime)
{
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
// The CLI wires the tester through ConfigureFromMaxPrime before invoking query helpers, so the guard
// would never trigger in production runs.
// if (!_isConfigured)
// {
// throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
// }
return ComputeAllowedMaxDivisor(prime, _divisorLimit);
}
return ComputeAllowedMaxDivisor(prime, _divisorLimit);
}
public bool IsPrime(ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null)
{
ulong allowedMax;
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
// The tester is configured once during startup, so runtime IsPrime calls always observe a ready
// configuration and the legacy guard can remain commented out.
// if (!_isConfigured)
// {
// throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
// }
allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit);
}
ulong divisorLimit = _divisorLimit;
ulong allowedMax = ComputeAllowedMaxDivisor(prime, divisorLimit);
if (allowedMax < 3UL)
{
@@ -128,10 +107,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
if (processedCount > 0UL)
{
lock (_sync)
{
UpdateStatusUnsafe(processedCount);
}
UpdateStatusUnsafe(processedCount);
}
if (composite)
@@ -146,22 +122,21 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
public void PrepareCandidates(ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues)
{
if (allowedMaxValues.Length < primes.Length)
{
throw new ArgumentException("allowedMaxValues span must be at least as long as primes span.", nameof(allowedMaxValues));
}
// EvenPerfectBitScanner always provides equally sized spans here, so the historical guard never triggers.
// if (allowedMaxValues.Length < primes.Length)
// {
// throw new ArgumentException("allowedMaxValues span must be at least as long as primes span.", nameof(allowedMaxValues));
// }
ulong divisorLimit;
// PrepareCandidates is scheduled only after ConfigureFromMaxPrime succeeds, so the tester state is
// stable here and the defensive guard is no longer needed.
// if (!_isConfigured)
// {
// throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
// }
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
divisorLimit = _divisorLimit;
}
// Configuration is immutable after startup, so reading the cached limit does not require locking.
ulong divisorLimit = _divisorLimit;
for (int index = 0; index < primes.Length; index++)
{
@@ -171,21 +146,20 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession()
{
lock (_sync)
// The session pool is populated only after configuration, so consumers never request a session before
// ConfigureFromMaxPrime completes.
// if (!_isConfigured)
// {
// throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
// }
if (_sessionPool.TryTake(out DivisorScanSession? session))
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
if (_sessionPool.TryTake(out DivisorScanSession? session))
{
session.Reset();
return session;
}
return new DivisorScanSession(this);
session.Reset();
return session;
}
return new DivisorScanSession(this);
}
private void ReturnSession(DivisorScanSession session)
@@ -236,11 +210,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
UInt128 step = (UInt128)prime << 1;
if (step == UInt128.Zero)
{
processedAll = true;
return false;
}
// EvenPerfectBitScanner only supplies primes greater than one here, so the shifted step never becomes zero.
// if (step == UInt128.Zero)
// {
// processedAll = true;
// return false;
// }
UInt128 limit = allowedMax;
UInt128 divisor = step + UInt128.One;
@@ -323,6 +298,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
if (divisorCycle == 0UL)
{
// Keep this guard active to surface mismatches between generated candidates and the snapshot cache.
Console.WriteLine($"Divisor cycle was not calculated for {prime}");
}
}
@@ -369,11 +345,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
GpuUInt128 step128 = new GpuUInt128(prime);
step128.ShiftLeft(1);
if (step128.IsZero)
{
processedAll = true;
return false;
}
// Primes above one yield a non-zero 2p value even after widening to GpuUInt128, so the defensive
// zero guard remains commented out to document the invariant without branching.
// if (step128.IsZero)
// {
// processedAll = true;
// return false;
// }
GpuUInt128 limit128 = new GpuUInt128(allowedMax);
GpuUInt128 divisor128 = step128;
@@ -449,10 +427,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
int chunkCount = remainingCandidates > (ulong)batchCapacity
? batchCapacity
: (int)remainingCandidates;
if (chunkCount <= 0)
{
break;
}
// batchCapacity and remainingCandidates are strictly positive in this path, so the zero-count
// shortcut never executes on production scans.
// if (chunkCount <= 0)
// {
// break;
// }
Span<ulong> candidateSpan = candidateBuffer.AsSpan(0, chunkCount);
Span<byte> maskSpan = maskBuffer.AsSpan(0, chunkCount);
@@ -603,17 +583,34 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
return;
}
ulong interval = PerfectNumberConstants.ConsoleInterval;
if (interval == 0UL)
int interval = PerfectNumberConstants.ConsoleInterval;
// PerfectNumberConstants.ConsoleInterval is a positive constant, so this guard never triggers on EvenPerfectBitScanner.
// if (interval <= 0)
// {
// Volatile.Write(ref _lastStatusDivisor, 0);
// return;
// }
int processedDelta = (int)(processedCount % (uint)interval);
if (processedDelta == 0)
{
_lastStatusDivisor = 0UL;
return;
}
ulong total = _lastStatusDivisor + processedCount;
// TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
// status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly.
_lastStatusDivisor = total % interval;
while (true)
{
int previous = Volatile.Read(ref _lastStatusDivisor);
int total = previous + processedDelta;
if (total >= interval)
{
total -= interval;
}
if (previous == Interlocked.CompareExchange(ref _lastStatusDivisor, total, previous))
{
return;
}
}
}
private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
@@ -633,17 +630,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
@@ -674,10 +677,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)
@@ -688,11 +693,14 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
ulong modulus = effectiveDivisorData.Modulus;
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
hits.Clear();
return;
}
// EvenPerfectBitScanner only schedules odd divisors greater than one, and the fallback
// above rebuilds Montgomery data when the modulus mismatches, so this guard never
// triggers for production scans.
// if (modulus <= 1UL || (modulus & 1UL) == 0UL)
// {
// hits.Clear();
// return;
// }
if (divisorCycle == 0UL)
{
@@ -795,12 +803,15 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
private void EnsureCapacity(int requiredLength)
{
if (requiredLength <= 0)
{
return;
}
// Prime batches are never empty on the by-divisor scanner path, so the zero-length guard stays commented
// out to remove the redundant branch.
// if (requiredLength <= 0)
// {
// return;
// }
if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _residues 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;
}
@@ -808,6 +819,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
int newCapacity = _bufferCapacity;
if (newCapacity == 0)
{
// EvenPerfectBitScanner hits this branch on the first batch of each session, so keep the seed allocation.
newCapacity = 1;
}
@@ -826,45 +838,36 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
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<ulong>.Shared.Return(_primeDeltas, clearArray: false);
_primeDeltas = null;
}
if (_cycleRemainders is not null)
{
ArrayPool<ulong>.Shared.Return(_cycleRemainders, clearArray: false);
_cycleRemainders = null;
}
if (_residues is not null)
{
ArrayPool<ulong>.Shared.Return(_residues, clearArray: false);
_residues = null;
return;
}
ArrayPool<ulong>.Shared.Return(_primeDeltas, clearArray: false);
ArrayPool<ulong>.Shared.Return(_cycleRemainders!, clearArray: false);
ArrayPool<ulong>.Shared.Return(_residues!, clearArray: false);
_primeDeltas = null;
_cycleRemainders = null;
_residues = null;
_bufferCapacity = 0;
}
private static void ComputePrimeDeltas(ReadOnlySpan<ulong> primes, Span<ulong> 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;
}
@@ -872,11 +875,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
private static void ComputeCycleRemainders(ulong cycleLength, ReadOnlySpan<ulong> primes, ReadOnlySpan<ulong> deltas, Span<ulong> remainders)
{
if (cycleLength == 0UL || primes.Length == 0)
{
remainders.Clear();
return;
}
// Divisor cycles are validated before reaching the GPU helper, and prime batches are non-empty, so
// the defensive zero checks remain commented out here.
// if (cycleLength == 0UL || primes.Length == 0)
// {
// remainders.Clear();
// return;
// }
ulong remainder = primes[0] % cycleLength;
remainders[0] = remainder;
@@ -892,6 +897,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
{
if (modulus == 0UL)
{
// A zero modulus indicates a missing cycle entry; preserve the fallback so snapshot regressions remain visible.
return 0UL;
}
+4 -4
View File
@@ -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)
+2 -2
View File
@@ -76,10 +76,10 @@ public readonly struct ResiduePrimeViews(
public readonly ArrayView1D<ulong, Stride1D.Dense> LastSevenPow2 = lastSevenPow2;
}
public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels)
public ref struct GpuKernelLease(GpuPrimeWorkLimiter.Lease limiter, GpuContextLease gpu, KernelContainer kernels)
{
private bool disposedValue;
private IDisposable? _limiter = limiter;
private GpuPrimeWorkLimiter.Lease? _limiter = limiter;
private GpuContextLease? _gpu = gpu;
private KernelContainer _kernels = kernels;
private AcceleratorStream? _stream;
@@ -6,13 +6,13 @@ public static class GpuPrimeWorkLimiter
private static readonly object _lock = new();
private static int _currentLimit = 1;
public static IDisposable Acquire()
public static Lease Acquire()
{
// TODO: Inline the pooled limiter guard from the GpuLimiterThroughputBenchmarks so prime-sieve
// GPU jobs stop allocating Releaser instances on every acquisition.
var sem = _semaphore;
sem.Wait();
return new Releaser(sem);
return new Lease(sem);
}
public static void SetLimit(int value)
@@ -36,12 +36,12 @@ public static class GpuPrimeWorkLimiter
}
}
private sealed class Releaser : IDisposable
public sealed class Lease
{
private readonly SemaphoreSlim _sem;
private bool _disposed;
public Releaser(SemaphoreSlim sem)
public Lease(SemaphoreSlim sem)
{
_sem = sem;
}
+4 -4
View File
@@ -6,14 +6,14 @@ public static class GpuWorkLimiter
private static readonly object _lock = new();
private static int _currentLimit = 1;
public static IDisposable Acquire()
public static Lease Acquire()
{
// TODO: Replace the per-call Releaser allocation with the pooled struct-based guard from the
// GpuLimiterThroughputBenchmarks so limiter acquisition does not allocate when we enter the
// GPU scanning hot path.
var sem = _semaphore;
sem.Wait();
return new Releaser(sem);
return new Lease(sem);
}
public static void SetLimit(int value)
@@ -40,12 +40,12 @@ public static class GpuWorkLimiter
}
}
private sealed class Releaser : IDisposable
public sealed class Lease
{
private readonly SemaphoreSlim _sem;
private bool _disposed;
public Releaser(SemaphoreSlim sem)
public Lease(SemaphoreSlim sem)
{
_sem = sem;
}
@@ -1,6 +1,7 @@
using System;
using System.Buffers;
using System.Collections.Concurrent;
using System.Collections.Generic;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
using System.Numerics;
@@ -19,7 +20,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
private const int RemainderSlotCount = 4;
private int _gpuBatchSize = GpuConstants.ScanBatchSize;
private readonly object _sync = new();
private ulong _divisorLimit;
private bool _isConfigured;
@@ -64,34 +64,30 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
public void ConfigureFromMaxPrime(ulong maxPrime)
{
lock (_sync)
{
_divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
_isConfigured = true;
}
// Configuration happens once during startup, so no synchronization is required to cache the divisor limit.
_divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
_isConfigured = true;
}
public bool IsPrime(ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null)
{
ulong allowedMax;
int batchCapacity;
// The CLI always configures the tester during startup, so the historical guard is commented out to avoid branching.
// if (!_isConfigured)
// {
// throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
// }
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
// Configuration is immutable, so the cached values can be read without locking.
ulong allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit);
int batchCapacity = _gpuBatchSize;
allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit);
batchCapacity = _gpuBatchSize;
}
if (allowedMax < 3UL)
{
divisorsExhausted = true;
return true;
}
// The GPU path always receives an allowed maximum well above the minimal divisor, so the legacy
// guard stays commented out to document the invariant without branching.
// if (allowedMax < 3UL)
// {
// divisorsExhausted = true;
// return true;
// }
bool composite;
bool coveredRange;
@@ -127,6 +123,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
resources.CycleCandidates,
resources.CycleLengths,
resources.CycleCapacity,
resources,
out lastProcessed,
out coveredRange,
out processedCount
@@ -155,18 +152,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
throw new ArgumentException("allowedMaxValues span must be at least as long as primes span.", nameof(allowedMaxValues));
}
ulong divisorLimit;
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
divisorLimit = _divisorLimit;
}
// The GPU pipeline calls ConfigureFromMaxPrime before preparing candidates, so the legacy guard
// remains commented out to document the invariant without branching.
// if (!_isConfigured)
// {
// throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
// }
ulong divisorLimit = _divisorLimit;
for (int index = 0; index < primes.Length; index++)
{
allowedMaxValues[index] = ComputeAllowedMaxDivisor(primes[index], divisorLimit);
@@ -195,6 +188,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
ulong[] cycleCandidates,
ulong[] cycleLengths,
int cycleCapacity,
BatchResources resources,
out ulong lastProcessed,
out bool coveredRange,
out ulong processedCount)
@@ -222,7 +216,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
ArrayView1D<ulong, Stride1D.Dense> exponentView;
ArrayView1D<byte, Stride1D.Dense> hitsView;
cycleCapacity = Math.Max(1, cycleCapacity);
// Batch resources expose a positive cycle capacity, so skip the redundant clamp here.
// cycleCapacity = Math.Max(1, cycleCapacity);
int chunkCapacity = Math.Min(batchCapacity, cycleCapacity);
int maxThreadsPerGroup = (int)accelerator.MaxNumThreadsPerGroup;
@@ -231,20 +226,22 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
chunkCapacity = Math.Min(chunkCapacity, maxThreadsPerGroup);
}
if (chunkCapacity <= 0)
{
chunkCapacity = cycleCapacity;
}
// cycleCapacity already bounds chunkCapacity above zero, so the fallback would never execute.
// if (chunkCapacity <= 0)
// {
// chunkCapacity = cycleCapacity;
// }
ulong[] filteredDivisors = ArrayPool<ulong>.Shared.Rent(chunkCapacity);
ulong[] divisorGaps = ArrayPool<ulong>.Shared.Rent(chunkCapacity);
UInt128 twoP128 = (UInt128)prime << 1;
if (twoP128 == UInt128.Zero)
{
coveredRange = true;
return false;
}
// Primes above one yield a non-zero 2p value, so the legacy zero guard remains commented out.
// if (twoP128 == UInt128.Zero)
// {
// coveredRange = true;
// return false;
// }
UInt128 allowedMax128 = allowedMax;
UInt128 firstDivisor128 = twoP128 + UInt128.One;
@@ -255,11 +252,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
}
UInt128 maxK128 = (allowedMax128 - UInt128.One) / twoP128;
if (maxK128 == UInt128.Zero)
{
coveredRange = true;
return false;
}
// Once the first divisor falls within the allowed range, at least one multiplier k is admissible,
// so keep the defensive zero check commented out.
// if (maxK128 == UInt128.Zero)
// {
// coveredRange = true;
// return false;
// }
ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128;
ulong currentK = 1UL;
@@ -301,10 +300,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
chunkCount = (int)remainingK;
}
if (chunkCount <= 0)
{
break;
}
// chunkCapacity and the remaining multiplier range are strictly positive here, so the zero
// chunk fallback never triggers on production scans.
// if (chunkCount <= 0)
// {
// break;
// }
Span<ulong> candidateSpan = cycleCandidates.AsSpan(0, chunkCount);
Span<ulong> gapSpan = divisorGaps.AsSpan(0, chunkCount);
@@ -342,10 +343,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
gapView.CopyFromCPU(ref MemoryMarshal.GetReference(gapSpan), chunkCount);
ComputeRemaindersOnGpu(chunkCount, remainder10, 10, gapView, deltaView, remainder10View, remainderDeltaKernel, remainderScanKernel);
ComputeRemaindersOnGpu(chunkCount, remainder8, 8, gapView, deltaView, remainder8View, remainderDeltaKernel, remainderScanKernel);
ComputeRemaindersOnGpu(chunkCount, remainder5, 5, gapView, deltaView, remainder5View, remainderDeltaKernel, remainderScanKernel);
ComputeRemaindersOnGpu(chunkCount, remainder3, 3, gapView, deltaView, remainder3View, remainderDeltaKernel, remainderScanKernel);
ComputeRemaindersOnGpu(chunkCount, remainder10, 10, gapView, deltaView, remainder10View, remainderDeltaKernel, remainderScanKernel, resources.GetRemainderScanConfig);
ComputeRemaindersOnGpu(chunkCount, remainder8, 8, gapView, deltaView, remainder8View, remainderDeltaKernel, remainderScanKernel, resources.GetRemainderScanConfig);
ComputeRemaindersOnGpu(chunkCount, remainder5, 5, gapView, deltaView, remainder5View, remainderDeltaKernel, remainderScanKernel, resources.GetRemainderScanConfig);
ComputeRemaindersOnGpu(chunkCount, remainder3, 3, gapView, deltaView, remainder3View, remainderDeltaKernel, remainderScanKernel, resources.GetRemainderScanConfig);
BuildCandidateMaskOnGpu(chunkCount, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView, maskSpan, candidateMaskKernel);
@@ -463,18 +464,18 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
ArrayView1D<byte, Stride1D.Dense> deltaView,
ArrayView1D<byte, Stride1D.Dense> remainderView,
Action<Index1D, ArrayView<ulong>, ArrayView<byte>, byte> deltaKernel,
Action<KernelConfig, Index1D, ArrayView<byte>, ArrayView<byte>, byte, byte> scanKernel)
Action<KernelConfig, Index1D, ArrayView<byte>, ArrayView<byte>, byte, byte> scanKernel,
Func<int, KernelConfig> configProvider)
{
if (count <= 0)
{
return;
}
// GPU batches always contain at least one divisor gap, so the empty-count guard remains commented
// to eliminate redundant branching.
// if (count <= 0)
// {
// return;
// }
deltaKernel(count, gapView, deltaView, modulus);
SharedMemoryConfig sharedMemoryConfig = SharedMemoryConfig.RequestDynamic<int>(count);
Index1D gridDim = new Index1D(1);
Index1D groupDim = new Index1D(count);
KernelConfig kernelConfig = new KernelConfig(gridDim, groupDim, sharedMemoryConfig);
KernelConfig kernelConfig = configProvider(count);
scanKernel(kernelConfig, new Index1D(count), deltaView, remainderView, baseRemainder, modulus);
}
@@ -491,11 +492,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
Span<byte> destination,
Action<Index1D, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, byte, ArrayView<byte>> maskKernel)
{
if (count <= 0)
{
destination.Clear();
return;
}
// Candidate batches always provide at least one entry, so the zero-count guard stays commented out
// to avoid clearing spans unnecessarily.
// if (count <= 0)
// {
// destination.Clear();
// return;
// }
maskKernel(count, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView);
maskView.CopyToCPU(ref MemoryMarshal.GetReference(destination), count);
@@ -503,10 +506,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
private static byte FetchLastRemainder(ArrayView1D<byte, Stride1D.Dense> remainderView, int count)
{
if (count <= 0)
{
return 0;
}
// Kernels always process at least one remainder in this path, so the empty-span guard remains
// commented to document the invariant.
// if (count <= 0)
// {
// return 0;
// }
byte value = 0;
remainderView.SubView(count - 1, 1).CopyToCPU(ref value, 1);
@@ -519,10 +524,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
if (result >= modulus)
{
result -= modulus;
if (result >= modulus)
{
result %= modulus;
}
// The step deltas remain below the modulus, so a single subtraction suffices on this path.
// if (result >= modulus)
// {
// result %= modulus;
// }
}
return (byte)result;
@@ -530,46 +536,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<MontgomeryDivisorData> divisors, ArrayView<ulong> exponents, ArrayView<byte> 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.Pow2ModBinaryGpu(modulus) == 1UL ? (byte)1 : (byte)0;
hits[index] = exponent.Pow2ModWindowedGpu(modulus) == 1UL ? (byte)1 : (byte)0;
}
public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -609,16 +629,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<ulong, Stride1D.Dense> exponentsView = _exponentsBuffer.View;
ArrayView1D<ulong, Stride1D.Dense> resultsView = _resultsBuffer.View;
@@ -780,30 +804,31 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
private static void ComputeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> 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.Pow2ModBinaryGpu(modulus);
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;
}
}
@@ -820,15 +845,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
private void ReturnBatchResources(Accelerator accelerator, BatchResources resources) => _resourcePools.GetOrAdd(accelerator, static _ => []).Add(resources);
private sealed class BatchResources : IDisposable
private sealed class BatchResources
{
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);
@@ -880,10 +907,24 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
internal readonly Action<KernelConfig, Index1D, ArrayView<byte>, ArrayView<byte>, byte, byte> RemainderScanKernel;
internal readonly Action<Index1D, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, byte, ArrayView<byte>> CandidateMaskKernel;
private readonly Dictionary<int, KernelConfig> _remainderScanConfigs = new();
internal readonly int Capacity;
internal readonly int CycleCapacity;
internal KernelConfig GetRemainderScanConfig(int count)
{
if (!_remainderScanConfigs.TryGetValue(count, out KernelConfig config))
{
SharedMemoryConfig sharedMemoryConfig = SharedMemoryConfig.RequestDynamic<int>(count);
config = new KernelConfig(new Index1D(1), new Index1D(count), sharedMemoryConfig);
_remainderScanConfigs[count] = config;
}
return config;
}
public void Dispose()
{
DivisorsBuffer.Dispose();
@@ -912,37 +953,33 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
public ulong GetAllowedMaxDivisor(ulong prime)
{
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
// The CLI configures the GPU tester once at startup, so callers never request divisor limits before
// initialization completes.
// if (!_isConfigured)
// {
// throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
// }
return ComputeAllowedMaxDivisor(prime, _divisorLimit);
}
return ComputeAllowedMaxDivisor(prime, _divisorLimit);
}
public DivisorScanSession CreateDivisorSession()
public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession()
{
lock (_sync)
// Divisor sessions are handed out only after configuration, so this guard would never trigger on the
// production path.
// if (!_isConfigured)
// {
// throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
// }
if (_sessionPool.TryTake(out DivisorScanSession? session))
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
if (_sessionPool.TryTake(out DivisorScanSession? session))
{
session.Reset();
return session;
}
return new DivisorScanSession(this);
session.Reset();
return session;
}
}
IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() => CreateDivisorSession();
return new DivisorScanSession(this);
}
}
@@ -6,7 +6,7 @@ using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable
internal sealed class MersenneNumberDivisorCandidateGpuEvaluator
{
private const int RemainderCount = 6;
@@ -68,12 +68,17 @@ internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable
}
int count = candidates.Length;
if (count == 0)
{
return;
}
// Candidate batches always carry at least one divisor on the by-divisor path, so the empty-count
// guard remains commented out to eliminate the branch.
// if (count == 0)
// {
// return;
// }
EnsureCapacity(count);
if (count > _capacity)
{
EnsureCapacity(count);
}
Span<ulong> hostCandidateSpan = _hostCandidates.AsSpan(0, count);
Span<byte> hostMaskSpan = _hostMask.AsSpan(0, count);
@@ -115,7 +120,12 @@ internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable
return;
}
int newCapacity = _capacity == 0 ? 1 : _capacity;
int newCapacity = _capacity;
// The evaluator always starts with a positive capacity, so the zero-capacity fallback remains commented.
// if (newCapacity == 0)
// {
// newCapacity = 1;
// }
while (newCapacity < required)
{
newCapacity <<= 1;
@@ -188,10 +198,11 @@ internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable
private static int ComputeRemainder(byte remainder, byte step, int index, int modulus)
{
if (modulus == 0)
{
return 0;
}
// Modulus lookups only cover {10, 8, 5, 3, 7, 11}, so the zero-modulus guard stays commented out.
// if (modulus == 0)
// {
// return 0;
// }
long value = remainder;
value += (long)step * index;
@@ -7,7 +7,7 @@ using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable
internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder
{
private readonly GpuContextPool.GpuContextLease _lease;
private readonly Accelerator _accelerator;
@@ -54,10 +54,12 @@ internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable
}
int length = moduli.Length;
if (length == 0)
{
return;
}
// Montgomery batches always include at least one modulus on the by-divisor path, so the empty-length
// guard stays commented out to avoid branching.
// if (length == 0)
// {
// return;
// }
EnsureCapacity(length);
@@ -83,7 +85,12 @@ internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable
return;
}
int newCapacity = _capacity == 0 ? 1 : _capacity;
int newCapacity = _capacity;
// The builder always starts with a positive capacity, so the zero-capacity fallback remains commented.
// if (newCapacity == 0)
// {
// newCapacity = 1;
// }
while (newCapacity < required)
{
newCapacity <<= 1;
@@ -104,21 +111,25 @@ internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable
private static void BuildKernel(Index1D index, ArrayView<ulong> moduli, ArrayView<MontgomeryDivisorData> 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(
@@ -6,7 +6,7 @@ using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable
internal sealed class MersenneNumberDivisorRemainderGpuStepper
{
private const int RemainderCount = 6;
@@ -52,7 +52,7 @@ internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable
{
if (_disposed)
{
return;
throw new InvalidOperationException("MersenneNumberDivisorRemainderGpuStepper was disposed more than once.");
}
_disposed = true;
@@ -71,15 +71,18 @@ internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable
throw new ObjectDisposedException(nameof(MersenneNumberDivisorRemainderGpuStepper));
}
if (advanceCount <= 0 || remainders.Length == 0)
{
return;
}
// GPU remainder batches always advance by a positive chunk count and operate on fixed-size tables, so
// the defensive early-return remains commented out.
// if (advanceCount <= 0 || remainders.Length == 0)
// {
// return;
// }
if (remainders.Length != _length)
{
throw new ArgumentException("Remainder buffer length must match the configured length.");
}
// Production code always calls the stepper with a fixed-length table, so the historical length guard stays commented.
// if (remainders.Length != _length)
// {
// throw new ArgumentException("Remainder buffer length must match the configured length.");
// }
Span<byte> hostRemainderSpan = _hostRemainders.AsSpan(0, _length);
remainders.CopyTo(hostRemainderSpan);
@@ -99,17 +102,20 @@ internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable
private static void AdvanceRemaindersKernel(Index1D index, int advanceCount, ArrayView<byte> remainders, ArrayView<byte> steps, ArrayView<byte> results)
{
int idx = index;
if (idx >= remainders.Length)
{
return;
}
// Kernel launches match the remainder table length, so the bounds guard stays commented to avoid redundant branching.
// if (idx >= remainders.Length)
// {
// return;
// }
int modulus = GetModulus(idx);
if (modulus == 0)
{
results[idx] = remainders[idx];
return;
}
// The modulus selector only emits positive bases (10, 8, 5, 3, 7, 11) for the configured indices, so the zero-modulus
// fallback remains commented out.
// if (modulus == 0)
// {
// results[idx] = remainders[idx];
// return;
// }
int remainder = remainders[idx] % modulus;
int step = steps[idx] % modulus;
@@ -7,7 +7,7 @@ using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable
internal sealed class MersenneNumberDivisorResidueGpuEvaluator
{
private readonly GpuContextPool.GpuContextLease _lease;
private readonly Accelerator _accelerator;
@@ -33,7 +33,7 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable
{
if (_disposed)
{
return;
throw new InvalidOperationException("MersenneNumberDivisorResidueGpuEvaluator was disposed more than once.");
}
_disposed = true;
@@ -51,10 +51,12 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable
}
int length = exponents.Length;
if (length == 0)
{
return;
}
// Residue batches always include at least one exponent when scheduled from the by-divisor scanner, so
// the empty-length guard remains commented out to avoid branching.
// if (length == 0)
// {
// return;
// }
EnsureCapacity(length);
@@ -79,10 +81,11 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable
}
int newCapacity = _capacity;
if (newCapacity == 0)
{
newCapacity = 1;
}
// The evaluator always starts with a positive capacity, so the zero-capacity fallback stays commented.
// if (newCapacity == 0)
// {
// newCapacity = 1;
// }
while (newCapacity < required)
{
@@ -103,19 +106,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];
ulong modulus = divisorData.Modulus;
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
results[idx] = 0UL;
return;
}
// 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.Pow2ModBinaryGpu(modulus);
results[idx] = exponent.Pow2ModWindowedGpu(modulus);
}
}
@@ -327,8 +327,9 @@ public class MersenneNumberLucasLehmerGpuTester
// TODO: Replace StreamWriter with the pooled TextFileWriter pipeline so persisting NTT
// parameters reuses the zero-allocation buffered writes highlighted in the scanner I/O
// benchmarks instead of allocating a new encoder per append.
using var writer = new StreamWriter(ParameterFilePath, append: true);
var writer = new StreamWriter(ParameterFilePath, append: true);
writer.WriteLine($"{length} {modulus.High} {modulus.Low} {primitiveRoot.High} {primitiveRoot.Low}");
writer.Dispose();
}
}
@@ -338,7 +339,7 @@ public class MersenneNumberLucasLehmerGpuTester
ulong foundRoot = 0UL;
bool success = false;
var processorCount = Environment.ProcessorCount;
using var cts = new System.Threading.CancellationTokenSource();
var cts = new System.Threading.CancellationTokenSource();
object sync = new();
// TODO: Move this Parallel.For to the shared low-overhead work scheduler once the NTT parameter
@@ -381,11 +382,13 @@ public class MersenneNumberLucasLehmerGpuTester
{
modulus = new GpuUInt128(0UL, foundCandidate);
primitiveRoot = new GpuUInt128(0UL, foundRoot);
cts.Dispose();
return true;
}
modulus = new GpuUInt128(0UL, 18446744069414584321UL);
primitiveRoot = new GpuUInt128(0UL, 7UL);
cts.Dispose();
return false;
}
@@ -64,31 +64,25 @@ internal static class MontgomeryOddPowerGpu
var accelerator = lease.Accelerator;
var stream = lease.Stream;
using var divisorBuffer = accelerator.Allocate1D<MontgomeryDivisorData>(divisorCount);
using var countBuffer = accelerator.Allocate1D<int>(divisorCount);
using var resultBuffer = accelerator.Allocate1D<ulong>(divisorCount * stride);
var divisorBuffer = accelerator.Allocate1D<MontgomeryDivisorData>(divisorCount);
var countBuffer = accelerator.Allocate1D<int>(divisorCount);
var resultBuffer = accelerator.Allocate1D<ulong>(divisorCount * stride);
try
{
divisorBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(divisors), divisorCount);
countBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(oddPowerCounts), divisorCount);
resultBuffer.MemSetToZero();
divisorBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(divisors), divisorCount);
countBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(oddPowerCounts), divisorCount);
resultBuffer.MemSetToZero();
var kernel = lease.MontgomeryOddPowerKernel;
kernel(stream, divisorCount, divisorBuffer.View, countBuffer.View, resultBuffer.View, stride);
stream.Synchronize();
var kernel = lease.MontgomeryOddPowerKernel;
kernel(stream, divisorCount, divisorBuffer.View, countBuffer.View, resultBuffer.View, stride);
stream.Synchronize();
resultBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(results), divisorCount * stride);
return true;
}
catch
{
return false;
}
finally
{
execution.Dispose();
lease.Dispose();
}
resultBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(results), divisorCount * stride);
resultBuffer.Dispose();
countBuffer.Dispose();
divisorBuffer.Dispose();
execution.Dispose();
lease.Dispose();
return true;
}
}
+1 -1
View File
@@ -47,7 +47,7 @@ public static class NttGpuMath
private static readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, int>> BitReverseKernelCache = new(); // TODO: Precompute bit-reversal kernels during boot so we can remove the concurrent lookup overhead flagged in the benchmarks.
private sealed class SquareCacheEntry : IDisposable
private sealed class SquareCacheEntry
{
public MemoryBuffer1D<GpuUInt128, Stride1D.Dense> Buffer { get; }
public GpuUInt128 Root { get; }
@@ -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.Pow2ModBinaryGpu(divisor.Modulus) == 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<ulong, Stride1D.Dense> compositeStack,
ArrayView1D<byte, Stride1D.Dense> 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<ulong, Stride1D.Dense> compositeStack,
ArrayView1D<byte, Stride1D.Dense> 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,15 +1608,10 @@ internal static partial class PrimeOrderGpuHeuristics
}
}
private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
{
return exponent.Pow2ModBinaryGpu(divisor.Modulus) == 1UL;
}
private static ulong CalculateByDoublingKernel(ulong prime)
{
ulong order = 1UL;
ulong value = 2UL % prime;
ulong value = 2UL; // Candidate primes exceed two, so the modulo would always return two.
while (value != 1UL)
{
@@ -1619,10 +1629,11 @@ internal static partial class PrimeOrderGpuHeuristics
private static GpuPow2ModStatus TryPow2ModBatchInternal(ReadOnlySpan<UInt128> exponents, UInt128 prime, Span<UInt128> 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<ulong> exponents, ulong prime, in MontgomeryDivisorData divisorData, Span<ulong> 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<UInt128> exponents, UInt128 prime, Span<UInt128> 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<GpuUInt128>.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<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> 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<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> 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<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)
@@ -18,8 +18,10 @@ public interface IMersenneNumberDivisorByDivisorTester
void PrepareCandidates(ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues);
public interface IDivisorScanSession : IDisposable
public interface IDivisorScanSession
{
void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits);
void Dispose();
}
}
+83 -40
View File
@@ -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<ulong, FactorCacheEntry>? 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<ulong, FactorCacheEntry>(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;
@@ -156,7 +196,7 @@ public class MersenneDivisorCycles
public static IEnumerable<(ulong divisor, ulong cycleLength)> EnumerateStream(Stream compressor)
{
// Binary pairs: (ulong divisor, ulong cycle)
using var reader = new BinaryReader(compressor, Encoding.UTF8, leaveOpen: true);
var reader = new BinaryReader(compressor, Encoding.UTF8, leaveOpen: true);
while (true)
{
ulong d, c;
@@ -167,6 +207,7 @@ public class MersenneDivisorCycles
}
catch (EndOfStreamException)
{
reader.Dispose();
yield break;
}
@@ -789,7 +830,7 @@ public class MersenneDivisorCycles
// Start fresh to avoid mixing formats
outputStream.SetLength(0L);
outputStream.Position = 0L;
using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false);
var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false);
for (var taskIndex = 0; taskIndex < threads; taskIndex++)
{
ulong threadStart = start + (ulong)taskIndex * blockSize;
@@ -871,6 +912,7 @@ public class MersenneDivisorCycles
}
Task.WaitAll(tasks);
writer.Dispose();
}
public static void GenerateGpu(string path, ulong maxDivisor, int batchSize = 1_000_000, long skipCount = 0L, long nextPosition = 0L)
@@ -907,7 +949,7 @@ public class MersenneDivisorCycles
outputStream.Position = 0L;
}
using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false);
var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false);
try
{
while (skipCount > 0L && start <= maxDivisor)
@@ -982,6 +1024,7 @@ public class MersenneDivisorCycles
pool.Return(outCycles, clearArray: false);
}
writer.Dispose();
stream.Dispose();
lease.Dispose();
@@ -180,11 +180,13 @@ public static class MersenneNumberDivisorByDivisorTester
return;
}
if (maxPrime <= 1UL)
{
Console.WriteLine("The filter specified by --filter-p must contain at least one prime exponent greater than 1 for --mersenne=bydivisor.");
return;
}
// EvenPerfectBitScanner removes candidates <= 1 before reaching this point,
// so this guard never triggers during production runs.
// if (maxPrime <= 1UL)
// {
// Console.WriteLine("The filter specified by --filter-p must contain at least one prime exponent greater than 1 for --mersenne=bydivisor.");
// return;
// }
tester.ConfigureFromMaxPrime(maxPrime);
@@ -3,7 +3,7 @@ using System.Buffers;
namespace PerfectNumbers.Core
{
internal sealed class PrimeFactorPartialResult : IDisposable
internal sealed class PrimeFactorPartialResult
{
[ThreadStatic]
private static PrimeFactorPartialResult? s_poolHead;
+1 -1
View File
@@ -176,7 +176,7 @@ public sealed class PrimeTester(bool useInternal = false)
DevicePrimes.View.CopyFromCPU(primes);
}
internal sealed class ScratchBuffers : IDisposable
internal sealed class ScratchBuffers
{
public MemoryBuffer1D<ulong, Stride1D.Dense> Input { get; private set; }
public MemoryBuffer1D<byte, Stride1D.Dense> Output { get; private set; }
+1 -1
View File
@@ -3,7 +3,7 @@ namespace PerfectNumbers.Core;
/// <summary>
/// Provides thread-safe file writing for both text and binary data.
/// </summary>
public sealed class TextFileWriter : IDisposable
public sealed class TextFileWriter
{
private readonly FileStream _stream;
private StreamWriter? _textWriter;
+392 -34
View File
@@ -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,45 +292,112 @@ 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;
// }
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2ModBinaryGpu(this ulong exponent, ulong modulus)
{
if (modulus <= 1UL)
// 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 0UL;
return Pow2ModBinaryFallback(exponent, modulus);
}
int bitLength = GetPortableBitLength(exponent);
int windowSize = GetWindowSize(bitLength);
int oddPowerCount = 1 << (windowSize - 1);
Span<ulong> 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;
ulong result = 1UL % modulus;
ulong baseValue = 2UL % modulus;
// 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.MulMod64GpuCompatible(baseValue, modulus);
result = result.MulMod64(baseValue, modulus);
}
remainingExponent >>= 1;
@@ -335,7 +406,103 @@ public static class ULongExtensions
break;
}
baseValue = baseValue.MulMod64GpuCompatible(baseValue, modulus);
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;
@@ -363,7 +530,9 @@ public static class ULongExtensions
ulong modulus = divisor.Modulus;
if (exponent == 0UL)
{
return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus;
return keepMontgomery
? divisor.MontgomeryOne
: 1UL; // Montgomery moduli stay above one on the scanning path, so the modulo would always yield one.
}
if (exponent <= Pow2WindowFallbackThreshold)
@@ -378,10 +547,16 @@ public static class ULongExtensions
ulong result = divisor.MontgomeryOne;
ulong nPrime = divisor.NPrime;
Span<ulong> oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount] : default;
Span<ulong> 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];
}
@@ -449,7 +624,9 @@ public static class ULongExtensions
ulong modulus = divisor.Modulus;
if (exponent == 0UL)
{
return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus;
return keepMontgomery
? divisor.MontgomeryOne
: 1UL; // Montgomery moduli stay above one on the scanning path, so the modulo would always yield one.
}
if (exponent <= Pow2WindowFallbackThreshold)
@@ -579,6 +756,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<ulong> 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; // Moduli arrive as odd primes greater than two, so the modulo would always return two.
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<ulong> destination)
{
int oddPowerCount = destination.Length;
+7 -4
View File
@@ -335,10 +335,13 @@ if (patchResults.Any(static r => !r.Success))
static List<PatchFile> ParsePatchFromFile(string patchFilePath)
{
using var stream = File.OpenRead(patchFilePath);
using var reader = new StreamReader(stream, Encoding.UTF8, detectEncodingFromByteOrderMarks: true, bufferSize: 16_384, leaveOpen: false);
using var lineReader = new PatchLineReader(reader);
var stream = File.OpenRead(patchFilePath);
var reader = new StreamReader(stream, Encoding.UTF8, detectEncodingFromByteOrderMarks: true, bufferSize: 16_384, leaveOpen: false);
var lineReader = new PatchLineReader(reader);
var files = ParsePatch(lineReader);
lineReader.Dispose();
reader.Dispose();
stream.Dispose();
return files;
}
@@ -500,7 +503,7 @@ static void FinalizePatchFile(
files.Add(new PatchFile(originalPath, modifiedPath, hunks ?? new List<PatchHunk>(), type));
}
private sealed class PatchLineReader : IDisposable
private sealed class PatchLineReader
{
private const int ReadBufferSize = 4096;
private const int InitialLineBufferSize = 256;
-18
View File
@@ -1,18 +0,0 @@
# GPU Modular Arithmetic in the `--mersenne=bydivisor` Path
The GPU implementation now routes single-use exponent checks through standard modular arithmetic while keeping Montgomery reduction only where repeated multiplications reuse the same modulus.
## Kernels using standard modular arithmetic
- `CheckKernel` in `MersenneNumberDivisorByDivisorGpuTester` evaluates each divisor candidate with `Pow2ModBinaryGpu` and compares the result directly against 1. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L561-L645】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
- `ComputeExponentKernel` in the same tester precomputes batched residues with the plain helper so CPU lookups receive standard-domain results. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L780-L791】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
- `ComputeResiduesKernel` inside `MersenneNumberDivisorResidueGpuEvaluator` now mirrors that logic when building residue tables for mixed CPU/GPU scans. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs†L102-L114】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
- The order-search helpers in `PrimeOrderGpuHeuristics` (`TrySpecialMaxKernel`, `ExponentLoweringKernel`, and `Pow2EqualsOneKernel`) also rely on `Pow2ModBinaryGpu` for their equality checks. 【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L798-L866】【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L1596-L1598】
## Kernels still using Montgomery arithmetic
- `MontgomeryOddPowerKernelScan` in `GpuKernelPool` continues to build windowed odd-power tables with `MontgomeryMultiply`, which amortizes the Montgomery setup over many multiplications per modulus. 【F:PerfectNumbers.Core/Gpu/GpuKernelPool.cs†L200-L237】【F:PerfectNumbers.Core/ULongExtensions.cs†L532-L559】
## Shared helper
- `Pow2ModBinaryGpu` in `ULongExtensions` implements the GPU-friendly binary exponentiation loop that powers the standard modular path across these kernels. 【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
-13
View File
@@ -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.
+8 -4
View File
@@ -54,10 +54,11 @@ static ProbeResult ProbeBigInteger(Accelerator accelerator)
try
{
var kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<int, Stride1D.Dense>>(BigIntegerKernel);
using var buffer = accelerator.Allocate1D<int>(1);
var buffer = accelerator.Allocate1D<int>(1);
var extent = new Index1D((int)buffer.Length);
kernel(extent, buffer.View);
accelerator.Synchronize();
buffer.Dispose();
return ProbeResult.Success(
"System.Numerics.BigInteger",
"Addition, subtraction, multiplication, division, and modulo compile on the CPU accelerator.");
@@ -73,9 +74,10 @@ static ProbeResult ProbeEInteger(Accelerator accelerator)
try
{
var kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<int, Stride1D.Dense>>(EIntegerKernel);
using var buffer = accelerator.Allocate1D<int>(1);
var buffer = accelerator.Allocate1D<int>(1);
kernel(new Index1D((int)buffer.Length), buffer.View);
accelerator.Synchronize();
buffer.Dispose();
return ProbeResult.Success(
"PeterO.Numbers.EInteger",
"Addition, subtraction, multiplication, division, and remainder compile on the CPU accelerator.");
@@ -93,7 +95,7 @@ static ProbeResult ProbeERational(Accelerator accelerator)
try
{
var kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<int, Stride1D.Dense>>(ERationalKernel);
using var buffer = accelerator.Allocate1D<int>(1);
var buffer = accelerator.Allocate1D<int>(1);
kernel(new Index1D((int)buffer.Length), buffer.View);
accelerator.Synchronize();
@@ -101,6 +103,7 @@ static ProbeResult ProbeERational(Accelerator accelerator)
details += hasRemainder
? " ERational exposes a remainder API that requires further investigation."
: " ERational does not expose a remainder/modulo API.";
buffer.Dispose();
return ProbeResult.Success("PeterO.Numbers.ERational", details);
}
catch (Exception ex)
@@ -120,9 +123,10 @@ static ProbeResult ProbeOpenNumericPrimesIsPrime(Accelerator accelerator)
try
{
var kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<int, Stride1D.Dense>>(OpenNumericPrimesIsPrimeKernel);
using var buffer = accelerator.Allocate1D<int>(4);
var buffer = accelerator.Allocate1D<int>(4);
kernel(new Index1D((int)buffer.Length), buffer.View);
accelerator.Synchronize();
buffer.Dispose();
return ProbeResult.Success(
"Open.Numeric.Primes.Prime.Numbers.IsPrime",
"IsPrime compiles and runs on the CPU accelerator kernel.");