From 16181903a1e0d328d2f9d5b893d06382db7d004a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 5 Oct 2025 21:32:40 +0200 Subject: [PATCH] Streamline by-divisor cycle equality checks --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 115 +++------- ...MersenneNumberDivisorByDivisorGpuTester.cs | 215 +++--------------- PerfectNumbers.Core/MersenneDivisorCycles.cs | 94 ++++++++ 3 files changed, 149 insertions(+), 275 deletions(-) diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 476fb4a..8cec408 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,6 +1,5 @@ using System; using System.Collections.Concurrent; -using System.Collections.Generic; using System.Numerics; using PerfectNumbers.Core; @@ -174,70 +173,40 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } - Dictionary? factorCache = null; - ulong divisor = 3UL; - while (divisor <= allowedMax) + UInt128 step = (UInt128)prime << 1; + if (step == UInt128.Zero) { + processedAll = true; + return false; + } + + UInt128 limit = allowedMax; + UInt128 divisor = step + UInt128.One; + if (divisor > limit) + { + processedAll = true; + return false; + } + + while (divisor <= limit) + { + ulong candidate = (ulong)divisor; processedCount++; - lastProcessed = divisor; + lastProcessed = candidate; - MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); - // TODO: Hoist MontgomeryDivisorData acquisition into the divisor-cycle cache so we reuse the - // staged ProcessEightBitWindows-ready metadata instead of reloading the slower Montgomery - // structs for every divisor. - ulong divisorCycle; - if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) - { - divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); - } - else - { - factorCache ??= new Dictionary(8); - if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(divisor, prime, factorCache, out divisorCycle) || divisorCycle == 0UL) - { - divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); - } - } - if (divisorCycle != 0UL) - { - byte hit = CheckDivisor(prime, divisorCycle, divisorData); - if (hit != 0) - { - processedAll = true; - return true; - } - } - - if (divisor >= allowedMax - 1UL) + if (MersenneDivisorCycles.CycleEqualsExponentForMersenneCandidate(candidate, prime)) { processedAll = true; - break; + return true; } - if (divisor > ulong.MaxValue - 2UL) - { - break; - } - - divisor += 2UL; - // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the - // CPU by-divisor benchmarks so we advance directly to the next cached divisor candidate - // instead of testing every odd integer. + divisor += step; } - if (!processedAll) - { - processedAll = divisor > allowedMax; - } + processedAll = divisor > limit; return false; } - private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) - { - ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData); - return residue == 1UL ? (byte)1 : (byte)0; - } - private void UpdateStatusUnsafe(ulong lastProcessed, ulong processedCount) { if (processedCount == 0UL) @@ -316,44 +285,14 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return; } - if (divisorCycle == 0UL) + for (int i = 0; i < length; i++) { - divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); - if (divisorCycle == 0UL) - { - hits.Clear(); - return; - } - } - - MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); - - var exponentStepper = new ExponentRemainderStepper(divisorData); - if (!exponentStepper.IsValidModulus) - { - hits.Clear(); - return; - } - - var stepper = new CycleRemainderStepper(divisorCycle); - - ulong remainder = stepper.Initialize(primes[0]); - hits[0] = remainder == 0UL - ? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0) - : (byte)0; - - for (int i = 1; i < length; i++) - { - remainder = stepper.ComputeNext(primes[i]); - if (remainder != 0UL) - { - hits[i] = 0; - continue; - } - - hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0; + ulong prime = primes[i]; + bool divides = MersenneDivisorCycles.CycleEqualsExponentForMersenneCandidate(divisor, prime); + hits[i] = divides ? (byte)1 : (byte)0; } } + public void Dispose() { if (_disposed) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 17132e0..cb3eec2 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -1,7 +1,6 @@ using System; using System.Buffers; using System.Collections.Concurrent; -using System.Collections.Generic; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using System.Numerics; @@ -197,27 +196,19 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayView1D exponentView; ArrayView1D hitsView; - DivisorCycleCache cycleCache = DivisorCycleCache.Shared; - int cycleBatchCapacity = Math.Min(batchCapacity, cycleCache.PreferredBatchSize); - int maxKernelGroupSize = (int)accelerator.MaxNumThreadsPerGroup; - if (cycleBatchCapacity > maxKernelGroupSize) + int chunkCapacity = Math.Min(batchCapacity, (int)accelerator.MaxNumThreadsPerGroup); + if (chunkCapacity <= 0) { - cycleBatchCapacity = maxKernelGroupSize; + chunkCapacity = batchCapacity; } - if (cycleBatchCapacity <= 0) - { - cycleBatchCapacity = batchCapacity; - } - - ulong[] candidateDivisors = ArrayPool.Shared.Rent(cycleBatchCapacity); - ulong[] filteredDivisors = ArrayPool.Shared.Rent(cycleBatchCapacity); - ulong[] filteredCycles = ArrayPool.Shared.Rent(cycleBatchCapacity); - ulong[] divisorGaps = ArrayPool.Shared.Rent(cycleBatchCapacity); - byte[] remainder10Buffer = ArrayPool.Shared.Rent(cycleBatchCapacity); - byte[] remainder8Buffer = ArrayPool.Shared.Rent(cycleBatchCapacity); - byte[] remainder5Buffer = ArrayPool.Shared.Rent(cycleBatchCapacity); - byte[] remainder3Buffer = ArrayPool.Shared.Rent(cycleBatchCapacity); + ulong[] candidateDivisors = ArrayPool.Shared.Rent(chunkCapacity); + ulong[] filteredDivisors = ArrayPool.Shared.Rent(chunkCapacity); + ulong[] divisorGaps = ArrayPool.Shared.Rent(chunkCapacity); + byte[] remainder10Buffer = ArrayPool.Shared.Rent(chunkCapacity); + byte[] remainder8Buffer = ArrayPool.Shared.Rent(chunkCapacity); + byte[] remainder5Buffer = ArrayPool.Shared.Rent(chunkCapacity); + byte[] remainder3Buffer = ArrayPool.Shared.Rent(chunkCapacity); UInt128 twoP128 = (UInt128)prime << 1; if (twoP128 == UInt128.Zero) @@ -260,7 +251,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { while (currentK <= maxK) { - int chunkCount = Math.Min(cycleBatchCapacity, batchCapacity); + int chunkCount = Math.Min(chunkCapacity, batchCapacity); ulong remainingK = maxK - currentK + 1UL; if ((ulong)chunkCount > remainingK) { @@ -318,7 +309,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv int filteredCount = 0; Span filteredDivisorsSpan = filteredDivisors.AsSpan(0, chunkCount); - Span filteredCyclesSpan = filteredCycles.AsSpan(0, chunkCount); for (int i = 0; i < chunkCount; i++) { @@ -360,42 +350,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv continue; } - cycleCache.GetCycleLengths(filteredDivisorsSpan[..filteredCount], filteredCyclesSpan[..filteredCount]); - int batchSize = 0; for (int i = 0; i < filteredCount; i++) { ulong divisorValue = filteredDivisorsSpan[i]; - ulong cycleLength = filteredCyclesSpan[i]; - - bool includeDivisor = true; - ulong exponent = prime; - - if (cycleLength > 0UL) - { - ulong remainder = prime; - if (remainder >= cycleLength) - { - remainder %= cycleLength; - } - - if (remainder != 0UL) - { - includeDivisor = false; - } - else - { - exponent = 0UL; - } - } - - if (!includeDivisor) - { - continue; - } - divisors[batchSize] = divisorValue; - exponents[batchSize] = exponent; + exponents[batchSize] = prime; batchSize++; } @@ -450,7 +410,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { ArrayPool.Shared.Return(candidateDivisors, clearArray: false); ArrayPool.Shared.Return(filteredDivisors, clearArray: false); - ArrayPool.Shared.Return(filteredCycles, clearArray: false); ArrayPool.Shared.Return(divisorGaps, clearArray: false); ArrayPool.Shared.Return(remainder10Buffer, clearArray: false); ArrayPool.Shared.Return(remainder8Buffer, clearArray: false); @@ -559,8 +518,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private MemoryBuffer1D _exponentsBuffer; private MemoryBuffer1D _resultsBuffer; private ulong[] _hostBuffer; - private int[]? _positionBuffer; - private Dictionary? _factorCache; private int _capacity; private bool _disposed; @@ -574,13 +531,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv _exponentsBuffer = _accelerator.Allocate1D(_capacity); _resultsBuffer = _accelerator.Allocate1D(_capacity); _hostBuffer = ArrayPool.Shared.Rent(_capacity); - _positionBuffer = null; } internal void Reset() { _disposed = false; - _factorCache?.Clear(); } public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits) @@ -590,153 +545,41 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv throw new ObjectDisposedException(nameof(DivisorScanSession)); } - int primesLength = primes.Length; - if (primesLength == 0) + int length = primes.Length; + if (length == 0) { return; } MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); - var exponentStepper = new ExponentRemainderStepper(divisorData); - if (!exponentStepper.IsValidModulus) + if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL) { hits.Clear(); return; } - if (divisorCycle == 0UL) - { - Dictionary cache = _factorCache ??= new Dictionary(16); - if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(divisor, primes[0], cache, out divisorCycle) || divisorCycle == 0UL) - { - divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); - if (divisorCycle == 0UL) - { - hits.Clear(); - return; - } - } - } - - var cycleStepper = new CycleRemainderStepper(divisorCycle); - ulong initialCycleRemainder = cycleStepper.Initialize(primes[0]); - ProcessWithCycle(ref exponentStepper, ref cycleStepper, initialCycleRemainder, divisorCycle, divisorData, primes, hits); - } - - - private void ProcessWithCycle( - ref ExponentRemainderStepper exponentStepper, - ref CycleRemainderStepper cycleStepper, - ulong initialCycleRemainder, - ulong divisorCycle, - in MontgomeryDivisorData divisorData, - ReadOnlySpan primes, - Span hits) - { - Action, ArrayView> kernel = _kernel; ArrayView1D exponentsView = _exponentsBuffer.View; ArrayView1D resultsView = _resultsBuffer.View; - int capacity = _capacity; - Span hostSpan = _hostBuffer.AsSpan(0, capacity); - int[]? positionBuffer = _positionBuffer; - if (positionBuffer is null) - { - // Session construction already enforces the batch capacity, so we only - // need to rent the position buffer on the first cycle-enabled invocation. - positionBuffer = ArrayPool.Shared.Rent(capacity); - _positionBuffer = positionBuffer; - } + Span hostSpan = _hostBuffer.AsSpan(0, _capacity); - Span positionSpan = positionBuffer.AsSpan(0, capacity); - - int gpuBatchSize = _owner._gpuBatchSize; - int primesLength = primes.Length; int offset = 0; - bool cycleRemainderPending = true; - - while (offset < primesLength) + while (offset < length) { - int batchSize = Math.Min(gpuBatchSize, primesLength - offset); - ReadOnlySpan primesSlice = primes.Slice(offset, batchSize); - Span hitsSlice = hits.Slice(offset, batchSize); + int batchSize = Math.Min(_capacity, length - offset); + ReadOnlySpan batch = primes.Slice(offset, batchSize); + batch.CopyTo(hostSpan); - int computeCount = 0; - bool hasState = exponentStepper.HasState; - ulong previousExponent = exponentStepper.PreviousExponent; - bool useInitialRemainder = cycleRemainderPending; + ArrayView1D exponentSlice = exponentsView.SubView(0, batchSize); + exponentSlice.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), batchSize); + ArrayView1D resultSlice = resultsView.SubView(0, batchSize); + _kernel(batchSize, divisorData, exponentSlice, resultSlice); + resultSlice.CopyToCPU(ref MemoryMarshal.GetReference(hostSpan), batchSize); + + Span hitSlice = hits.Slice(offset, batchSize); for (int i = 0; i < batchSize; i++) { - hitsSlice[i] = 0; - ulong primeValue = primesSlice[i]; - - ulong remainder = useInitialRemainder - ? initialCycleRemainder - : cycleStepper.ComputeNext(primeValue); - useInitialRemainder = false; - - if (remainder != 0UL) - { - continue; - } - - ulong workValue; - if (!hasState) - { - workValue = primeValue; - hasState = true; - } - else - { - workValue = primeValue - previousExponent; - } - - previousExponent = primeValue; - - if (workValue >= divisorCycle) - { - // TODO: Replace this `%` with the divisor-cycle remainder helper so we reuse cached - // modulo results instead of recomputing them in the GPU staging loop. - workValue %= divisorCycle; - } - - hostSpan[computeCount] = workValue; - positionSpan[computeCount] = i; - computeCount++; - } - - cycleRemainderPending = useInitialRemainder; - - if (computeCount > 0) - { - Span exponentSlice = hostSpan[..computeCount]; - ArrayView1D exponentView = exponentsView.SubView(0, computeCount); - exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); - - ArrayView1D resultView = resultsView.SubView(0, computeCount); - kernel(computeCount, divisorData, exponentView, resultView); - resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); - - bool requiresInitialization = !exponentStepper.HasState; - for (int i = 0; i < computeCount; i++) - { - int position = positionSpan[i]; - ulong primeValue = primesSlice[position]; - ulong montgomeryResult = exponentSlice[i]; - bool isUnity; - - if (requiresInitialization) - { - exponentStepper.TryInitializeFromMontgomeryResult(primeValue, montgomeryResult, out isUnity); - requiresInitialization = false; - } - else - { - exponentStepper.TryAdvanceWithMontgomeryDelta(primeValue, montgomeryResult, out isUnity); - } - - hitsSlice[position] = isUnity ? (byte)1 : (byte)0; - } + hitSlice[i] = hostSpan[i] == 1UL ? (byte)1 : (byte)0; } offset += batchSize; @@ -751,11 +594,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } _disposed = true; - _factorCache?.Clear(); _owner._sessionPool.Add(this); } } - private static void ComputeRemainderDeltasKernel(Index1D index, ArrayView gaps, ArrayView deltas, byte modulus) { deltas[index] = (byte)(gaps[index] % modulus); diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 2e01d4f..c585dcc 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -54,6 +54,100 @@ public class MersenneDivisorCycles return snapshot; } + public static bool CycleEqualsExponent(ulong divisor, ulong exponent) + { + Dictionary? cache = null; + return CycleEqualsExponent(divisor, exponent, ref cache); + } + + public static bool CycleEqualsExponent(ulong divisor, ulong exponent, ref Dictionary? factorCache) + { + if (CycleEqualsExponentForMersenneCandidate(divisor, exponent)) + { + return true; + } + + 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; + } + } + } + + if (TryCalculateCycleLengthHeuristic(divisor, out ulong heuristicCycle) && heuristicCycle != 0UL) + { + return heuristicCycle == exponent; + } + + Dictionary? cache = factorCache; + if (cache is null) + { + cache = new Dictionary(8); + factorCache = cache; + } + + if (TryCalculateCycleLengthForExponent(divisor, exponent, cache, out ulong cycleLength) && cycleLength != 0UL) + { + return cycleLength == exponent; + } + + MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); + return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL; + } + + public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, ulong exponent) + { + if (!IsValidMersenneDivisorCandidate(divisor, exponent)) + { + return false; + } + + 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; + } + } + } + + MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); + return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL; + } + + private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) + { + if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) + { + return false; + } + + UInt128 step = (UInt128)exponent << 1; + if (step == UInt128.Zero) + { + return false; + } + + UInt128 adjusted = (UInt128)divisor - UInt128.One; + return adjusted % step == UInt128.Zero; + } + public static IEnumerable<(ulong divisor, ulong cycleLength)> EnumerateStream(Stream compressor) { // Binary pairs: (ulong divisor, ulong cycle)