diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 5254d7a..b356ebe 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -7,10 +7,7 @@ namespace PerfectNumbers.Core.Cpu; public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester { - private readonly object _sync = new(); - private ulong _divisorLimit; - private ulong _lastStatusDivisor; private int _batchSize = 1_024; [ThreadStatic] @@ -32,7 +29,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv public void ConfigureFromMaxPrime(ulong maxPrime) { _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); - _lastStatusDivisor = 0UL; } public ulong DivisorLimit @@ -60,22 +56,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv // return true; // } - ulong processedCount; bool processedAll; bool composite = CheckDivisors( prime, allowedMax, - out processedAll, - out processedCount); - - if (processedCount > 0UL) - { - lock (_sync) - { - UpdateStatusUnsafe(processedCount); - } - } + out processedAll); if (composite) { @@ -111,10 +97,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private bool CheckDivisors( ulong prime, ulong allowedMax, - out bool processedAll, - out ulong processedCount) + out bool processedAll) { - processedCount = 0UL; processedAll = false; // The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs. @@ -170,24 +154,20 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv // 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); + ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); while (divisor.CompareTo(limit) <= 0) { ulong candidate = divisor.Low; - processedCount++; - ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; - bool admissible = ((decimalMask >> remainder10) & 1) != 0; - - if (admissible + bool admissible = ((decimalMask >> remainder10) & 1) != 0 && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 - && remainder11 != 0) + && remainder11 != 0; + + if (admissible) { MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); ulong divisorCycle; @@ -234,7 +214,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv divisor.Add(step); remainder10 = AddMod(remainder10, step10, (byte)10); - remainder8 = AddMod(remainder8, step8, (byte)8); + remainder8 = AddMod8(remainder8, step8); remainder5 = AddMod(remainder5, step5, (byte)5); remainder3 = AddMod(remainder3, step3, (byte)3); remainder7 = AddMod(remainder7, step7, (byte)7); @@ -251,9 +231,11 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return residue == 1UL ? (byte)1 : (byte)0; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static byte AddMod(byte value, byte delta, byte modulus) { int sum = value + delta; + if (sum >= modulus) { sum -= modulus; @@ -262,24 +244,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return (byte)sum; } - private void UpdateStatusUnsafe(ulong processedCount) + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod8(byte value, byte delta) { - if (processedCount == 0UL) - { - return; - } - - ulong interval = PerfectNumberConstants.ConsoleInterval; - if (interval == 0UL) - { - _lastStatusDivisor = 0UL; - return; - } - - ulong total = _lastStatusDivisor + processedCount; - // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI - // status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly. - _lastStatusDivisor = total % interval; + return (byte)((value + delta) & 7); } private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime) diff --git a/PerfectNumbers.Core/DivisorGenerator.cs b/PerfectNumbers.Core/DivisorGenerator.cs index 0001267..44235fc 100644 --- a/PerfectNumbers.Core/DivisorGenerator.cs +++ b/PerfectNumbers.Core/DivisorGenerator.cs @@ -7,6 +7,12 @@ internal static class DivisorGenerator private const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); private const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ushort GetDecimalMask(LastDigit lastDigit) + { + return lastDigit == LastDigit.Seven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static bool IsValidDivisor( byte remainder10, @@ -15,9 +21,8 @@ internal static class DivisorGenerator byte remainder5, byte remainder7, byte remainder11, - LastDigit lastDigit) + ushort decimalMask) { - ushort decimalMask = lastDigit == LastDigit.Seven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; if (((decimalMask >> remainder10) & 1) == 0) { return false; diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index 327d25b..370b641 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -36,28 +36,30 @@ internal static class DivisorByDivisorKernels byte firstHitValue; ulong previousExponent = firstExponent; ulong cycleRemainder = 0UL; - if (cycleLength == 0UL) - { - firstHitValue = firstUnity ? (byte)1 : (byte)0; - } - else + bool hasCycle = cycleLength != 0UL; + if (hasCycle) { cycleRemainder = firstExponent % cycleLength; firstHitValue = cycleRemainder == 0UL && firstUnity ? (byte)1 : (byte)0; } + else + { + firstHitValue = firstUnity ? (byte)1 : (byte)0; + } hits[offset] = firstHitValue; byte anyHit = firstHitValue; - for (int i = 1; i < count; i++) + int endIndex = offset + count; + if (hasCycle) { - ulong exponent = exponents[offset + i]; - ulong delta = exponent - previousExponent; - previousExponent = exponent; - - if (cycleLength != 0UL) + GpuUInt128 remainderValue = new(cycleRemainder); + for (int exponentIndex = offset + 1; exponentIndex < endIndex; exponentIndex++) { - var remainderValue = new GpuUInt128(cycleRemainder); + ulong exponent = exponents[exponentIndex]; + ulong delta = exponent - previousExponent; + previousExponent = exponent; + remainderValue.AddMod(delta, cycleLength); cycleRemainder = remainderValue.Low; @@ -67,13 +69,17 @@ internal static class DivisorByDivisorKernels hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; } - hits[offset + i] = hit; + hits[exponentIndex] = hit; anyHit |= hit; } - else + } + else + { + for (int exponentIndex = offset + 1; exponentIndex < endIndex; exponentIndex++) { + ulong exponent = exponents[exponentIndex]; byte hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; - hits[offset + i] = hit; + hits[exponentIndex] = hit; anyHit |= hit; } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 8d83ee2..1a9a31f 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -192,7 +192,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv processedCount = 0UL; lastProcessed = 0UL; - int chunkCapacity = Math.Max(1, batchCapacity); + // EvenPerfectBitScanner always configures non-empty GPU divisor buffers, so the defensive guard stays + // commented out. Tests and benchmarks that rely on zero-capacity buffers must be updated instead of re-enabling + // this branch. + // if (batchCapacity <= 0) + // { + // throw new InvalidOperationException("Divisor buffers must be non-empty."); + // } + + int chunkCountBaseline = batchCapacity; UInt128 twoP128 = (UInt128)prime << 1; UInt128 allowedMax128 = allowedMax; @@ -209,7 +217,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128; - ulong currentK = 1UL; + ulong remainingCount = maxK; byte step10 = (byte)(twoP128 % 10UL); byte step8 = (byte)(twoP128 % 8UL); @@ -226,6 +234,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv byte remainder7 = (byte)((ulong)(currentDivisor128 % 7UL)); byte remainder11 = (byte)((ulong)(currentDivisor128 % 11UL)); LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; + ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); Span filteredStorage = filteredDivisors.AsSpan(); Span divisorStorage = divisors.AsSpan(); @@ -234,19 +243,32 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv 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; - } + Span divisorSpan = divisorStorage; + Span exponentSpan = exponentStorage; + Span divisorDataSpan = divisorDataStorage; + Span offsetSpan = offsetStorage; + Span countSpan = countStorage; + Span cycleSpan = cycleStorage; + ref GpuDivisorPartialData divisorDataRef = ref MemoryMarshal.GetReference(divisorDataSpan); + ref int offsetRef = ref MemoryMarshal.GetReference(offsetSpan); + ref int countRef = ref MemoryMarshal.GetReference(countSpan); + ref ulong cycleRef = ref MemoryMarshal.GetReference(cycleSpan); + ref ulong exponentRef = ref MemoryMarshal.GetReference(exponentSpan); - if (chunkCount <= 0) + ArrayView1D divisorDataView = divisorDataBuffer.View; + ArrayView1D offsetView = offsetBuffer.View; + ArrayView1D countView = countBuffer.View; + ArrayView1D cycleView = cycleBuffer.View; + ArrayView1D exponentViewDevice = exponentBuffer.View; + ArrayView1D hitsView = hitsBuffer.View; + ArrayView1D hitIndexView = hitIndexBuffer.View; + + while (remainingCount > 0UL && !composite) + { + int chunkCount = chunkCountBaseline; + if ((ulong)chunkCount > remainingCount) { - processedAll = true; - break; + chunkCount = (int)remainingCount; } int filteredCount = 0; @@ -263,14 +285,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { ulong candidate = (ulong)nextDivisor128; - if (DivisorGenerator.IsValidDivisor(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, lastDigit)) + if (DivisorGenerator.IsValidDivisor(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, decimalMask)) { filteredStorage[filteredCount++] = candidate; } nextDivisor128 += twoP128; localRemainder10 = AddMod(localRemainder10, step10, 10); - localRemainder8 = AddMod(localRemainder8, step8, 8); + localRemainder8 = AddMod8(localRemainder8, step8); localRemainder5 = AddMod(localRemainder5, step5, 5); localRemainder3 = AddMod(localRemainder3, step3, 3); localRemainder7 = AddMod(localRemainder7, step7, 7); @@ -289,28 +311,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv UInt128 lastDivisor128 = nextDivisor128 - twoP128; lastProcessed = (ulong)lastDivisor128; currentDivisor128 = nextDivisor128; - - currentK += (ulong)chunkCount; - processedAll = currentK > maxK; + remainingCount -= (ulong)chunkCount; 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; for (int i = 0; i < filteredCount; i++) { - ulong divisorValue = filteredDivisorsSpan[i]; + ulong divisorValue = filteredStorage[i]; MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue); ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData); if (divisorCycle != prime) @@ -319,14 +330,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } divisorSpan[admissibleCount] = divisorValue; - exponentSpan[exponentIndex] = prime; + exponentSpan[admissibleCount] = prime; divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue); - offsetSpan[admissibleCount] = exponentIndex; + offsetSpan[admissibleCount] = admissibleCount; countSpan[admissibleCount] = 1; cycleSpan[admissibleCount] = divisorCycle; admissibleCount++; - exponentIndex++; } if (admissibleCount == 0) @@ -334,24 +344,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv continue; } - ArrayView1D divisorView = divisorDataBuffer.View.SubView(0, admissibleCount); - ArrayView1D offsetView = offsetBuffer.View.SubView(0, admissibleCount); - ArrayView1D countView = countBuffer.View.SubView(0, admissibleCount); - ArrayView1D cycleView = cycleBuffer.View.SubView(0, admissibleCount); - ArrayView1D exponentView = exponentBuffer.View.SubView(0, exponentIndex); - ArrayView1D hitsView = hitsBuffer.View.SubView(0, exponentIndex); - ArrayView1D hitIndexView = hitIndexBuffer.View.SubView(0, 1); - - divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), admissibleCount); - offsetView.CopyFromCPU(ref MemoryMarshal.GetReference(offsetSpan), admissibleCount); - countView.CopyFromCPU(ref MemoryMarshal.GetReference(countSpan), admissibleCount); - cycleView.CopyFromCPU(ref MemoryMarshal.GetReference(cycleSpan), admissibleCount); - exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSpan), exponentIndex); + divisorDataView.CopyFromCPU(ref divisorDataRef, admissibleCount); + offsetView.CopyFromCPU(ref offsetRef, admissibleCount); + countView.CopyFromCPU(ref countRef, admissibleCount); + cycleView.CopyFromCPU(ref cycleRef, admissibleCount); + exponentViewDevice.CopyFromCPU(ref exponentRef, admissibleCount); int sentinel = int.MaxValue; hitIndexView.CopyFromCPU(ref sentinel, 1); - kernel(admissibleCount, divisorView, offsetView, countView, exponentView, cycleView, hitsView, hitIndexView); + kernel(admissibleCount, divisorDataView, offsetView, countView, exponentViewDevice, cycleView, hitsView, hitIndexView); int firstHit = sentinel; hitIndexView.CopyToCPU(ref firstHit, 1); @@ -362,6 +364,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv lastProcessed = hitFound ? divisorSpan[hitIndex] : divisorSpan[lastIndex]; } + processedAll = remainingCount == 0UL; coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128); return composite; } @@ -382,19 +385,23 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv 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; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod8(byte value, byte delta) + { + return (byte)((value + delta) & 7); + } + + @@ -557,27 +564,25 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayView1D hitView = _hitBuffer.View; GpuDivisorPartialData partialData = new GpuDivisorPartialData(divisor); - divisorView.SubView(0, 1).CopyFromCPU(ref partialData, 1); + divisorView.CopyFromCPU(ref partialData, 1); int offsetValue = 0; - offsetView.SubView(0, 1).CopyFromCPU(ref offsetValue, 1); + offsetView.CopyFromCPU(ref offsetValue, 1); int countValue = length; - countView.SubView(0, 1).CopyFromCPU(ref countValue, 1); + countView.CopyFromCPU(ref countValue, 1); ulong cycleValue = cycle; - cycleView.SubView(0, 1).CopyFromCPU(ref cycleValue, 1); + cycleView.CopyFromCPU(ref cycleValue, 1); Span hostSpan = _hostBuffer.AsSpan(0, length); primes.CopyTo(hostSpan); + ref ulong hostRef = ref MemoryMarshal.GetReference(hostSpan); - ArrayView1D exponentSlice = exponentView.SubView(0, length); - exponentSlice.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), length); - - ArrayView1D hitSliceView = hitView.SubView(0, length); + exponentView.CopyFromCPU(ref hostRef, length); int sentinel = int.MaxValue; - firstHitView.SubView(0, 1).CopyFromCPU(ref sentinel, 1); + firstHitView.CopyFromCPU(ref sentinel, 1); var kernel = _kernel; if (kernel is null) @@ -586,10 +591,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv throw new InvalidOperationException("GPU divisor kernel was not initialized."); } - kernel(1, divisorView, offsetView, countView, exponentSlice, cycleView, hitSliceView, firstHitView); + kernel(1, divisorView, offsetView, countView, exponentView, cycleView, hitView, firstHitView); Span hitSlice = hits.Slice(0, length); - hitSliceView.CopyToCPU(ref MemoryMarshal.GetReference(hitSlice), length); + ref byte hitRef = ref MemoryMarshal.GetReference(hitSlice); + hitView.CopyToCPU(ref hitRef, length); Monitor.Exit(_lease.ExecutionLock); } diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index 09c4435..131df86 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -93,9 +93,15 @@ public static class MersenneNumberDivisorByDivisorTester } else { - while (startIndex > 0 && candidates[startIndex - 1] >= startPrime) + while (startIndex > 0) { - startIndex--; + int previousIndex = startIndex - 1; + if (candidates[previousIndex] < startPrime) + { + break; + } + + startIndex = previousIndex; } } } diff --git a/docs/HeuristicIsPrimePlan.md b/docs/HeuristicIsPrimePlan.md index fd6671e..e2157c9 100644 --- a/docs/HeuristicIsPrimePlan.md +++ b/docs/HeuristicIsPrimePlan.md @@ -136,7 +136,7 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels, * [done] Kept GPU divisibility in standard modular arithmetic while reserving Montgomery transforms for CPU confirmation logic. * [done] When the temporary Group B gate is disabled, delegate GPU calls to the CPU fallback so Group A coverage remains in place while the accelerator path stays available for future re-enablement. 8. [done] **Refactor by-divisor CPU/GPU testers to reuse heuristics** - * [done] Drove `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` from `HeuristicPrimeTester.CreateMersenneDivisorEnumerator`, eliminating the bespoke residue loop while preserving the existing status counters. + * [done] Drove `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` from `HeuristicPrimeTester.CreateMersenneDivisorEnumerator`, eliminating the bespoke residue loop and retiring the legacy status counter. * [done] Streamed GPU batches from the same enumerator and Montgomery preparation data, removing the filtered-divisor scratch arrays while continuing to feed the existing kernel and hit bookkeeping. 9. [done] **Route legacy prime checks through heuristics** * [done] Introduce `HeuristicPrimeTester` so heuristic CPU/GPU checks live alongside the legacy `PrimeTester` implementations while callers opt in explicitly. diff --git a/todo_list.txt b/todo_list.txt index 72016e5..9ceb2df 100644 --- a/todo_list.txt +++ b/todo_list.txt @@ -41,7 +41,6 @@ Montgomery Primitives Cycle-Aware Modular Arithmetic Status & Progress Reporting -- PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs:296: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly. Cycle Remainder Steppers - PerfectNumbers.Core/CycleRemainderStepper.cs:28: Inline this reset at the call sites so the hot loops reuse struct reinitialization measured fastest in MersenneDivisorCycleLengthGpuBenchmarks, avoiding the extra method call when scanners need to restart stepping.