diff --git a/EvenPerfectBitScanner.Tests/ProgramTests.cs b/EvenPerfectBitScanner.Tests/ProgramTests.cs index 7858749..02aae50 100644 --- a/EvenPerfectBitScanner.Tests/ProgramTests.cs +++ b/EvenPerfectBitScanner.Tests/ProgramTests.cs @@ -2,6 +2,8 @@ using FluentAssertions; using Xunit; using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; +using System; +using System.Numerics; using System.Reflection; namespace EvenPerfectBitScanner.Tests; @@ -79,7 +81,7 @@ public class ProgramTests var forceCpuProp = typeof(GpuContextPool).GetProperty("ForceCpu"); useDivisorField.SetValue(null, true); - divisorField.SetValue(null, 0UL); + divisorField.SetValue(null, UInt128.Zero); testerField.SetValue(null, new MersenneNumberDivisorGpuTester()); primeField.SetValue(null, new ThreadLocal(() => new PrimeTester(), trackAllValues: true)); residueField.SetValue(null, new ThreadLocal(() => new ModResidueTracker(ResidueModel.Identity, 2UL, true), trackAllValues: true)); @@ -88,9 +90,9 @@ public class ProgramTests try { - Program.IsEvenPerfectCandidate(11UL, 32).Should().BeFalse(); - Program.IsEvenPerfectCandidate(5UL, 32).Should().BeTrue(); - Program.IsEvenPerfectCandidate(127UL, 32).Should().BeTrue(); + Program.IsEvenPerfectCandidate(11UL, 32UL).Should().BeFalse(); + Program.IsEvenPerfectCandidate(5UL, 32UL).Should().BeTrue(); + Program.IsEvenPerfectCandidate(127UL, 32UL).Should().BeTrue(); } 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(() => new PrimeTester(), trackAllValues: true)); + residueField.SetValue(null, new ThreadLocal(() => 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(() => new PrimeTester(), trackAllValues: true)); + residueField.SetValue(null, new ThreadLocal(() => 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] [Trait("Category", "Fast")] public void Main_displays_help_when_requested() @@ -127,3 +198,4 @@ public class ProgramTests output.Should().Contain("--primes-device=cpu|gpu"); } } + diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 7064b9f..0b7452d 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -25,8 +25,8 @@ internal static class Program private static int _primeCount; private static bool _primeFoundAfterInit; private static bool _useGcdFilter; - private static bool _useDivisor; - private static ulong _divisor; + private static bool _useDivisor; + private static UInt128 _divisor; private static MersenneNumberDivisorGpuTester? _divisorTester; private static ulong? _orderWarmupLimitOverride; private static unsafe delegate* _transformP; @@ -58,9 +58,9 @@ internal static class Program bool showHelp = false; bool useBitTransform = false; bool useLucas = false; - bool useResidue = false; // M_p test via residue divisors (replaces LL/incremental) - bool useDivisor = false; // M_p divisibility by specific divisor - ulong divisor = 0UL; + bool useResidue = false; // M_p test via residue divisors (replaces LL/incremental) + bool useDivisor = false; // M_p divisibility by specific divisor + UInt128 divisor = UInt128.Zero; // Device routing bool useGpuCycles = true; bool mersenneOnGpu = true; // controls Lucas/incremental/pow2mod device @@ -69,9 +69,9 @@ internal static class Program ulong residueKMax = 5_000_000UL; string filterFile = string.Empty; string cyclesPath = DefaultCyclesPath; - int cyclesBatchSize = 512; - bool continueCyclesGeneration = false; - int divisorCyclesSearchLimit = PerfectNumberConstants.ExtraDivisorCycleSearchLimit; + int cyclesBatchSize = 512; + bool continueCyclesGeneration = false; + ulong divisorCyclesSearchLimit = PerfectNumberConstants.ExtraDivisorCycleSearchLimit; // NTT backend selection (GPU): reference vs staged for (int i = 0; i < args.Length; i++) @@ -128,16 +128,16 @@ internal static class Program kernelType = GpuKernelType.Incremental; } } - else if (arg.StartsWith("--divisor=", StringComparison.OrdinalIgnoreCase)) - { - int eq = arg.IndexOf('='); - divisor = ulong.Parse(arg.AsSpan(eq + 1)); - } - else if (arg.StartsWith("--divisor-cycles-limit=", StringComparison.OrdinalIgnoreCase)) - { - int eq = arg.IndexOf('='); - divisorCyclesSearchLimit = int.Parse(arg.AsSpan(eq + 1)); - } + else if (arg.StartsWith("--divisor=", StringComparison.OrdinalIgnoreCase)) + { + int eq = arg.IndexOf('='); + divisor = UInt128.Parse(arg.AsSpan(eq + 1)); + } + else if (arg.StartsWith("--divisor-cycles-limit=", StringComparison.OrdinalIgnoreCase)) + { + int eq = arg.IndexOf('='); + divisorCyclesSearchLimit = ulong.Parse(arg.AsSpan(eq + 1)); + } else if (arg.StartsWith("--residue-max-k=", StringComparison.OrdinalIgnoreCase)) { int eq = arg.IndexOf('='); @@ -438,15 +438,15 @@ internal static class Program return tester; }, trackAllValues: true); } - else - { - if (divisor == 0UL) - { - MersenneNumberDivisorGpuTester.BuildDivisorCandidates(); - } + else + { + if (divisor == UInt128.Zero) + { + MersenneNumberDivisorGpuTester.BuildDivisorCandidates(); + } - _divisorTester = new MersenneNumberDivisorGpuTester(); - } + _divisorTester = new MersenneNumberDivisorGpuTester(); + } // Load RLE blacklist (optional) if (!string.IsNullOrEmpty(_rleBlacklistPath)) @@ -906,10 +906,10 @@ internal static class Program return next + value; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static bool IsEvenPerfectCandidate(ulong p, int divisorCyclesSearchLimit) => IsEvenPerfectCandidate(p, divisorCyclesSearchLimit, out _, out _); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + 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; detailedCheck = false; @@ -972,15 +972,15 @@ internal static class Program return false; } - searchedMersenne = true; - if (_useDivisor) - { - return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck); - } + searchedMersenne = true; + if (_useDivisor) + { + return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck); + } - detailedCheck = MersenneTesters.Value!.IsMersennePrime(p); - return detailedCheck; - } + detailedCheck = MersenneTesters.Value!.IsMersennePrime(p); + return detailedCheck; + } // Use ModResidueTracker with a small set of primes to pre-filter composite p. private static bool IsCompositeByResidues(ulong p) diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index d50a48a..18f1258 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -1,5 +1,7 @@ using FluentAssertions; using PerfectNumbers.Core.Gpu; +using System.Numerics; +using System.Reflection; using Xunit; namespace PerfectNumbers.Core.Tests; @@ -17,4 +19,46 @@ public class MersenneNumberDivisorGpuTesterTests var tester = new MersenneNumberDivisorGpuTester(); 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(); + } } + diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs index 13819e8..46d02b4 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs @@ -1,4 +1,5 @@ using System.Collections.Concurrent; +using System.Numerics; using ILGPU; using ILGPU.Algorithms; using ILGPU.Runtime; @@ -7,11 +8,11 @@ namespace PerfectNumbers.Core.Gpu; public sealed class MersenneNumberDivisorGpuTester { - private readonly ConcurrentDictionary>> _kernelCache = new(); + private readonly ConcurrentDictionary>> _kernelCache = new(); private readonly ConcurrentDictionary> _resultBuffers = new(); - private Action> GetKernel(Accelerator accelerator) => - _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(Kernel)); + private Action> GetKernel(Accelerator accelerator) => + _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(Kernel)); public static void BuildDivisorCandidates() { @@ -33,13 +34,13 @@ public sealed class MersenneNumberDivisorGpuTester _divisorCandidates = list[..count].ToArray(); } - public bool IsDivisible(ulong exponent, ulong divisor) - { - var gpu = GpuContextPool.RentPreferred(preferCpu: false); - var accelerator = gpu.Accelerator; + public bool IsDivisible(ulong exponent, UInt128 divisor) + { + var gpu = GpuContextPool.RentPreferred(preferCpu: false); + var accelerator = gpu.Accelerator; var kernel = GetKernel(accelerator); var resultBuffer = _resultBuffers.GetOrAdd(accelerator, acc => acc.Allocate1D(1)); - kernel(1, exponent, divisor, resultBuffer.View); + kernel(1, exponent, (GpuUInt128)divisor, resultBuffer.View); accelerator.Synchronize(); bool divisible = resultBuffer.GetAsArray1D()[0] != 0; gpu.Dispose(); @@ -48,98 +49,109 @@ public sealed class MersenneNumberDivisorGpuTester private static (ulong divisor, uint cycle)[]? _divisorCandidates; - public bool IsPrime(ulong p, ulong d, int divisorCyclesSearchLimit, out bool detailedCheck) - { - if (d != 0UL) - { - if (IsDivisible(p, d)) - { - detailedCheck = true; - return false; - } + public bool IsPrime(ulong p, UInt128 d, ulong divisorCyclesSearchLimit, out bool divisorsExhausted) + { + if (d != UInt128.Zero) + { + if (IsDivisible(p, d)) + { + divisorsExhausted = true; + return false; + } - detailedCheck = false; - return true; - } + divisorsExhausted = false; + return true; + } - var candidates = _divisorCandidates!; - int k, len = candidates.Length; - uint cycle; + var candidates = _divisorCandidates!; + int k, len = candidates.Length; + uint cycle; - // k is simple iterator here - for (k = 0; k < len; k++) - { - (d, cycle) = candidates[k]; - 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)) + for (k = 0; k < len; k++) + { + (ulong dSmall, cycle) = candidates[k]; + if (p % cycle != 0UL) { 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)) - { - detailedCheck = true; - return false; - } - } + if (IsDivisible(p, dSmall)) + { + divisorsExhausted = true; + return false; + } + } - detailedCheck = false; - return true; - } + UInt128 kMul2; + 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 result) - { - if (divisor <= 1UL) - { - result[0] = 0; - return; - } + ulong k2; + for (k2 = 1UL; k2 <= limit; k2++) + { + kMul2 = (UInt128)k2 << 1; + d = kMul2 * p + UInt128.One; + if (p < 64UL && d == ((UInt128)1 << (int)p) - UInt128.One) + { + continue; + } - GpuUInt128 baseVal, mod; - int x = 63 - XMath.LeadingZeroCount(divisor); - if (x == 0) - { - 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; - } + UInt128 cycle128 = MersenneDivisorCycles.GetCycle(d); + if ((UInt128)p % cycle128 != UInt128.Zero) + { + continue; + } - ulong ux = (ulong)x; - mod = new GpuUInt128(0UL, divisor); - baseVal = new GpuUInt128(0UL, 1UL << x); // 2^x ≡ -y (mod divisor) - 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; - } + if (IsDivisible(p, d)) + { + divisorsExhausted = true; + return false; + } + } + + divisorsExhausted = (UInt128)divisorCyclesSearchLimit >= maxK2; + return true; + } + + private static void Kernel(Index1D _, ulong exponent, GpuUInt128 divisor, ArrayView 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; + } } + diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index 6c4f3cf..92fafa9 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -5,6 +5,6 @@ public static class PerfectNumberConstants public const ulong BiggestKnownEvenPerfectP = 136279841UL; public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; public const int MaxQForDivisorCycles = 4_000_000; - public const int ExtraDivisorCycleSearchLimit = 64; + public const ulong ExtraDivisorCycleSearchLimit = 64UL; -} \ No newline at end of file +}