diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index c061da1..bfe2808 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,5 +1,6 @@ using System.Collections.Generic; using System.Runtime.CompilerServices; +using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core.Cpu; @@ -12,6 +13,16 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private ulong _lastStatusDivisor; private int _batchSize = 1_024; + [ThreadStatic] + private static GpuUInt128WorkSet _divisorScanGpuWorkSet; + + private struct GpuUInt128WorkSet + { + public GpuUInt128 Step; + public GpuUInt128 Divisor; + public GpuUInt128 Limit; + } + public int BatchSize { get => _batchSize; @@ -103,41 +114,152 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv out bool processedAll, out ulong processedCount) { - processedAll = false; processedCount = 0UL; + processedAll = false; - // var enumerator = HeuristicPrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); - var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); + // The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs. + // Keeping the guard commented out documents the reasoning for benchmarks and tests. + // if (allowedMax < 3UL) + // { + // return false; + // } - // while (enumerator.TryGetNext(out HeuristicPrimeTester.HeuristicDivisorCandidate candidate)) - while (enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate)) + ref GpuUInt128WorkSet workSet = ref _divisorScanGpuWorkSet; + + ref GpuUInt128 step = ref workSet.Step; + step.High = 0UL; + step.Low = prime; + step.ShiftLeft(1); + + ref GpuUInt128 limit = ref workSet.Limit; + limit.High = 0UL; + limit.Low = allowedMax; + + ref GpuUInt128 divisor = ref workSet.Divisor; + divisor.High = step.High; + divisor.Low = step.Low; + divisor.Add(1UL); + if (divisor.CompareTo(limit) > 0) { - processedCount = enumerator.ProcessedCount; - - // HeuristicPrimeTester.HeuristicDivisorPreparation preparation = HeuristicPrimeTester.PrepareHeuristicDivisor(in candidate); - PrimeTester.HeuristicDivisorPreparation preparation = PrimeTester.PrepareHeuristicDivisor(in candidate); - // ulong divisorCycle = HeuristicPrimeTester.ResolveHeuristicCycleLength( - ulong divisorCycle = PrimeTester.ResolveHeuristicCycleLength( - prime, - in preparation, - out _, - out _, - out _); - - if (divisorCycle == prime) - { - processedAll = true; - return true; - } - - if (divisorCycle == 0UL) - { - Console.WriteLine($"Divisor cycle was not calculated for {prime}"); - } + processedAll = true; + divisor.High = 0UL; + divisor.Low = 0UL; + step.High = 0UL; + step.Low = 0UL; + limit.High = 0UL; + limit.Low = 0UL; + return false; } - processedCount = enumerator.ProcessedCount; - processedAll = enumerator.Exhausted; + // Intentionally recomputes factorizations without a per-thread cache. + // The previous factor cache recorded virtually no hits and only slowed down the scan. + DivisorCycleCache cycleCache = DivisorCycleCache.Shared; + + ulong stepHigh = step.High; + ulong stepLow = step.Low; + byte step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL); + byte step8 = (byte)(stepLow % 8UL); + byte step5 = (byte)(((stepHigh % 5UL) + (stepLow % 5UL)) % 5UL); + byte step3 = (byte)(((stepHigh % 3UL) + (stepLow % 3UL)) % 3UL); + byte step7 = (byte)((((stepHigh % 7UL) * 2UL) + (stepLow % 7UL)) % 7UL); + byte step11 = (byte)((((stepHigh % 11UL) * 5UL) + (stepLow % 11UL)) % 11UL); + + ulong divisorHigh = divisor.High; + ulong divisorLow = divisor.Low; + byte remainder10 = (byte)((((divisorHigh % 10UL) * 6UL) + (divisorLow % 10UL)) % 10UL); + byte remainder8 = (byte)(divisorLow % 8UL); + byte remainder5 = (byte)(((divisorHigh % 5UL) + (divisorLow % 5UL)) % 5UL); + byte remainder3 = (byte)(((divisorHigh % 3UL) + (divisorLow % 3UL)) % 3UL); + byte remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL); + byte remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL); + + // Keep the divisibility filters aligned with the divisor-cycle generator so the + // CPU path never requests cycles that were skipped during cache creation. + LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; + bool lastIsSeven = lastDigit == LastDigit.Seven; + const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); + const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); + + while (divisor.CompareTo(limit) <= 0) + { + ulong candidate = divisor.Low; + processedCount++; + + ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; + bool admissible = ((decimalMask >> remainder10) & 1) != 0; + + if (admissible + && (remainder8 == 1 || remainder8 == 7) + && remainder3 != 0 + && remainder5 != 0 + && remainder7 != 0 + && remainder11 != 0) + { + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); + ulong divisorCycle; + // Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs. + // if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles) + // { + // divisorCycle = cycleCache.GetCycleLength(candidate); + // } + // else + { + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( + candidate, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) + { + // Divisors produced by 2 * k * p + 1 always exceed PerfectNumberConstants.MaxQForDivisorCycles + // for the exponents scanned here, so skip the unused cache fallback and compute directly. + divisorCycle = MersenneDivisorCycles.CalculateCycleLength( + candidate, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + } + + if (divisorCycle == prime) + { + // A cycle equal to the tested exponent (which is prime in this path) guarantees that the candidate divides + // the corresponding Mersenne number because the order of 2 modulo the divisor is exactly p. + divisor.High = 0UL; + divisor.Low = 0UL; + step.High = 0UL; + step.Low = 0UL; + limit.High = 0UL; + limit.Low = 0UL; + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + divisor.Add(step); + remainder10 = AddMod(remainder10, step10, (byte)10); + remainder8 = AddMod(remainder8, step8, (byte)8); + remainder5 = AddMod(remainder5, step5, (byte)5); + remainder3 = AddMod(remainder3, step3, (byte)3); + remainder7 = AddMod(remainder7, step7, (byte)7); + remainder11 = AddMod(remainder11, step11, (byte)11); + } + + processedAll = divisor.CompareTo(limit) > 0; + divisor.High = 0UL; + divisor.Low = 0UL; + step.High = 0UL; + step.Low = 0UL; + limit.High = 0UL; + limit.Low = 0UL; return false; } @@ -147,6 +269,17 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return residue == 1UL ? (byte)1 : (byte)0; } + private static byte AddMod(byte value, byte delta, byte modulus) + { + int sum = value + delta; + if (sum >= modulus) + { + sum -= modulus; + } + + return (byte)sum; + } + private void UpdateStatusUnsafe(ulong processedCount) { if (processedCount == 0UL) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index e1bc1de..e1b1fe5 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -112,6 +112,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv resources.HitIndexBuffer, resources.Divisors, resources.Exponents, + resources.FilteredDivisors, resources.DivisorData, resources.Offsets, resources.Counts, @@ -176,6 +177,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv MemoryBuffer1D hitIndexBuffer, ulong[] divisors, ulong[] exponents, + ulong[] filteredDivisors, GpuDivisorPartialData[] divisorData, int[] offsets, int[] counts, @@ -186,40 +188,131 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { int batchCapacity = (int)divisorDataBuffer.Length; bool composite = false; + bool processedAll = false; processedCount = 0UL; lastProcessed = 0UL; - // var enumerator = HeuristicPrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); - var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); + int chunkCapacity = Math.Max(1, batchCapacity); - Span divisorSpan = divisors.AsSpan(); - Span exponentSpan = exponents.AsSpan(); - Span divisorDataSpan = divisorData.AsSpan(); - Span offsetSpan = offsets.AsSpan(); - Span countSpan = counts.AsSpan(); - Span cycleSpan = cycles.AsSpan(); - - while (!composite) + UInt128 twoP128 = (UInt128)prime << 1; + UInt128 allowedMax128 = allowedMax; + UInt128 firstDivisor128 = twoP128 + UInt128.One; + bool invalidStride = twoP128 == UInt128.Zero; + bool outOfRange = firstDivisor128 > allowedMax128; + UInt128 numerator = allowedMax128 - UInt128.One; + UInt128 maxK128 = invalidStride || outOfRange ? UInt128.Zero : numerator / twoP128; + bool hasCandidates = !invalidStride && !outOfRange && maxK128 != UInt128.Zero; + if (!hasCandidates) { + coveredRange = true; + return false; + } + + ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128; + ulong currentK = 1UL; + + byte step10 = (byte)(twoP128 % 10UL); + byte step8 = (byte)(twoP128 % 8UL); + byte step5 = (byte)(twoP128 % 5UL); + byte step3 = (byte)(twoP128 % 3UL); + byte step7 = (byte)(twoP128 % 7UL); + byte step11 = (byte)(twoP128 % 11UL); + + UInt128 currentDivisor128 = firstDivisor128; + byte remainder10 = (byte)((ulong)(currentDivisor128 % 10UL)); + byte remainder8 = (byte)((ulong)(currentDivisor128 & 7UL)); + byte remainder5 = (byte)((ulong)(currentDivisor128 % 5UL)); + byte remainder3 = (byte)((ulong)(currentDivisor128 % 3UL)); + byte remainder7 = (byte)((ulong)(currentDivisor128 % 7UL)); + byte remainder11 = (byte)((ulong)(currentDivisor128 % 11UL)); + LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; + + Span filteredStorage = filteredDivisors.AsSpan(); + Span divisorStorage = divisors.AsSpan(); + Span exponentStorage = exponents.AsSpan(); + Span divisorDataStorage = divisorData.AsSpan(); + Span offsetStorage = offsets.AsSpan(); + Span countStorage = counts.AsSpan(); + Span cycleStorage = cycles.AsSpan(); + while (currentK <= maxK && !composite) + { + int chunkCount = Math.Min(chunkCapacity, batchCapacity); + ulong remainingK = maxK - currentK + 1UL; + if ((ulong)chunkCount > remainingK) + { + chunkCount = (int)remainingK; + } + + if (chunkCount <= 0) + { + processedAll = true; + break; + } + + int filteredCount = 0; + UInt128 nextDivisor128 = currentDivisor128; + + byte localRemainder10 = remainder10; + byte localRemainder8 = remainder8; + byte localRemainder5 = remainder5; + byte localRemainder3 = remainder3; + byte localRemainder7 = remainder7; + byte localRemainder11 = remainder11; + + for (int i = 0; i < chunkCount; i++) + { + ulong candidate = (ulong)nextDivisor128; + + if (CandidatePassesHeuristics(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, lastDigit)) + { + filteredStorage[filteredCount++] = candidate; + } + + nextDivisor128 += twoP128; + localRemainder10 = AddMod(localRemainder10, step10, 10); + localRemainder8 = AddMod(localRemainder8, step8, 8); + localRemainder5 = AddMod(localRemainder5, step5, 5); + localRemainder3 = AddMod(localRemainder3, step3, 3); + localRemainder7 = AddMod(localRemainder7, step7, 7); + localRemainder11 = AddMod(localRemainder11, step11, 11); + } + + remainder10 = localRemainder10; + remainder8 = localRemainder8; + remainder5 = localRemainder5; + remainder3 = localRemainder3; + remainder7 = localRemainder7; + remainder11 = localRemainder11; + + processedCount += (ulong)chunkCount; + + UInt128 lastDivisor128 = nextDivisor128 - twoP128; + lastProcessed = (ulong)lastDivisor128; + currentDivisor128 = nextDivisor128; + + currentK += (ulong)chunkCount; + processedAll = currentK > maxK; + + if (filteredCount == 0) + { + continue; + } + + Span filteredDivisorsSpan = filteredStorage[..filteredCount]; + + Span divisorSpan = divisorStorage; + Span exponentSpan = exponentStorage; + Span divisorDataSpan = divisorDataStorage; + Span offsetSpan = offsetStorage; + Span countSpan = countStorage; + Span cycleSpan = cycleStorage; int admissibleCount = 0; int exponentIndex = 0; - - // while (admissibleCount < batchCapacity && enumerator.TryGetNext(out HeuristicPrimeTester.HeuristicDivisorCandidate candidate)) - while (admissibleCount < batchCapacity && enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate)) + for (int i = 0; i < filteredCount; i++) { - processedCount = enumerator.ProcessedCount; - ulong divisorValue = candidate.Value; - - // HeuristicPrimeTester.HeuristicDivisorPreparation preparation = HeuristicPrimeTester.PrepareHeuristicDivisor(in candidate); - PrimeTester.HeuristicDivisorPreparation preparation = PrimeTester.PrepareHeuristicDivisor(in candidate); - // ulong divisorCycle = HeuristicPrimeTester.ResolveHeuristicCycleLength( - ulong divisorCycle = PrimeTester.ResolveHeuristicCycleLength( - prime, - in preparation, - out _, - out _, - out _); - + ulong divisorValue = filteredDivisorsSpan[i]; + MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue); + ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData); if (divisorCycle != prime) { continue; @@ -236,16 +329,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv exponentIndex++; } - processedCount = enumerator.ProcessedCount; - lastProcessed = enumerator.LastDivisor; - if (admissibleCount == 0) { - if (enumerator.Exhausted) - { - break; - } - continue; } @@ -277,10 +362,74 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv lastProcessed = hitFound ? divisorSpan[hitIndex] : divisorSpan[lastIndex]; } - coveredRange = composite || enumerator.Exhausted; + coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128); return composite; } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ResolveDivisorCycle(ulong divisor, ulong prime, in MontgomeryDivisorData divisorData) + { + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(divisor, prime, divisorData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL) + { + return MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); + } + + return computedCycle; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool CandidatePassesHeuristics( + byte remainder10, + byte remainder8, + byte remainder3, + byte remainder5, + byte remainder7, + byte remainder11, + LastDigit lastDigit) + { + const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); + const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); + + ushort decimalMask = lastDigit == LastDigit.Seven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; + if (((decimalMask >> remainder10) & 1) == 0) + { + return false; + } + + if (remainder8 != 1 && remainder8 != 7) + { + return false; + } + + if (remainder3 == 0 || remainder5 == 0 || remainder7 == 0 || remainder11 == 0) + { + return false; + } + + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod(byte value, byte delta, byte modulus) + { + int result = value + delta; + if (result >= modulus) + { + result -= modulus; + if (result >= modulus) + { + result -= modulus; + } + } + + return (byte)result; + } + + + + + private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime) { // The --mersenne=bydivisor flow in EvenPerfectBitScanner only calls ConfigureFromMaxPrime with primes greater than 1, @@ -571,6 +720,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal ulong[] Exponents { get; private set; } = null!; + internal ulong[] FilteredDivisors { get; private set; } = null!; + internal GpuDivisorPartialData[] DivisorData { get; private set; } = null!; internal int[] Offsets { get; private set; } = null!; @@ -645,6 +796,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv Divisors = ulongPool.Rent(capacity); Exponents = ulongPool.Rent(capacity); + FilteredDivisors = ulongPool.Rent(capacity); DivisorData = partialPool.Rent(capacity); Offsets = intPool.Rent(capacity); Counts = intPool.Rent(capacity); @@ -683,12 +835,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayPool partialPool = ThreadStaticPools.GpuDivisorPartialDataPool; ulongPool.Return(Divisors, clearArray: false); ulongPool.Return(Exponents, clearArray: false); + ulongPool.Return(FilteredDivisors, clearArray: false); partialPool.Return(DivisorData, clearArray: false); intPool.Return(Offsets, clearArray: false); intPool.Return(Counts, clearArray: false); ulongPool.Return(Cycles, clearArray: false); Divisors = Array.Empty(); Exponents = Array.Empty(); + FilteredDivisors = Array.Empty(); DivisorData = Array.Empty(); Offsets = Array.Empty(); Counts = Array.Empty();