diff --git a/EvenPerfectBitScanner.Benchmarks/Mod6ULongBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Mod6ULongBenchmarks.cs index 0b5c8a0..f75a5e2 100644 --- a/EvenPerfectBitScanner.Benchmarks/Mod6ULongBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Mod6ULongBenchmarks.cs @@ -1,6 +1,7 @@ using System.Runtime.CompilerServices; using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core; namespace EvenPerfectBitScanner.Benchmarks; diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 0d1209a..c411c60 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -14,1624 +14,1594 @@ namespace EvenPerfectBitScanner; internal static class Program { - private static ThreadLocal PrimeTesters = null!; - private static ThreadLocal MersenneTesters = null!; - private static ThreadLocal PResidue = null!; // p mod d tracker (per-thread) - private const ulong InitialP = PerfectNumberConstants.BiggestKnownEvenPerfectP; - private const int WriteBatchSize = 100; - private const string DefaultCyclesPath = "divisor_cycles.bin"; - private static string ResultsFileName = "even_perfect_bit_scan_results.csv"; - private const string PrimeFoundSuffix = " (FOUND VALID CANDIDATES)"; - private static StringBuilder? _outputBuilder; - private static readonly object Sync = new(); - private static readonly object FileWriteSync = new(); - private static int _consoleCounter; - private static int _writeIndex; - private static int _primeCount; - private static bool _primeFoundAfterInit; - private static bool _useGcdFilter; - private static bool _useDivisor; - private static bool _useResidueMode; - private static bool _useByDivisorMode; - private static UInt128 _divisor; - private static MersenneNumberDivisorGpuTester? _divisorTester; - private static IMersenneNumberDivisorByDivisorTester? _byDivisorTester; - private static Dictionary? _byDivisorPreviousResults; - 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 ulong _byDivisorStartPrime; + private static ThreadLocal PrimeTesters = null!; + private static ThreadLocal MersenneTesters = null!; + private static ThreadLocal PResidue = null!; // p mod d tracker (per-thread) + private const ulong InitialP = PerfectNumberConstants.BiggestKnownEvenPerfectP; + private const int WriteBatchSize = 100; + private const string DefaultCyclesPath = "divisor_cycles.bin"; + private static string ResultsFileName = "even_perfect_bit_scan_results.csv"; + private const string PrimeFoundSuffix = " (FOUND VALID CANDIDATES)"; + private static StringBuilder? _outputBuilder; + private static readonly object Sync = new(); + private static readonly object FileWriteSync = new(); + private static int _consoleCounter; + private static int _writeIndex; + private static int _primeCount; + private static bool _primeFoundAfterInit; + private static bool _useGcdFilter; + private static bool _useDivisor; + private static bool _useResidueMode; + private static bool _useByDivisorMode; + private static MersenneNumberDivisorGpuTester? _divisorTester; + private static IMersenneNumberDivisorByDivisorTester? _byDivisorTester; + private static Dictionary? _byDivisorPreviousResults; + 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 ulong _byDivisorStartPrime; - [ThreadStatic] - private static bool _lastCompositeP; + [ThreadStatic] + private static bool _lastCompositeP; - // RLE/bit-pattern filtering options - private static string? _rleBlacklistPath; - private static ulong _rleHardMaxP = ulong.MaxValue; - private static bool _rleOnlyLast7 = true; - private static int _writeBatchSize = WriteBatchSize; - private static double _zeroFracHard = -1.0; // disabled when < 0 - private static double _zeroFracConj = -1.0; // disabled when < 0 - private static int _maxZeroConj = -1; // disabled when < 0 + // RLE/bit-pattern filtering options + private static string? _rleBlacklistPath; + private static ulong _rleHardMaxP = ulong.MaxValue; + private static bool _rleOnlyLast7 = true; + private static int _writeBatchSize = WriteBatchSize; + private static double _zeroFracHard = -1.0; // disabled when < 0 + private static double _zeroFracConj = -1.0; // disabled when < 0 + private static int _maxZeroConj = -1; // disabled when < 0 - private static unsafe void Main(string[] args) - { - ulong currentP = InitialP; - int threadCount = Environment.ProcessorCount; - int blockSize = 1; - int gpuPrimeThreads = 1; - int gpuPrimeBatch = 262_144; - GpuKernelType kernelType = GpuKernelType.Incremental; - bool useModuloWorkaround = false; // TODO: Remove once the runtime defaults to the ImmediateModulo path measured fastest in GpuUInt128NativeModuloBenchmarks. - // removed: useModAutomaton - bool useOrder = false; - bool showHelp = false; - bool useBitTransform = false; - bool useLucas = false; - bool useResidue = false; // M_p test via residue divisors (replaces LL/incremental) - bool useDivisor = false; // M_p divisibility by specific divisor - bool useByDivisor = false; // Iterative divisor scan across primes - bool useDivisorCycles = false; // TODO: Default this flag to true and remove the opt-out once divisor cycle lookups are mandatory per the optimized by-divisor benchmarks. - UInt128 divisor = UInt128.Zero; - // Device routing - bool useGpuCycles = true; - bool mersenneOnGpu = true; // controls Lucas/incremental/pow2mod device - bool orderOnGpu = true; // controls order computations device - int scanBatchSize = 2_097_152, sliceSize = 32; - ulong residueKMax = 5_000_000UL; - string filterFile = string.Empty; - string cyclesPath = DefaultCyclesPath; - int cyclesBatchSize = 512; - bool continueCyclesGeneration = false; - ulong divisorCyclesSearchLimit = PerfectNumberConstants.ExtraDivisorCycleSearchLimit; - int argIndex = 0; - string arg = string.Empty; - ReadOnlySpan mersenneOption = default; - ulong parsedLimit = 0UL; - ulong parsedResidueMax = 0UL; - bool startPrimeProvided = false; + private static unsafe void Main(string[] args) + { + ulong currentP = InitialP; + int threadCount = Environment.ProcessorCount; + int blockSize = 1; + int gpuPrimeThreads = 1; + int gpuPrimeBatch = 262_144; + GpuKernelType kernelType = GpuKernelType.Incremental; + bool useModuloWorkaround = false; // TODO: Remove once the runtime defaults to the ImmediateModulo path measured fastest in GpuUInt128NativeModuloBenchmarks. + // removed: useModAutomaton + bool useOrder = false; + bool showHelp = false; + bool useBitTransform = false; + bool useLucas = false; + bool useResidue = false; // M_p test via residue divisors (replaces LL/incremental) + bool useDivisor = false; // M_p divisibility by specific divisor + bool useByDivisor = false; // Iterative divisor scan across primes UInt128 divisor = UInt128.Zero; + // Device routing + bool useGpuCycles = true; + bool mersenneOnGpu = true; // controls Lucas/incremental/pow2mod device + bool orderOnGpu = true; // controls order computations device + int scanBatchSize = 2_097_152, sliceSize = 32; + ulong residueKMax = 5_000_000UL; + string filterFile = string.Empty; + string cyclesPath = DefaultCyclesPath; + int cyclesBatchSize = 512; + bool continueCyclesGeneration = false; + ulong divisorCyclesSearchLimit = PerfectNumberConstants.ExtraDivisorCycleSearchLimit; + int argIndex = 0; + string arg = string.Empty; + ReadOnlySpan mersenneOption = default; + ulong parsedLimit = 0UL; + ulong parsedResidueMax = 0UL; + ulong remainder = 0UL; + bool startPrimeProvided = false; - // NTT backend selection (GPU): reference vs staged - for (; argIndex < args.Length; argIndex++) - { - arg = args[argIndex]; - if (arg.Equals("--?", StringComparison.OrdinalIgnoreCase) || arg.Equals("-?", StringComparison.OrdinalIgnoreCase) || arg.Equals("--help", StringComparison.OrdinalIgnoreCase) || arg.Equals("-help", StringComparison.OrdinalIgnoreCase) || arg.Equals("/?", StringComparison.OrdinalIgnoreCase)) - { - showHelp = true; - break; - } - else if (arg.StartsWith("--prime=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Replace this ulong.Parse with the Utf8Parser-based span helper once CLI option parsing - // adopts the zero-allocation fast path proven fastest in the parser benchmarks. - currentP = ulong.Parse(arg.AsSpan(arg.IndexOf('=') + 1)); - startPrimeProvided = true; - } - else if (arg.Equals("--increment=bit", StringComparison.OrdinalIgnoreCase)) - { - useBitTransform = true; - } - else if (arg.Equals("--increment=add", StringComparison.OrdinalIgnoreCase)) - { - useBitTransform = false; - } - else if (arg.StartsWith("--threads=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Swap int.Parse for the span-based Utf8Parser fast path to avoid transient strings when - // processing CLI numeric options. - threadCount = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); - } - else if (arg.StartsWith("--mersenne=", StringComparison.OrdinalIgnoreCase)) - { - mersenneOption = arg.AsSpan(arg.IndexOf('=') + 1); - if (mersenneOption.Equals("pow2mod", StringComparison.OrdinalIgnoreCase)) - { - kernelType = GpuKernelType.Pow2Mod; // TODO: Switch this mode to the ProcessEightBitWindows kernel once Pow2Minus1Mod adopts the benchmarked faster windowed ladder. - } - else if (mersenneOption.Equals("lucas", StringComparison.OrdinalIgnoreCase)) - { - useLucas = true; - } - else if (mersenneOption.Equals("residue", StringComparison.OrdinalIgnoreCase)) - { - useResidue = true; - useLucas = false; - } - else if (mersenneOption.Equals("divisor", StringComparison.OrdinalIgnoreCase)) - { - useDivisor = true; - useLucas = false; - useResidue = false; - } - else if (mersenneOption.Equals("bydivisor", StringComparison.OrdinalIgnoreCase)) - { - useByDivisor = true; - useLucas = false; - useResidue = false; - useDivisor = false; - } - else - { - kernelType = GpuKernelType.Incremental; - } - } - else if (arg.StartsWith("--divisor=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Replace UInt128.Parse with the benchmarked fast-path parser once we expose the - // TryParse-style span helper for UInt128 command-line arguments. - divisor = UInt128.Parse(arg.AsSpan(arg.IndexOf('=') + 1)); - } - else if (arg.StartsWith("--use-divisor-cycles=", StringComparison.OrdinalIgnoreCase)) - { - if (bool.TryParse(arg.AsSpan(arg.IndexOf('=') + 1), out bool parsedBool)) - { - useDivisorCycles = parsedBool; - } - } - else if (arg.StartsWith("--divisor-cycles-limit=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Replace ulong.TryParse with the Utf8Parser fast-path once CLI parsing utilities expose - // span-based helpers for optional arguments. - if (ulong.TryParse(arg.AsSpan(arg.IndexOf('=') + 1), out parsedLimit)) - { - divisorCyclesSearchLimit = parsedLimit; - } - } - else if (arg.StartsWith("--residue-max-k=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Swap this TryParse for the zero-allocation Utf8Parser helper to keep residue CLI option - // parsing in the fast path. - if (ulong.TryParse(arg.AsSpan(arg.IndexOf('=') + 1), out parsedResidueMax)) - { - residueKMax = parsedResidueMax; - } - } - // Replaces --lucas=cpu|gpu; controls device for Lucas and for - // incremental/pow2mod scanning - else if (arg.StartsWith("--mersenne-device=", StringComparison.OrdinalIgnoreCase)) - { - mersenneOnGpu = !arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); - } - else if (arg.Equals("--workaround-mod", StringComparison.OrdinalIgnoreCase)) - { - useModuloWorkaround = true; - } - else if (arg.Equals("--use-order", StringComparison.OrdinalIgnoreCase)) - { - useOrder = true; - } - else if (arg.StartsWith("--order-warmup-limit=", StringComparison.OrdinalIgnoreCase)) - { - // stored and used below when initializing testers - // Reusing parsedLimit to capture the --order-warmup-limit argument. - // TODO: Use the Utf8Parser-based reader here to eliminate the temporary substring allocation - // when parsing warm-up limits. - if (ulong.TryParse(arg.AsSpan(arg.IndexOf('=') + 1), out parsedLimit)) - { - _orderWarmupLimitOverride = parsedLimit; - } - } - else if (arg.Equals("--gcd-filter", StringComparison.OrdinalIgnoreCase)) - { - _useGcdFilter = true; - } - else if (arg.StartsWith("--ntt=", StringComparison.OrdinalIgnoreCase)) - { - ReadOnlySpan value = arg.AsSpan(arg.IndexOf('=') + 1); - if (value.Equals("staged", StringComparison.OrdinalIgnoreCase)) - { - NttGpuMath.GpuTransformBackend = NttBackend.Staged; - } - else - { - NttGpuMath.GpuTransformBackend = NttBackend.Reference; - } - } - else if (arg.StartsWith("--mod-reduction=", StringComparison.OrdinalIgnoreCase)) - { - ReadOnlySpan value = arg.AsSpan(arg.IndexOf('=') + 1); - if (value.Equals("mont64", StringComparison.OrdinalIgnoreCase)) - { - NttGpuMath.ReductionMode = ModReductionMode.Mont64; - } - else if (value.Equals("barrett128", StringComparison.OrdinalIgnoreCase)) - { - NttGpuMath.ReductionMode = ModReductionMode.Barrett128; - } - else if (value.Equals("uint128", StringComparison.OrdinalIgnoreCase)) - { - NttGpuMath.ReductionMode = ModReductionMode.GpuUInt128; - } - else - { - NttGpuMath.ReductionMode = ModReductionMode.Auto; - } - } - else if (arg.StartsWith("--primes-device=", StringComparison.OrdinalIgnoreCase)) - { - // Controls default device for library prime-related GPU kernels - // that are not explicitly parameterized (backward compatibility). - GpuContextPool.ForceCpu = arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); - } - // New: choose device for order computations (warm-ups and order scans) - else if (arg.StartsWith("--order-device=", StringComparison.OrdinalIgnoreCase)) - { - orderOnGpu = !arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); - } - // Optional: RLE blacklist for p - else if (arg.StartsWith("--rle-blacklist=", StringComparison.OrdinalIgnoreCase)) - { - _rleBlacklistPath = arg[(arg.IndexOf('=') + 1)..]; - } - else if (arg.StartsWith("--rle-hard-max=", StringComparison.OrdinalIgnoreCase)) - { - ulong rleMaxP; - // TODO: Replace this TryParse with the Utf8Parser span helper once shared CLI parsing utilities - // land, mirroring the other numeric option optimizations. - if (ulong.TryParse(arg.AsSpan(arg.IndexOf('=') + 1), out rleMaxP)) - { - _rleHardMaxP = rleMaxP; - } - } - else if (arg.StartsWith("--rle-only-last7=", StringComparison.OrdinalIgnoreCase)) - { - ReadOnlySpan value = arg.AsSpan(arg.IndexOf('=') + 1); - _rleOnlyLast7 = !value.Equals("false", StringComparison.OrdinalIgnoreCase) && !value.Equals("0", StringComparison.OrdinalIgnoreCase); - } - else if (arg.StartsWith("--zero-hard=", StringComparison.OrdinalIgnoreCase)) - { - double zeroFraction; - // TODO: Swap to the Utf8Parser double fast-path so zero-hard parsing avoids culture-dependent - // conversions in the hot CLI path. - if (double.TryParse(arg.AsSpan(arg.IndexOf('=') + 1), out zeroFraction)) - { - _zeroFracHard = zeroFraction; - } - } - else if (arg.StartsWith("--zero-conj=", StringComparison.OrdinalIgnoreCase)) - { - // format: : - ReadOnlySpan value = arg.AsSpan(arg.IndexOf('=') + 1); - int colon = value.IndexOf(':'); - if (colon > 0) - { - double zeroFrac; - int maxZero; - // TODO: Replace these double/int TryParse calls with span-based Utf8Parser helpers once - // available so parsing zero-conj stays allocation-free. - if (double.TryParse(value[..colon], out zeroFrac) && int.TryParse(value[(colon + 1)..], out maxZero)) - { - _zeroFracConj = zeroFrac; - _maxZeroConj = maxZero; - } - } - } - // Backward-compat: accept deprecated device flags and map to new ones - else if (arg.StartsWith("--lucas=", StringComparison.OrdinalIgnoreCase)) - { - mersenneOnGpu = !arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); - } - else if (arg.StartsWith("--primes=", StringComparison.OrdinalIgnoreCase)) - { - GpuContextPool.ForceCpu = arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); - } - else if (arg.StartsWith("--gpu-kernels=", StringComparison.OrdinalIgnoreCase) || arg.StartsWith("--accelerator=", StringComparison.OrdinalIgnoreCase)) - { - // Legacy alias: treat as primes-device for backward compat - GpuContextPool.ForceCpu = arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); - } - else if (arg.StartsWith("--results-dir=", StringComparison.OrdinalIgnoreCase)) - { - _resultsDir = arg[(arg.IndexOf('=') + 1)..]; - } - else if (arg.StartsWith("--results-prefix=", StringComparison.OrdinalIgnoreCase)) - { - _resultsPrefix = arg[(arg.IndexOf('=') + 1)..]; - } - - else if (arg.StartsWith("--gpu-prime-threads=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Convert this int.Parse to the Utf8Parser-based helper to keep CLI option parsing - // allocation-free in the hot startup path. - gpuPrimeThreads = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); - } - else if (arg.StartsWith("--gpu-prime-batch=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Swap int.Parse for Utf8Parser to align with the faster CLI numeric parsing path. - gpuPrimeBatch = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); - } - else if (arg.StartsWith("--ll-slice=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Replace int.Parse with Utf8Parser-based parsing to eliminate temporary strings when - // reading slice configuration. - sliceSize = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); - } - else if (arg.StartsWith("--gpu-scan-batch=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Use the Utf8Parser-based fast path here to match the other CLI numeric parsing - // optimizations we plan to adopt. - scanBatchSize = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); - } - else if (arg.StartsWith("--block-size=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Inline the Utf8Parser-based int reader here once shared CLI parsing utilities land so - // we avoid repeated string allocations. - blockSize = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); - } - else if (arg.StartsWith("--filter-p=", StringComparison.OrdinalIgnoreCase)) - { - filterFile = arg[(arg.IndexOf('=') + 1)..]; - } - else if (arg.StartsWith("--write-batch-size=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Replace int.Parse with the shared Utf8Parser fast-path to keep hot CLI option parsing - // allocation-free. - _writeBatchSize = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); - } - if (arg.StartsWith("--divisor-cycles=")) - { - cyclesPath = arg[(arg.IndexOf('=') + 1)..]; - } - else if (arg.StartsWith("--divisor-cycles-device=", StringComparison.OrdinalIgnoreCase)) - { - useGpuCycles = !arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); - } - else if (arg.StartsWith("--divisor-cycles-batch=", StringComparison.OrdinalIgnoreCase)) - { - // TODO: Convert this int.Parse to Utf8Parser once the optimized CLI parsing helper is shared - // across tools. - cyclesBatchSize = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); - } - else if (arg.Equals("--divisor-cycles-continue", StringComparison.OrdinalIgnoreCase)) - { - continueCyclesGeneration = true; - } - } - - if (useByDivisor) - { - _byDivisorStartPrime = startPrimeProvided ? currentP : 0UL; - } - else - { - _byDivisorStartPrime = 0UL; - } - - if (useByDivisor && string.IsNullOrEmpty(filterFile)) - { - Console.WriteLine("--mersenne=bydivisor requires --filter-p=."); - return; - } - - if (showHelp) - { - PrintHelp(); - return; - } - - // Apply GPU prime sieve runtime configuration - GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads); - PrimeTester.GpuBatchSize = Math.Max(1, gpuPrimeBatch); - - MersenneDivisorCycles mersenneDivisorCycles = new(); - if (!File.Exists(cyclesPath) || continueCyclesGeneration) - { - Console.WriteLine($"Generating divisor cycles '{cyclesPath}' for all p <= {PerfectNumberConstants.MaxQForDivisorCycles}..."); - - long nextPosition = 0L, completeCount = 0L; - if (continueCyclesGeneration && File.Exists(cyclesPath)) - { - Console.WriteLine("Finding last divisor position..."); - (nextPosition, completeCount) = MersenneDivisorCycles.FindLast(cyclesPath); - } - - if (useGpuCycles) - { - if (!continueCyclesGeneration || completeCount == 0L) - { - Console.WriteLine("Starting generation..."); - } - else - { - Console.WriteLine("Resuming generation..."); - } - - // TODO: Wire GenerateGpu to the unrolled-hex kernel that led the MersenneDivisorCycleLengthGpuBenchmarks once it lands. - MersenneDivisorCycles.GenerateGpu(cyclesPath, PerfectNumberConstants.MaxQForDivisorCycles, cyclesBatchSize, skipCount: completeCount, nextPosition: nextPosition); - } - else - { - MersenneDivisorCycles.Generate(cyclesPath, PerfectNumberConstants.MaxQForDivisorCycles); - } - } - - Console.WriteLine($"Loading divisor cycles into memory..."); - if (string.Equals(Path.GetExtension(cyclesPath), ".csv", StringComparison.OrdinalIgnoreCase)) - { - MersenneDivisorCycles.Shared.LoadFrom(cyclesPath); - } - else - { - MersenneDivisorCycles.Shared.LoadFrom(cyclesPath); - } - - DivisorCycleCache.SetDivisorCyclesBatchSize(cyclesBatchSize); - DivisorCycleCache.Shared.ConfigureGeneratorDevice(orderOnGpu); - // TODO: Keep a single cached block loaded from disk and honor the configured device when - // computing ad-hoc cycles for divisors that fall outside that snapshot instead of queuing - // generation of additional blocks. - DivisorCycleCache.Shared.ReloadFromCurrentSnapshot(); - // TODO: Stop reloading the full snapshot once the ad-hoc path streams results straight from - // the configured device without persisting them, so memory stays bounded while preserving - // the single cached block strategy. - - Console.WriteLine("Divisor cycles are ready"); - - if (useByDivisor) - { - threadCount = 1; - blockSize = 1; - } - - _useDivisor = useDivisor; - _useResidueMode = useResidue; - _useByDivisorMode = useByDivisor; - _divisor = divisor; - if (useBitTransform) - { - _transformP = &TransformPBit; - } - else if (useResidue) - { - _transformP = &TransformPAdd; - } - else - { - _transformP = &TransformPAddPrimes; - } - PResidue = new ThreadLocal(() => new ModResidueTracker(ResidueModel.Identity, initialNumber: currentP, initialized: true), trackAllValues: true); - PrimeTesters = new ThreadLocal(() => new PrimeTester(), trackAllValues: true); - // Note: --primes-device controls default device for library kernels; p primality remains CPU here. - // Initialize per-thread p residue tracker (Identity model) at currentP - if (!useDivisor && !useByDivisor) - { - MersenneTesters = new ThreadLocal(() => - { - // TODO: Swap the underlying pow2mod kernels to ProcessEightBitWindows once Pow2Minus1Mod migrates to the windowed helper highlighted in GpuPow2ModBenchmarks. - var tester = new MersenneNumberTester( - useIncremental: !useLucas, - useOrderCache: false, - kernelType: kernelType, - useModuloWorkaround: useModuloWorkaround, - useOrder: useOrder, - useGpuLucas: mersenneOnGpu, - useGpuScan: mersenneOnGpu, - useGpuOrder: orderOnGpu, - useResidue: useResidue, - maxK: residueKMax); - if (!useLucas) - { - tester.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL); - } - - return tester; - }, trackAllValues: true); - } - else if (useDivisor) - { - if (divisor == UInt128.Zero) - { - MersenneNumberDivisorGpuTester.BuildDivisorCandidates(); - } - - _divisorTester = new MersenneNumberDivisorGpuTester(); - } - else if (useByDivisor) - { - _byDivisorTester = mersenneOnGpu - ? new MersenneNumberDivisorByDivisorGpuTester() - : new MersenneNumberDivisorByDivisorCpuTester(); - _byDivisorTester.BatchSize = scanBatchSize; - _byDivisorTester.UseDivisorCycles = useDivisorCycles; // TODO: Collapse this assignment once the runtime enforces divisor-cycle acceleration for every by-divisor scan path. - } - - // Load RLE blacklist (optional) - if (!string.IsNullOrEmpty(_rleBlacklistPath)) - { - RleBlacklist.Load(_rleBlacklistPath!); - } - - // Mod6 lookup turned out slower (see Mod6ComparisonBenchmarks), so keep `%` here for large candidates. - if (currentP == InitialP && string.IsNullOrEmpty(filterFile)) - { - // bool passedAllTests = IsEvenPerfectCandidate(InitialP, out bool searchedMersenne, out bool detailedCheck); - // Skip the already processed range below 138 million - while (currentP < 138_000_000UL && !Volatile.Read(ref _limitReached)) - { - currentP = _transformP(currentP, ref remainder); - } - - if (Volatile.Read(ref _limitReached)) - { - return; - } - } - - Console.WriteLine("Initialization..."); - // Compose a results file name that encodes configuration (before opening file) - var builtName = BuildResultsFileName( - useBitTransform, - threadCount, - blockSize, - kernelType, - useLucas, - useDivisor, - useByDivisor, - mersenneOnGpu, - useOrder, - useModuloWorkaround, - _useGcdFilter, - NttGpuMath.GpuTransformBackend, - gpuPrimeThreads, - sliceSize, - scanBatchSize, - _orderWarmupLimitOverride ?? 5_000_000UL, - NttGpuMath.ReductionMode, - mersenneOnGpu ? "gpu" : "cpu", - (GpuContextPool.ForceCpu ? "cpu" : "gpu"), - orderOnGpu ? "gpu" : "cpu"); - - if (!string.IsNullOrEmpty(_resultsPrefix)) - { - builtName = _resultsPrefix + "_" + builtName; - } - ResultsFileName = string.IsNullOrEmpty(_resultsDir) ? builtName : Path.Combine(_resultsDir!, builtName); - var dir = Path.GetDirectoryName(ResultsFileName); - if (!string.IsNullOrEmpty(dir)) - { - Directory.CreateDirectory(dir!); - } - - - Dictionary? previousResults = useByDivisor ? new Dictionary() : null; - _byDivisorPreviousResults = previousResults; - - if (File.Exists(ResultsFileName)) + // NTT backend selection (GPU): reference vs staged + for (; argIndex < args.Length; argIndex++) + { + arg = args[argIndex]; + if (arg.Equals("--?", StringComparison.OrdinalIgnoreCase) || arg.Equals("-?", StringComparison.OrdinalIgnoreCase) || arg.Equals("--help", StringComparison.OrdinalIgnoreCase) || arg.Equals("-help", StringComparison.OrdinalIgnoreCase) || arg.Equals("/?", StringComparison.OrdinalIgnoreCase)) + { + showHelp = true; + break; + } + else if (arg.StartsWith("--prime=", StringComparison.OrdinalIgnoreCase)) + { + // TODO: Replace this ulong.Parse with the Utf8Parser-based span helper once CLI option parsing + // adopts the zero-allocation fast path proven fastest in the parser benchmarks. + currentP = ulong.Parse(arg.AsSpan(arg.IndexOf('=') + 1)); + startPrimeProvided = true; + } + else if (arg.Equals("--increment=bit", StringComparison.OrdinalIgnoreCase)) + { + useBitTransform = true; + } + else if (arg.Equals("--increment=add", StringComparison.OrdinalIgnoreCase)) + { + useBitTransform = false; + } + else if (arg.StartsWith("--threads=", StringComparison.OrdinalIgnoreCase)) + { + // TODO: Swap int.Parse for the span-based Utf8Parser fast path to avoid transient strings when + // processing CLI numeric options. + threadCount = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); + } + else if (arg.StartsWith("--mersenne=", StringComparison.OrdinalIgnoreCase)) + { + mersenneOption = arg.AsSpan(arg.IndexOf('=') + 1); + if (mersenneOption.Equals("pow2mod", StringComparison.OrdinalIgnoreCase)) { - Console.WriteLine("Processing previous results..."); - LoadResultsFile(ResultsFileName, (p, detailedCheck, passedAllTests) => - { - if (detailedCheck && passedAllTests) - { - _primeCount++; - } - - if (previousResults is not null) - { - previousResults[p] = (detailedCheck, passedAllTests); - } - }); + kernelType = GpuKernelType.Pow2Mod; // TODO: Switch this mode to the ProcessEightBitWindows kernel once Pow2Minus1Mod adopts the benchmarked faster windowed ladder. + } + else if (mersenneOption.Equals("lucas", StringComparison.OrdinalIgnoreCase)) + { + useLucas = true; + } + else if (mersenneOption.Equals("residue", StringComparison.OrdinalIgnoreCase)) + { + useResidue = true; + useLucas = false; + } + else if (mersenneOption.Equals("divisor", StringComparison.OrdinalIgnoreCase)) + { + useDivisor = true; + useLucas = false; + useResidue = false; + } + else if (mersenneOption.Equals("bydivisor", StringComparison.OrdinalIgnoreCase)) + { + useByDivisor = true; + useLucas = false; + useResidue = false; + useDivisor = false; } else { - // TODO: Replace this File.WriteAllText call with the pooled TextFileWriter pipeline once the - // scanner's output flush adopts the benchmarked buffered writer so we avoid allocating the - // entire header string and opening the file twice per run. - File.WriteAllText(ResultsFileName, $"p,searchedMersenne,detailedCheck,passedAllTests{Environment.NewLine}"); + kernelType = GpuKernelType.Incremental; } - - bool useFilter = !string.IsNullOrEmpty(filterFile); - HashSet filter = []; - List byDivisorCandidates = []; - ulong maxP = 0UL; - ulong[] localFilter = Array.Empty(); - int filterCount = 0; - if (useFilter) - { - Console.WriteLine("Loading filter..."); - if (_useByDivisorMode) - { - byDivisorCandidates = LoadByDivisorCandidates(filterFile); - } - else - { - // TODO: Rent this filter buffer from ArrayPool and reuse it across reload batches so - // we do not allocate fresh arrays while replaying large result filters. - localFilter = new ulong[1024]; - filterCount = 0; - LoadResultsFile(filterFile, (p, detailedCheck, passedAllTests) => - { - if (passedAllTests) - { - localFilter[filterCount++] = p; - if (p > maxP) - { - maxP = p; - } - - if (filterCount == 1024) - { - filter.AddRange(localFilter[..filterCount]); - filterCount = 0; - Console.WriteLine($"Added {p}"); - } - } - }); - - if (filterCount > 0) - { - filter.AddRange(localFilter[..filterCount]); - } - } - } - - if (_primeCount >= 2) - { - _primeFoundAfterInit = true; - } - - _outputBuilder = StringBuilderPool.Rent(); - - // Limit GPU concurrency only for prime checks (LL/NTT & GPU order scans). - GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads); - // Configure batch size for GPU primality sieve - PrimeTester.GpuBatchSize = gpuPrimeBatch; - - Console.WriteLine("Starting scan..."); - - if (_useByDivisorMode) - { - RunByDivisorMode(byDivisorCandidates, threadCount); - FlushBuffer(); - StringBuilderPool.Return(_outputBuilder!); - return; - } - - if (!useDivisor && !useByDivisor) - { - Console.WriteLine("Warming up orders..."); - - threadCount = Math.Max(1, threadCount); - _ = MersenneTesters.Value; - if (!useLucas) - { - MersenneTesters.Value!.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL); - } - } - - _state = ((long)currentP << 3) | (long)remainder; - Task[] tasks = new Task[threadCount]; - - int taskIndex = 0; - for (; taskIndex < threadCount; taskIndex++) - { - tasks[taskIndex] = Task.Run(() => - { - int count, j; - ulong p; - bool passedAllTests, searchedMersenne, detailedCheck; - // TODO: Rent this buffer from ArrayPool once the pooled prime-block allocator lands so worker - // threads stop allocating fresh arrays for every reservation. - ulong[] buffer = new ulong[blockSize]; - bool reachedMax = false; - while (!Volatile.Read(ref _limitReached)) - { - count = ReserveBlock(buffer, blockSize); - if (count == 0) - { - break; - } - - if (!useFilter) - { - for (j = 0; j < count && !Volatile.Read(ref _limitReached); j++) - { - p = buffer[j]; - passedAllTests = IsEvenPerfectCandidate(p, divisorCyclesSearchLimit, out searchedMersenne, out detailedCheck); - PrintResult(p, searchedMersenne, detailedCheck, passedAllTests); - } - } - else - { - reachedMax = false; - - for (j = 0; j < count; j++) - { - p = buffer[j]; - - if (Volatile.Read(ref _limitReached) && p > maxP) - { - break; - } - - if (useFilter && !filter.Contains(p)) - { - continue; - } - - passedAllTests = IsEvenPerfectCandidate(p, divisorCyclesSearchLimit, out searchedMersenne, out detailedCheck); - PrintResult(p, searchedMersenne, detailedCheck, passedAllTests); - - if (p == maxP) - { - reachedMax = true; - } - } - - if (reachedMax) - { - Volatile.Write(ref _limitReached, true); - break; - } - } - } - }); - } - - Task.WaitAll(tasks); - FlushBuffer(); - StringBuilderPool.Return(_outputBuilder!); - } - - private static List LoadByDivisorCandidates(string candidateFile) - { - List candidates = []; - using FileStream readStream = new(candidateFile, FileMode.Open, FileAccess.Read, FileShare.ReadWrite); - using StreamReader reader = new(readStream); - string? line; - ReadOnlySpan span; - while ((line = reader.ReadLine()) is not null) - { - if (string.IsNullOrWhiteSpace(line)) - { - continue; - } - - span = line.AsSpan().Trim(); - if (span.IsEmpty || span[0] == '#') - { - continue; - } - - // TODO: Swap this char-by-char parser for the benchmarked Utf8Parser-based fast-path once - // the CLI plumbing accepts spans to avoid redundant range checks while collecting primes. - int index = 0; - while (index < span.Length) - { - if (!char.IsDigit(span[index])) - { - index++; - continue; - } - - int start = index; - index++; - while (index < span.Length && char.IsDigit(span[index])) - { - index++; - } - - // TODO: Swap ulong.TryParse for the Utf8Parser-based fast path once the shared span helpers - // arrive so by-divisor candidate loading avoids per-value string allocations. - if (ulong.TryParse(span[start..index], NumberStyles.None, CultureInfo.InvariantCulture, out ulong parsed)) - { - candidates.Add(parsed); - } - } - } - - return candidates; - } - - private static void RunByDivisorMode(List candidates, int threadCount) - { - if (_byDivisorTester is null) - { - throw new InvalidOperationException("By-divisor tester has not been initialized."); - } - - // TODO: Inline this helper by invoking MersenneNumberDivisorByDivisorTester.Run directly at call sites so - // the CLI avoids this extra frame during the hot composite-reporting path. - MersenneNumberDivisorByDivisorTester.Run( - candidates, - _byDivisorTester, - _byDivisorPreviousResults, - _byDivisorStartPrime, - static () => _lastCompositeP = true, - static () => _lastCompositeP = false, - PrintResult, - threadCount); - } - - private static unsafe void LoadResultsFile(string resultsFileName, Action lineProcessorAction) - { - using FileStream readStream = new(resultsFileName, FileMode.Open, FileAccess.Read, FileShare.ReadWrite); - using StreamReader reader = new(readStream); - string? line; - bool headerSkipped = false; - ReadOnlySpan span; - int first = 0; - int second = 0; - int third = 0; - ulong parsedP = 0UL; - ReadOnlySpan detailedSpan = default; - ReadOnlySpan passedAllTestsSpan = default; - bool detailed = false; - bool passedAllTests = false; - while ((line = reader.ReadLine()) != null) - { - if (!headerSkipped) - { - headerSkipped = true; - continue; - } - - if (string.IsNullOrWhiteSpace(line)) - { - continue; - } - - span = line.AsSpan(); - first = span.IndexOf(','); - if (first < 0) - { - continue; - } - - // TODO: Replace this ulong.Parse call with the Utf8Parser-based span fast-path we benchmarked so results - // reload skips the slower string allocation and culture-aware parsing. - parsedP = ulong.Parse(span[..first]); - span = span[(first + 1)..]; - second = span.IndexOf(','); - if (second < 0) - { - continue; - } - - span = span[(second + 1)..]; - third = span.IndexOf(','); - if (third < 0) - { - continue; - } - - detailedSpan = span[..third]; - passedAllTestsSpan = span[(third + 1)..]; - - // TODO: Replace these bool.TryParse calls with the span-based fast path once we expose the - // Utf8Parser-powered helpers so results reload avoids per-line string allocations for booleans. - if (bool.TryParse(detailedSpan, out detailed) && bool.TryParse(passedAllTestsSpan, out passedAllTests)) - { - lineProcessorAction(parsedP, detailed, passedAllTests); - } - } - } - - private static void PrintHelp() - { - Console.WriteLine("Usage: EvenPerfectBitScanner [options]"); - Console.WriteLine(); - Console.WriteLine("Options:"); - Console.WriteLine(" --prime= starting exponent (p)"); - Console.WriteLine(" --increment=bit|add exponent increment method"); - Console.WriteLine(" --threads= number of worker threads"); - Console.WriteLine(" --block-size= values processed per thread batch"); - Console.WriteLine(" --mersenne=pow2mod|incremental|lucas|residue|divisor|bydivisor Mersenne test method"); - Console.WriteLine(" --divisor= optional divisor for --mersenne=divisor mode"); - Console.WriteLine(" --residue-max-k= max k for residue Mersenne test (q = 2*p*k + 1)"); - Console.WriteLine(" --mersenne-device=cpu|gpu Device for Mersenne method (default gpu)"); - Console.WriteLine(" --primes-device=cpu|gpu Device for prime-scan kernels (default gpu)"); - Console.WriteLine(" --gpu-prime-batch= Batch size for GPU primality sieve (default 262144)"); - Console.WriteLine(" --order-device=cpu|gpu Device for order computations (default gpu)"); - Console.WriteLine(" --ntt=reference|staged GPU NTT backend (default staged)"); - Console.WriteLine(" --mod-reduction=auto|uint128|mont64|barrett128 staged NTT reduction (default auto)"); - Console.WriteLine(" --gpu-prime-threads= max concurrent GPU prime checks (default 1)"); - Console.WriteLine(" --ll-slice= Lucas–Lehmer iterations per slice (default 32)"); - Console.WriteLine(" --gpu-scan-batch= GPU q-scan batch size (default 2_097_152)"); - Console.WriteLine(" --order-warmup-limit= Warm-up order candidates (default 5_000_000)"); - Console.WriteLine(" --rle-blacklist= enable RLE blacklist for p (hard filter up to --rle-hard-max)"); - Console.WriteLine(" --rle-hard-max=

apply RLE blacklist only for p <= this (default ulong.MaxValue = no limit)"); - Console.WriteLine(" --rle-only-last7=true|false apply RLE only when p % 10 == 7 (default true)"); - Console.WriteLine(" --zero-hard= hard reject if zero_fraction(p) > f (default off)"); - Console.WriteLine(" --zero-conj=: hard reject if zero_fraction(p) > f AND max_zero_block >= r (default off)"); - Console.WriteLine(" --results-dir= directory for results file"); - Console.WriteLine(" --results-prefix= prefix to prepend to results filename"); - Console.WriteLine(" --divisor-cycles= divisor cycles data file"); - Console.WriteLine(" --divisor-cycles-device=cpu|gpu device for cycles generation (default gpu)"); - Console.WriteLine(" --divisor-cycles-batch= 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"); - Console.WriteLine(" --use-divisor-cycles=true|false enable divisor cycle filtering in --mersenne=bydivisor (default false)"); - - Console.WriteLine(" --use-order test primality via q order"); - Console.WriteLine(" --workaround-mod avoid '%' operator on the GPU"); - // mod-automaton removed - Console.WriteLine(" --filter-p= process only p from previous run results (required for --mersenne=bydivisor)"); - Console.WriteLine(" --write-batch-size= overwrite frequency of disk writes (default 100 lines)"); - Console.WriteLine(" --gcd-filter enable early sieve based on GCD"); - Console.WriteLine(" --help, -help, --?, -?, /? show this help message"); - } - - private static string BuildResultsFileName(bool bitInc, int threads, int block, GpuKernelType kernelType, bool useLucasFlag, bool useDivisorFlag, bool useByDivisorFlag, bool mersenneOnGpu, bool useOrder, bool useModWorkaround, bool useGcd, NttBackend nttBackend, int gpuPrimeThreads, int llSlice, int gpuScanBatch, ulong warmupLimit, ModReductionMode reduction, string mersenneDevice, string primesDevice, string orderDevice) - { - string inc = bitInc ? "bit" : "add"; - string mers = useDivisorFlag - ? "divisor" - : (useLucasFlag - ? "lucas" - : (useByDivisorFlag - ? "bydivisor" - : (kernelType == GpuKernelType.Pow2Mod ? "pow2mod" : "incremental"))); - string ntt = nttBackend == NttBackend.Staged ? "staged" : "reference"; - string red = reduction switch { ModReductionMode.Mont64 => "mont64", ModReductionMode.Barrett128 => "barrett128", ModReductionMode.GpuUInt128 => "uint128", _ => "auto" }; - string order = useOrder ? "order-on" : "order-off"; - string gcd = useGcd ? "gcd-on" : "gcd-off"; - string work = useModWorkaround ? "modfix-on" : "modfix-off"; - // TODO: Replace this string interpolation with the pooled ValueStringBuilder pipeline from the - // results-writer benchmarks so filename generation reuses the zero-allocation formatter once - // we consolidate output paths. - return $"even_perfect_bit_scan_inc-{inc}_thr-{threads}_blk-{block}_mers-{mers}_mersdev-{mersenneDevice}_ntt-{ntt}_red-{red}_primesdev-{primesDevice}_{order}_orderdev-{orderDevice}_{gcd}_{work}_gputh-{gpuPrimeThreads}_llslice-{llSlice}_scanb-{gpuScanBatch}_warm-{warmupLimit}.csv"; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static unsafe int ReserveBlock(ulong[] buffer, int blockSize) - { - Span bufferSpan = new(buffer); - ulong p, remainder; - long state, newState, original; - int count; - while (true) + } + else if (arg.StartsWith("--divisor-cycles-limit=", StringComparison.OrdinalIgnoreCase)) { - state = Volatile.Read(ref _state); - p = (ulong)state >> 3; - remainder = ((ulong)state) & 7UL; - - count = 0; - while (count < blockSize && !Volatile.Read(ref _limitReached)) + // TODO: Replace ulong.TryParse with the Utf8Parser fast-path once CLI parsing utilities expose + // span-based helpers for optional arguments. + if (ulong.TryParse(arg.AsSpan(arg.IndexOf('=') + 1), out parsedLimit)) { - bufferSpan[count++] = p; - // TODO: Inline the fast-path transform here once we collapse the delegate* indirection so - // block reservations stop paying the extra jump highlighted in the prime-stepping benchmarks. - p = _transformP(p, ref remainder); + divisorCyclesSearchLimit = parsedLimit; + } + } + else if (arg.StartsWith("--residue-max-k=", StringComparison.OrdinalIgnoreCase)) + { + // TODO: Swap this TryParse for the zero-allocation Utf8Parser helper to keep residue CLI option + // parsing in the fast path. + if (ulong.TryParse(arg.AsSpan(arg.IndexOf('=') + 1), out parsedResidueMax)) + { + residueKMax = parsedResidueMax; + } + } + // Replaces --lucas=cpu|gpu; controls device for Lucas and for + // incremental/pow2mod scanning + else if (arg.StartsWith("--mersenne-device=", StringComparison.OrdinalIgnoreCase)) + { + mersenneOnGpu = !arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + } + else if (arg.Equals("--workaround-mod", StringComparison.OrdinalIgnoreCase)) + { + useModuloWorkaround = true; + } + else if (arg.Equals("--use-order", StringComparison.OrdinalIgnoreCase)) + { + useOrder = true; + } + else if (arg.StartsWith("--order-warmup-limit=", StringComparison.OrdinalIgnoreCase)) + { + // stored and used below when initializing testers + // Reusing parsedLimit to capture the --order-warmup-limit argument. + // TODO: Use the Utf8Parser-based reader here to eliminate the temporary substring allocation + // when parsing warm-up limits. + if (ulong.TryParse(arg.AsSpan(arg.IndexOf('=') + 1), out parsedLimit)) + { + _orderWarmupLimitOverride = parsedLimit; + } + } + else if (arg.Equals("--gcd-filter", StringComparison.OrdinalIgnoreCase)) + { + _useGcdFilter = true; + } + else if (arg.StartsWith("--ntt=", StringComparison.OrdinalIgnoreCase)) + { + ReadOnlySpan value = arg.AsSpan(arg.IndexOf('=') + 1); + if (value.Equals("staged", StringComparison.OrdinalIgnoreCase)) + { + NttGpuMath.GpuTransformBackend = NttBackend.Staged; + } + else + { + NttGpuMath.GpuTransformBackend = NttBackend.Reference; + } + } + else if (arg.StartsWith("--mod-reduction=", StringComparison.OrdinalIgnoreCase)) + { + ReadOnlySpan value = arg.AsSpan(arg.IndexOf('=') + 1); + if (value.Equals("mont64", StringComparison.OrdinalIgnoreCase)) + { + NttGpuMath.ReductionMode = ModReductionMode.Mont64; + } + else if (value.Equals("barrett128", StringComparison.OrdinalIgnoreCase)) + { + NttGpuMath.ReductionMode = ModReductionMode.Barrett128; + } + else if (value.Equals("uint128", StringComparison.OrdinalIgnoreCase)) + { + NttGpuMath.ReductionMode = ModReductionMode.GpuUInt128; + } + else + { + NttGpuMath.ReductionMode = ModReductionMode.Auto; + } + } + else if (arg.StartsWith("--primes-device=", StringComparison.OrdinalIgnoreCase)) + { + // Controls default device for library prime-related GPU kernels + // that are not explicitly parameterized (backward compatibility). + GpuContextPool.ForceCpu = arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + } + // New: choose device for order computations (warm-ups and order scans) + else if (arg.StartsWith("--order-device=", StringComparison.OrdinalIgnoreCase)) + { + orderOnGpu = !arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + } + // Optional: RLE blacklist for p + else if (arg.StartsWith("--rle-blacklist=", StringComparison.OrdinalIgnoreCase)) + { + _rleBlacklistPath = arg[(arg.IndexOf('=') + 1)..]; + } + else if (arg.StartsWith("--rle-hard-max=", StringComparison.OrdinalIgnoreCase)) + { + ulong rleMaxP; + // TODO: Replace this TryParse with the Utf8Parser span helper once shared CLI parsing utilities + // land, mirroring the other numeric option optimizations. + if (ulong.TryParse(arg.AsSpan(arg.IndexOf('=') + 1), out rleMaxP)) + { + _rleHardMaxP = rleMaxP; + } + } + else if (arg.StartsWith("--rle-only-last7=", StringComparison.OrdinalIgnoreCase)) + { + ReadOnlySpan value = arg.AsSpan(arg.IndexOf('=') + 1); + _rleOnlyLast7 = !value.Equals("false", StringComparison.OrdinalIgnoreCase) && !value.Equals("0", StringComparison.OrdinalIgnoreCase); + } + else if (arg.StartsWith("--zero-hard=", StringComparison.OrdinalIgnoreCase)) + { + double zeroFraction; + // TODO: Swap to the Utf8Parser double fast-path so zero-hard parsing avoids culture-dependent + // conversions in the hot CLI path. + if (double.TryParse(arg.AsSpan(arg.IndexOf('=') + 1), out zeroFraction)) + { + _zeroFracHard = zeroFraction; + } + } + else if (arg.StartsWith("--zero-conj=", StringComparison.OrdinalIgnoreCase)) + { + // format: : + ReadOnlySpan value = arg.AsSpan(arg.IndexOf('=') + 1); + int colon = value.IndexOf(':'); + if (colon > 0) + { + double zeroFrac; + int maxZero; + // TODO: Replace these double/int TryParse calls with span-based Utf8Parser helpers once + // available so parsing zero-conj stays allocation-free. + if (double.TryParse(value[..colon], out zeroFrac) && int.TryParse(value[(colon + 1)..], out maxZero)) + { + _zeroFracConj = zeroFrac; + _maxZeroConj = maxZero; + } + } + } + // Backward-compat: accept deprecated device flags and map to new ones + else if (arg.StartsWith("--lucas=", StringComparison.OrdinalIgnoreCase)) + { + mersenneOnGpu = !arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + } + else if (arg.StartsWith("--primes=", StringComparison.OrdinalIgnoreCase)) + { + GpuContextPool.ForceCpu = arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + } + else if (arg.StartsWith("--gpu-kernels=", StringComparison.OrdinalIgnoreCase) || arg.StartsWith("--accelerator=", StringComparison.OrdinalIgnoreCase)) + { + // Legacy alias: treat as primes-device for backward compat + GpuContextPool.ForceCpu = arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + } + else if (arg.StartsWith("--results-dir=", StringComparison.OrdinalIgnoreCase)) + { + _resultsDir = arg[(arg.IndexOf('=') + 1)..]; + } + else if (arg.StartsWith("--results-prefix=", StringComparison.OrdinalIgnoreCase)) + { + _resultsPrefix = arg[(arg.IndexOf('=') + 1)..]; + } + + else if (arg.StartsWith("--gpu-prime-threads=", StringComparison.OrdinalIgnoreCase)) + { + // TODO: Convert this int.Parse to the Utf8Parser-based helper to keep CLI option parsing + // allocation-free in the hot startup path. + gpuPrimeThreads = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); + } + else if (arg.StartsWith("--gpu-prime-batch=", StringComparison.OrdinalIgnoreCase)) + { + // TODO: Swap int.Parse for Utf8Parser to align with the faster CLI numeric parsing path. + gpuPrimeBatch = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); + } + else if (arg.StartsWith("--ll-slice=", StringComparison.OrdinalIgnoreCase)) + { + // TODO: Replace int.Parse with Utf8Parser-based parsing to eliminate temporary strings when + // reading slice configuration. + sliceSize = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); + } + else if (arg.StartsWith("--gpu-scan-batch=", StringComparison.OrdinalIgnoreCase)) + { + // TODO: Use the Utf8Parser-based fast path here to match the other CLI numeric parsing + // optimizations we plan to adopt. + scanBatchSize = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); + } + else if (arg.StartsWith("--block-size=", StringComparison.OrdinalIgnoreCase)) + { + // TODO: Inline the Utf8Parser-based int reader here once shared CLI parsing utilities land so + // we avoid repeated string allocations. + blockSize = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); + } + else if (arg.StartsWith("--filter-p=", StringComparison.OrdinalIgnoreCase)) + { + filterFile = arg[(arg.IndexOf('=') + 1)..]; + } + else if (arg.StartsWith("--write-batch-size=", StringComparison.OrdinalIgnoreCase)) + { + // TODO: Replace int.Parse with the shared Utf8Parser fast-path to keep hot CLI option parsing + // allocation-free. + _writeBatchSize = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); + } + if (arg.StartsWith("--divisor-cycles=")) + { + cyclesPath = arg[(arg.IndexOf('=') + 1)..]; + } + else if (arg.StartsWith("--divisor-cycles-device=", StringComparison.OrdinalIgnoreCase)) + { + useGpuCycles = !arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); + } + else if (arg.StartsWith("--divisor-cycles-batch=", StringComparison.OrdinalIgnoreCase)) + { + // TODO: Convert this int.Parse to Utf8Parser once the optimized CLI parsing helper is shared + // across tools. + cyclesBatchSize = Math.Max(1, int.Parse(arg.AsSpan(arg.IndexOf('=') + 1))); + } + else if (arg.Equals("--divisor-cycles-continue", StringComparison.OrdinalIgnoreCase)) + { + continueCyclesGeneration = true; + } + } + + if (useByDivisor) + { + _byDivisorStartPrime = startPrimeProvided ? currentP : 0UL; + } + else + { + _byDivisorStartPrime = 0UL; + } + + if (useByDivisor && string.IsNullOrEmpty(filterFile)) + { + Console.WriteLine("--mersenne=bydivisor requires --filter-p=."); + return; + } + + if (showHelp) + { + PrintHelp(); + return; + } + + // Apply GPU prime sieve runtime configuration + GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads); + PrimeTester.GpuBatchSize = Math.Max(1, gpuPrimeBatch); + + MersenneDivisorCycles mersenneDivisorCycles = new(); + if (!File.Exists(cyclesPath) || continueCyclesGeneration) + { + Console.WriteLine($"Generating divisor cycles '{cyclesPath}' for all p <= {PerfectNumberConstants.MaxQForDivisorCycles}..."); + + long nextPosition = 0L, completeCount = 0L; + if (continueCyclesGeneration && File.Exists(cyclesPath)) + { + Console.WriteLine("Finding last divisor position..."); + (nextPosition, completeCount) = MersenneDivisorCycles.FindLast(cyclesPath); + } + + if (useGpuCycles) + { + if (!continueCyclesGeneration || completeCount == 0L) + { + Console.WriteLine("Starting generation..."); + } + else + { + Console.WriteLine("Resuming generation..."); } - if (count == 0) + // TODO: Wire GenerateGpu to the unrolled-hex kernel that led the MersenneDivisorCycleLengthGpuBenchmarks once it lands. + MersenneDivisorCycles.GenerateGpu(cyclesPath, PerfectNumberConstants.MaxQForDivisorCycles, cyclesBatchSize, skipCount: completeCount, nextPosition: nextPosition); + } + else + { + MersenneDivisorCycles.Generate(cyclesPath, PerfectNumberConstants.MaxQForDivisorCycles); + } + } + + Console.WriteLine($"Loading divisor cycles into memory..."); + if (string.Equals(Path.GetExtension(cyclesPath), ".csv", StringComparison.OrdinalIgnoreCase)) + { + MersenneDivisorCycles.Shared.LoadFrom(cyclesPath); + } + else + { + MersenneDivisorCycles.Shared.LoadFrom(cyclesPath); + } + + DivisorCycleCache.SetDivisorCyclesBatchSize(cyclesBatchSize); + DivisorCycleCache.Shared.ConfigureGeneratorDevice(useGpuCycles); + // TODO: Keep a single cached block loaded from disk and honor the configured device when + // computing ad-hoc cycles for divisors that fall outside that snapshot instead of queuing + // generation of additional blocks. + DivisorCycleCache.Shared.RefreshSnapshot(); + // TODO: Stop reloading the full snapshot once the ad-hoc path streams results straight from + // the configured device without persisting them, so memory stays bounded while preserving + // the single cached block strategy. + + Console.WriteLine("Divisor cycles are ready"); + + if (useByDivisor) + { + threadCount = 1; + blockSize = 1; + } + + _useDivisor = useDivisor; + _useResidueMode = useResidue; + _useByDivisorMode = useByDivisor; + if (useBitTransform) + { + _transformP = &TransformPBit; + } + else if (useResidue) + { + _transformP = &TransformPAdd; + } + else + { + _transformP = &TransformPAddPrimes; + } + PResidue = new ThreadLocal(() => new ModResidueTracker(ResidueModel.Identity, initialNumber: currentP, initialized: true), trackAllValues: true); + PrimeTesters = new ThreadLocal(() => new PrimeTester(), trackAllValues: true); + // Note: --primes-device controls default device for library kernels; p primality remains CPU here. + // Initialize per-thread p residue tracker (Identity model) at currentP + if (!useDivisor && !useByDivisor) + { + MersenneTesters = new ThreadLocal(() => + { + // TODO: Swap the underlying pow2mod kernels to ProcessEightBitWindows once Pow2Minus1Mod migrates to the windowed helper highlighted in GpuPow2ModBenchmarks. + var tester = new MersenneNumberTester( + useIncremental: !useLucas, + useOrderCache: false, + kernelType: kernelType, + useModuloWorkaround: useModuloWorkaround, + useOrder: useOrder, + useGpuLucas: mersenneOnGpu, + useGpuScan: mersenneOnGpu, + useGpuOrder: orderOnGpu, + useResidue: useResidue, + maxK: residueKMax); + if (!useLucas) { - return 0; + tester.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL); } - newState = ((long)p << 3) | (long)remainder; - original = Interlocked.CompareExchange(ref _state, newState, state); - if (original == state) + return tester; + }, trackAllValues: true); + } + else if (useDivisor) + { + MersenneNumberDivisorGpuTester.BuildDivisorCandidates(); + + _divisorTester = new MersenneNumberDivisorGpuTester(); + } + else if (useByDivisor) + { + _byDivisorTester = mersenneOnGpu + ? new MersenneNumberDivisorByDivisorGpuTester() + : new MersenneNumberDivisorByDivisorCpuTester(); + _byDivisorTester.BatchSize = scanBatchSize; + } + + // Load RLE blacklist (optional) + if (!string.IsNullOrEmpty(_rleBlacklistPath)) + { + RleBlacklist.Load(_rleBlacklistPath!); + } + + // Mod6 lookup turned out slower (see Mod6ComparisonBenchmarks), so keep `%` here for large candidates. + if (currentP == InitialP && string.IsNullOrEmpty(filterFile)) + { + // bool passedAllTests = IsEvenPerfectCandidate(InitialP, out bool searchedMersenne, out bool detailedCheck); + // Skip the already processed range below 138 million + while (currentP < 138_000_000UL && !Volatile.Read(ref _limitReached)) + { + currentP = _transformP(currentP, ref remainder); + } + + if (Volatile.Read(ref _limitReached)) + { + return; + } + } + + Console.WriteLine("Initialization..."); + // Compose a results file name that encodes configuration (before opening file) + var builtName = BuildResultsFileName( + useBitTransform, + threadCount, + blockSize, + kernelType, + useLucas, + useDivisor, + useByDivisor, + mersenneOnGpu, + useOrder, + useModuloWorkaround, + _useGcdFilter, + NttGpuMath.GpuTransformBackend, + gpuPrimeThreads, + sliceSize, + scanBatchSize, + _orderWarmupLimitOverride ?? 5_000_000UL, + NttGpuMath.ReductionMode, + mersenneOnGpu ? "gpu" : "cpu", + (GpuContextPool.ForceCpu ? "cpu" : "gpu"), + orderOnGpu ? "gpu" : "cpu"); + + if (!string.IsNullOrEmpty(_resultsPrefix)) + { + builtName = _resultsPrefix + "_" + builtName; + } + ResultsFileName = string.IsNullOrEmpty(_resultsDir) ? builtName : Path.Combine(_resultsDir!, builtName); + var dir = Path.GetDirectoryName(ResultsFileName); + if (!string.IsNullOrEmpty(dir)) + { + Directory.CreateDirectory(dir!); + } + + + Dictionary? previousResults = useByDivisor ? new Dictionary() : null; + _byDivisorPreviousResults = previousResults; + + if (File.Exists(ResultsFileName)) + { + Console.WriteLine("Processing previous results..."); + LoadResultsFile(ResultsFileName, (p, detailedCheck, passedAllTests) => + { + if (detailedCheck && passedAllTests) { - return count; + _primeCount++; + } + + if (previousResults is not null) + { + previousResults[p] = (detailedCheck, passedAllTests); + } + }); + } + else + { + // TODO: Replace this File.WriteAllText call with the pooled TextFileWriter pipeline once the + // scanner's output flush adopts the benchmarked buffered writer so we avoid allocating the + // entire header string and opening the file twice per run. + File.WriteAllText(ResultsFileName, $"p,searchedMersenne,detailedCheck,passedAllTests{Environment.NewLine}"); + } + + bool useFilter = !string.IsNullOrEmpty(filterFile); + HashSet filter = []; + List byDivisorCandidates = []; + ulong maxP = 0UL; + ulong[] localFilter = Array.Empty(); + int filterCount = 0; + if (useFilter) + { + Console.WriteLine("Loading filter..."); + if (_useByDivisorMode) + { + byDivisorCandidates = LoadByDivisorCandidates(filterFile); + } + else + { + // TODO: Rent this filter buffer from ArrayPool and reuse it across reload batches so + // we do not allocate fresh arrays while replaying large result filters. + localFilter = new ulong[1024]; + filterCount = 0; + LoadResultsFile(filterFile, (p, detailedCheck, passedAllTests) => + { + if (passedAllTests) + { + localFilter[filterCount++] = p; + if (p > maxP) + { + maxP = p; + } + + if (filterCount == 1024) + { + filter.AddRange(localFilter[..filterCount]); + filterCount = 0; + Console.WriteLine($"Added {p}"); + } + } + }); + + if (filterCount > 0) + { + filter.AddRange(localFilter[..filterCount]); } } } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool TryAcquireConsoleSlot(bool primeCandidate, out bool primeFlag) - { - int consoleInterval = PerfectNumberConstants.ConsoleInterval; - int counterValue = Interlocked.Increment(ref _consoleCounter); - if (counterValue >= consoleInterval) - { - int previous = Interlocked.Exchange(ref _consoleCounter, 0); - if (previous >= consoleInterval) - { - primeFlag = primeCandidate; - return true; - } - } - - primeFlag = false; - return false; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static StringBuilder? AppendToOutputBuffer(StringBuilder source) - { - StringBuilder? builderToFlush = null; - - lock (Sync) - { - StringBuilder outputBuilder = _outputBuilder!; - _ = outputBuilder.Append(source); - - int nextIndex = _writeIndex + 1; - if (nextIndex >= _writeBatchSize) - { - builderToFlush = outputBuilder; - _outputBuilder = StringBuilderPool.Rent(); - _writeIndex = 0; - } - else - { - _writeIndex = nextIndex; - } - } - - return builderToFlush; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void PrintResult(ulong currentP, bool searchedMersenne, bool detailedCheck, bool passedAllTests) - { - if (passedAllTests) - { - int newCount = Interlocked.Increment(ref _primeCount); - if (newCount >= 2) - { - Volatile.Write(ref _primeFoundAfterInit, true); - } - } - - bool lastWasComposite = _lastCompositeP; - bool primeCandidate = passedAllTests && !lastWasComposite && Volatile.Read(ref _primeFoundAfterInit); - - // TODO: Switch this StringBuilder-based formatter to a span-friendly Utf8Formatter pipeline so - // console/file output avoids intermediate builders in the hot logging loop. - StringBuilder localBuilder = StringBuilderPool.Rent(); - _ = localBuilder - .Append(currentP).Append(',') - .Append(searchedMersenne).Append(',') - .Append(detailedCheck).Append(',') - .Append(passedAllTests).Append('\n'); - - bool printToConsole = TryAcquireConsoleSlot(primeCandidate, out bool primeFlag); - StringBuilder? builderToFlush = AppendToOutputBuffer(localBuilder); - - if (builderToFlush is not null) - { - FlushBuffer(builderToFlush); - } - - if (printToConsole) - { - if (primeFlag) - { - _ = localBuilder - .Remove(localBuilder.Length - 1, 1) - .Append(PrimeFoundSuffix) - .Append('\n'); - } - - Console.WriteLine(localBuilder.Remove(localBuilder.Length - 1, 1).ToString()); - } - - localBuilder.Clear(); - StringBuilderPool.Return(localBuilder); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void FlushBuffer() - { - StringBuilder? builderToFlush = null; - - lock (Sync) - { - if (_writeIndex == 0 || _outputBuilder is null || _outputBuilder.Length == 0) - { - return; - } - - builderToFlush = _outputBuilder; - _outputBuilder = StringBuilderPool.Rent(); - _writeIndex = 0; - } - - FlushBuffer(builderToFlush!); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void FlushBuffer(StringBuilder builderToFlush) - { - try - { - lock (FileWriteSync) - { - // TODO: Replace these per-flush FileStream/StreamWriter allocations with a pooled TextFileWriter-style - // helper so batched result writes keep the hot path on the benchmarked persistent-handle pipeline. - using FileStream stream = new(ResultsFileName, FileMode.Append, FileAccess.Write, FileShare.Read); - using StreamWriter writer = new(stream) { AutoFlush = false }; - writer.Write(builderToFlush); - writer.Flush(); - } - } - finally - { - builderToFlush.Clear(); - StringBuilderPool.Return(builderToFlush); - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong CountOnes(ulong value) + if (_primeCount >= 2) { - // TODO: Remove this wrapper once callers can invoke BitOperations.PopCount directly so the hot bit-stat path - // avoids the extra call layer. - return (ulong)BitOperations.PopCount(value); + _primeFoundAfterInit = true; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong GetNextAddDiff(ulong remainder) - { - return remainder switch - { - 0UL => 1UL, - 1UL => 4UL, - 2UL => 3UL, - 3UL => 2UL, - 4UL => 1UL, - _ => 2UL, - }; - } + _outputBuilder = StringBuilderPool.Rent(); - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong TransformPBit(ulong value, ref ulong remainder) - { - ulong original = value; - if (value > (ulong.MaxValue >> 1)) - { - Volatile.Write(ref _limitReached, true); - return value; - } + // Limit GPU concurrency only for prime checks (LL/NTT & GPU order scans). + GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads); + // Configure batch size for GPU primality sieve + PrimeTester.GpuBatchSize = gpuPrimeBatch; - ulong next = (value << 1) | 1UL; - remainder = (remainder << 1) + 1UL; - // Mod6 lookup loses to `%` in the benches; stick with subtraction + modulo for correctness and speed. - while (remainder >= 6UL) + Console.WriteLine("Starting scan..."); + + if (_useByDivisorMode) + { + RunByDivisorMode(byDivisorCandidates, threadCount); + FlushBuffer(); + StringBuilderPool.Return(_outputBuilder!); + return; + } + + if (!useDivisor && !useByDivisor) + { + Console.WriteLine("Warming up orders..."); + + threadCount = Math.Max(1, threadCount); + _ = MersenneTesters.Value; + if (!useLucas) + { + MersenneTesters.Value!.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL); + } + } + + _state = ((long)currentP << 3) | (long)remainder; + Task[] tasks = new Task[threadCount]; + + int taskIndex = 0; + for (; taskIndex < threadCount; taskIndex++) + { + tasks[taskIndex] = Task.Run(() => + { + int count, j; + ulong p; + bool passedAllTests, searchedMersenne, detailedCheck; + // TODO: Rent this buffer from ArrayPool once the pooled prime-block allocator lands so worker + // threads stop allocating fresh arrays for every reservation. + ulong[] buffer = new ulong[blockSize]; + bool reachedMax = false; + while (!Volatile.Read(ref _limitReached)) { - remainder -= 6UL; - } + count = ReserveBlock(buffer, blockSize); + if (count == 0) + { + break; + } - value = remainder switch - { - 0UL => 1UL, - 1UL => 0UL, - 2UL => 5UL, - 3UL => 4UL, - 4UL => 3UL, - _ => 2UL, - }; // 'value' now holds diff - - if (next > ulong.MaxValue - value) - { - Volatile.Write(ref _limitReached, true); - return original; - } - - remainder += value; - // `% 6` remains faster per Mod6ComparisonBenchmarks; keep this modulo fold and continue using subtraction. - while (remainder >= 6UL) - { - remainder -= 6UL; - } - - return next + value; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong TransformPAdd(ulong value, ref ulong remainder) - { - ulong next = value; - ulong diff = GetNextAddDiff(remainder); - - if (next > ulong.MaxValue - diff) - { - Volatile.Write(ref _limitReached, true); - return next; - } - - remainder += diff; - // Retain direct modulo because the Mod6 helper underperforms on 64-bit operands. - while (remainder >= 6UL) - { - remainder -= 6UL; - } - - return next + diff; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong TransformPAddPrimes(ulong value, ref ulong remainder) - { - ulong originalRemainder = remainder; - ulong addRemainder = remainder; - ulong candidate = value; - ulong primeCandidate = value; - bool advanceAdd = true; - bool advancePrime = true; - ulong diff = 0UL; - ulong nextPrime = 0UL; - - while (true) - { - if (advanceAdd) - { - diff = GetNextAddDiff(addRemainder); - if (candidate > ulong.MaxValue - diff) - { - Volatile.Write(ref _limitReached, true); - remainder = originalRemainder; - return candidate; - } - - candidate += diff; - addRemainder += diff; - // Skip the Mod6 lookup: benchmarked `%` stays ahead for these prime increments. - // benchmarked fastest remainder updates instead of looping subtraction. - while (addRemainder >= 6UL) - { - addRemainder -= 6UL; - } - } - - if (advancePrime) - { - try - { - nextPrime = PrimeIterator.Next(in primeCandidate); - } - catch (InvalidOperationException) - { - Volatile.Write(ref _limitReached, true); - remainder = originalRemainder; - return value; - } - - if (nextPrime <= primeCandidate) - { - Volatile.Write(ref _limitReached, true); - remainder = originalRemainder; - return value; - } - - primeCandidate = nextPrime; - } - - if (candidate == primeCandidate) - { - remainder = addRemainder; - return candidate; - } - - advanceAdd = candidate < primeCandidate; - advancePrime = candidate > primeCandidate; - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static bool IsEvenPerfectCandidate(ulong p, ulong divisorCyclesSearchLimit) => - // TODO: Remove this pass-through once callers can provide the full out parameters - // directly; shaving this hop keeps the hot candidate filter on the shortest path. - 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; - - if (_useGcdFilter && IsCompositeByGcd(p)) - { - _lastCompositeP = true; - return false; - } - - // Optional: RLE blacklist and binary-threshold filters on p (safe only when configured) - if (p <= _rleHardMaxP) - { - if (!_rleOnlyLast7 || p.Mod10() == 7UL) - { - if (RleBlacklist.IsLoaded() && RleBlacklist.Matches(p)) - { - return false; - } - } - } - - if (_zeroFracHard >= 0 || (_zeroFracConj >= 0 && _maxZeroConj >= 0)) - { - double zf; - int bitLength = 0; - int zeroCountValue = 0; - int maxZeroBlockValue = 0; - ComputeBitStats(p, out bitLength, out zeroCountValue, out maxZeroBlockValue); - if (_zeroFracHard >= 0) - { - zf = (double)zeroCountValue / bitLength; - if (zf > _zeroFracHard) - { - return false; - } - } - - if (_zeroFracConj >= 0 && _maxZeroConj >= 0) - { - zf = (double)zeroCountValue / bitLength; - if (zf > _zeroFracConj && maxZeroBlockValue >= _maxZeroConj) - { - return false; - } - } - } - - // Fast residue-based composite check for p using small primes - if (!_useDivisor && !_useByDivisorMode && IsCompositeByResidues(p)) - { - searchedMersenne = false; - detailedCheck = false; - _lastCompositeP = true; - return false; - } - - // If primes-device=gpu, route p primality through GPU-assisted sieve with deterministic MR validation. - if (!_useByDivisorMode) - { - if (!(GpuContextPool.ForceCpu - ? PrimeTesters.Value!.IsPrime(p, CancellationToken.None) - : PrimeTesters.Value!.IsPrimeGpu(p, CancellationToken.None))) - { - _lastCompositeP = true; - return false; - } - } - - searchedMersenne = true; - if (_useByDivisorMode) - { - // TODO: Route by-divisor scans through the ProcessEightBitWindows pow2mod kernel once the windowed Pow2Minus1Mod helper replaces the slower single-bit ladder. - return _byDivisorTester!.IsPrime(p, out detailedCheck); - } - - if (_useDivisor) - { - // TODO: Adopt the windowed Pow2Minus1Mod helper here so divisor checks benefit from the ProcessEightBitWindows speedups. - return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck); - } - - detailedCheck = MersenneTesters.Value!.IsMersennePrime(p); - return detailedCheck; - } - - // Use ModResidueTracker with a small set of primes to pre-filter composite p. - private static bool IsCompositeByResidues(ulong p) - { - var tracker = PResidue.Value!; - // TODO: Integrate the divisor-cycle cache here so the small-prime sweep reuses precomputed remainders instead of - // running MergeOrAppend for every candidate and missing the cycle-accelerated early exits. - tracker.BeginMerge(p); - // Use the small prime list from PerfectNumbers.Core to the extent of sqrt(p) - var primes = PrimesGenerator.SmallPrimes; - var primesPow2 = PrimesGenerator.SmallPrimesPow2; - int len = primes.Length; - int primeIndex = 0; - bool divisible = false; - for (; primeIndex < len; primeIndex++) - { - if (primesPow2[primeIndex] > p) - { - break; - } - - if (tracker.MergeOrAppend(p, primes[primeIndex], out divisible) && divisible) - { - return true; - } - } - - return false; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void ComputeBitStats(ulong value, out int bitLen, out int zeroCount, out int maxZeroBlock) - { - bitLen = 64 - int.CreateChecked(ulong.LeadingZeroCount(value)); - zeroCount = 0; - maxZeroBlock = 0; - if (bitLen <= 0) - { - return; - } - - int msbIndex = (bitLen - 1) >> 3; // 0..7 - int bitsInTopByte = ((bitLen - 1) & 7) + 1; // 1..8 - int currentRun = 0; - byte inspectedByte = 0; - int zeroCountInByte = 0; - int prefixZeros = 0; - int suffixZeros = 0; - int maxZeroRunInByte = 0; - int candidate = 0; - int byteIndex = msbIndex; - - for (; byteIndex >= 0; byteIndex--) - { - inspectedByte = (byte)(value >> (byteIndex * 8)); - if (byteIndex == msbIndex && bitsInTopByte < 8) + if (!useFilter) + { + for (j = 0; j < count && !Volatile.Read(ref _limitReached); j++) { - // Mask off leading unused bits by setting them to 1 (ignore) - inspectedByte |= (byte)(0xFF << bitsInTopByte); + p = buffer[j]; + passedAllTests = IsEvenPerfectCandidate(p, divisorCyclesSearchLimit, out searchedMersenne, out detailedCheck); + PrintResult(p, searchedMersenne, detailedCheck, passedAllTests); + } + } + else + { + reachedMax = false; + + for (j = 0; j < count; j++) + { + p = buffer[j]; + + if (Volatile.Read(ref _limitReached) && p > maxP) + { + break; + } + + if (useFilter && !filter.Contains(p)) + { + continue; + } + + passedAllTests = IsEvenPerfectCandidate(p, divisorCyclesSearchLimit, out searchedMersenne, out detailedCheck); + PrintResult(p, searchedMersenne, detailedCheck, passedAllTests); + + if (p == maxP) + { + reachedMax = true; + } } - // TODO: Replace this byte-by-byte scan with the lookup-table based statistics collector validated in the - // BitStats benchmarks so zero runs leverage cached results instead of recomputing per bit. + if (reachedMax) + { + Volatile.Write(ref _limitReached, true); + break; + } + } + } + }); + } - zeroCountInByte = ByteZeroCount[inspectedByte]; - zeroCount += zeroCountInByte; + Task.WaitAll(tasks); + FlushBuffer(); + StringBuilderPool.Return(_outputBuilder!); + } - if (zeroCountInByte == 8) - { - currentRun += 8; - if (currentRun > maxZeroBlock) - { - maxZeroBlock = currentRun; - } + private static List LoadByDivisorCandidates(string candidateFile) + { + List candidates = []; + using FileStream readStream = new(candidateFile, FileMode.Open, FileAccess.Read, FileShare.ReadWrite); + using StreamReader reader = new(readStream); + string? line; + ReadOnlySpan span; + while ((line = reader.ReadLine()) is not null) + { + if (string.IsNullOrWhiteSpace(line)) + { + continue; + } - continue; - } + span = line.AsSpan().Trim(); + if (span.IsEmpty || span[0] == '#') + { + continue; + } - prefixZeros = BytePrefixZero[inspectedByte]; - suffixZeros = ByteSuffixZero[inspectedByte]; - maxZeroRunInByte = ByteMaxZeroRun[inspectedByte]; + // TODO: Swap this char-by-char parser for the benchmarked Utf8Parser-based fast-path once + // the CLI plumbing accepts spans to avoid redundant range checks while collecting primes. + int index = 0; + while (index < span.Length) + { + if (!char.IsDigit(span[index])) + { + index++; + continue; + } - candidate = currentRun + prefixZeros; - if (candidate > maxZeroBlock) - { - maxZeroBlock = candidate; - } + int start = index; + index++; + while (index < span.Length && char.IsDigit(span[index])) + { + index++; + } - if (maxZeroRunInByte > maxZeroBlock) - { - maxZeroBlock = maxZeroRunInByte; - } + // TODO: Swap ulong.TryParse for the Utf8Parser-based fast path once the shared span helpers + // arrive so by-divisor candidate loading avoids per-value string allocations. + if (ulong.TryParse(span[start..index], NumberStyles.None, CultureInfo.InvariantCulture, out ulong parsed)) + { + candidates.Add(parsed); + } + } + } - currentRun = suffixZeros; - } + return candidates; + } - if (currentRun > maxZeroBlock) - { - maxZeroBlock = currentRun; - } - } + private static void RunByDivisorMode(List candidates, int threadCount) + { + if (_byDivisorTester is null) + { + throw new InvalidOperationException("By-divisor tester has not been initialized."); + } - // LUTs for fast per-byte zero stats (MSB-first in each byte) - private static readonly byte[] ByteZeroCount = new byte[256]; - private static readonly byte[] BytePrefixZero = new byte[256]; - private static readonly byte[] ByteSuffixZero = new byte[256]; - private static readonly byte[] ByteMaxZeroRun = new byte[256]; + // TODO: Inline this helper by invoking MersenneNumberDivisorByDivisorTester.Run directly at call sites so + // the CLI avoids this extra frame during the hot composite-reporting path. + MersenneNumberDivisorByDivisorTester.Run( + candidates, + _byDivisorTester, + _byDivisorPreviousResults, + _byDivisorStartPrime, + static () => _lastCompositeP = true, + static () => _lastCompositeP = false, + PrintResult, + threadCount); + } - static Program() - { - int valueIndex = 0; - int zeros = 0; - int pref = 0; - int suff = 0; - int maxIn = 0; - int run = 0; - int bit = 0; - bool isZero = false; - int boundaryCounter = 0; - int boundaryBitIndex = 0; - for (; valueIndex < 256; valueIndex++) - { - zeros = 0; - pref = 0; - suff = 0; - maxIn = 0; - run = 0; + private static unsafe void LoadResultsFile(string resultsFileName, Action lineProcessorAction) + { + using FileStream readStream = new(resultsFileName, FileMode.Open, FileAccess.Read, FileShare.ReadWrite); + using StreamReader reader = new(readStream); + string? line; + bool headerSkipped = false; + ReadOnlySpan span; + int first = 0; + int second = 0; + int third = 0; + ulong parsedP = 0UL; + ReadOnlySpan detailedSpan = default; + ReadOnlySpan passedAllTestsSpan = default; + bool detailed = false; + bool passedAllTests = false; + while ((line = reader.ReadLine()) != null) + { + if (!headerSkipped) + { + headerSkipped = true; + continue; + } - // MSB-first within byte: bit 7 .. bit 0 - for (bit = 7; bit >= 0; bit--) - { - isZero = ((valueIndex >> bit) & 1) == 0; - if (isZero) - { - zeros++; - run++; - if (run > maxIn) - { - maxIn = run; - } - } - else - { - run = 0; - } + if (string.IsNullOrWhiteSpace(line)) + { + continue; + } - if (bit == 7) - { - // leading zeros (prefix) - boundaryCounter = 0; - for (boundaryBitIndex = 7; boundaryBitIndex >= 0; boundaryBitIndex--) - { - if (((valueIndex >> boundaryBitIndex) & 1) == 0) - { - boundaryCounter++; - } - else - { - break; - } - } - pref = boundaryCounter; - } + span = line.AsSpan(); + first = span.IndexOf(','); + if (first < 0) + { + continue; + } - if (bit == 0) - { - // trailing zeros (suffix) - // Reusing boundaryCounter to count trailing zeros after prefix detection above. - boundaryCounter = 0; - for (boundaryBitIndex = 0; boundaryBitIndex < 8; boundaryBitIndex++) - { - if (((valueIndex >> boundaryBitIndex) & 1) == 0) - { - boundaryCounter++; - } - else - { - break; - } - } - suff = boundaryCounter; - } - } + // TODO: Replace this ulong.Parse call with the Utf8Parser-based span fast-path we benchmarked so results + // reload skips the slower string allocation and culture-aware parsing. + parsedP = ulong.Parse(span[..first]); + span = span[(first + 1)..]; + second = span.IndexOf(','); + if (second < 0) + { + continue; + } - ByteZeroCount[valueIndex] = (byte)zeros; - BytePrefixZero[valueIndex] = (byte)pref; - ByteSuffixZero[valueIndex] = (byte)suff; - ByteMaxZeroRun[valueIndex] = (byte)maxIn; - } - } + span = span[(second + 1)..]; + third = span.IndexOf(','); + if (third < 0) + { + continue; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static bool IsCompositeByGcd(ulong p) - { - if (p < 2) - { - return true; - } + detailedSpan = span[..third]; + passedAllTestsSpan = span[(third + 1)..]; - ulong m = (ulong)BitOperations.Log2(p); - return BinaryGcd(p, m) > 1UL; - } + // TODO: Replace these bool.TryParse calls with the span-based fast path once we expose the + // Utf8Parser-powered helpers so results reload avoids per-line string allocations for booleans. + if (bool.TryParse(detailedSpan, out detailed) && bool.TryParse(passedAllTestsSpan, out passedAllTests)) + { + lineProcessorAction(parsedP, detailed, passedAllTests); + } + } + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong BinaryGcd(ulong a, ulong b) - { - if (a == 0UL) - { - return b; - } + private static void PrintHelp() + { + Console.WriteLine("Usage: EvenPerfectBitScanner [options]"); + Console.WriteLine(); + Console.WriteLine("Options:"); + Console.WriteLine(" --prime= starting exponent (p)"); + Console.WriteLine(" --increment=bit|add exponent increment method"); + Console.WriteLine(" --threads= number of worker threads"); + Console.WriteLine(" --block-size= values processed per thread batch"); + Console.WriteLine(" --mersenne=pow2mod|incremental|lucas|residue|divisor|bydivisor Mersenne test method"); + Console.WriteLine(" --residue-max-k= max k for residue Mersenne test (q = 2*p*k + 1)"); + Console.WriteLine(" --mersenne-device=cpu|gpu Device for Mersenne method (default gpu)"); + Console.WriteLine(" --primes-device=cpu|gpu Device for prime-scan kernels (default gpu)"); + Console.WriteLine(" --gpu-prime-batch= Batch size for GPU primality sieve (default 262144)"); + Console.WriteLine(" --order-device=cpu|gpu Device for order computations (default gpu)"); + Console.WriteLine(" --ntt=reference|staged GPU NTT backend (default staged)"); + Console.WriteLine(" --mod-reduction=auto|uint128|mont64|barrett128 staged NTT reduction (default auto)"); + Console.WriteLine(" --gpu-prime-threads= max concurrent GPU prime checks (default 1)"); + Console.WriteLine(" --ll-slice= Lucas–Lehmer iterations per slice (default 32)"); + Console.WriteLine(" --gpu-scan-batch= GPU q-scan batch size (default 2_097_152)"); + Console.WriteLine(" --order-warmup-limit= Warm-up order candidates (default 5_000_000)"); + Console.WriteLine(" --rle-blacklist= enable RLE blacklist for p (hard filter up to --rle-hard-max)"); + Console.WriteLine(" --rle-hard-max=

apply RLE blacklist only for p <= this (default ulong.MaxValue = no limit)"); + Console.WriteLine(" --rle-only-last7=true|false apply RLE only when p % 10 == 7 (default true)"); + Console.WriteLine(" --zero-hard= hard reject if zero_fraction(p) > f (default off)"); + Console.WriteLine(" --zero-conj=: hard reject if zero_fraction(p) > f AND max_zero_block >= r (default off)"); + Console.WriteLine(" --results-dir= directory for results file"); + Console.WriteLine(" --results-prefix= prefix to prepend to results filename"); + Console.WriteLine(" --divisor-cycles= divisor cycles data file"); + Console.WriteLine(" --divisor-cycles-device=cpu|gpu device for cycles generation (default gpu)"); + Console.WriteLine(" --divisor-cycles-batch= 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"); + Console.WriteLine(" --use-order test primality via q order"); + Console.WriteLine(" --workaround-mod avoid '%' operator on the GPU"); + // mod-automaton removed + Console.WriteLine(" --filter-p= process only p from previous run results (required for --mersenne=bydivisor)"); + Console.WriteLine(" --write-batch-size= overwrite frequency of disk writes (default 100 lines)"); + Console.WriteLine(" --gcd-filter enable early sieve based on GCD"); + Console.WriteLine(" --help, -help, --?, -?, /? show this help message"); + } - if (b == 0UL) - { - return a; - } + private static string BuildResultsFileName(bool bitInc, int threads, int block, GpuKernelType kernelType, bool useLucasFlag, bool useDivisorFlag, bool useByDivisorFlag, bool mersenneOnGpu, bool useOrder, bool useModWorkaround, bool useGcd, NttBackend nttBackend, int gpuPrimeThreads, int llSlice, int gpuScanBatch, ulong warmupLimit, ModReductionMode reduction, string mersenneDevice, string primesDevice, string orderDevice) + { + string inc = bitInc ? "bit" : "add"; + string mers = useDivisorFlag + ? "divisor" + : (useLucasFlag + ? "lucas" + : (useByDivisorFlag + ? "bydivisor" + : (kernelType == GpuKernelType.Pow2Mod ? "pow2mod" : "incremental"))); + string ntt = nttBackend == NttBackend.Staged ? "staged" : "reference"; + string red = reduction switch { ModReductionMode.Mont64 => "mont64", ModReductionMode.Barrett128 => "barrett128", ModReductionMode.GpuUInt128 => "uint128", _ => "auto" }; + string order = useOrder ? "order-on" : "order-off"; + string gcd = useGcd ? "gcd-on" : "gcd-off"; + string work = useModWorkaround ? "modfix-on" : "modfix-off"; + // TODO: Replace this string interpolation with the pooled ValueStringBuilder pipeline from the + // results-writer benchmarks so filename generation reuses the zero-allocation formatter once + // we consolidate output paths. + return $"even_perfect_bit_scan_inc-{inc}_thr-{threads}_blk-{block}_mers-{mers}_mersdev-{mersenneDevice}_ntt-{ntt}_red-{red}_primesdev-{primesDevice}_{order}_orderdev-{orderDevice}_{gcd}_{work}_gputh-{gpuPrimeThreads}_llslice-{llSlice}_scanb-{gpuScanBatch}_warm-{warmupLimit}.csv"; + } - int shift = BitOperations.TrailingZeroCount(a | b); - a >>= BitOperations.TrailingZeroCount(a); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static unsafe int ReserveBlock(ulong[] buffer, int blockSize) + { + Span bufferSpan = new(buffer); + ulong p, remainder; + long state, newState, original; + int count; + while (true) + { + state = Volatile.Read(ref _state); + p = (ulong)state >> 3; + remainder = ((ulong)state) & 7UL; - do - { - b >>= BitOperations.TrailingZeroCount(b); + count = 0; + while (count < blockSize && !Volatile.Read(ref _limitReached)) + { + bufferSpan[count++] = p; + // TODO: Inline the fast-path transform here once we collapse the delegate* indirection so + // block reservations stop paying the extra jump highlighted in the prime-stepping benchmarks. + p = _transformP(p, ref remainder); + } - if (a > b) - { - (a, b) = (b, a); - } + if (count == 0) + { + return 0; + } - b -= a; - } - while (b != 0UL); + newState = ((long)p << 3) | (long)remainder; + original = Interlocked.CompareExchange(ref _state, newState, state); + if (original == state) + { + return count; + } + } + } - return a << shift; - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool TryAcquireConsoleSlot(bool primeCandidate, out bool primeFlag) + { + int consoleInterval = PerfectNumberConstants.ConsoleInterval; + int counterValue = Interlocked.Increment(ref _consoleCounter); + if (counterValue >= consoleInterval) + { + int previous = Interlocked.Exchange(ref _consoleCounter, 0); + if (previous >= consoleInterval) + { + primeFlag = primeCandidate; + return true; + } + } + primeFlag = false; + return false; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static StringBuilder? AppendToOutputBuffer(StringBuilder source) + { + StringBuilder? builderToFlush = null; + + lock (Sync) + { + StringBuilder outputBuilder = _outputBuilder!; + _ = outputBuilder.Append(source); + + int nextIndex = _writeIndex + 1; + if (nextIndex >= _writeBatchSize) + { + builderToFlush = outputBuilder; + _outputBuilder = StringBuilderPool.Rent(); + _writeIndex = 0; + } + else + { + _writeIndex = nextIndex; + } + } + + return builderToFlush; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void PrintResult(ulong currentP, bool searchedMersenne, bool detailedCheck, bool passedAllTests) + { + if (passedAllTests) + { + int newCount = Interlocked.Increment(ref _primeCount); + if (newCount >= 2) + { + Volatile.Write(ref _primeFoundAfterInit, true); + } + } + + bool lastWasComposite = _lastCompositeP; + bool primeCandidate = passedAllTests && !lastWasComposite && Volatile.Read(ref _primeFoundAfterInit); + + // TODO: Switch this StringBuilder-based formatter to a span-friendly Utf8Formatter pipeline so + // console/file output avoids intermediate builders in the hot logging loop. + StringBuilder localBuilder = StringBuilderPool.Rent(); + _ = localBuilder + .Append(currentP).Append(',') + .Append(searchedMersenne).Append(',') + .Append(detailedCheck).Append(',') + .Append(passedAllTests).Append('\n'); + + bool printToConsole = TryAcquireConsoleSlot(primeCandidate, out bool primeFlag); + StringBuilder? builderToFlush = AppendToOutputBuffer(localBuilder); + + if (builderToFlush is not null) + { + FlushBuffer(builderToFlush); + } + + if (printToConsole) + { + if (primeFlag) + { + _ = localBuilder + .Remove(localBuilder.Length - 1, 1) + .Append(PrimeFoundSuffix) + .Append('\n'); + } + + Console.WriteLine(localBuilder.Remove(localBuilder.Length - 1, 1).ToString()); + } + + localBuilder.Clear(); + StringBuilderPool.Return(localBuilder); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void FlushBuffer() + { + StringBuilder? builderToFlush = null; + + lock (Sync) + { + if (_writeIndex == 0 || _outputBuilder is null || _outputBuilder.Length == 0) + { + return; + } + + builderToFlush = _outputBuilder; + _outputBuilder = StringBuilderPool.Rent(); + _writeIndex = 0; + } + + FlushBuffer(builderToFlush!); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void FlushBuffer(StringBuilder builderToFlush) + { + try + { + lock (FileWriteSync) + { + // TODO: Replace these per-flush FileStream/StreamWriter allocations with a pooled TextFileWriter-style + // helper so batched result writes keep the hot path on the benchmarked persistent-handle pipeline. + using FileStream stream = new(ResultsFileName, FileMode.Append, FileAccess.Write, FileShare.Read); + using StreamWriter writer = new(stream) { AutoFlush = false }; + writer.Write(builderToFlush); + writer.Flush(); + } + } + finally + { + builderToFlush.Clear(); + StringBuilderPool.Return(builderToFlush); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong CountOnes(ulong value) + { + // TODO: Remove this wrapper once callers can invoke BitOperations.PopCount directly so the hot bit-stat path + // avoids the extra call layer. + return (ulong)BitOperations.PopCount(value); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong GetNextAddDiff(ulong remainder) + { + return remainder switch + { + 0UL => 1UL, + 1UL => 4UL, + 2UL => 3UL, + 3UL => 2UL, + 4UL => 1UL, + _ => 2UL, + }; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong TransformPBit(ulong value, ref ulong remainder) + { + ulong original = value; + if (value > (ulong.MaxValue >> 1)) + { + Volatile.Write(ref _limitReached, true); + return value; + } + + ulong next = (value << 1) | 1UL; + remainder = (remainder << 1) + 1UL; + // Mod6 lookup loses to `%` in the benches; stick with subtraction + modulo for correctness and speed. + while (remainder >= 6UL) + { + remainder -= 6UL; + } + + value = remainder switch + { + 0UL => 1UL, + 1UL => 0UL, + 2UL => 5UL, + 3UL => 4UL, + 4UL => 3UL, + _ => 2UL, + }; // 'value' now holds diff + + if (next > ulong.MaxValue - value) + { + Volatile.Write(ref _limitReached, true); + return original; + } + + remainder += value; + // `% 6` remains faster per Mod6ComparisonBenchmarks; keep this modulo fold and continue using subtraction. + while (remainder >= 6UL) + { + remainder -= 6UL; + } + + return next + value; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong TransformPAdd(ulong value, ref ulong remainder) + { + ulong next = value; + ulong diff = GetNextAddDiff(remainder); + + if (next > ulong.MaxValue - diff) + { + Volatile.Write(ref _limitReached, true); + return next; + } + + remainder += diff; + // Retain direct modulo because the Mod6 helper underperforms on 64-bit operands. + while (remainder >= 6UL) + { + remainder -= 6UL; + } + + return next + diff; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong TransformPAddPrimes(ulong value, ref ulong remainder) + { + ulong originalRemainder = remainder; + ulong addRemainder = remainder; + ulong candidate = value; + ulong primeCandidate = value; + bool advanceAdd = true; + bool advancePrime = true; + ulong diff = 0UL; + ulong nextPrime = 0UL; + + while (true) + { + if (advanceAdd) + { + diff = GetNextAddDiff(addRemainder); + if (candidate > ulong.MaxValue - diff) + { + Volatile.Write(ref _limitReached, true); + remainder = originalRemainder; + return candidate; + } + + candidate += diff; + addRemainder += diff; + // Skip the Mod6 lookup: benchmarked `%` stays ahead for these prime increments. + // benchmarked fastest remainder updates instead of looping subtraction. + while (addRemainder >= 6UL) + { + addRemainder -= 6UL; + } + } + + if (advancePrime) + { + try + { + nextPrime = PrimeIterator.Next(in primeCandidate); + } + catch (InvalidOperationException) + { + Volatile.Write(ref _limitReached, true); + remainder = originalRemainder; + return value; + } + + if (nextPrime <= primeCandidate) + { + Volatile.Write(ref _limitReached, true); + remainder = originalRemainder; + return value; + } + + primeCandidate = nextPrime; + } + + if (candidate == primeCandidate) + { + remainder = addRemainder; + return candidate; + } + + advanceAdd = candidate < primeCandidate; + advancePrime = candidate > primeCandidate; + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static bool IsEvenPerfectCandidate(ulong p, ulong divisorCyclesSearchLimit) => + // TODO: Remove this pass-through once callers can provide the full out parameters + // directly; shaving this hop keeps the hot candidate filter on the shortest path. + 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; + + if (_useGcdFilter && IsCompositeByGcd(p)) + { + _lastCompositeP = true; + return false; + } + + // Optional: RLE blacklist and binary-threshold filters on p (safe only when configured) + if (p <= _rleHardMaxP) + { + if (!_rleOnlyLast7 || p.Mod10() == 7UL) + { + if (RleBlacklist.IsLoaded() && RleBlacklist.Matches(p)) + { + return false; + } + } + } + + if (_zeroFracHard >= 0 || (_zeroFracConj >= 0 && _maxZeroConj >= 0)) + { + double zf; + int bitLength = 0; + int zeroCountValue = 0; + int maxZeroBlockValue = 0; + ComputeBitStats(p, out bitLength, out zeroCountValue, out maxZeroBlockValue); + if (_zeroFracHard >= 0) + { + zf = (double)zeroCountValue / bitLength; + if (zf > _zeroFracHard) + { + return false; + } + } + + if (_zeroFracConj >= 0 && _maxZeroConj >= 0) + { + zf = (double)zeroCountValue / bitLength; + if (zf > _zeroFracConj && maxZeroBlockValue >= _maxZeroConj) + { + return false; + } + } + } + + // Fast residue-based composite check for p using small primes + if (!_useDivisor && !_useByDivisorMode && IsCompositeByResidues(p)) + { + searchedMersenne = false; + detailedCheck = false; + _lastCompositeP = true; + return false; + } + + // If primes-device=gpu, route p primality through GPU-assisted sieve with deterministic MR validation. + if (!_useByDivisorMode) + { + if (!(GpuContextPool.ForceCpu + ? PrimeTesters.Value!.IsPrime(p, CancellationToken.None) + : PrimeTesters.Value!.IsPrimeGpu(p, CancellationToken.None))) + { + _lastCompositeP = true; + return false; + } + } + + searchedMersenne = true; + if (_useByDivisorMode) + { + // TODO: Route by-divisor scans through the ProcessEightBitWindows pow2mod kernel once the windowed Pow2Minus1Mod helper replaces the slower single-bit ladder. + return _byDivisorTester!.IsPrime(p, out detailedCheck); + } + + detailedCheck = MersenneTesters.Value!.IsMersennePrime(p); + return detailedCheck; + } + + // Use ModResidueTracker with a small set of primes to pre-filter composite p. + private static bool IsCompositeByResidues(ulong p) + { + var tracker = PResidue.Value!; + // TODO: Integrate the divisor-cycle cache here so the small-prime sweep reuses precomputed remainders instead of + // running MergeOrAppend for every candidate and missing the cycle-accelerated early exits. + tracker.BeginMerge(p); + // Use the small prime list from PerfectNumbers.Core to the extent of sqrt(p) + var primes = PrimesGenerator.SmallPrimes; + var primesPow2 = PrimesGenerator.SmallPrimesPow2; + int len = primes.Length; + int primeIndex = 0; + bool divisible = false; + for (; primeIndex < len; primeIndex++) + { + if (primesPow2[primeIndex] > p) + { + break; + } + + if (tracker.MergeOrAppend(p, primes[primeIndex], out divisible) && divisible) + { + return true; + } + } + + return false; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void ComputeBitStats(ulong value, out int bitLen, out int zeroCount, out int maxZeroBlock) + { + bitLen = 64 - int.CreateChecked(ulong.LeadingZeroCount(value)); + zeroCount = 0; + maxZeroBlock = 0; + if (bitLen <= 0) + { + return; + } + + int msbIndex = (bitLen - 1) >> 3; // 0..7 + int bitsInTopByte = ((bitLen - 1) & 7) + 1; // 1..8 + int currentRun = 0; + byte inspectedByte = 0; + int zeroCountInByte = 0; + int prefixZeros = 0; + int suffixZeros = 0; + int maxZeroRunInByte = 0; + int candidate = 0; + int byteIndex = msbIndex; + + for (; byteIndex >= 0; byteIndex--) + { + inspectedByte = (byte)(value >> (byteIndex * 8)); + if (byteIndex == msbIndex && bitsInTopByte < 8) + { + // Mask off leading unused bits by setting them to 1 (ignore) + inspectedByte |= (byte)(0xFF << bitsInTopByte); + } + + // TODO: Replace this byte-by-byte scan with the lookup-table based statistics collector validated in the + // BitStats benchmarks so zero runs leverage cached results instead of recomputing per bit. + + zeroCountInByte = ByteZeroCount[inspectedByte]; + zeroCount += zeroCountInByte; + + if (zeroCountInByte == 8) + { + currentRun += 8; + if (currentRun > maxZeroBlock) + { + maxZeroBlock = currentRun; + } + + continue; + } + + prefixZeros = BytePrefixZero[inspectedByte]; + suffixZeros = ByteSuffixZero[inspectedByte]; + maxZeroRunInByte = ByteMaxZeroRun[inspectedByte]; + + candidate = currentRun + prefixZeros; + if (candidate > maxZeroBlock) + { + maxZeroBlock = candidate; + } + + if (maxZeroRunInByte > maxZeroBlock) + { + maxZeroBlock = maxZeroRunInByte; + } + + currentRun = suffixZeros; + } + + if (currentRun > maxZeroBlock) + { + maxZeroBlock = currentRun; + } + } + + // LUTs for fast per-byte zero stats (MSB-first in each byte) + private static readonly byte[] ByteZeroCount = new byte[256]; + private static readonly byte[] BytePrefixZero = new byte[256]; + private static readonly byte[] ByteSuffixZero = new byte[256]; + private static readonly byte[] ByteMaxZeroRun = new byte[256]; + + static Program() + { + int valueIndex = 0; + int zeros = 0; + int pref = 0; + int suff = 0; + int maxIn = 0; + int run = 0; + int bit = 0; + bool isZero = false; + int boundaryCounter = 0; + int boundaryBitIndex = 0; + for (; valueIndex < 256; valueIndex++) + { + zeros = 0; + pref = 0; + suff = 0; + maxIn = 0; + run = 0; + + // MSB-first within byte: bit 7 .. bit 0 + for (bit = 7; bit >= 0; bit--) + { + isZero = ((valueIndex >> bit) & 1) == 0; + if (isZero) + { + zeros++; + run++; + if (run > maxIn) + { + maxIn = run; + } + } + else + { + run = 0; + } + + if (bit == 7) + { + // leading zeros (prefix) + boundaryCounter = 0; + for (boundaryBitIndex = 7; boundaryBitIndex >= 0; boundaryBitIndex--) + { + if (((valueIndex >> boundaryBitIndex) & 1) == 0) + { + boundaryCounter++; + } + else + { + break; + } + } + pref = boundaryCounter; + } + + if (bit == 0) + { + // trailing zeros (suffix) + // Reusing boundaryCounter to count trailing zeros after prefix detection above. + boundaryCounter = 0; + for (boundaryBitIndex = 0; boundaryBitIndex < 8; boundaryBitIndex++) + { + if (((valueIndex >> boundaryBitIndex) & 1) == 0) + { + boundaryCounter++; + } + else + { + break; + } + } + suff = boundaryCounter; + } + } + + ByteZeroCount[valueIndex] = (byte)zeros; + BytePrefixZero[valueIndex] = (byte)pref; + ByteSuffixZero[valueIndex] = (byte)suff; + ByteMaxZeroRun[valueIndex] = (byte)maxIn; + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static bool IsCompositeByGcd(ulong p) + { + if (p < 2) + { + return true; + } + + ulong m = (ulong)BitOperations.Log2(p); + return BinaryGcd(p, m) > 1UL; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong BinaryGcd(ulong a, ulong b) + { + if (a == 0UL) + { + return b; + } + + if (b == 0UL) + { + return a; + } + + int shift = BitOperations.TrailingZeroCount(a | b); + a >>= BitOperations.TrailingZeroCount(a); + + do + { + b >>= BitOperations.TrailingZeroCount(b); + + if (a > b) + { + (a, b) = (b, a); + } + + b -= a; + } + while (b != 0UL); + + return a << shift; + } } + diff --git a/PerfectNumbers.Core.Tests/EulerPrimeTesterTests.cs b/PerfectNumbers.Core.Tests/EulerPrimeTesterTests.cs deleted file mode 100644 index 20f587b..0000000 --- a/PerfectNumbers.Core.Tests/EulerPrimeTesterTests.cs +++ /dev/null @@ -1,22 +0,0 @@ -using FluentAssertions; -using Xunit; - -namespace PerfectNumbers.Core.Tests; - -[Trait("Category", "Fast")] -public class EulerPrimeTesterTests -{ - private readonly EulerPrimeTester _tester = new(new PrimeTester()); - - [Theory] - [InlineData(13UL, true)] - [InlineData(29UL, true)] - [InlineData(3UL, false)] - [InlineData(21UL, false)] - [InlineData(1UL, false)] - public void IsEulerPrime_returns_expected(ulong n, bool expected) - { - _tester.IsEulerPrime(n).Should().Be(expected); - } -} - diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index 8dc9343..8aaefb1 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -110,7 +110,6 @@ public class MersenneNumberDivisorGpuTesterTests { var tester = new MersenneNumberDivisorByDivisorGpuTester { - UseDivisorCycles = true, GpuBatchSize = 8, }; @@ -147,7 +146,6 @@ public class MersenneNumberDivisorGpuTesterTests { var tester = new MersenneNumberDivisorByDivisorGpuTester { - UseDivisorCycles = true, GpuBatchSize = 2, }; @@ -169,7 +167,6 @@ public class MersenneNumberDivisorGpuTesterTests { var tester = new MersenneNumberDivisorByDivisorCpuTester { - UseDivisorCycles = true, BatchSize = 2, }; diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index c830d8e..9370406 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -12,14 +12,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private ulong _divisorLimit; private ulong _lastStatusDivisor; private bool _isConfigured; - private bool _useDivisorCycles; private int _batchSize = 1_024; - public bool UseDivisorCycles - { - get => _useDivisorCycles; - set => _useDivisorCycles = value; // TODO: Delete the toggle once divisor cycle data is always consulted so CPU scans never fall back to the slower pure-Montgomery path. - } public int BatchSize { @@ -69,8 +63,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv public bool IsPrime(ulong prime, out bool divisorsExhausted) { ulong allowedMax; - bool useCycles; - lock (_sync) { if (!_isConfigured) @@ -79,7 +71,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); - useCycles = _useDivisorCycles; } if (allowedMax < 3UL) @@ -93,12 +84,11 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv bool processedAll = false; bool composite = CheckDivisors( - prime, - allowedMax, - useCycles, - out lastProcessed, - out processedAll, - out processedCount); + prime, + allowedMax, + out lastProcessed, + out processedAll, + out processedCount); if (processedCount > 0UL) { @@ -168,12 +158,11 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } private bool CheckDivisors( - ulong prime, - ulong allowedMax, - bool useCycles, - out ulong lastProcessed, - out bool processedAll, - out ulong processedCount) + ulong prime, + ulong allowedMax, + out ulong lastProcessed, + out bool processedAll, + out ulong processedCount) { lastProcessed = 0UL; processedCount = 0UL; @@ -194,35 +183,15 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv // 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. - byte hit; - ulong divisorCycle = 0UL; - DivisorCycleCache.CycleBlock? cycleBlock = null; - if (useCycles) + ulong divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); + if (divisorCycle != 0UL) { - cycleBlock = DivisorCycleCache.Shared.Acquire(divisor); - // TODO: Swap this block lease for the direct single-cycle lookup once the cache exposes it - // so the CPU scanner mirrors the benchmarked single-block policy without retaining mutable - // blocks beyond the startup snapshot. - divisorCycle = cycleBlock.GetCycle(divisor); - // TODO: When divisorCycle is zero (cache miss), compute only that single cycle on the device - // selected by the current settings, skip inserting the result into the shared cache, and keep - // operating with the single on-disk block without requesting additional blocks. - hit = CheckDivisor(prime, divisorCycle != 0UL, divisorCycle, divisorData); - } - else - { - // TODO: Remove this no-cycle branch and require on-demand cycle computation (without - // mutating the cache) so CPU scans never fall back to the slower Montgomery-only path - // identified in the benchmarks. - hit = CheckDivisor(prime, false, 0UL, divisorData); - } - - cycleBlock?.Dispose(); - - if (hit != 0) - { - processedAll = true; - return true; + byte hit = CheckDivisor(prime, divisorCycle, divisorData); + if (hit != 0) + { + processedAll = true; + return true; + } } if (divisor >= allowedMax - 1UL) @@ -249,22 +218,18 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } - private static byte CheckDivisor(ulong prime, bool useCycles, ulong divisorCycle, in MontgomeryDivisorData divisorData) + private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) + { + ulong modulus = divisorData.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL || divisorCycle == 0UL) { - ulong modulus = divisorData.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return 0; - } - - ulong residue = useCycles && divisorCycle != 0UL - ? prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData) - : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the - // scalar Pow2Minus1Mod adopts it so CPU scans match the GPU speedups. - - return residue == 1UL ? (byte)1 : (byte)0; + return 0; } + ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData); + return residue == 1UL ? (byte)1 : (byte)0; + } + private void UpdateStatusUnsafe(ulong lastProcessed, ulong processedCount) { if (processedCount == 0UL) @@ -343,23 +308,22 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return; } - ulong modulus = divisor; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) + if (divisor <= 1UL || (divisor & 1UL) == 0UL) { hits.Clear(); return; } - bool useDivisorCycles = _owner._useDivisorCycles; - if (useDivisorCycles && divisorCycle == 0UL) + if (divisorCycle == 0UL) { - // TODO: Replace this guard with an on-demand CPU/GPU cycle computation path that skips - // cache insertion and continues operating with the single cached block instead of - // requesting additional cycle batches when a lookup misses. - throw new InvalidOperationException($"Missing divisor cycle for divisor {divisor}."); + divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); + if (divisorCycle == 0UL) + { + hits.Clear(); + return; + } } - bool cycleEnabled = useDivisorCycles; MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); var exponentStepper = new ExponentRemainderStepper(divisorData); @@ -369,39 +333,25 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return; } - if (cycleEnabled) + 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++) { - 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) { - remainder = stepper.ComputeNext(primes[i]); - if (remainder != 0UL) - { - hits[i] = 0; - continue; - } - - hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0; + hits[i] = 0; + continue; } - return; - } - - // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here - // keeps the slower Montgomery stepping alive even though the benchmarks showed the cached cycle path is far - // faster for large divisor sets. - for (int i = 0; i < length; i++) - { hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0; } } - public void Dispose() { if (_disposed) @@ -415,3 +365,4 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } + diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index 4599100..16db1f0 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -1,10 +1,10 @@ using System; using System.Buffers; using System.Collections.Concurrent; +using System.Collections.Generic; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using System.Threading; -using System.Threading.Tasks; using ILGPU; using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; @@ -13,535 +13,232 @@ namespace PerfectNumbers.Core; public sealed class DivisorCycleCache { - public sealed class CycleBlock : IDisposable - { - private readonly DivisorCycleCache _owner; - private int _referenceCount; - - internal CycleBlock(DivisorCycleCache owner, int index, ulong start, in ulong[] cycles) - { - _owner = owner; - Index = index; - Start = start; - End = start + (ulong)cycles.Length - 1UL; - Cycles = cycles; - } - - internal int Index { get; } - - internal ulong Start { get; } - - internal ulong End { get; } - - internal ulong[] Cycles { get; } - - internal int ReferenceCount => Volatile.Read(ref _referenceCount); - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public ulong GetCycle(ulong divisor) => Cycles[(int)(divisor - Start)]; - - internal void Retain() - { - Interlocked.Increment(ref _referenceCount); - } - - public void Dispose() - { - int remaining = Interlocked.Decrement(ref _referenceCount); - if (remaining < 0) - { - throw new InvalidOperationException("Cycle block released more times than it was acquired."); - } - - if (remaining == 0) - { - _owner.OnCycleBlockReleased(this); - } - } - } - - private const int CycleGenerationBatchSize = 262_144; // TODO: Remove the block-sized batch once the single-cycle generator - // replaces block hydration so we stop materializing large buffers - // just to satisfy one missing divisor lookup. private const byte ByteZero = 0; private const byte ByteOne = 1; + private const int StackBufferThreshold = 128; - private readonly object _sync = new(); // TODO: Delete the synchronization object once the cache exposes only snapshot reads and single-cycle computation so callers no longer contend on locks. - private readonly ConcurrentDictionary> _pending = new(); // TODO: Drop the concurrent dictionary entirely when the block pipeline is replaced with single-cycle generation; the new plan never queues background work. - private readonly ConcurrentDictionary, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>> _gpuKernelCache = new(AcceleratorReferenceComparer.Instance); // TODO: Swap this ConcurrentDictionary for a simple accelerator-indexed array once the single-cycle helper precompiles kernels during startup, eliminating the need for thread-safe lookups. - private CycleBlock _baseBlock = null!; - private CycleBlock? _activeBlock; // TODO: Remove dynamic blocks altogether once lookups route to the snapshot or compute an individual cycle without mutating shared state. - private CycleBlock? _prefetchedBlock; // TODO: Delete prefetch support so the cache never materializes extra blocks and callers always compute the missing cycle inline. - private bool _useGpuGeneration = true; - private readonly int _divisorCyclesBatchSize; - private static int _sharedDivisorCyclesBatchSize = GpuConstants.GpuCycleStepsPerInvocation; + private readonly ConcurrentDictionary, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>> _gpuKernelCache = new(AcceleratorReferenceComparer.Instance); + private volatile ulong[] _snapshot; + private volatile bool _useGpuGeneration = true; + private readonly int _divisorCyclesBatchSize; + private static int _sharedDivisorCyclesBatchSize = GpuConstants.GpuCycleStepsPerInvocation; + + public static void SetDivisorCyclesBatchSize(int divisorCyclesBatchSize) + { + _sharedDivisorCyclesBatchSize = Math.Max(1, divisorCyclesBatchSize); + } - public static void SetDivisorCyclesBatchSize(int divisorCyclesBatchSize) => _sharedDivisorCyclesBatchSize = divisorCyclesBatchSize; public static DivisorCycleCache Shared { get; } = new DivisorCycleCache(_sharedDivisorCyclesBatchSize); + public int PreferredBatchSize => _divisorCyclesBatchSize; + private DivisorCycleCache(int divisorCyclesBatchSize) { - _divisorCyclesBatchSize = divisorCyclesBatchSize; - ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); - Initialize(snapshot); - } - - // TODO: Remove the ability to reload at runtime once the cache becomes a read-only snapshot plus - // transient single-cycle generator; future callers should load the snapshot during startup only. - public void ReloadFromCurrentSnapshot() - { - ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); - Initialize(snapshot); + _divisorCyclesBatchSize = Math.Max(1, divisorCyclesBatchSize); + _snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); } public void ConfigureGeneratorDevice(bool useGpu) { - lock (_sync) - { - _useGpuGeneration = useGpu; - } - // TODO: Remove the lock once single-cycle generation eliminates shared mutable state; a simple - // volatile write will be enough when the cache becomes read-only. + _useGpuGeneration = useGpu; } - private void Initialize(ulong[] snapshot) + public void RefreshSnapshot() { - CycleBlock? activeToDispose = null; - CycleBlock? prefetchedToDispose = null; + _snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); + } - lock (_sync) + public ulong GetCycleLength(ulong divisor) + { + Span result = stackalloc ulong[1]; + ReadOnlySpan singleDivisor = stackalloc ulong[1] { divisor }; + GetCycleLengths(singleDivisor, result); + return result[0]; + } + + public void GetCycleLengths(ReadOnlySpan divisors, Span cycles) + { + ulong[] snapshot = _snapshot; + int length = divisors.Length; + int[]? rentedMissing = null; + Span missingBuffer = length <= StackBufferThreshold + ? stackalloc int[length] + : new Span(rentedMissing = ArrayPool.Shared.Rent(length), 0, length); + + int missingCount = 0; + + for (int i = 0; i < length; i++) { - if (_activeBlock is not null && _activeBlock.Index != 0) + ulong divisor = divisors[i]; + if (divisor <= 1UL) { - activeToDispose = _activeBlock; + throw new InvalidDataException("Divisor must be > 1"); } - if (_prefetchedBlock is not null && _prefetchedBlock.Index != 0) + if (divisor < (ulong)snapshot.Length) { - prefetchedToDispose = _prefetchedBlock; - } - - _pending.Clear(); - // TODO: Remove the dictionary once prefetching and background block tracking disappear; with single-cycle generation we will keep only the snapshot and skip additional bookkeeping entirely. - _baseBlock = new CycleBlock(this, index: 0, start: 0UL, snapshot); - _activeBlock = null; - _prefetchedBlock = null; - StartPrefetchLocked(1); - } - - activeToDispose?.Dispose(); - prefetchedToDispose?.Dispose(); - } - - public CycleBlock Acquire(ulong divisor) - { - CycleBlock block = divisor <= _baseBlock.End ? _baseBlock : AcquireDynamicBlock(divisor); - // TODO: Replace this block-based acquisition with a direct GetCycle API that computes the missing cycle synchronously on the configured device so the cache never hands out mutable blocks or requires disposal. - block.Retain(); - return block; - } - - private CycleBlock AcquireDynamicBlock(ulong divisor) - { - int blockIndex = GetBlockIndex(divisor); - - while (true) - { - Task? pendingTask = null; - lock (_sync) - { - if (TryGetCachedBlockLocked(blockIndex, out CycleBlock? cached)) + ulong cached = snapshot[divisor]; + if (cached == 0UL) { - return cached!; + throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); } - if (!_pending.TryGetValue(blockIndex, out pendingTask)) - { - // TODO: Replace Task.Run with an inline single-cycle computation that runs on the configured - // device so we materialize exactly one missing cycle and stop scheduling thread-pool work for - // entire blocks. - pendingTask = Task.Run(() => GenerateBlock(blockIndex)); - _pending[blockIndex] = pendingTask; - } + cycles[i] = cached; } - - CycleBlock block = pendingTask.GetAwaiter().GetResult(); - lock (_sync) + else { - _pending.TryRemove(blockIndex, out _); - - if (TryGetCachedBlockLocked(blockIndex, out CycleBlock? cachedAfterAwait)) - { - return cachedAfterAwait!; - } - - SetActiveBlockLocked(block); - // TODO: Stop promoting freshly generated blocks into the shared cache and, instead of - // returning a block, surface the computed cycle directly to the caller so no extra cache - // state survives beyond the single snapshot block. - return block; + missingBuffer[missingCount++] = i; } } - } - private bool TryGetCachedBlockLocked(int blockIndex, out CycleBlock? block) - { - CycleBlock? active = _activeBlock; - if (active is not null && active.Index == blockIndex) + if (missingCount != 0) { - block = active; - return true; - } - - CycleBlock? prefetched = _prefetchedBlock; - if (prefetched is not null && prefetched.Index == blockIndex) - { - PromotePrefetchedLocked(); - block = _activeBlock!; - return true; - } - - block = null; - return false; - } - - private void PromotePrefetchedLocked() - { - CycleBlock? prefetched = _prefetchedBlock; - if (prefetched is null) - { - return; - } - - CycleBlock? previousActive = _activeBlock; - _activeBlock = prefetched; - _prefetchedBlock = null; - prefetched.Retain(); - prefetched.Dispose(); - - if (previousActive is not null && previousActive.Index != 0 && previousActive.Index != prefetched.Index) - { - previousActive.Dispose(); - } - - StartPrefetchLocked(prefetched.Index + 1); // TODO: Remove this call once prefetching is dropped - // so we strictly keep a single active block in memory. - } - - private void SetActiveBlockLocked(CycleBlock block) - { - CycleBlock? previousActive = _activeBlock; - _activeBlock = block; - block.Retain(); - - CycleBlock? prefetched = _prefetchedBlock; - if (prefetched is not null && prefetched.Index != block.Index + 1) - { - if (prefetched.Index != 0) + ReadOnlySpan missingIndices = missingBuffer[..missingCount]; + if (_useGpuGeneration) { - prefetched.Dispose(); + ComputeCyclesGpu(divisors, cycles, missingIndices); } - - _prefetchedBlock = null; - } - - if (previousActive is not null && previousActive.Index != 0 && previousActive.Index != block.Index) - { - previousActive.Dispose(); - } - - StartPrefetchLocked(block.Index + 1); // TODO: Remove once prefetching is removed and the cache - // no longer speculatively requests the next block after a miss. - } - - private void StartPrefetchLocked(int blockIndex) - { - // TODO: Remove prefetching entirely so the cache maintains a single snapshot block and computes - // individual cycles on demand instead of triggering background generation for future indices. - if (blockIndex <= 0) - { - return; - } - - CycleBlock? prefetched = _prefetchedBlock; - if (prefetched is not null && prefetched.Index == blockIndex) - { - return; - } - - if (_pending.TryGetValue(blockIndex, out Task? existing)) - { - if (existing.IsCompletedSuccessfully) + else { - _pending.TryRemove(blockIndex, out _); - AssignPrefetchedBlockLocked(blockIndex, existing.Result); + ComputeCyclesCpu(divisors, cycles, missingIndices); } - - return; } - // TODO: Delete this once prefetching is removed; single-cycle generation will happen inline and - // no asynchronous work should be scheduled. - Task task = Task.Run(() => GenerateBlock(blockIndex)); - _pending[blockIndex] = task; - task.ContinueWith(t => OnPrefetchCompleted(blockIndex, t), TaskScheduler.Default); + if (rentedMissing is not null) + { + ArrayPool.Shared.Return(rentedMissing, clearArray: false); + } } - - private void OnPrefetchCompleted(int blockIndex, Task task) + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void ComputeCyclesCpu(ReadOnlySpan divisors, Span cycles, ReadOnlySpan indices) { - // TODO: Delete this callback when we remove prefetching. Single-cycle generation will run inline, - // eliminating the need for asynchronous continuation handlers and concurrent dictionaries. - if (!task.IsCompletedSuccessfully) + for (int i = 0; i < indices.Length; i++) { - _pending.TryRemove(blockIndex, out _); - return; - } - - lock (_sync) - { - _pending.TryRemove(blockIndex, out _); - AssignPrefetchedBlockLocked(blockIndex, task.Result); + int targetIndex = indices[i]; + cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLength(divisors[targetIndex]); } } - private void AssignPrefetchedBlockLocked(int blockIndex, CycleBlock block) + private void ComputeCyclesGpu(ReadOnlySpan divisors, Span cycles, ReadOnlySpan indices) { - // TODO: Remove this entire method when prefetching goes away. The single-cycle flow should never - // allocate or retain extra blocks beyond the startup snapshot, so prefetched block bookkeeping - // becomes unnecessary. - CycleBlock? active = _activeBlock; - if (active is not null && blockIndex != active.Index + 1) + int count = indices.Length; + ulong[]? rentedDivisors = null; + Span workDivisors = count <= StackBufferThreshold + ? stackalloc ulong[count] + : new Span(rentedDivisors = ArrayPool.Shared.Rent(count), 0, count); + + ulong[]? rentedResults = null; + Span workResults = count <= StackBufferThreshold + ? stackalloc ulong[count] + : new Span(rentedResults = ArrayPool.Shared.Rent(count), 0, count); + + for (int i = 0; i < count; i++) { - MontgomeryDivisorDataCache.ReleaseBlock(block); - return; + workDivisors[i] = divisors[indices[i]]; } - CycleBlock? previousPrefetched = _prefetchedBlock; - if (ReferenceEquals(previousPrefetched, block)) + ComputeCyclesGpuCore(workDivisors, workResults); + + for (int i = 0; i < count; i++) { - return; + cycles[indices[i]] = workResults[i]; } - if (_prefetchedBlock is null || _prefetchedBlock.Index <= blockIndex) + if (rentedDivisors is not null) { - if (previousPrefetched is not null && previousPrefetched.Index != 0) - { - previousPrefetched.Dispose(); - } - - _prefetchedBlock = block; - block.Retain(); - } - else - { - MontgomeryDivisorDataCache.ReleaseBlock(block); + ArrayPool.Shared.Return(rentedDivisors, clearArray: false); + ArrayPool.Shared.Return(rentedResults!, clearArray: false); } } - private void OnCycleBlockReleased(CycleBlock block) + private void ComputeCyclesGpuCore(ReadOnlySpan divisors, Span destination) { - if (block.Index == 0) - { - return; - } - - for (int attempt = 0; attempt < 8; attempt++) - { - if (block.ReferenceCount != 0) - { - return; - } - - Thread.SpinWait(1 << attempt); - } - - bool stillTracked; - - lock (_sync) - { - stillTracked = ReferenceEquals(_activeBlock, block) || ReferenceEquals(_prefetchedBlock, block); - } - - if (!stillTracked && block.ReferenceCount == 0) - { - MontgomeryDivisorDataCache.ReleaseBlock(block); - } - } - - private CycleBlock GenerateBlock(int blockIndex) - { - if (blockIndex == 0) - { - return _baseBlock; - } - - ulong start = GetBlockStart(blockIndex); - int length = _baseBlock.Cycles.Length; - ulong[] cycles = new ulong[length]; - - if (_useGpuGeneration) - { - ComputeCyclesGpu(start, cycles, _divisorCyclesBatchSize); - } - else - { - ComputeCyclesCpu(start, cycles); - } - - // TODO: Replace this block generator with a single-cycle helper that reuses shared buffers, - // returns just the requested divisor's cycle, and avoids mutating cache state beyond the - // preloaded snapshot. - return new CycleBlock(this, blockIndex, start, cycles); - } - - private int GetBlockIndex(ulong divisor) - { - if (divisor <= _baseBlock.End) - { - return 0; - } - - ulong offset = divisor - (_baseBlock.End + 1UL); - return (int)(offset / (ulong)PerfectNumberConstants.MaxQForDivisorCycles) + 1; - } - - private ulong GetBlockStart(int blockIndex) - { - if (blockIndex == 0) - { - return 0UL; - } - - ulong baseEnd = _baseBlock.End; - return baseEnd + 1UL + (ulong)(blockIndex - 1) * (ulong)PerfectNumberConstants.MaxQForDivisorCycles; - } - - private void ComputeCyclesCpu(ulong start, ulong[] destination) - { - ulong divisor = start; - for (int i = 0; i < destination.Length; i++) - { - // TODO: Port this CPU fallback to the unrolled-hex cycle generator once it is shared so misses outside the snapshot - // stop looping through CalculateCycleLength, which benchmarks showed is far slower than the GPU-assisted pipeline. - destination[i] = MersenneDivisorCycles.CalculateCycleLength(divisor); - divisor++; - } - } - - private void ComputeCyclesGpu(ulong start, ulong[] destination, int divisorCyclesBatchSize) - { - var gpuLease = GpuKernelPool.GetKernel(_useGpuGeneration); + int length = divisors.Length; + var gpuLease = GpuKernelPool.GetKernel(useGpuOrder: true); var execution = gpuLease.EnterExecutionScope(); - try + Accelerator accelerator = gpuLease.Accelerator; + var kernel = _gpuKernelCache.GetOrAdd(accelerator, LoadKernel); + + using MemoryBuffer1D divisorBuffer = accelerator.Allocate1D(length); + using MemoryBuffer1D powBuffer = accelerator.Allocate1D(length); + using MemoryBuffer1D orderBuffer = accelerator.Allocate1D(length); + using MemoryBuffer1D resultBuffer = accelerator.Allocate1D(length); + using MemoryBuffer1D statusBuffer = accelerator.Allocate1D(length); + + ulong[]? rentedPow = null; + Span powSpan = length <= StackBufferThreshold + ? stackalloc ulong[length] + : new Span(rentedPow = ArrayPool.Shared.Rent(length), 0, length); + + ulong[]? rentedOrder = null; + Span orderSpan = length <= StackBufferThreshold + ? stackalloc ulong[length] + : new Span(rentedOrder = ArrayPool.Shared.Rent(length), 0, length); + + ulong[]? rentedResult = null; + Span resultSpan = length <= StackBufferThreshold + ? stackalloc ulong[length] + : new Span(rentedResult = ArrayPool.Shared.Rent(length), 0, length); + + byte[]? rentedStatus = null; + Span statusSpan = length <= StackBufferThreshold + ? stackalloc byte[length] + : new Span(rentedStatus = ArrayPool.Shared.Rent(length), 0, length); + + ref ulong divisorRef = ref MemoryMarshal.GetReference(divisors); + divisorBuffer.View.CopyFromCPU(ref divisorRef, length); + + for (int i = 0; i < length; i++) { - Accelerator accelerator = gpuLease.Accelerator; - Action, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> kernel = _gpuKernelCache.GetOrAdd(accelerator, LoadKernel); - // TODO: Once single-cycle generation is in place, bypass this block-sized kernel entirely and - // call the ProcessEightBitWindows-based GPU helper that materializes only the requested - // divisor, reusing a pinned upload buffer instead of allocating per block. - - int length = destination.Length; - int chunkCapacity = Math.Min(length, CycleGenerationBatchSize); - ulong[] divisors = ArrayPool.Shared.Rent(chunkCapacity); - ulong[] powScratch = ArrayPool.Shared.Rent(chunkCapacity); - ulong[] orderScratch = ArrayPool.Shared.Rent(chunkCapacity); - byte[] statusScratch = ArrayPool.Shared.Rent(chunkCapacity); - - try + ulong divisor = divisors[i]; + ulong initialPow = 2UL; + if (initialPow >= divisor) { - using MemoryBuffer1D divisorBuffer = accelerator.Allocate1D(chunkCapacity); - using MemoryBuffer1D powBuffer = accelerator.Allocate1D(chunkCapacity); - using MemoryBuffer1D orderBuffer = accelerator.Allocate1D(chunkCapacity); - using MemoryBuffer1D resultBuffer = accelerator.Allocate1D(chunkCapacity); - using MemoryBuffer1D statusBuffer = accelerator.Allocate1D(chunkCapacity); + initialPow -= divisor; + } - int offset = 0; - while (offset < length) + powSpan[i] = initialPow; + orderSpan[i] = 1UL; + resultSpan[i] = 0UL; + statusSpan[i] = ByteZero; + } + + powBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(powSpan), length); + orderBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(orderSpan), length); + resultBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(resultSpan), length); + statusBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(statusSpan), length); + + int pending; + do + { + kernel(length, _divisorCyclesBatchSize, divisorBuffer.View, powBuffer.View, orderBuffer.View, resultBuffer.View, statusBuffer.View); + accelerator.Synchronize(); + + statusBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(statusSpan), length); + pending = 0; + for (int i = 0; i < length; i++) + { + if (statusSpan[i] == ByteZero) { - int batchSize = Math.Min(chunkCapacity, length - offset); - Span divisorsSpan = divisors.AsSpan(0, batchSize); - Span destinationSpan = destination.AsSpan(offset, batchSize); - Span powSpan = powScratch.AsSpan(0, batchSize); - Span orderSpan = orderScratch.AsSpan(0, batchSize); - Span statusSpan = statusScratch.AsSpan(0, batchSize); - - int pending = 0; - for (int i = 0; i < batchSize; i++) - { - ulong divisorValue = start + (ulong)(offset + i); - divisorsSpan[i] = divisorValue; - - if ((divisorValue & (divisorValue - 1UL)) == 0UL) - { - destinationSpan[i] = 1UL; - powSpan[i] = 0UL; - orderSpan[i] = 1UL; - statusSpan[i] = ByteOne; - } - else - { - ulong initialPow = 2UL; - if (initialPow >= divisorValue) - { - initialPow -= divisorValue; - } - - destinationSpan[i] = 0UL; - powSpan[i] = initialPow; - orderSpan[i] = 1UL; - statusSpan[i] = ByteZero; - pending++; - } - } - - ArrayView1D divisorView = divisorBuffer.View.SubView(0, batchSize); - ArrayView1D powView = powBuffer.View.SubView(0, batchSize); - ArrayView1D orderView = orderBuffer.View.SubView(0, batchSize); - ArrayView1D resultView = resultBuffer.View.SubView(0, batchSize); - ArrayView1D statusView = statusBuffer.View.SubView(0, batchSize); - - divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorsSpan), batchSize); - powView.CopyFromCPU(ref MemoryMarshal.GetReference(powSpan), batchSize); - orderView.CopyFromCPU(ref MemoryMarshal.GetReference(orderSpan), batchSize); - resultView.CopyFromCPU(ref MemoryMarshal.GetReference(destinationSpan), batchSize); - statusView.CopyFromCPU(ref MemoryMarshal.GetReference(statusSpan), batchSize); - - while (pending > 0) - { - kernel(batchSize, divisorCyclesBatchSize, divisorView, powView, orderView, resultView, statusView); - accelerator.Synchronize(); - - statusView.CopyToCPU(ref MemoryMarshal.GetReference(statusSpan), batchSize); - - pending = 0; - for (int i = 0; i < batchSize; i++) - { - if (statusSpan[i] == ByteZero) - { - pending++; - } - } - } - - resultView.CopyToCPU(ref MemoryMarshal.GetReference(destinationSpan), batchSize); - - offset += batchSize; + pending++; } } - finally - { - ArrayPool.Shared.Return(divisors, clearArray: false); - ArrayPool.Shared.Return(powScratch, clearArray: false); - ArrayPool.Shared.Return(orderScratch, clearArray: false); - ArrayPool.Shared.Return(statusScratch, clearArray: false); - } } - finally + while (pending > 0); + + resultBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(resultSpan), length); + resultSpan.CopyTo(destination); + if (rentedPow is not null) { - execution.Dispose(); - gpuLease.Dispose(); + ArrayPool.Shared.Return(rentedPow, clearArray: false); + ArrayPool.Shared.Return(rentedOrder!, clearArray: false); + ArrayPool.Shared.Return(rentedResult!, clearArray: false); + ArrayPool.Shared.Return(rentedStatus!, clearArray: false); } + + execution.Dispose(); + gpuLease.Dispose(); } private static Action, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> LoadKernel(Accelerator accelerator) @@ -551,7 +248,7 @@ public sealed class DivisorCycleCache private static void GpuAdvanceDivisorCyclesKernel( Index1D index, - int steps, + int steps, ArrayView1D divisors, ArrayView1D pow, ArrayView1D order, @@ -602,3 +299,6 @@ public sealed class DivisorCycleCache } } + + + diff --git a/PerfectNumbers.Core/EulerPrimeTester.cs b/PerfectNumbers.Core/EulerPrimeTester.cs deleted file mode 100644 index cd3e36b..0000000 --- a/PerfectNumbers.Core/EulerPrimeTester.cs +++ /dev/null @@ -1,22 +0,0 @@ -using System.Runtime.CompilerServices; - -namespace PerfectNumbers.Core; - -public sealed class EulerPrimeTester -{ - private readonly PrimeTester _primeTester; - - public EulerPrimeTester(PrimeTester primeTester) - { - _primeTester = primeTester; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public bool IsEulerPrime(ulong n) - { - // TODO: Inline this helper into callers so we skip the wrapper and pass the divisibility check directly to the optimized - // prime tester path highlighted in the CPU benchmarks. - return (n & 3UL) == 1UL && _primeTester.IsPrime(n, CancellationToken.None); - } -} - diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 2242764..6d451df 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -18,7 +18,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private readonly object _sync = new(); private ulong _divisorLimit; private bool _isConfigured; - private bool _useDivisorCycles; private readonly ConcurrentDictionary, ArrayView, ArrayView>> _kernelCache = new(); private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); @@ -43,12 +42,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv set => GpuBatchSize = value; } - public bool UseDivisorCycles - { - get => _useDivisorCycles; - set => _useDivisorCycles = value; // TODO: Remove this flag once GPU runs unconditionally fetch cycle lengths so every divisor check benefits from the measured speedups. - } - public void ConfigureFromMaxPrime(ulong maxPrime) { lock (_sync) @@ -61,7 +54,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv public bool IsPrime(ulong prime, out bool divisorsExhausted) { ulong allowedMax; - bool useCycles; int batchCapacity; lock (_sync) @@ -72,7 +64,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); - useCycles = _useDivisorCycles; batchCapacity = _gpuBatchSize; } @@ -97,7 +88,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv composite = CheckDivisors( prime, allowedMax, - useCycles, accelerator, kernel, resources.DivisorsBuffer, @@ -156,7 +146,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private static bool CheckDivisors( ulong prime, ulong allowedMax, - bool useCycles, Accelerator accelerator, Action, ArrayView, ArrayView> kernel, MemoryBuffer1D divisorsBuffer, @@ -170,68 +159,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv out bool coveredRange, out ulong processedCount) { - lastProcessed = 0UL; + int batchCapacity = (int)divisorsBuffer.Length; + bool composite = false; + bool processedAll = false; + ulong currentDivisor = 3UL; processedCount = 0UL; - - if (allowedMax < 3UL) - { - coveredRange = true; - return false; - } - - return useCycles - ? CheckDivisorsWithCycles( - prime, - allowedMax, - accelerator, - kernel, - divisorsBuffer, - exponentBuffer, - hitsBuffer, - divisors, - exponents, - hits, - divisorData, - ref lastProcessed, - ref processedCount, - out coveredRange) - : CheckDivisorsWithoutCycles( - prime, - allowedMax, - accelerator, - kernel, - divisorsBuffer, - exponentBuffer, - hitsBuffer, - divisors, - exponents, - hits, - divisorData, - ref lastProcessed, - ref processedCount, - out coveredRange); - } - - private static bool CheckDivisorsWithoutCycles( - ulong prime, - ulong allowedMax, - Accelerator accelerator, - Action, ArrayView, ArrayView> kernel, - MemoryBuffer1D divisorsBuffer, - MemoryBuffer1D exponentBuffer, - MemoryBuffer1D hitsBuffer, - in ulong[] divisors, - ulong[] exponents, - byte[] hits, - MontgomeryDivisorData[] divisorData, - ref ulong lastProcessed, - ref ulong processedCount, - out bool coveredRange) - { - int batchCapacity = (int)divisorsBuffer.Length; - bool composite = false; - bool processedAll = false; - ulong currentDivisor = 3UL; + lastProcessed = 0UL; Span divisorDataSpan; Span divisorSpan; @@ -241,138 +174,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayView1D exponentView; ArrayView1D hitsView; - while (currentDivisor <= allowedMax) + DivisorCycleCache cycleCache = DivisorCycleCache.Shared; + int cycleBatchCapacity = Math.Min(batchCapacity, cycleCache.PreferredBatchSize); + if (cycleBatchCapacity <= 0) { - int batchSize = 0; - bool reachedEndInBatch = false; - - // The GPU buffers are sized exactly to the configured batch capacity, so we - // must stop once the current batch is full to avoid writing past the spans. - while (batchSize < batchCapacity && currentDivisor <= allowedMax) - { - ulong nextDivisor = currentDivisor + 2UL; - processedCount++; - - divisors[batchSize] = currentDivisor; - exponents[batchSize] = prime; - batchSize++; - lastProcessed = currentDivisor; - - if (nextDivisor <= currentDivisor) - { - reachedEndInBatch = true; - break; - } - - if (nextDivisor > allowedMax) - { - currentDivisor = nextDivisor; - reachedEndInBatch = true; - break; - } - - currentDivisor = nextDivisor; - } - - if (batchSize == 0) - { - if (reachedEndInBatch) - { - processedAll = true; - break; - } - - continue; - } - - divisorDataSpan = divisorData.AsSpan(0, batchSize); - divisorSpan = divisors.AsSpan(0, batchSize); - exponentSpan = exponents.AsSpan(0, batchSize); - hitsSpan = hits.AsSpan(0, batchSize); - divisorView = divisorsBuffer.View.SubView(0, batchSize); - exponentView = exponentBuffer.View.SubView(0, batchSize); - hitsView = hitsBuffer.View.SubView(0, batchSize); - - for (int i = 0; i < batchSize; i++) - { - // TODO: Stage MontgomeryDivisorData in a GPU-ready cache keyed by divisor so we can copy - // precomputed ProcessEightBitWindows metadata directly instead of rehydrating the slower - // Montgomery structs for each batch item. - divisorDataSpan[i] = MontgomeryDivisorDataCache.Get(divisorSpan[i]); - } - - divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize); - exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSpan), batchSize); - kernel(batchSize, divisorView, exponentView, hitsView); - hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize); - - for (int i = 0; i < batchSize; i++) - { - if (hitsSpan[i] != 0) - { - composite = true; - lastProcessed = divisorSpan[i]; - break; - } - } - - if (!composite) - { - lastProcessed = divisorSpan[batchSize - 1]; - } - - if (composite) - { - break; - } - - if (reachedEndInBatch) - { - processedAll = true; - break; - } + cycleBatchCapacity = batchCapacity; } - coveredRange = composite || processedAll || currentDivisor > allowedMax; - return composite; - } - - private static bool CheckDivisorsWithCycles( - ulong prime, - ulong allowedMax, - Accelerator accelerator, - Action, ArrayView, ArrayView> kernel, - MemoryBuffer1D divisorsBuffer, - MemoryBuffer1D exponentBuffer, - MemoryBuffer1D hitsBuffer, - in ulong[] divisors, - ulong[] exponents, - byte[] hits, - MontgomeryDivisorData[] divisorData, - ref ulong lastProcessed, - ref ulong processedCount, - out bool coveredRange) - { - int batchCapacity = (int)divisorsBuffer.Length; - bool composite = false; - bool processedAll = false; - ulong currentDivisor = 3UL; - ulong cycleLeaseEnd; - - Span divisorDataSpan; - Span divisorSpan; - Span exponentSpan; - Span hitsSpan; - ArrayView1D divisorView; - ArrayView1D exponentView; - ArrayView1D hitsView; - - DivisorCycleCache divisorCyclesCache = DivisorCycleCache.Shared; - DivisorCycleCache.CycleBlock cycleLease = divisorCyclesCache.Acquire(currentDivisor); - // TODO: Replace the block lease with a direct single-cycle lookup so the GPU path - // follows the benchmarked single-block policy and computes the missing cycle inline on the - // configured device without mutating or swapping cache blocks after startup. - cycleLeaseEnd = cycleLease.End; + ulong[] cycleDivisors = ArrayPool.Shared.Rent(cycleBatchCapacity); + ulong[] cycleLengths = ArrayPool.Shared.Rent(cycleBatchCapacity); try { @@ -381,75 +191,77 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv int batchSize = 0; bool reachedEndInBatch = false; - // Batch writes share the same buffers as the kernel launch, so we guard - // against exceeding the rented capacity on each iteration. while (batchSize < batchCapacity && currentDivisor <= allowedMax) { - ulong nextDivisor = currentDivisor + 2UL; - processedCount++; - - if (currentDivisor > cycleLeaseEnd) + int chunkCount = Math.Min(cycleBatchCapacity, batchCapacity - batchSize); + ulong remainingWidth = allowedMax - currentDivisor; + ulong maxOddCount = (remainingWidth >> 1) + 1UL; + if (maxOddCount < (ulong)chunkCount) { - cycleLease.Dispose(); - cycleLease = divisorCyclesCache.Acquire(currentDivisor); - // TODO: Drop this block hand-off once the cache exposes a direct - // GetCycle helper that materializes only the requested divisor length - // without touching additional blocks or scheduling background work. - cycleLeaseEnd = cycleLease.End; + chunkCount = (int)maxOddCount; } - bool includeDivisor = true; - ulong exponent = prime; - ulong cycleLength = cycleLease.GetCycle(currentDivisor); - if (cycleLength > 0UL) - { - ulong remainder = prime; - if (remainder >= cycleLength) - { - // TODO: Replace this `%` with the cycle-cache helper that subtracts precomputed - // remainders (or uses a lookup table) so GPU batching avoids modulo in this inner loop. - remainder %= cycleLength; - } - - if (remainder != 0UL) - { - includeDivisor = false; - } - else - { - exponent = 0UL; - } - } - else - { - // TODO: Trigger an immediate GPU computation for the single missing cycle without storing - // the result back or requesting extra cache blocks so very large divisors stay within the - // single-block constraint while still reusing the snapshot for other divisors. - } - - if (includeDivisor) - { - divisors[batchSize] = currentDivisor; - exponents[batchSize] = exponent; - batchSize++; - } - - lastProcessed = currentDivisor; - - if (nextDivisor <= currentDivisor) + if (chunkCount <= 0) { reachedEndInBatch = true; break; } - if (nextDivisor > allowedMax) + Span chunkDivisors = cycleDivisors.AsSpan(0, chunkCount); + Span chunkCycles = cycleLengths.AsSpan(0, chunkCount); + + ulong fillDivisor = currentDivisor; + for (int i = 0; i < chunkCount; i++) { - currentDivisor = nextDivisor; - reachedEndInBatch = true; - break; + chunkDivisors[i] = fillDivisor; + fillDivisor += 2UL; } - currentDivisor = nextDivisor; + cycleCache.GetCycleLengths(chunkDivisors, chunkCycles); + + for (int i = 0; i < chunkCount && batchSize < batchCapacity; i++) + { + ulong divisorValue = chunkDivisors[i]; + ulong cycleLength = chunkCycles[i]; + + processedCount++; + 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) + { + divisors[batchSize] = divisorValue; + exponents[batchSize] = exponent; + batchSize++; + } + + lastProcessed = divisorValue; + currentDivisor = divisorValue + 2UL; + + if (currentDivisor > allowedMax) + { + reachedEndInBatch = true; + break; + } + } } if (batchSize == 0) @@ -473,8 +285,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv for (int i = 0; i < batchSize; i++) { - // TODO: Preload ProcessEightBitWindows-friendly divisor metadata for GPU batches so this - // hot loop stops rebuilding MontgomeryDivisorData structures for every scan chunk. divisorDataSpan[i] = MontgomeryDivisorDataCache.Get(divisorSpan[i]); } @@ -512,7 +322,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } finally { - cycleLease.Dispose(); + ArrayPool.Shared.Return(cycleDivisors, clearArray: false); + ArrayPool.Shared.Return(cycleLengths, clearArray: false); } coveredRange = composite || processedAll || currentDivisor > allowedMax; @@ -601,17 +412,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv throw new ObjectDisposedException(nameof(DivisorScanSession)); } - MersenneNumberDivisorByDivisorGpuTester owner = _owner; int primesLength = primes.Length; if (primesLength == 0) { return; } - // The production pipeline never adjusts the GPU batch size after a session - // is initialized, so the buffers allocated in the constructor remain - // large enough for every launch. - MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); var exponentStepper = new ExponentRemainderStepper(divisorData); if (!exponentStepper.IsValidModulus) @@ -620,19 +426,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return; } - bool useDivisorCycles = owner._useDivisorCycles; - if (!useDivisorCycles) - { - ProcessWithoutCycle(ref exponentStepper, divisorData, primes, hits); - return; - } - if (divisorCycle == 0UL) { - // TODO: Replace this guard with an on-demand cycle computation that runs on the device selected - // for the session, keeps the result transient instead of caching it, and continues operating with the - // single shared block loaded from disk. - throw new InvalidOperationException($"Missing divisor cycle for divisor {divisor}."); + divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); + if (divisorCycle == 0UL) + { + hits.Clear(); + return; + } } var cycleStepper = new CycleRemainderStepper(divisorCycle); @@ -640,78 +441,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ProcessWithCycle(ref exponentStepper, ref cycleStepper, initialCycleRemainder, divisorCycle, divisorData, primes, hits); } - private void ProcessWithoutCycle( - ref ExponentRemainderStepper exponentStepper, - in MontgomeryDivisorData divisorData, - ReadOnlySpan primes, - Span hits) - { - Action, ArrayView> kernel = _kernel; - ArrayView1D exponentsView = _exponentsBuffer.View; - ArrayView1D resultsView = _resultsBuffer.View; - Span hostSpan = _hostBuffer.AsSpan(0, _capacity); - - int primesLength = primes.Length; - int gpuBatchSize = _owner._gpuBatchSize; - int offset = 0; - - while (offset < primesLength) - { - int batchSize = Math.Min(gpuBatchSize, primesLength - offset); - ReadOnlySpan primesSlice = primes.Slice(offset, batchSize); - Span hitsSlice = hits.Slice(offset, batchSize); - Span exponentSlice = hostSpan[..batchSize]; - - bool hasState = exponentStepper.HasState; - ulong previousExponent = exponentStepper.PreviousExponent; - - int startIndex = 0; - if (!hasState) - { - ulong primeValue = primesSlice[0]; - exponentSlice[0] = primeValue; - previousExponent = primeValue; - hasState = true; - startIndex = 1; - } - - for (int i = startIndex; i < batchSize; i++) - { - ulong primeValue = primesSlice[i]; - exponentSlice[i] = primeValue - previousExponent; - previousExponent = primeValue; - } - - ArrayView1D exponentView = exponentsView.SubView(0, batchSize); - exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), batchSize); - - ArrayView1D resultView = resultsView.SubView(0, batchSize); - kernel(batchSize, divisorData, exponentView, resultView); - resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), batchSize); - - bool requiresInitialization = !exponentStepper.HasState; - for (int i = 0; i < batchSize; i++) - { - ulong primeValue = primesSlice[i]; - ulong montgomeryResult = exponentSlice[i]; - bool isUnity; - - if (requiresInitialization) - { - exponentStepper.TryInitializeFromMontgomeryResult(primeValue, montgomeryResult, out isUnity); - requiresInitialization = false; - } - else - { - exponentStepper.TryAdvanceWithMontgomeryDelta(primeValue, montgomeryResult, out isUnity); - } - - hitsSlice[i] = isUnity ? (byte)1 : (byte)0; - } - - offset += batchSize; - } - } private void ProcessWithCycle( ref ExponentRemainderStepper exponentStepper, @@ -973,3 +702,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() => CreateDivisorSession(); } + + + diff --git a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs index 03cd04b..cfdbbd6 100644 --- a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs @@ -4,8 +4,6 @@ namespace PerfectNumbers.Core; public interface IMersenneNumberDivisorByDivisorTester { - bool UseDivisorCycles { get; set; } // TODO: Remove the setter once divisor cycle usage becomes mandatory so all implementations always leverage the cached cycles. - int BatchSize { get; set; } void ConfigureFromMaxPrime(ulong maxPrime); @@ -25,4 +23,3 @@ public interface IMersenneNumberDivisorByDivisorTester void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits); } } - diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index bbbcaab..1a2291f 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -237,7 +237,6 @@ public static class MersenneNumberDivisorByDivisorTester Span hitsSpan = default; int hitIndex = 0; int index; - bool useDivisorCycles = tester.UseDivisorCycles; // TODO: Remove the conditional path once all testers always enable divisor cycles to keep every divisor scan on the faster cached-length track. ulong divisorCycle = 0UL; try @@ -255,23 +254,7 @@ public static class MersenneNumberDivisorByDivisorTester divisor = AcquireNextDivisor(ref nextDivisor, divisorLimit, ref divisorsExhaustedFlag, ref finalDivisorBits, out exhausted, ref localDivisorCursor, ref localDivisorsRemaining); if (divisor != 0UL) { - DivisorCycleCache.CycleBlock? cycleBlock = null; - if (useDivisorCycles) - { - cycleBlock = DivisorCycleCache.Shared.Acquire(divisor); - // TODO: Remove this block leasing once the cache exposes a direct - // single-cycle retrieval API so coordinator threads never mutate the - // cache state beyond the snapshot loaded at startup. - divisorCycle = cycleBlock.GetCycle(divisor); - // TODO: When divisorCycle returns zero, compute only that single cycle on the - // device selected by the current configuration, keep the result transient, and - // continue operating with the single cached block instead of scheduling additional - // batches. - } - else - { - divisorCycle = 0UL; - } + divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); activeCount = BuildPrimeBuffer( divisor, @@ -296,7 +279,6 @@ public static class MersenneNumberDivisorByDivisorTester if (activeCount == 0) { - cycleBlock?.Dispose(); continue; } @@ -330,11 +312,8 @@ public static class MersenneNumberDivisorByDivisorTester FlushPendingResults(compositesBuffer, ref compositesCount, markComposite, clearComposite, printResult); } - cycleBlock?.Dispose(); - continue; } - if (!exhausted) { if (Volatile.Read(ref remainingStates) == 0) @@ -738,3 +717,6 @@ public static class MersenneNumberDivisorByDivisorTester internal bool PassedAllTests { get; } } } + + + diff --git a/PerfectNumbers.Core/MontgomeryDivisorData.cs b/PerfectNumbers.Core/MontgomeryDivisorData.cs index cde8ffe..fa73db7 100644 --- a/PerfectNumbers.Core/MontgomeryDivisorData.cs +++ b/PerfectNumbers.Core/MontgomeryDivisorData.cs @@ -17,13 +17,7 @@ public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong internal static class MontgomeryDivisorDataCache { - private const ulong BaseBlockEnd = PerfectNumberConstants.MaxQForDivisorCycles; - - private static readonly ConcurrentDictionary Cache = new(); // TODO: Replace this with a plain Dictionary that never mutates after startup once the divisor-cycle snapshot becomes read-only so Montgomery lookups avoid any concurrent collections. - // TODO: Delete the block-membership tracking entirely when the divisor-cycle cache stops producing - // dynamic blocks. The single-cycle plan keeps the snapshot immutable, so we should not retain block - // indices or concurrent maps at all. - private static readonly ConcurrentDictionary> BlockMembership = new(); + private static readonly ConcurrentDictionary Cache = new(); public static MontgomeryDivisorData Get(ulong modulus) { @@ -32,97 +26,7 @@ internal static class MontgomeryDivisorDataCache return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL); } - int blockIndex = ComputeBlockIndex(modulus); - // TODO: Short-circuit block indices greater than zero once the cache stops prefetching future - // blocks; when only the startup snapshot remains we should bypass block arithmetic entirely and - // compute single cycles inline on cache misses without updating shared state. - return GetManaged(modulus, blockIndex); - } - - internal static void ReleaseBlock(DivisorCycleCache.CycleBlock block) - { - // TODO: Remove this entire release hook once the cycle cache no longer hands out dynamic blocks; - // the single-cycle pipeline will never transfer ownership back to Montgomery metadata. - if (block.Index <= 0) - { - return; - } - - if (BlockMembership.TryRemove(block.Index, out ConcurrentDictionary? members)) - { - foreach (ulong modulus in members.Keys) - { - Cache.TryRemove(modulus, out _); - } - } - } - - private static MontgomeryDivisorData GetManaged(ulong modulus, int blockIndex) - { - while (true) - { - if (Cache.TryGetValue(modulus, out CacheEntry entry)) - { - if (entry.BlockIndex != blockIndex) - { - if (blockIndex > 0 && entry.BlockIndex <= 0) - { - CacheEntry updated = new(entry.Data, blockIndex); - if (Cache.TryUpdate(modulus, updated, entry)) - { - RegisterBlockMembership(modulus, blockIndex); - return entry.Data; - } - - continue; - } - - if (blockIndex > 0) - { - RegisterBlockMembership(modulus, blockIndex); - } - - return entry.Data; - } - return entry.Data; - } - - MontgomeryDivisorData created = Create(modulus); - if (Cache.TryAdd(modulus, new CacheEntry(created, blockIndex))) - { - RegisterBlockMembership(modulus, blockIndex); - return created; - } - } - // TODO: Collapse this retry loop once the cache becomes a simple snapshot dictionary. With the - // single-cycle policy we can compute the requested cycle, populate the dictionary exactly once, - // and return without spinning through concurrent retries. - } - - private static void RegisterBlockMembership(ulong modulus, int blockIndex) - { - // TODO: Delete this method once block membership goes away; the single snapshot model must avoid - // tracking additional indices or allocating concurrent dictionaries for Montgomery lookups. - if (blockIndex <= 0) - { - return; - } - - ConcurrentDictionary set = BlockMembership.GetOrAdd(blockIndex, static _ => new ConcurrentDictionary()); - set.TryAdd(modulus, 0); - } - - private static int ComputeBlockIndex(ulong modulus) - { - if (modulus <= BaseBlockEnd) - { - return 0; - } - - // TODO: Return 0 here once the divisor-cycle cache no longer tracks additional blocks so - // MontgomeryDivisorData lookups never attempt to follow prefetched ranges or trigger new cycle - // computations for entire blocks. - return (int)((modulus - (BaseBlockEnd + 1UL)) / (ulong)PerfectNumberConstants.MaxQForDivisorCycles) + 1; + return Cache.GetOrAdd(modulus, static m => Create(m)); } private static MontgomeryDivisorData Create(ulong modulus) @@ -135,8 +39,6 @@ internal static class MontgomeryDivisorDataCache } [MethodImpl(MethodImplOptions.AggressiveInlining)] - // TODO: Replace this generic `%` reduction with the UInt128 Montgomery folding helper measured faster in - // MontgomeryMultiplyBenchmarks so CPU paths avoid the slow BigInteger-style remainder in hot divisor cache lookups. private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus); private static ulong ComputeMontgomeryNPrime(ulong modulus) @@ -148,19 +50,7 @@ internal static class MontgomeryDivisorDataCache inv *= unchecked(2UL - modulus * inv); inv *= unchecked(2UL - modulus * inv); inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); return unchecked(0UL - inv); } - - private readonly struct CacheEntry - { - internal CacheEntry(MontgomeryDivisorData data, int blockIndex) - { - Data = data; - BlockIndex = blockIndex; - } - - internal MontgomeryDivisorData Data { get; } - - internal int BlockIndex { get; } - } }