diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index f9ddf0d..64bb07e 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -26,20 +26,20 @@ internal static class Program private static int _primeCount; private static bool _primeFoundAfterInit; private static bool _useGcdFilter; - private static bool _useDivisor; - private static bool _useResidueMode; - private static UInt128 _divisor; - private static MersenneNumberDivisorGpuTester? _divisorTester; - private static ulong? _orderWarmupLimitOverride; - private static unsafe delegate* _transformP; - private static long _state; - private static bool _limitReached; - private static string? _resultsDir; - private static string? _resultsPrefix; - private static readonly Optimized PrimeIterator = new(); + private static bool _useDivisor; + private static bool _useResidueMode; + private static UInt128 _divisor; + private static MersenneNumberDivisorGpuTester? _divisorTester; + private static ulong? _orderWarmupLimitOverride; + private static unsafe delegate* _transformP; + private static long _state; + private static bool _limitReached; + private static string? _resultsDir; + private static string? _resultsPrefix; + private static readonly Optimized PrimeIterator = new(); - [ThreadStatic] - private static bool _lastCompositeP; + [ThreadStatic] + private static bool _lastCompositeP; // RLE/bit-pattern filtering options private static string? _rleBlacklistPath; @@ -414,20 +414,20 @@ internal static class Program Console.WriteLine("Divisor cycles are ready"); - _useDivisor = useDivisor; - _useResidueMode = useResidue; - _divisor = divisor; - if (useBitTransform) - { - _transformP = &TransformPBit; - } + _useDivisor = useDivisor; + _useResidueMode = useResidue; + _divisor = divisor; + if (useBitTransform) + { + _transformP = &TransformPBit; + } else if (useResidue) { - _transformP = &TransformPAdd; + _transformP = &TransformPAdd; } else { - _transformP = &TransformPAddPrimes; + _transformP = &TransformPAddPrimes; } PResidue = new ThreadLocal(() => new ModResidueTracker(ResidueModel.Identity, initialNumber: currentP, initialized: true), trackAllValues: true); PrimeTesters = new ThreadLocal(() => new PrimeTester(), trackAllValues: true); @@ -435,20 +435,20 @@ internal static class Program // Initialize per-thread p residue tracker (Identity model) at currentP if (!useDivisor) { - MersenneTesters = new ThreadLocal(() => - { - var tester = new MersenneNumberTester( - useIncremental: !useLucas, - useOrderCache: false, - kernelType: kernelType, - useModuloWorkaround: useModuloWorkaround, - useOrder: useOrder, - useGpuLucas: mersenneOnGpu, - useGpuScan: mersenneOnGpu, - useGpuOrder: orderOnGpu, - useResidue: useResidue, - maxK: residueKMax, - residueDivisorSets: divisorCyclesSearchLimit); + MersenneTesters = new ThreadLocal(() => + { + var tester = new MersenneNumberTester( + useIncremental: !useLucas, + useOrderCache: false, + kernelType: kernelType, + useModuloWorkaround: useModuloWorkaround, + useOrder: useOrder, + useGpuLucas: mersenneOnGpu, + useGpuScan: mersenneOnGpu, + useGpuOrder: orderOnGpu, + useResidue: useResidue, + maxK: residueKMax, + residueDivisorSets: divisorCyclesSearchLimit); if (!useLucas) { tester.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL); @@ -732,7 +732,7 @@ internal static class Program Console.WriteLine(" --divisor-cycles-device=cpu|gpu device for cycles generation (default gpu)"); Console.WriteLine(" --divisor-cycles-batch-size= batch size for cycles generation (default 512)"); Console.WriteLine(" --divisor-cycles-continue continue divisor cycles generation"); - Console.WriteLine(" --divisor-cycles-limit= cycle search iterations when --mersenne=divisor or residue"); + Console.WriteLine(" --divisor-cycles-limit= cycle search iterations when --mersenne=divisor or residue"); Console.WriteLine(" --use-order test primality via q order"); Console.WriteLine(" --workaround-mod avoid '%' operator on the GPU"); @@ -790,30 +790,30 @@ internal static class Program } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void PrintResult(ulong currentP, bool searchedMersenne, bool detailedCheck, bool isPerfect) - { - if (isPerfect) - { - int newCount = Interlocked.Increment(ref _primeCount); - if (newCount >= 2) - { - Volatile.Write(ref _primeFoundAfterInit, true); - } - } + private static void PrintResult(ulong currentP, bool searchedMersenne, bool detailedCheck, bool isPerfect) + { + if (isPerfect) + { + int newCount = Interlocked.Increment(ref _primeCount); + if (newCount >= 2) + { + Volatile.Write(ref _primeFoundAfterInit, true); + } + } - if (_useResidueMode && _lastCompositeP) - { - return; - } + if (_useResidueMode && _lastCompositeP) + { + return; + } - bool primeFlag = false, printToConsole = false; + bool primeFlag = false, printToConsole = false; - StringBuilder localBuilder = StringBuilderPool.Rent(); - _ = localBuilder - .Append(currentP).Append(',') - .Append(searchedMersenne).Append(',') - .Append(detailedCheck).Append(',') - .Append(isPerfect).Append('\n'); + StringBuilder localBuilder = StringBuilderPool.Rent(); + _ = localBuilder + .Append(currentP).Append(',') + .Append(searchedMersenne).Append(',') + .Append(detailedCheck).Append(',') + .Append(isPerfect).Append('\n'); lock (Sync) { @@ -1011,20 +1011,20 @@ internal static class Program [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static bool IsEvenPerfectCandidate(ulong p, ulong divisorCyclesSearchLimit) => IsEvenPerfectCandidate(p, divisorCyclesSearchLimit, out _, out _); - internal static bool IsEvenPerfectCandidate(ulong p, ulong divisorCyclesSearchLimit, out bool searchedMersenne, out bool detailedCheck) - { - searchedMersenne = false; - detailedCheck = false; - _lastCompositeP = false; + internal static bool IsEvenPerfectCandidate(ulong p, ulong divisorCyclesSearchLimit, out bool searchedMersenne, out bool detailedCheck) + { + searchedMersenne = false; + detailedCheck = false; + _lastCompositeP = false; - if (_useGcdFilter && IsCompositeByGcd(p)) - { - _lastCompositeP = true; - return false; - } + if (_useGcdFilter && IsCompositeByGcd(p)) + { + _lastCompositeP = true; + return false; + } - // Optional: RLE blacklist and binary-threshold filters on p (safe only when configured) - if (p <= _rleHardMaxP) + // Optional: RLE blacklist and binary-threshold filters on p (safe only when configured) + if (p <= _rleHardMaxP) { if (!_rleOnlyLast7 || (p % 10UL) == 7UL) { @@ -1058,36 +1058,36 @@ internal static class Program } } - // Fast residue-based composite check for p using small primes - if (!_useDivisor && IsCompositeByResidues(p)) - { - searchedMersenne = false; - detailedCheck = false; - _lastCompositeP = true; - return false; - } + // Fast residue-based composite check for p using small primes + if (!_useDivisor && IsCompositeByResidues(p)) + { + searchedMersenne = false; + detailedCheck = false; + _lastCompositeP = true; + return false; + } - // If primes-device=gpu, route p primality through GPU-assisted sieve (optional MR later) - // TODO: Add deterministic Miller–Rabin rounds in GPU path for 64-bit range. - bool isPrimeP = GpuContextPool.ForceCpu - ? PrimeTesters.Value!.IsPrime(p, CancellationToken.None) - : PrimeTesters.Value!.IsPrimeGpu(p, CancellationToken.None); - if (!isPrimeP) - { - _lastCompositeP = true; - return false; - } + // If primes-device=gpu, route p primality through GPU-assisted sieve (optional MR later) + // TODO: Add deterministic Miller–Rabin rounds in GPU path for 64-bit range. + bool isPrimeP = GpuContextPool.ForceCpu + ? PrimeTesters.Value!.IsPrime(p, CancellationToken.None) + : PrimeTesters.Value!.IsPrimeGpu(p, CancellationToken.None); + if (!isPrimeP) + { + _lastCompositeP = true; + 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); + } - bool mersennePrime = MersenneTesters.Value!.IsMersennePrime(p, out bool residueExhausted); - detailedCheck = residueExhausted; - return mersennePrime; - } + bool mersennePrime = MersenneTesters.Value!.IsMersennePrime(p, out bool residueExhausted); + detailedCheck = residueExhausted; + return mersennePrime; + } // Use ModResidueTracker with a small set of primes to pre-filter composite p. private static bool IsCompositeByResidues(ulong p) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs index 598de60..32f0ee4 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs @@ -11,103 +11,103 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) // candidates without launching heavy kernels. Consider a device-visible constant buffer per-accelerator // via GpuKernelPool to avoid host round-trips. - // GPU residue variant: check 2^p % q == 1 for q = 2*p*k + 1. - public void Scan( - ulong exponent, - UInt128 twoP, - bool lastIsSeven, - UInt128 perSetLimit, - ulong setCount, - UInt128 overallLimit, - ref bool isPrime, - ref bool divisorsExhausted) - { - if (setCount == 0UL || perSetLimit == UInt128.Zero || overallLimit == UInt128.Zero) - { - return; - } + // GPU residue variant: check 2^p % q == 1 for q = 2*p*k + 1. + public void Scan( + ulong exponent, + UInt128 twoP, + bool lastIsSeven, + UInt128 perSetLimit, + ulong setCount, + UInt128 overallLimit, + ref bool isPrime, + ref bool divisorsExhausted) + { + if (setCount == 0UL || perSetLimit == UInt128.Zero || overallLimit == UInt128.Zero) + { + return; + } - var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder); - var accelerator = gpuLease.Accelerator; - // Ensure device has small cycles and primes tables for in-kernel lookup - var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator); - ResiduePrimeViews primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator); - int batchSize = GpuConstants.ScanBatchSize; - byte last = lastIsSeven ? (byte)1 : (byte)0; - var kernel = gpuLease.Pow2ModKernel; - ulong step10 = ((exponent % 10UL) << 1) % 10UL; - ulong step8 = ((exponent & 7UL) << 1) & 7UL; - ulong step3 = ((exponent % 3UL) << 1) % 3UL; - ulong step5 = ((exponent % 5UL) << 1) % 5UL; - GpuUInt128 twoPGpu = (GpuUInt128)twoP; + var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder); + var accelerator = gpuLease.Accelerator; + // Ensure device has small cycles and primes tables for in-kernel lookup + var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator); + ResiduePrimeViews primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator); + int batchSize = GpuConstants.ScanBatchSize; + byte last = lastIsSeven ? (byte)1 : (byte)0; + var kernel = gpuLease.Pow2ModKernel; + ulong step10 = ((exponent % 10UL) << 1) % 10UL; + ulong step8 = ((exponent & 7UL) << 1) & 7UL; + ulong step3 = ((exponent % 3UL) << 1) % 3UL; + ulong step5 = ((exponent % 5UL) << 1) % 5UL; + GpuUInt128 twoPGpu = (GpuUInt128)twoP; - var orderBuffer = accelerator.Allocate1D(batchSize); - ulong[] orderArray = ArrayPool.Shared.Rent(batchSize); - UInt128 batchSize128 = (UInt128)batchSize; - UInt128 limitInclusive = overallLimit + UInt128.One; - Span orders = orderArray.AsSpan(0, batchSize); - try - { - for (ulong setIndex = 0; setIndex < setCount && Volatile.Read(ref isPrime); setIndex++) - { - UInt128 setOffset = perSetLimit * (UInt128)setIndex; - UInt128 setStart = setOffset + UInt128.One; - if (setStart >= limitInclusive) - { - break; - } + var orderBuffer = accelerator.Allocate1D(batchSize); + ulong[] orderArray = ArrayPool.Shared.Rent(batchSize); + UInt128 batchSize128 = (UInt128)batchSize; + UInt128 limitInclusive = overallLimit + UInt128.One; + Span orders = orderArray.AsSpan(0, batchSize); + try + { + for (ulong setIndex = 0; setIndex < setCount && Volatile.Read(ref isPrime); setIndex++) + { + UInt128 setOffset = perSetLimit * (UInt128)setIndex; + UInt128 setStart = setOffset + UInt128.One; + if (setStart >= limitInclusive) + { + break; + } - UInt128 setLimitExclusive = setOffset + perSetLimit + UInt128.One; - if (setLimitExclusive > limitInclusive) - { - setLimitExclusive = limitInclusive; - } + UInt128 setLimitExclusive = setOffset + perSetLimit + UInt128.One; + if (setLimitExclusive > limitInclusive) + { + setLimitExclusive = limitInclusive; + } - UInt128 kStart = setStart; - while (kStart < setLimitExclusive && Volatile.Read(ref isPrime)) - { - UInt128 remaining = setLimitExclusive - kStart; - int currentSize = remaining > batchSize128 ? batchSize : (int)remaining; - if (currentSize <= 0) - { - break; - } - orders = orderArray.AsSpan(0, currentSize); + UInt128 kStart = setStart; + while (kStart < setLimitExclusive && Volatile.Read(ref isPrime)) + { + UInt128 remaining = setLimitExclusive - kStart; + int currentSize = remaining > batchSize128 ? batchSize : (int)remaining; + if (currentSize <= 0) + { + break; + } + orders = orderArray.AsSpan(0, currentSize); - UInt128 q = twoP * kStart + UInt128.One; - q.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3); - var ra = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5); + UInt128 q = twoP * kStart + UInt128.One; + q.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3); + var ra = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5); - kernel(currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, 0UL, - ra, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven, primeViews.LastOnePow2, primeViews.LastSevenPow2); + kernel(currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, 0UL, + ra, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven, primeViews.LastOnePow2, primeViews.LastSevenPow2); - accelerator.Synchronize(); - orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize); - if (!Volatile.Read(ref isPrime)) - { - break; - } + accelerator.Synchronize(); + orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize); + if (!Volatile.Read(ref isPrime)) + { + break; + } - for (int i = 0; i < currentSize; i++) - { - if (orders[i] != 0UL) - { - Volatile.Write(ref isPrime, false); - break; - } - } + for (int i = 0; i < currentSize; i++) + { + if (orders[i] != 0UL) + { + Volatile.Write(ref isPrime, false); + break; + } + } - kStart += batchSize128; - } - } - } - finally - { - ArrayPool.Shared.Return(orderArray); - orderBuffer.Dispose(); - gpuLease.Dispose(); - } + kStart += batchSize128; + } + } + } + finally + { + ArrayPool.Shared.Return(orderArray); + orderBuffer.Dispose(); + gpuLease.Dispose(); + } - divisorsExhausted = true; - } + divisorsExhausted = true; + } } diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index 92fafa9..eb53ca1 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 ulong ExtraDivisorCycleSearchLimit = 64UL; + public const ulong ExtraDivisorCycleSearchLimit = 64UL; }