diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 81f5df7..65f0a9c 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -363,21 +363,29 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } - UInt128 step = (UInt128)prime << 1; - if (step == UInt128.Zero) + UInt128 step128 = (UInt128)prime << 1; + if (step128 == UInt128.Zero) { processedAll = true; return false; } + if (step128 > ulong.MaxValue) + { + return CheckDivisorsCpu(prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit); + } + UInt128 limit = allowedMax; - UInt128 divisor = step + UInt128.One; - if (divisor > limit) + UInt128 divisor128 = step128 + UInt128.One; + if (divisor128 > limit) { processedAll = true; return false; } + ulong step = (ulong)step128; + ulong limit64 = allowedMax; + int batchCapacity = Math.Max(1, _batchSize); DivisorCycleCache cycleCache = DivisorCycleCache.Shared; int preferredCycleBatch = cycleCache.PreferredBatchSize; @@ -386,17 +394,16 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv batchCapacity = preferredCycleBatch; } - ulong[] candidateBuffer = ArrayPool.Shared.Rent(batchCapacity); - ulong[] cycleBuffer = ArrayPool.Shared.Rent(batchCapacity); - MontgomeryDivisorData[] divisorDataBuffer = ArrayPool.Shared.Rent(batchCapacity); - byte[] remainderBuffer = ArrayPool.Shared.Rent(batchCapacity * 6); + var candidateEvaluator = new MersenneNumberDivisorCandidateGpuEvaluator(batchCapacity); - Span remainder10Span = remainderBuffer.AsSpan(0, batchCapacity); - Span remainder8Span = remainderBuffer.AsSpan(batchCapacity, batchCapacity); - Span remainder5Span = remainderBuffer.AsSpan(batchCapacity * 2, batchCapacity); - Span remainder3Span = remainderBuffer.AsSpan(batchCapacity * 3, batchCapacity); - Span remainder7Span = remainderBuffer.AsSpan(batchCapacity * 4, batchCapacity); - Span remainder11Span = remainderBuffer.AsSpan(batchCapacity * 5, batchCapacity); + ulong[] candidateBuffer = ArrayPool.Shared.Rent(batchCapacity); + byte[] maskBuffer = ArrayPool.Shared.Rent(batchCapacity); + + Dictionary? factorCache = null; + + ulong cachedModulus = 0UL; + MontgomeryDivisorData cachedDivisorData = default; + bool hasCachedDivisorData = false; byte step10 = (byte)(step % 10UL); byte step8 = (byte)(step % 8UL); @@ -405,31 +412,49 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv byte step7 = (byte)(step % 7UL); byte step11 = (byte)(step % 11UL); - byte remainder10 = (byte)(divisor % 10UL); - byte remainder8 = (byte)(divisor % 8UL); - byte remainder5 = (byte)(divisor % 5UL); - byte remainder3 = (byte)(divisor % 3UL); - byte remainder7 = (byte)(divisor % 7UL); - byte remainder11 = (byte)(divisor % 11UL); + byte remainder10 = (byte)((ulong)divisor128 % 10UL); + byte remainder8 = (byte)((ulong)divisor128 % 8UL); + byte remainder5 = (byte)((ulong)divisor128 % 5UL); + byte remainder3 = (byte)((ulong)divisor128 % 3UL); + byte remainder7 = (byte)((ulong)divisor128 % 7UL); + byte remainder11 = (byte)((ulong)divisor128 % 11UL); bool lastIsSeven = (prime & 3UL) == 3UL; try { - while (divisor <= limit) + while (divisor128 <= limit) { - Span candidateSpan = candidateBuffer.AsSpan(0, batchCapacity); - int chunkCount = 0; - UInt128 localDivisor = divisor; + int chunkCount = ComputeChunkCount(divisor128, limit, step128, batchCapacity); + if (chunkCount <= 0) + { + break; + } - byte localRemainder10 = remainder10; - byte localRemainder8 = remainder8; - byte localRemainder5 = remainder5; - byte localRemainder3 = remainder3; - byte localRemainder7 = remainder7; - byte localRemainder11 = remainder11; + Span candidateSpan = candidateBuffer.AsSpan(0, chunkCount); + Span maskSpan = maskBuffer.AsSpan(0, chunkCount); - while (chunkCount < batchCapacity && localDivisor <= limit) + candidateEvaluator.EvaluateCandidates( + (ulong)divisor128, + step, + limit64, + remainder10, + remainder8, + remainder5, + remainder3, + remainder7, + remainder11, + step10, + step8, + step5, + step3, + step7, + step11, + lastIsSeven, + candidateSpan, + maskSpan); + + for (int i = 0; i < chunkCount; i++) { if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired()) { @@ -437,126 +462,116 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } - ulong candidate = (ulong)localDivisor; - candidateSpan[chunkCount] = candidate; - remainder10Span[chunkCount] = localRemainder10; - remainder8Span[chunkCount] = localRemainder8; - remainder5Span[chunkCount] = localRemainder5; - remainder3Span[chunkCount] = localRemainder3; - remainder7Span[chunkCount] = localRemainder7; - remainder11Span[chunkCount] = localRemainder11; + iterationCounter++; + ulong candidate = candidateSpan[i]; processedCount++; lastProcessed = candidate; - iterationCounter++; - chunkCount++; - localDivisor += step; - localRemainder10 = AddMod(localRemainder10, step10, (byte)10); - localRemainder8 = AddMod(localRemainder8, step8, (byte)8); - localRemainder5 = AddMod(localRemainder5, step5, (byte)5); - localRemainder3 = AddMod(localRemainder3, step3, (byte)3); - localRemainder7 = AddMod(localRemainder7, step7, (byte)7); - localRemainder11 = AddMod(localRemainder11, step11, (byte)11); - } - - if (chunkCount == 0) - { - break; - } - - divisor = localDivisor; - remainder10 = localRemainder10; - remainder8 = localRemainder8; - remainder5 = localRemainder5; - remainder3 = localRemainder3; - remainder7 = localRemainder7; - remainder11 = localRemainder11; - - // Reuse the candidate buffer in-place so the filtered results stay contiguous without additional arrays. - Span chunkCandidates = candidateSpan[..chunkCount]; - int filteredCount = 0; - - for (int i = 0; i < chunkCount; i++) - { - bool admissible = lastIsSeven - ? (remainder10Span[i] == 3 || remainder10Span[i] == 7 || remainder10Span[i] == 9) - : (remainder10Span[i] == 1 || remainder10Span[i] == 3 || remainder10Span[i] == 9); - - if (!admissible) + if (maskSpan[i] == 0) { continue; } - if ((remainder8Span[i] != 1 && remainder8Span[i] != 7) - || remainder3Span[i] == 0 - || remainder5Span[i] == 0 - || remainder7Span[i] == 0 - || remainder11Span[i] == 0) + if (!hasCachedDivisorData || cachedModulus != candidate) { - continue; + cachedDivisorData = MontgomeryDivisorData.FromModulus(candidate); + cachedModulus = candidate; + hasCachedDivisorData = true; } - // Store the accepted candidate at the next free slot so the front of the span becomes the filtered set. - chunkCandidates[filteredCount++] = chunkCandidates[i]; - } + MontgomeryDivisorData divisorData = cachedDivisorData; + ulong divisorCycle; - if (filteredCount == 0) - { - continue; - } - - Span filteredSpan = chunkCandidates[..filteredCount]; - Span cycleSpan = cycleBuffer.AsSpan(0, filteredCount); - cycleCache.GetCycleLengths(filteredSpan, cycleSpan); - - int admissibleCount = 0; - Span divisorDataSpan = divisorDataBuffer.AsSpan(0, filteredCount); - - for (int i = 0; i < filteredCount; i++) - { - if (cycleSpan[i] != prime) + if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles) { - continue; + divisorCycle = cycleCache.GetCycleLength(candidate); + } + else + { + factorCache ??= new Dictionary(8); + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL) + { + divisorCycle = cycleCache.GetCycleLength(candidate); + } } - ulong candidate = filteredSpan[i]; - // Reuse the filtered span's leading segment to accumulate the admissible divisors for Montgomery checks. - filteredSpan[admissibleCount] = candidate; - divisorDataSpan[admissibleCount] = MontgomeryDivisorData.FromModulus(candidate); - admissibleCount++; - } - - if (admissibleCount == 0) - { - continue; - } - - Span admissibleSpan = filteredSpan[..admissibleCount]; - Span admissibleDataSpan = divisorDataSpan[..admissibleCount]; - - for (int i = 0; i < admissibleCount; i++) - { - if (CheckDivisor(prime, prime, in admissibleDataSpan[i]) != 0) + if (divisorCycle == prime) { - lastProcessed = admissibleSpan[i]; processedAll = true; return true; } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } } + + divisor128 += step128 * (UInt128)chunkCount; + remainder10 = AdvanceRemainder(remainder10, step10, 10, chunkCount); + remainder8 = AdvanceRemainder(remainder8, step8, 8, chunkCount); + remainder5 = AdvanceRemainder(remainder5, step5, 5, chunkCount); + remainder3 = AdvanceRemainder(remainder3, step3, 3, chunkCount); + remainder7 = AdvanceRemainder(remainder7, step7, 7, chunkCount); + remainder11 = AdvanceRemainder(remainder11, step11, 11, chunkCount); } } finally { ArrayPool.Shared.Return(candidateBuffer, clearArray: false); - ArrayPool.Shared.Return(cycleBuffer, clearArray: false); - ArrayPool.Shared.Return(divisorDataBuffer, clearArray: false); - ArrayPool.Shared.Return(remainderBuffer, clearArray: false); + ArrayPool.Shared.Return(maskBuffer, clearArray: false); + candidateEvaluator.Dispose(); } - processedAll = divisor > limit; + processedAll = divisor128 > limit; return false; } + + private static int ComputeChunkCount(UInt128 currentDivisor, UInt128 limit, UInt128 step, int maxCount) + { + if (maxCount <= 0 || step == UInt128.Zero || currentDivisor > limit) + { + return 0; + } + + UInt128 distance = limit - currentDivisor; + UInt128 stepsAvailable = distance / step; + UInt128 total = stepsAvailable + UInt128.One; + if (total > (UInt128)maxCount) + { + total = (UInt128)maxCount; + } + + if (total > (UInt128)int.MaxValue) + { + return int.MaxValue; + } + + return (int)total; + } + + private static byte AdvanceRemainder(byte remainder, byte step, byte modulus, int count) + { + if (modulus == 0 || count <= 0) + { + return remainder; + } + + int modulusValue = modulus; + int remainderValue = remainder % modulusValue; + int stepValue = step % modulusValue; + + long delta = (long)stepValue * count; + int adjusted = (int)((remainderValue + delta) % modulusValue); + if (adjusted < 0) + { + adjusted += modulusValue; + } + + return (byte)adjusted; + } + private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) { ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData); diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs new file mode 100644 index 0000000..9a75acb --- /dev/null +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs @@ -0,0 +1,265 @@ +using System; +using System.Buffers; +using System.Runtime.InteropServices; +using ILGPU; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable +{ + private readonly GpuContextPool.GpuContextLease _lease; + private readonly Accelerator _accelerator; + private readonly Action, ArrayView> _kernel; + private MemoryBuffer1D _candidateBuffer; + private MemoryBuffer1D _maskBuffer; + private ulong[] _hostCandidates; + private byte[] _hostMask; + private int _capacity; + private bool _disposed; + + internal MersenneNumberDivisorCandidateGpuEvaluator(int batchSize) + { + _lease = GpuContextPool.RentPreferred(preferCpu: false); + _accelerator = _lease.Accelerator; + _kernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(EvaluateCandidatesKernel); + _capacity = Math.Max(1, batchSize); + _candidateBuffer = _accelerator.Allocate1D(_capacity); + _maskBuffer = _accelerator.Allocate1D(_capacity); + _hostCandidates = ArrayPool.Shared.Rent(_capacity); + _hostMask = ArrayPool.Shared.Rent(_capacity); + } + + public void Dispose() + { + if (_disposed) + { + return; + } + + _disposed = true; + _candidateBuffer.Dispose(); + _maskBuffer.Dispose(); + ArrayPool.Shared.Return(_hostCandidates, clearArray: false); + ArrayPool.Shared.Return(_hostMask, clearArray: false); + _lease.Dispose(); + } + + internal void EvaluateCandidates( + ulong startDivisor, + ulong step, + ulong limit, + byte remainder10, + byte remainder8, + byte remainder5, + byte remainder3, + byte remainder7, + byte remainder11, + byte step10, + byte step8, + byte step5, + byte step3, + byte step7, + byte step11, + bool lastIsSeven, + Span candidates, + Span mask) + { + if (_disposed) + { + throw new ObjectDisposedException(nameof(MersenneNumberDivisorCandidateGpuEvaluator)); + } + + int count = candidates.Length; + if (count == 0) + { + return; + } + + EnsureCapacity(count); + + Span hostCandidateSpan = _hostCandidates.AsSpan(0, count); + Span hostMaskSpan = _hostMask.AsSpan(0, count); + + ArrayView1D candidateView = _candidateBuffer.View.SubView(0, count); + ArrayView1D maskView = _maskBuffer.View.SubView(0, count); + + var args = new CandidateFilterArgs( + startDivisor, + step, + limit, + remainder10, + remainder8, + remainder5, + remainder3, + remainder7, + remainder11, + step10, + step8, + step5, + step3, + step7, + step11, + lastIsSeven ? (byte)1 : (byte)0); + + _kernel(count, args, candidateView, maskView); + + candidateView.CopyToCPU(ref MemoryMarshal.GetReference(hostCandidateSpan), count); + maskView.CopyToCPU(ref MemoryMarshal.GetReference(hostMaskSpan), count); + + hostCandidateSpan.CopyTo(candidates); + hostMaskSpan.CopyTo(mask); + } + + private void EnsureCapacity(int required) + { + if (required <= _capacity) + { + return; + } + + int newCapacity = _capacity == 0 ? 1 : _capacity; + while (newCapacity < required) + { + newCapacity <<= 1; + } + + _candidateBuffer.Dispose(); + _maskBuffer.Dispose(); + ArrayPool.Shared.Return(_hostCandidates, clearArray: false); + ArrayPool.Shared.Return(_hostMask, clearArray: false); + + _capacity = newCapacity; + _candidateBuffer = _accelerator.Allocate1D(_capacity); + _maskBuffer = _accelerator.Allocate1D(_capacity); + _hostCandidates = ArrayPool.Shared.Rent(_capacity); + _hostMask = ArrayPool.Shared.Rent(_capacity); + } + + private static void EvaluateCandidatesKernel( + Index1D index, + CandidateFilterArgs args, + ArrayView candidates, + ArrayView mask) + { + int idx = index; + if (idx >= mask.Length) + { + return; + } + + ulong offset = (ulong)idx; + ulong candidate = args.StartDivisor + args.Step * offset; + candidates[idx] = candidate; + + if (candidate > args.Limit || (args.Step != 0UL && candidate < args.StartDivisor)) + { + mask[idx] = 0; + return; + } + + int remainderValue10 = ComputeRemainder(args.Remainder10, args.Step10, idx, 10); + int remainderValue8 = ComputeRemainder(args.Remainder8, args.Step8, idx, 8); + int remainderValue5 = ComputeRemainder(args.Remainder5, args.Step5, idx, 5); + int remainderValue3 = ComputeRemainder(args.Remainder3, args.Step3, idx, 3); + int remainderValue7 = ComputeRemainder(args.Remainder7, args.Step7, idx, 7); + int remainderValue11 = ComputeRemainder(args.Remainder11, args.Step11, idx, 11); + + bool lastIsSeven = args.LastIsSevenFlag != 0; + bool admissible = lastIsSeven + ? (remainderValue10 == 3 || remainderValue10 == 7 || remainderValue10 == 9) + : (remainderValue10 == 1 || remainderValue10 == 3 || remainderValue10 == 9); + + if (!admissible) + { + mask[idx] = 0; + return; + } + + if ((remainderValue8 != 1 && remainderValue8 != 7) + || remainderValue3 == 0 + || remainderValue5 == 0 + || remainderValue7 == 0 + || remainderValue11 == 0) + { + mask[idx] = 0; + return; + } + + mask[idx] = 1; + } + + private static int ComputeRemainder(byte remainder, byte step, int index, int modulus) + { + if (modulus == 0) + { + return 0; + } + + long value = remainder; + value += (long)step * index; + int result = (int)(value % modulus); + if (result < 0) + { + result += modulus; + } + + return result; + } + + private readonly struct CandidateFilterArgs + { + public CandidateFilterArgs( + ulong startDivisor, + ulong step, + ulong limit, + byte remainder10, + byte remainder8, + byte remainder5, + byte remainder3, + byte remainder7, + byte remainder11, + byte step10, + byte step8, + byte step5, + byte step3, + byte step7, + byte step11, + byte lastIsSevenFlag) + { + StartDivisor = startDivisor; + Step = step; + Limit = limit; + Remainder10 = remainder10; + Remainder8 = remainder8; + Remainder5 = remainder5; + Remainder3 = remainder3; + Remainder7 = remainder7; + Remainder11 = remainder11; + Step10 = step10; + Step8 = step8; + Step5 = step5; + Step3 = step3; + Step7 = step7; + Step11 = step11; + LastIsSevenFlag = lastIsSevenFlag; + } + + public ulong StartDivisor { get; } + public ulong Step { get; } + public ulong Limit { get; } + public byte Remainder10 { get; } + public byte Remainder8 { get; } + public byte Remainder5 { get; } + public byte Remainder3 { get; } + public byte Remainder7 { get; } + public byte Remainder11 { get; } + public byte Step10 { get; } + public byte Step8 { get; } + public byte Step5 { get; } + public byte Step3 { get; } + public byte Step7 { get; } + public byte Step11 { get; } + public byte LastIsSevenFlag { get; } + } +}