Merge pull request #18 from mlemanczyk/codex/modify-evenperfectbitscanner-mersenne-behavior

Handle 128-bit divisors in Mersenne scans
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-09-15 14:16:07 +02:00
commit 87c5ade66e
5 files changed
+261 -133

No files matched your search

+76 -4
View File
@@ -2,6 +2,8 @@ using FluentAssertions;
using Xunit; using Xunit;
using PerfectNumbers.Core; using PerfectNumbers.Core;
using PerfectNumbers.Core.Gpu; using PerfectNumbers.Core.Gpu;
using System;
using System.Numerics;
using System.Reflection; using System.Reflection;
namespace EvenPerfectBitScanner.Tests; namespace EvenPerfectBitScanner.Tests;
@@ -79,7 +81,7 @@ public class ProgramTests
var forceCpuProp = typeof(GpuContextPool).GetProperty("ForceCpu"); var forceCpuProp = typeof(GpuContextPool).GetProperty("ForceCpu");
useDivisorField.SetValue(null, true); useDivisorField.SetValue(null, true);
divisorField.SetValue(null, 0UL); divisorField.SetValue(null, UInt128.Zero);
testerField.SetValue(null, new MersenneNumberDivisorGpuTester()); testerField.SetValue(null, new MersenneNumberDivisorGpuTester());
primeField.SetValue(null, new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true)); primeField.SetValue(null, new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true));
residueField.SetValue(null, new ThreadLocal<ModResidueTracker>(() => new ModResidueTracker(ResidueModel.Identity, 2UL, true), trackAllValues: true)); residueField.SetValue(null, new ThreadLocal<ModResidueTracker>(() => new ModResidueTracker(ResidueModel.Identity, 2UL, true), trackAllValues: true));
@@ -88,9 +90,9 @@ public class ProgramTests
try try
{ {
Program.IsEvenPerfectCandidate(11UL, 32).Should().BeFalse(); Program.IsEvenPerfectCandidate(11UL, 32UL).Should().BeFalse();
Program.IsEvenPerfectCandidate(5UL, 32).Should().BeTrue(); Program.IsEvenPerfectCandidate(5UL, 32UL).Should().BeTrue();
Program.IsEvenPerfectCandidate(127UL, 32).Should().BeTrue(); Program.IsEvenPerfectCandidate(127UL, 32UL).Should().BeTrue();
} }
finally finally
{ {
@@ -103,6 +105,75 @@ public class ProgramTests
} }
} }
[Fact]
public void Divisor_mode_sets_detailedCheck_based_on_exhaustion()
{
var useDivisorField = typeof(Program).GetField("_useDivisor", BindingFlags.NonPublic | BindingFlags.Static)!;
var divisorField = typeof(Program).GetField("_divisor", BindingFlags.NonPublic | BindingFlags.Static)!;
var testerField = typeof(Program).GetField("_divisorTester", BindingFlags.NonPublic | BindingFlags.Static)!;
var candidatesField = typeof(MersenneNumberDivisorGpuTester).GetField("_divisorCandidates", BindingFlags.NonPublic | BindingFlags.Static)!;
var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!;
var residueField = typeof(Program).GetField("PResidue", BindingFlags.NonPublic | BindingFlags.Static)!;
var forceCpuProp = typeof(GpuContextPool).GetProperty("ForceCpu");
useDivisorField.SetValue(null, true);
divisorField.SetValue(null, UInt128.Zero);
testerField.SetValue(null, new MersenneNumberDivisorGpuTester());
primeField.SetValue(null, new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true));
residueField.SetValue(null, new ThreadLocal<ModResidueTracker>(() => new ModResidueTracker(ResidueModel.Identity, 2UL, true), trackAllValues: true));
candidatesField.SetValue(null, Array.Empty<(ulong, uint)>());
forceCpuProp!.SetValue(null, true);
try
{
Program.IsEvenPerfectCandidate(11UL, 0UL, out bool searched, out bool detailedCheck).Should().BeTrue();
searched.Should().BeTrue();
detailedCheck.Should().BeFalse();
}
finally
{
useDivisorField.SetValue(null, false);
testerField.SetValue(null, null);
candidatesField.SetValue(null, null);
primeField.SetValue(null, null);
residueField.SetValue(null, null);
forceCpuProp!.SetValue(null, false);
}
}
[Fact]
public void Divisor_mode_accepts_large_search_limit()
{
var useDivisorField = typeof(Program).GetField("_useDivisor", BindingFlags.NonPublic | BindingFlags.Static)!;
var divisorField = typeof(Program).GetField("_divisor", BindingFlags.NonPublic | BindingFlags.Static)!;
var testerField = typeof(Program).GetField("_divisorTester", BindingFlags.NonPublic | BindingFlags.Static)!;
var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!;
var residueField = typeof(Program).GetField("PResidue", BindingFlags.NonPublic | BindingFlags.Static)!;
var forceCpuProp = typeof(GpuContextPool).GetProperty("ForceCpu");
useDivisorField.SetValue(null, true);
divisorField.SetValue(null, (UInt128)7);
testerField.SetValue(null, new MersenneNumberDivisorGpuTester());
primeField.SetValue(null, new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true));
residueField.SetValue(null, new ThreadLocal<ModResidueTracker>(() => new ModResidueTracker(ResidueModel.Identity, 2UL, true), trackAllValues: true));
forceCpuProp!.SetValue(null, true);
try
{
Program.IsEvenPerfectCandidate(3UL, 4_000_000_000UL, out bool searched, out bool detailed).Should().BeFalse();
searched.Should().BeTrue();
detailed.Should().BeTrue();
}
finally
{
useDivisorField.SetValue(null, false);
testerField.SetValue(null, null);
primeField.SetValue(null, null);
residueField.SetValue(null, null);
forceCpuProp!.SetValue(null, false);
}
}
[Fact] [Fact]
[Trait("Category", "Fast")] [Trait("Category", "Fast")]
public void Main_displays_help_when_requested() public void Main_displays_help_when_requested()
@@ -127,3 +198,4 @@ public class ProgramTests
output.Should().Contain("--primes-device=cpu|gpu"); output.Should().Contain("--primes-device=cpu|gpu");
} }
} }
+37 -37
View File
@@ -25,8 +25,8 @@ internal static class Program
private static int _primeCount; private static int _primeCount;
private static bool _primeFoundAfterInit; private static bool _primeFoundAfterInit;
private static bool _useGcdFilter; private static bool _useGcdFilter;
private static bool _useDivisor; private static bool _useDivisor;
private static ulong _divisor; private static UInt128 _divisor;
private static MersenneNumberDivisorGpuTester? _divisorTester; private static MersenneNumberDivisorGpuTester? _divisorTester;
private static ulong? _orderWarmupLimitOverride; private static ulong? _orderWarmupLimitOverride;
private static unsafe delegate*<ulong, ref ulong, ulong> _transformP; private static unsafe delegate*<ulong, ref ulong, ulong> _transformP;
@@ -58,9 +58,9 @@ internal static class Program
bool showHelp = false; bool showHelp = false;
bool useBitTransform = false; bool useBitTransform = false;
bool useLucas = false; bool useLucas = false;
bool useResidue = false; // M_p test via residue divisors (replaces LL/incremental) bool useResidue = false; // M_p test via residue divisors (replaces LL/incremental)
bool useDivisor = false; // M_p divisibility by specific divisor bool useDivisor = false; // M_p divisibility by specific divisor
ulong divisor = 0UL; UInt128 divisor = UInt128.Zero;
// Device routing // Device routing
bool useGpuCycles = true; bool useGpuCycles = true;
bool mersenneOnGpu = true; // controls Lucas/incremental/pow2mod device bool mersenneOnGpu = true; // controls Lucas/incremental/pow2mod device
@@ -69,9 +69,9 @@ internal static class Program
ulong residueKMax = 5_000_000UL; ulong residueKMax = 5_000_000UL;
string filterFile = string.Empty; string filterFile = string.Empty;
string cyclesPath = DefaultCyclesPath; string cyclesPath = DefaultCyclesPath;
int cyclesBatchSize = 512; int cyclesBatchSize = 512;
bool continueCyclesGeneration = false; bool continueCyclesGeneration = false;
int divisorCyclesSearchLimit = PerfectNumberConstants.ExtraDivisorCycleSearchLimit; ulong divisorCyclesSearchLimit = PerfectNumberConstants.ExtraDivisorCycleSearchLimit;
// NTT backend selection (GPU): reference vs staged // NTT backend selection (GPU): reference vs staged
for (int i = 0; i < args.Length; i++) for (int i = 0; i < args.Length; i++)
@@ -128,16 +128,16 @@ internal static class Program
kernelType = GpuKernelType.Incremental; kernelType = GpuKernelType.Incremental;
} }
} }
else if (arg.StartsWith("--divisor=", StringComparison.OrdinalIgnoreCase)) else if (arg.StartsWith("--divisor=", StringComparison.OrdinalIgnoreCase))
{ {
int eq = arg.IndexOf('='); int eq = arg.IndexOf('=');
divisor = ulong.Parse(arg.AsSpan(eq + 1)); divisor = UInt128.Parse(arg.AsSpan(eq + 1));
} }
else if (arg.StartsWith("--divisor-cycles-limit=", StringComparison.OrdinalIgnoreCase)) else if (arg.StartsWith("--divisor-cycles-limit=", StringComparison.OrdinalIgnoreCase))
{ {
int eq = arg.IndexOf('='); int eq = arg.IndexOf('=');
divisorCyclesSearchLimit = int.Parse(arg.AsSpan(eq + 1)); divisorCyclesSearchLimit = ulong.Parse(arg.AsSpan(eq + 1));
} }
else if (arg.StartsWith("--residue-max-k=", StringComparison.OrdinalIgnoreCase)) else if (arg.StartsWith("--residue-max-k=", StringComparison.OrdinalIgnoreCase))
{ {
int eq = arg.IndexOf('='); int eq = arg.IndexOf('=');
@@ -438,15 +438,15 @@ internal static class Program
return tester; return tester;
}, trackAllValues: true); }, trackAllValues: true);
} }
else else
{ {
if (divisor == 0UL) if (divisor == UInt128.Zero)
{ {
MersenneNumberDivisorGpuTester.BuildDivisorCandidates(); MersenneNumberDivisorGpuTester.BuildDivisorCandidates();
} }
_divisorTester = new MersenneNumberDivisorGpuTester(); _divisorTester = new MersenneNumberDivisorGpuTester();
} }
// Load RLE blacklist (optional) // Load RLE blacklist (optional)
if (!string.IsNullOrEmpty(_rleBlacklistPath)) if (!string.IsNullOrEmpty(_rleBlacklistPath))
@@ -906,10 +906,10 @@ internal static class Program
return next + value; return next + value;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static bool IsEvenPerfectCandidate(ulong p, int divisorCyclesSearchLimit) => IsEvenPerfectCandidate(p, divisorCyclesSearchLimit, out _, out _); internal static bool IsEvenPerfectCandidate(ulong p, ulong divisorCyclesSearchLimit) => IsEvenPerfectCandidate(p, divisorCyclesSearchLimit, out _, out _);
internal static bool IsEvenPerfectCandidate(ulong p, int divisorCyclesSearchLimit, out bool searchedMersenne, out bool detailedCheck) internal static bool IsEvenPerfectCandidate(ulong p, ulong divisorCyclesSearchLimit, out bool searchedMersenne, out bool detailedCheck)
{ {
searchedMersenne = false; searchedMersenne = false;
detailedCheck = false; detailedCheck = false;
@@ -972,15 +972,15 @@ internal static class Program
return false; return false;
} }
searchedMersenne = true; searchedMersenne = true;
if (_useDivisor) if (_useDivisor)
{ {
return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck); return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck);
} }
detailedCheck = MersenneTesters.Value!.IsMersennePrime(p); detailedCheck = MersenneTesters.Value!.IsMersennePrime(p);
return detailedCheck; return detailedCheck;
} }
// Use ModResidueTracker with a small set of primes to pre-filter composite p. // Use ModResidueTracker with a small set of primes to pre-filter composite p.
private static bool IsCompositeByResidues(ulong p) private static bool IsCompositeByResidues(ulong p)
@@ -1,5 +1,7 @@
using FluentAssertions; using FluentAssertions;
using PerfectNumbers.Core.Gpu; using PerfectNumbers.Core.Gpu;
using System.Numerics;
using System.Reflection;
using Xunit; using Xunit;
namespace PerfectNumbers.Core.Tests; namespace PerfectNumbers.Core.Tests;
@@ -17,4 +19,46 @@ public class MersenneNumberDivisorGpuTesterTests
var tester = new MersenneNumberDivisorGpuTester(); var tester = new MersenneNumberDivisorGpuTester();
tester.IsDivisible(exponent, divisor).Should().Be(expected); tester.IsDivisible(exponent, divisor).Should().Be(expected);
} }
[Fact]
[Trait("Category", "Fast")]
public void IsDivisible_handles_large_divisors()
{
var tester = new MersenneNumberDivisorGpuTester();
UInt128 divisor = (UInt128.One << 65) - UInt128.One;
tester.IsDivisible(65UL, divisor).Should().BeTrue();
tester.IsDivisible(65UL, divisor + 2).Should().BeFalse();
}
[Fact]
[Trait("Category", "Fast")]
public void IsPrime_sets_divisorsExhausted_false_when_search_range_not_exhausted()
{
var tester = new MersenneNumberDivisorGpuTester();
typeof(MersenneNumberDivisorGpuTester)
.GetField("_divisorCandidates", BindingFlags.NonPublic | BindingFlags.Static)!
.SetValue(null, Array.Empty<(ulong, uint)>());
tester.IsPrime(11UL, UInt128.Zero, 0UL, out bool divisorsExhausted).Should().BeTrue();
divisorsExhausted.Should().BeFalse();
}
[Fact]
[Trait("Category", "Fast")]
public void IsPrime_sets_divisorsExhausted_true_when_divisible_by_specified_divisor()
{
var tester = new MersenneNumberDivisorGpuTester();
tester.IsPrime(3UL, 7UL, 0UL, out bool divisorsExhausted).Should().BeFalse();
divisorsExhausted.Should().BeTrue();
}
[Fact]
[Trait("Category", "Fast")]
public void IsPrime_accepts_large_search_limits()
{
var tester = new MersenneNumberDivisorGpuTester();
tester.IsPrime(3UL, 7UL, ulong.MaxValue, out bool exhausted).Should().BeFalse();
exhausted.Should().BeTrue();
}
} }
@@ -1,4 +1,5 @@
using System.Collections.Concurrent; using System.Collections.Concurrent;
using System.Numerics;
using ILGPU; using ILGPU;
using ILGPU.Algorithms; using ILGPU.Algorithms;
using ILGPU.Runtime; using ILGPU.Runtime;
@@ -7,11 +8,11 @@ namespace PerfectNumbers.Core.Gpu;
public sealed class MersenneNumberDivisorGpuTester public sealed class MersenneNumberDivisorGpuTester
{ {
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ulong, ArrayView<byte>>> _kernelCache = new(); private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, GpuUInt128, ArrayView<byte>>> _kernelCache = new();
private readonly ConcurrentDictionary<Accelerator, MemoryBuffer1D<byte, Stride1D.Dense>> _resultBuffers = new(); private readonly ConcurrentDictionary<Accelerator, MemoryBuffer1D<byte, Stride1D.Dense>> _resultBuffers = new();
private Action<Index1D, ulong, ulong, ArrayView<byte>> GetKernel(Accelerator accelerator) => private Action<Index1D, ulong, GpuUInt128, ArrayView<byte>> GetKernel(Accelerator accelerator) =>
_kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ulong, ArrayView<byte>>(Kernel)); _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, ArrayView<byte>>(Kernel));
public static void BuildDivisorCandidates() public static void BuildDivisorCandidates()
{ {
@@ -33,13 +34,13 @@ public sealed class MersenneNumberDivisorGpuTester
_divisorCandidates = list[..count].ToArray(); _divisorCandidates = list[..count].ToArray();
} }
public bool IsDivisible(ulong exponent, ulong divisor) public bool IsDivisible(ulong exponent, UInt128 divisor)
{ {
var gpu = GpuContextPool.RentPreferred(preferCpu: false); var gpu = GpuContextPool.RentPreferred(preferCpu: false);
var accelerator = gpu.Accelerator; var accelerator = gpu.Accelerator;
var kernel = GetKernel(accelerator); var kernel = GetKernel(accelerator);
var resultBuffer = _resultBuffers.GetOrAdd(accelerator, acc => acc.Allocate1D<byte>(1)); var resultBuffer = _resultBuffers.GetOrAdd(accelerator, acc => acc.Allocate1D<byte>(1));
kernel(1, exponent, divisor, resultBuffer.View); kernel(1, exponent, (GpuUInt128)divisor, resultBuffer.View);
accelerator.Synchronize(); accelerator.Synchronize();
bool divisible = resultBuffer.GetAsArray1D()[0] != 0; bool divisible = resultBuffer.GetAsArray1D()[0] != 0;
gpu.Dispose(); gpu.Dispose();
@@ -48,98 +49,109 @@ public sealed class MersenneNumberDivisorGpuTester
private static (ulong divisor, uint cycle)[]? _divisorCandidates; private static (ulong divisor, uint cycle)[]? _divisorCandidates;
public bool IsPrime(ulong p, ulong d, int divisorCyclesSearchLimit, out bool detailedCheck) public bool IsPrime(ulong p, UInt128 d, ulong divisorCyclesSearchLimit, out bool divisorsExhausted)
{ {
if (d != 0UL) if (d != UInt128.Zero)
{ {
if (IsDivisible(p, d)) if (IsDivisible(p, d))
{ {
detailedCheck = true; divisorsExhausted = true;
return false; return false;
} }
detailedCheck = false; divisorsExhausted = false;
return true; return true;
} }
var candidates = _divisorCandidates!; var candidates = _divisorCandidates!;
int k, len = candidates.Length; int k, len = candidates.Length;
uint cycle; uint cycle;
// k is simple iterator here for (k = 0; k < len; k++)
for (k = 0; k < len; k++) {
{ (ulong dSmall, cycle) = candidates[k];
(d, cycle) = candidates[k]; if (p % cycle != 0UL)
if (p % cycle != 0UL)
{
continue;
}
if (IsDivisible(p, d))
{
detailedCheck = true;
return false;
}
}
ulong kMul2;
var divisorCycles = MersenneDivisorCycles.Shared;
for (k = 1; k <= divisorCyclesSearchLimit; k++)
{
kMul2 = (ulong)k << 1;
if (p > ulong.MaxValue / kMul2)
{
break;
}
d = kMul2 * p + 1UL;
if (p < 64UL && d == ((1UL << (int)p) - 1UL))
{ {
continue; continue;
} }
// kkMul2 now becomes cycle. We're reusing existing variable for the best performance
kMul2 = divisorCycles.GetCycle(d);
if (p % kMul2 != 0UL)
{
continue;
}
if (IsDivisible(p, d)) if (IsDivisible(p, dSmall))
{ {
detailedCheck = true; divisorsExhausted = true;
return false; return false;
} }
} }
detailedCheck = false; UInt128 kMul2;
return true; var divisorCycles = MersenneDivisorCycles.Shared;
} UInt128 maxK2 = UInt128.MaxValue / ((UInt128)p << 1);
ulong limit = divisorCyclesSearchLimit;
if ((UInt128)limit > maxK2)
{
limit = (ulong)maxK2;
}
private static void Kernel(Index1D _, ulong exponent, ulong divisor, ArrayView<byte> result) ulong k2;
{ for (k2 = 1UL; k2 <= limit; k2++)
if (divisor <= 1UL) {
{ kMul2 = (UInt128)k2 << 1;
result[0] = 0; d = kMul2 * p + UInt128.One;
return; if (p < 64UL && d == ((UInt128)1 << (int)p) - UInt128.One)
} {
continue;
}
GpuUInt128 baseVal, mod; UInt128 cycle128 = MersenneDivisorCycles.GetCycle(d);
int x = 63 - XMath.LeadingZeroCount(divisor); if ((UInt128)p % cycle128 != UInt128.Zero)
if (x == 0) {
{ continue;
mod = new GpuUInt128(0UL, divisor); }
// baseVal is temp value here. We're reusing existing variable in a different meaning for the best performance.
baseVal = GpuUInt128.Pow2Mod(exponent, mod);
result[0] = baseVal.High == 0UL && baseVal.Low == 1UL ? (byte)1 : (byte)0;
return;
}
ulong ux = (ulong)x; if (IsDivisible(p, d))
mod = new GpuUInt128(0UL, divisor); {
baseVal = new GpuUInt128(0UL, 1UL << x); // 2^x ≡ -y (mod divisor) divisorsExhausted = true;
baseVal.ModPow(exponent / ux, mod); return false;
var part2 = GpuUInt128.Pow2Mod(exponent % ux, mod); }
baseVal.MulMod(part2, mod); }
result[0] = baseVal.High == 0UL && baseVal.Low == 1UL ? (byte)1 : (byte)0;
} divisorsExhausted = (UInt128)divisorCyclesSearchLimit >= maxK2;
return true;
}
private static void Kernel(Index1D _, ulong exponent, GpuUInt128 divisor, ArrayView<byte> result)
{
if (divisor.High == 0UL && divisor.Low <= 1UL)
{
result[0] = 0;
return;
}
GpuUInt128 baseVal, mod;
int x = divisor.High == 0UL
? 63 - XMath.LeadingZeroCount(divisor.Low)
: 127 - XMath.LeadingZeroCount(divisor.High);
mod = divisor;
if (x == 0)
{
baseVal = GpuUInt128.Pow2Mod(exponent, mod);
result[0] = baseVal.High == 0UL && baseVal.Low == 1UL ? (byte)1 : (byte)0;
return;
}
ulong ux = (ulong)x;
if (x >= 64)
{
baseVal = new GpuUInt128(1UL << (x - 64), 0UL);
}
else
{
baseVal = new GpuUInt128(0UL, 1UL << x);
}
baseVal.ModPow(exponent / ux, mod);
var part2 = GpuUInt128.Pow2Mod(exponent % ux, mod);
baseVal.MulMod(part2, mod);
result[0] = baseVal.High == 0UL && baseVal.Low == 1UL ? (byte)1 : (byte)0;
}
} }
@@ -5,6 +5,6 @@ public static class PerfectNumberConstants
public const ulong BiggestKnownEvenPerfectP = 136279841UL; public const ulong BiggestKnownEvenPerfectP = 136279841UL;
public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000;
public const int MaxQForDivisorCycles = 4_000_000; public const int MaxQForDivisorCycles = 4_000_000;
public const int ExtraDivisorCycleSearchLimit = 64; public const ulong ExtraDivisorCycleSearchLimit = 64UL;
} }