diff --git a/EvenPerfectBitScanner/CliArguments.cs b/EvenPerfectBitScanner/CliArguments.cs index af98eeb..5629ad2 100644 --- a/EvenPerfectBitScanner/CliArguments.cs +++ b/EvenPerfectBitScanner/CliArguments.cs @@ -5,19 +5,20 @@ using PerfectNumbers.Core.Gpu; namespace EvenPerfectBitScanner; - internal enum ByDivisorKIncrementMode - { - Sequential, - Pow2Groups, - Predictive, - Percentile, - Additive, - TopDown, - BitContradiction, +internal enum ByDivisorKIncrementMode +{ + Sequential, + Pow2Groups, + Predictive, + Percentile, + Additive, + TopDown, + BitContradiction, + BitTree, } - internal readonly struct CliArguments - { +internal readonly struct CliArguments +{ private const string DefaultCyclesPath = "divisor_cycles.bin"; internal readonly ulong StartPrime; @@ -35,19 +36,19 @@ namespace EvenPerfectBitScanner; internal readonly bool UseLucas; internal readonly bool UseResidue; internal readonly bool UseDivisor; - internal readonly bool UseByDivisor; - internal readonly ByDivisorKIncrementMode ByDivisorKIncrementMode; - internal readonly int ByDivisorSpecialRange; - internal readonly bool UseGcdFilter; - internal readonly bool TestMode; - internal readonly bool UseGpuCycles; - internal readonly ComputationDevice MersenneDevice; - internal readonly ComputationDevice OrderDevice; - internal readonly int ScanBatchSize; - internal readonly int SliceSize; - internal readonly ulong MaxKLimit; + internal readonly bool UseByDivisor; + internal readonly ByDivisorKIncrementMode ByDivisorKIncrementMode; + internal readonly int ByDivisorSpecialRange; + internal readonly bool UseGcdFilter; + internal readonly bool TestMode; + internal readonly bool UseGpuCycles; + internal readonly ComputationDevice MersenneDevice; + internal readonly ComputationDevice OrderDevice; + internal readonly int ScanBatchSize; + internal readonly int SliceSize; + internal readonly ulong MaxKLimit; internal readonly bool MaxKConfigured; - internal readonly string FilterFile; + internal readonly string FilterFile; internal readonly bool ForcePrimeKernelsOnCpu; internal readonly ModReductionMode ModReductionMode; internal readonly NttBackend NttBackend; @@ -66,56 +67,56 @@ namespace EvenPerfectBitScanner; internal readonly double ZeroFractionConjecture; internal readonly int MaxZeroConjecture; internal readonly BigInteger MinK; - internal readonly int GpuRatio; + internal readonly int GpuRatio; - private CliArguments( - ulong startPrime, - bool startPrimeProvided, - UInt128 maxPrimeLimit, - bool maxPrimeConfigured, - int threadCount, - int blockSize, - int gpuPrimeThreads, - int gpuPrimeBatch, - GpuKernelType kernelType, - bool useOrder, - bool showHelp, - bool useBitTransform, - bool useLucas, - bool useResidue, - bool useDivisor, - bool useByDivisor, - bool useGcdFilter, - ByDivisorKIncrementMode byDivisorKIncrementMode, - int byDivisorSpecialRange, - bool testMode, - bool useGpuCycles, - ComputationDevice mersenneDevice, - ComputationDevice orderDevice, - int scanBatchSize, - int sliceSize, - ulong maxKLimit, - bool maxKConfigured, - string filterFile, - bool forcePrimeKernelsOnCpu, - ModReductionMode modReductionMode, - NttBackend nttBackend, - string? resultsDirectory, - string? resultsPrefix, - int writeBatchSize, - string cyclesPath, - bool continueCyclesGeneration, - int cyclesBatchSize, - ulong divisorCyclesSearchLimit, - ulong? orderWarmupLimitOverride, - string? rleBlacklistPath, - ulong rleHardMaxP, - bool rleOnlyLast7, - double zeroFractionHard, - double zeroFractionConjecture, - int maxZeroConjecture, - BigInteger minK, - int gpuRatio) + private CliArguments( + ulong startPrime, + bool startPrimeProvided, + UInt128 maxPrimeLimit, + bool maxPrimeConfigured, + int threadCount, + int blockSize, + int gpuPrimeThreads, + int gpuPrimeBatch, + GpuKernelType kernelType, + bool useOrder, + bool showHelp, + bool useBitTransform, + bool useLucas, + bool useResidue, + bool useDivisor, + bool useByDivisor, + bool useGcdFilter, + ByDivisorKIncrementMode byDivisorKIncrementMode, + int byDivisorSpecialRange, + bool testMode, + bool useGpuCycles, + ComputationDevice mersenneDevice, + ComputationDevice orderDevice, + int scanBatchSize, + int sliceSize, + ulong maxKLimit, + bool maxKConfigured, + string filterFile, + bool forcePrimeKernelsOnCpu, + ModReductionMode modReductionMode, + NttBackend nttBackend, + string? resultsDirectory, + string? resultsPrefix, + int writeBatchSize, + string cyclesPath, + bool continueCyclesGeneration, + int cyclesBatchSize, + ulong divisorCyclesSearchLimit, + ulong? orderWarmupLimitOverride, + string? rleBlacklistPath, + ulong rleHardMaxP, + bool rleOnlyLast7, + double zeroFractionHard, + double zeroFractionConjecture, + int maxZeroConjecture, + BigInteger minK, + int gpuRatio) { StartPrime = startPrime; StartPrimeProvided = startPrimeProvided; @@ -163,7 +164,7 @@ namespace EvenPerfectBitScanner; ZeroFractionConjecture = zeroFractionConjecture; MaxZeroConjecture = maxZeroConjecture; MinK = minK; - GpuRatio = gpuRatio; + GpuRatio = gpuRatio; } internal bool UseFilter => !TestMode && !string.IsNullOrEmpty(FilterFile); @@ -204,7 +205,7 @@ namespace EvenPerfectBitScanner; string? resultsDirectory = null; string? resultsPrefix = null; int writeBatchSize = CalculationResultHandler.DefaultWriteBatchSize; - string cyclesPath = DefaultCyclesPath; + string cyclesPath = DefaultCyclesPath; bool continueCyclesGeneration = false; int cyclesBatchSize = 512; ulong divisorCyclesSearchLimit = PerfectNumberConstants.ExtraDivisorCycleSearchLimit; @@ -216,7 +217,7 @@ namespace EvenPerfectBitScanner; double zeroFractionConjecture = -1.0; int maxZeroConjecture = -1; BigInteger minK = BigInteger.One; - int gpuRatio = PerfectNumberConstants.GpuRatio; + int gpuRatio = PerfectNumberConstants.GpuRatio; foreach (string argument in args) { @@ -279,6 +280,19 @@ namespace EvenPerfectBitScanner; continue; } + if (argument.Equals("--use-order", StringComparison.OrdinalIgnoreCase)) + { + useOrder = true; + continue; + } + + if (argument.StartsWith("--max-k=", StringComparison.OrdinalIgnoreCase)) + { + maxKLimit = Utf8CliParser.ParseUInt64(argument.AsSpan("--max-k=".Length)); + maxKConfigured = true; + continue; + } + if (argument.StartsWith("--mersenne=", StringComparison.OrdinalIgnoreCase)) { ReadOnlySpan value = argument.AsSpan("--mersenne=".Length); @@ -318,138 +332,46 @@ namespace EvenPerfectBitScanner; if (argument.StartsWith("--bydivisor-k-increment=", StringComparison.OrdinalIgnoreCase)) { - string value = argument["--bydivisor-k-increment=".Length..].Trim(); + ReadOnlySpan value = argument.AsSpan("--bydivisor-k-increment=".Length); if (value.Equals("pow2groups", StringComparison.OrdinalIgnoreCase)) { byDivisorKIncrementMode = ByDivisorKIncrementMode.Pow2Groups; } - else if (value.Equals("sequential", StringComparison.OrdinalIgnoreCase)) - { - byDivisorKIncrementMode = ByDivisorKIncrementMode.Sequential; - } else if (value.Equals("predictive", StringComparison.OrdinalIgnoreCase)) { byDivisorKIncrementMode = ByDivisorKIncrementMode.Predictive; } - else if (value.Equals("percentile", StringComparison.OrdinalIgnoreCase)) - { - byDivisorKIncrementMode = ByDivisorKIncrementMode.Percentile; - } - else if (value.Equals("additive", StringComparison.OrdinalIgnoreCase)) - { - byDivisorKIncrementMode = ByDivisorKIncrementMode.Additive; - } - else if (value.Equals("topdown", StringComparison.OrdinalIgnoreCase)) - { - byDivisorKIncrementMode = ByDivisorKIncrementMode.TopDown; - } - else if (value.Equals("bitcontradiction", StringComparison.OrdinalIgnoreCase)) - { - byDivisorKIncrementMode = ByDivisorKIncrementMode.BitContradiction; - } - else - { - Console.WriteLine("Invalid value for --bydivisor-k-increment. Use sequential, pow2groups, predictive, percentile, additive, topdown, or bitcontradiction."); - return default; - } + else if (value.Equals("percentile", StringComparison.OrdinalIgnoreCase)) + { + byDivisorKIncrementMode = ByDivisorKIncrementMode.Percentile; + } + else if (value.Equals("additive", StringComparison.OrdinalIgnoreCase)) + { + byDivisorKIncrementMode = ByDivisorKIncrementMode.Additive; + } + else if (value.Equals("topdown", StringComparison.OrdinalIgnoreCase)) + { + byDivisorKIncrementMode = ByDivisorKIncrementMode.TopDown; + } + else if (value.Equals("bitcontradiction", StringComparison.OrdinalIgnoreCase)) + { + byDivisorKIncrementMode = ByDivisorKIncrementMode.BitContradiction; + } + else if (value.Equals("bittree", StringComparison.OrdinalIgnoreCase)) + { + byDivisorKIncrementMode = ByDivisorKIncrementMode.BitTree; + } + else + { + byDivisorKIncrementMode = ByDivisorKIncrementMode.Sequential; + } continue; } if (argument.StartsWith("--bydivisor-special-range=", StringComparison.OrdinalIgnoreCase)) { - if (Utf8CliParser.TryParseInt32(argument.AsSpan("--bydivisor-special-range=".Length), out int parsedRange) && parsedRange >= 0) - { - byDivisorSpecialRange = parsedRange; - } - else - { - Console.WriteLine("Invalid value for --bydivisor-special-range. Use a non-negative integer."); - return default; - } - - continue; - } - - if (argument.StartsWith("--divisor-cycles-limit=", StringComparison.OrdinalIgnoreCase)) - { - if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--divisor-cycles-limit=".Length), out ulong parsedLimit)) - { - divisorCyclesSearchLimit = parsedLimit; - } - - continue; - } - - if (argument.StartsWith("--max-k=", StringComparison.OrdinalIgnoreCase)) - { - if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--max-k=".Length), out ulong parsedMaxK)) - { - maxKLimit = parsedMaxK; - maxKConfigured = true; - } - - continue; - } - - if (argument.StartsWith("--min-k=", StringComparison.OrdinalIgnoreCase)) - { - ReadOnlySpan minKSpan = argument.AsSpan("--min-k=".Length); - if (Utf8CliParser.TryParseUInt64(minKSpan, out ulong parsedMinK) && parsedMinK > 0UL) - { - minK = parsedMinK; - } - else - { - try - { - BigInteger parsedMinKBig = BigInteger.Parse(minKSpan, CultureInfo.InvariantCulture); - if (parsedMinKBig > BigInteger.Zero) - { - minK = parsedMinKBig; - } - else - { - Console.WriteLine("Invalid value for --min-k."); - return default; - } - } - catch - { - Console.WriteLine("Invalid value for --min-k."); - return default; - } - } - - continue; - } - - if (argument.StartsWith("--mersenne-device=", StringComparison.OrdinalIgnoreCase)) - { - ReadOnlySpan value = argument.AsSpan("--mersenne-device=".Length); - mersenneDevice = ParseDevice(value); - continue; - } - - if (argument.Equals("--test", StringComparison.OrdinalIgnoreCase)) - { - testMode = true; - continue; - } - - if (argument.Equals("--use-order", StringComparison.OrdinalIgnoreCase)) - { - useOrder = true; - continue; - } - - if (argument.StartsWith("--order-warmup-limit=", StringComparison.OrdinalIgnoreCase)) - { - if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--order-warmup-limit=".Length), out ulong parsedLimit)) - { - orderWarmupLimitOverride = parsedLimit; - } - + byDivisorSpecialRange = Utf8CliParser.ParseInt32(argument.AsSpan("--bydivisor-special-range=".Length)); continue; } @@ -459,138 +381,30 @@ namespace EvenPerfectBitScanner; continue; } - if (argument.StartsWith("--ntt=", StringComparison.OrdinalIgnoreCase)) + if (argument.Equals("--test-mode", StringComparison.OrdinalIgnoreCase)) { - ReadOnlySpan value = argument.AsSpan("--ntt=".Length); - nttBackend = value.Equals("staged", StringComparison.OrdinalIgnoreCase) ? NttBackend.Staged : NttBackend.Reference; + testMode = true; continue; } - if (argument.StartsWith("--mod-reduction=", StringComparison.OrdinalIgnoreCase)) + if (argument.StartsWith("--mersenne-device=", StringComparison.OrdinalIgnoreCase)) { - ReadOnlySpan value = argument.AsSpan("--mod-reduction=".Length); - if (value.Equals("mont64", StringComparison.OrdinalIgnoreCase)) - { - modReductionMode = ModReductionMode.Mont64; - } - else if (value.Equals("barrett128", StringComparison.OrdinalIgnoreCase)) - { - modReductionMode = ModReductionMode.Barrett128; - } - else if (value.Equals("uint128", StringComparison.OrdinalIgnoreCase)) - { - modReductionMode = ModReductionMode.GpuUInt128; - } - else - { - modReductionMode = ModReductionMode.Auto; - } - + mersenneDevice = ParseDevice(argument.AsSpan("--mersenne-device=".Length)); continue; } - if (argument.StartsWith("--primes-device=", StringComparison.OrdinalIgnoreCase)) + if (argument.StartsWith("--primes=", StringComparison.OrdinalIgnoreCase) || + argument.StartsWith("--gpu-kernels=", StringComparison.OrdinalIgnoreCase) || + argument.StartsWith("--accelerator=", StringComparison.OrdinalIgnoreCase)) { - forcePrimeKernelsOnCpu = argument.AsSpan("--primes-device=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); + ReadOnlySpan value = argument.Contains('=') ? argument.AsSpan(argument.IndexOf('=') + 1) : argument.AsSpan(); + forcePrimeKernelsOnCpu = value.Equals("cpu", StringComparison.OrdinalIgnoreCase); continue; } if (argument.StartsWith("--order-device=", StringComparison.OrdinalIgnoreCase)) { - ReadOnlySpan value = argument.AsSpan("--order-device=".Length); - orderDevice = ParseDevice(value); - continue; - } - - if (argument.StartsWith("--rle-blacklist=", StringComparison.OrdinalIgnoreCase)) - { - rleBlacklistPath = argument["--rle-blacklist=".Length..]; - continue; - } - - if (argument.StartsWith("--rle-hard-max=", StringComparison.OrdinalIgnoreCase)) - { - if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--rle-hard-max=".Length), out ulong parsedMax)) - { - rleHardMaxP = parsedMax; - } - - continue; - } - - if (argument.StartsWith("--rle-only-last7=", StringComparison.OrdinalIgnoreCase)) - { - ReadOnlySpan value = argument.AsSpan("--rle-only-last7=".Length); - rleOnlyLast7 = !value.Equals("false", StringComparison.OrdinalIgnoreCase) && !value.Equals("0", StringComparison.OrdinalIgnoreCase); - continue; - } - - if (argument.StartsWith("--zero-hard=", StringComparison.OrdinalIgnoreCase)) - { - if (Utf8CliParser.TryParseDouble(argument.AsSpan("--zero-hard=".Length), out double zeroFraction)) - { - zeroFractionHard = zeroFraction; - } - - continue; - } - - if (argument.StartsWith("--zero-conj=", StringComparison.OrdinalIgnoreCase)) - { - ReadOnlySpan value = argument.AsSpan("--zero-conj=".Length); - int colon = value.IndexOf(':'); - if (colon > 0) - { - if (Utf8CliParser.TryParseDouble(value[..colon], out double zeroFrac) && - Utf8CliParser.TryParseInt32(value[(colon + 1)..], out int maxZero)) - { - zeroFractionConjecture = zeroFrac; - maxZeroConjecture = maxZero; - } - } - - continue; - } - - if (argument.StartsWith("--lucas=", StringComparison.OrdinalIgnoreCase)) - { - mersenneDevice = ParseDevice(argument.AsSpan("--lucas=".Length)); - continue; - } - - if (argument.StartsWith("--primes=", StringComparison.OrdinalIgnoreCase)) - { - forcePrimeKernelsOnCpu = argument.AsSpan("--primes=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); - continue; - } - - if (argument.StartsWith("--gpu-kernels=", StringComparison.OrdinalIgnoreCase)) - { - forcePrimeKernelsOnCpu = argument.AsSpan("--gpu-kernels=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); - continue; - } - - if (argument.StartsWith("--accelerator=", StringComparison.OrdinalIgnoreCase)) - { - forcePrimeKernelsOnCpu = argument.AsSpan("--accelerator=".Length).Equals("cpu", StringComparison.OrdinalIgnoreCase); - continue; - } - - if (argument.StartsWith("--results-dir=", StringComparison.OrdinalIgnoreCase)) - { - resultsDirectory = argument["--results-dir=".Length..]; - continue; - } - - if (argument.StartsWith("--results-prefix=", StringComparison.OrdinalIgnoreCase)) - { - resultsPrefix = argument["--results-prefix=".Length..]; - continue; - } - - if (argument.StartsWith("--gpu-prime-threads=", StringComparison.OrdinalIgnoreCase)) - { - gpuPrimeThreads = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--gpu-prime-threads=".Length))); + orderDevice = ParseDevice(argument.AsSpan("--order-device=".Length)); continue; } @@ -600,6 +414,12 @@ namespace EvenPerfectBitScanner; continue; } + if (argument.StartsWith("--gpu-prime-threads=", StringComparison.OrdinalIgnoreCase)) + { + gpuPrimeThreads = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--gpu-prime-threads=".Length))); + continue; + } + if (argument.StartsWith("--ll-slice=", StringComparison.OrdinalIgnoreCase)) { sliceSize = Math.Max(1, Utf8CliParser.ParseInt32(argument.AsSpan("--ll-slice=".Length))); @@ -618,6 +438,54 @@ namespace EvenPerfectBitScanner; continue; } + if (argument.StartsWith("--mod-reduction=", StringComparison.OrdinalIgnoreCase)) + { + ReadOnlySpan value = argument.AsSpan("--mod-reduction=".Length); + if (value.Equals("uint128", StringComparison.OrdinalIgnoreCase)) + { + modReductionMode = ModReductionMode.GpuUInt128; + } + else if (value.Equals("mont64", StringComparison.OrdinalIgnoreCase)) + { + modReductionMode = ModReductionMode.Mont64; + } + else if (value.Equals("barrett128", StringComparison.OrdinalIgnoreCase)) + { + modReductionMode = ModReductionMode.Barrett128; + } + else + { + modReductionMode = ModReductionMode.Auto; + } + + continue; + } + + if (argument.StartsWith("--ntt=", StringComparison.OrdinalIgnoreCase)) + { + ReadOnlySpan value = argument.AsSpan("--ntt=".Length); + nttBackend = value.Equals("staged", StringComparison.OrdinalIgnoreCase) ? NttBackend.Staged : NttBackend.Reference; + continue; + } + + if (argument.StartsWith("--order-warmup-limit=", StringComparison.OrdinalIgnoreCase)) + { + orderWarmupLimitOverride = Utf8CliParser.ParseUInt64(argument.AsSpan("--order-warmup-limit=".Length)); + continue; + } + + if (argument.StartsWith("--divisor-cycles-limit=", StringComparison.OrdinalIgnoreCase)) + { + divisorCyclesSearchLimit = Utf8CliParser.ParseUInt64(argument.AsSpan("--divisor-cycles-limit=".Length)); + continue; + } + + if (argument.StartsWith("--min-k=", StringComparison.OrdinalIgnoreCase)) + { + minK = BigInteger.Parse(argument.AsSpan("--min-k=".Length), NumberStyles.None, CultureInfo.InvariantCulture); + continue; + } + if (argument.StartsWith("--filter-p=", StringComparison.OrdinalIgnoreCase)) { filterFile = argument["--filter-p=".Length..]; @@ -651,6 +519,7 @@ namespace EvenPerfectBitScanner; if (argument.Equals("--divisor-cycles-continue", StringComparison.OrdinalIgnoreCase)) { continueCyclesGeneration = true; + continue; } } @@ -659,15 +528,15 @@ namespace EvenPerfectBitScanner; startPrimeProvided, maxPrimeLimit, maxPrimeConfigured, - threadCount, - blockSize, - gpuPrimeThreads, - gpuPrimeBatch, - kernelType, - useOrder, - showHelp, - useBitTransform, - useLucas, + threadCount, + blockSize, + gpuPrimeThreads, + gpuPrimeBatch, + kernelType, + useOrder, + showHelp, + useBitTransform, + useLucas, useResidue, useDivisor, useByDivisor, @@ -701,7 +570,7 @@ namespace EvenPerfectBitScanner; zeroFractionConjecture, maxZeroConjecture, minK, - gpuRatio); + gpuRatio); } internal static void PrintHelp() @@ -715,14 +584,14 @@ namespace EvenPerfectBitScanner; Console.WriteLine(" --increment=bit|add exponent increment method (default add)"); Console.WriteLine(" --threads= number of worker threads"); Console.WriteLine(" --mersenne=pow2mod|incremental|lucas|residue|divisor|bydivisor Mersenne test method"); - Console.WriteLine(" --bydivisor-k-increment=sequential|pow2groups|predictive|percentile|additive|topdown|bitcontradiction k iteration strategy for --mersenne=bydivisor (default sequential)"); + Console.WriteLine(" --bydivisor-k-increment=sequential|pow2groups|predictive|percentile|additive|topdown|bitcontradiction|bittree k iteration strategy for --mersenne=bydivisor (default sequential)"); Console.WriteLine(" --bydivisor-special-range= widen pow2groups special values to +/- n percent (default 0 = points)"); Console.WriteLine(" --max-k= max k for Mersenne tests (q = 2*p*k + 1), including --mersenne=bydivisor"); Console.WriteLine(" --mersenne-device=cpu|gpu|hybrid 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($" --gpu-ratio= ratio of GPU operations vs CPU in hybrid device mode (default {PerfectNumberConstants.GpuRatio})"); + Console.WriteLine($" --gpu-ratio= ratio of GPU operations vs CPU in hybrid device mode (default {PerfectNumberConstants.GpuRatio})"); Console.WriteLine(" --gpu-prime-threads= max concurrent GPU prime checks (default 1)"); + Console.WriteLine(" --gpu-prime-batch= batch size for GPU primality sieve (default 262144)"); Console.WriteLine(" --gpu-scan-batch= GPU q-scan batch size (default 2_097_152)"); Console.WriteLine(" --order-device=cpu|gpu|hybrid Device for order computations (default gpu)"); Console.WriteLine(" --ll-slice= Lucas–Lehmer iterations per slice (default 32)"); diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index a28a1fe..7d7c741 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -32,6 +32,9 @@ using ByDivisorGpuOrderTopDown = PerfectNumbers.Core.Cpu.MersenneNumberDivisorBy using ByDivisorCpuOrderBitContradiction = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitContradictionDivisorSetForCpuOrder; using ByDivisorHybridOrderBitContradiction = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitContradictionDivisorSetForHybridOrder; using ByDivisorGpuOrderBitContradiction = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitContradictionDivisorSetForGpuOrder; +using ByDivisorCpuOrderBitTree = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitTreeDivisorSetForCpuOrder; +using ByDivisorHybridOrderBitTree = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitTreeDivisorSetForHybridOrder; +using ByDivisorGpuOrderBitTree = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitTreeDivisorSetForGpuOrder; using System.Globalization; namespace EvenPerfectBitScanner; @@ -64,6 +67,9 @@ private static ThreadLocal? MersenneTesters; private static ByDivisorCpuOrderBitContradiction _byDivisorCpuOrderTesterBitContradiction; private static ByDivisorHybridOrderBitContradiction _byDivisorHybridOrderTesterBitContradiction; private static ByDivisorGpuOrderBitContradiction _byDivisorGpuOrderTesterBitContradiction; + private static ByDivisorCpuOrderBitTree _byDivisorCpuOrderTesterBitTree; + private static ByDivisorHybridOrderBitTree _byDivisorHybridOrderTesterBitTree; + private static ByDivisorGpuOrderBitTree _byDivisorGpuOrderTesterBitTree; private static Dictionary? _byDivisorPreviousResults; private static string? _byDivisorModelPath; private static bool _byDivisorModeUsesModel; @@ -335,6 +341,12 @@ private static ThreadLocal? MersenneTesters; case ByDivisorKIncrementMode.TopDown: _byDivisorCpuOrderTesterTopDown = new ByDivisorCpuOrderTopDown(); break; + case ByDivisorKIncrementMode.BitContradiction: + _byDivisorCpuOrderTesterBitContradiction = new ByDivisorCpuOrderBitContradiction(); + break; + case ByDivisorKIncrementMode.BitTree: + _byDivisorCpuOrderTesterBitTree = new ByDivisorCpuOrderBitTree(); + break; default: _byDivisorCpuOrderTester = new MersenneNumberDivisorByDivisorCpuTesterWithCpuOrder(); break; @@ -359,6 +371,12 @@ private static ThreadLocal? MersenneTesters; case ByDivisorKIncrementMode.TopDown: _byDivisorHybridOrderTesterTopDown = new ByDivisorHybridOrderTopDown(); break; + case ByDivisorKIncrementMode.BitContradiction: + _byDivisorHybridOrderTesterBitContradiction = new ByDivisorHybridOrderBitContradiction(); + break; + case ByDivisorKIncrementMode.BitTree: + _byDivisorHybridOrderTesterBitTree = new ByDivisorHybridOrderBitTree(); + break; default: _byDivisorHybridOrderTester = new MersenneNumberDivisorByDivisorCpuTesterWithHybridOrder(); break; @@ -383,6 +401,12 @@ private static ThreadLocal? MersenneTesters; case ByDivisorKIncrementMode.TopDown: _byDivisorGpuOrderTesterTopDown = new ByDivisorGpuOrderTopDown(); break; + case ByDivisorKIncrementMode.BitContradiction: + _byDivisorGpuOrderTesterBitContradiction = new ByDivisorGpuOrderBitContradiction(); + break; + case ByDivisorKIncrementMode.BitTree: + _byDivisorGpuOrderTesterBitTree = new ByDivisorGpuOrderBitTree(); + break; default: _byDivisorGpuOrderTester = new MersenneNumberDivisorByDivisorCpuTesterWithGpuOrder(); break; @@ -545,6 +569,7 @@ private static ThreadLocal? MersenneTesters; ByDivisorKIncrementMode.Additive => "additive.kstate.faster", ByDivisorKIncrementMode.TopDown => "topdown.kstate.faster", ByDivisorKIncrementMode.BitContradiction => "bitcontradiction.kstate.faster", + ByDivisorKIncrementMode.BitTree => "bittree.kstate.faster", _ => "sequential.kstate.faster", }; string lmdbPath = Path.Combine(checksDir, Path.ChangeExtension(kRepoFile, ".lmdb")); @@ -664,7 +689,7 @@ private static ThreadLocal? MersenneTesters; CalculationResultHandler.InitializeOutputBuffer(); // Limit GPU concurrency only for prime checks (LL/NTT & GPU order scans). - GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads); + GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads * PerfectNumberConstants.ThreadsByAccelerator); // Configure batch size for GPU primality sieve PrimeTester.GpuBatchSize = gpuPrimeBatch; @@ -751,6 +776,19 @@ private static ThreadLocal? MersenneTesters; threadCount, byDivisorPrimeRange); break; + case ByDivisorKIncrementMode.BitTree: + Console.WriteLine("Using by-divisor bit-tree increment (CPU order)."); + MersenneNumberDivisorByDivisorTester.Run( + byDivisorCandidates, + ref _byDivisorCpuOrderTesterBitTree, + _byDivisorPreviousResults, + _byDivisorStartPrime, + static () => _lastCompositeP = true, + static () => _lastCompositeP = false, + HandleByDivisorResult, + threadCount, + byDivisorPrimeRange); + break; default: MersenneNumberDivisorByDivisorTester.Run( byDivisorCandidates, @@ -834,6 +872,19 @@ private static ThreadLocal? MersenneTesters; threadCount, byDivisorPrimeRange); break; + case ByDivisorKIncrementMode.BitTree: + Console.WriteLine("Using by-divisor bit-tree increment (Hybrid order)."); + MersenneNumberDivisorByDivisorTester.Run( + byDivisorCandidates, + ref _byDivisorHybridOrderTesterBitTree, + _byDivisorPreviousResults, + _byDivisorStartPrime, + static () => _lastCompositeP = true, + static () => _lastCompositeP = false, + HandleByDivisorResult, + threadCount, + byDivisorPrimeRange); + break; default: MersenneNumberDivisorByDivisorTester.Run( byDivisorCandidates, @@ -917,6 +968,19 @@ private static ThreadLocal? MersenneTesters; threadCount, byDivisorPrimeRange); break; + case ByDivisorKIncrementMode.BitTree: + Console.WriteLine("Using by-divisor bit-tree increment (GPU order)."); + MersenneNumberDivisorByDivisorTester.Run( + byDivisorCandidates, + ref _byDivisorGpuOrderTesterBitTree, + _byDivisorPreviousResults, + _byDivisorStartPrime, + static () => _lastCompositeP = true, + static () => _lastCompositeP = false, + HandleByDivisorResult, + threadCount, + byDivisorPrimeRange); + break; default: MersenneNumberDivisorByDivisorTester.Run( byDivisorCandidates, @@ -1327,6 +1391,7 @@ private static ThreadLocal? MersenneTesters; ByDivisorKIncrementMode.Additive => _byDivisorCpuOrderTesterAdditive.IsPrime(p, out detailedCheck, out _), ByDivisorKIncrementMode.TopDown => _byDivisorCpuOrderTesterTopDown.IsPrime(p, out detailedCheck, out _), ByDivisorKIncrementMode.BitContradiction => _byDivisorCpuOrderTesterBitContradiction.IsPrime(p, out detailedCheck, out _), + ByDivisorKIncrementMode.BitTree => _byDivisorCpuOrderTesterBitTree.IsPrime(p, out detailedCheck, out _), _ => _byDivisorCpuOrderTester.IsPrime(p, out detailedCheck, out _), }; } @@ -1341,6 +1406,7 @@ private static ThreadLocal? MersenneTesters; ByDivisorKIncrementMode.Additive => _byDivisorHybridOrderTesterAdditive.IsPrime(p, out detailedCheck, out _), ByDivisorKIncrementMode.TopDown => _byDivisorHybridOrderTesterTopDown.IsPrime(p, out detailedCheck, out _), ByDivisorKIncrementMode.BitContradiction => _byDivisorHybridOrderTesterBitContradiction.IsPrime(p, out detailedCheck, out _), + ByDivisorKIncrementMode.BitTree => _byDivisorHybridOrderTesterBitTree.IsPrime(p, out detailedCheck, out _), _ => _byDivisorHybridOrderTester.IsPrime(p, out detailedCheck, out _), }; } @@ -1353,6 +1419,7 @@ private static ThreadLocal? MersenneTesters; ByDivisorKIncrementMode.Additive => _byDivisorGpuOrderTesterAdditive.IsPrime(p, out detailedCheck, out _), ByDivisorKIncrementMode.TopDown => _byDivisorGpuOrderTesterTopDown.IsPrime(p, out detailedCheck, out _), ByDivisorKIncrementMode.BitContradiction => _byDivisorGpuOrderTesterBitContradiction.IsPrime(p, out detailedCheck, out _), + ByDivisorKIncrementMode.BitTree => _byDivisorGpuOrderTesterBitTree.IsPrime(p, out detailedCheck, out _), _ => _byDivisorGpuOrderTester.IsPrime(p, out detailedCheck, out _), }; } diff --git a/PerfectNumbers.Core.Tests/BitTreeSolverTests.cs b/PerfectNumbers.Core.Tests/BitTreeSolverTests.cs new file mode 100644 index 0000000..796c5fe --- /dev/null +++ b/PerfectNumbers.Core.Tests/BitTreeSolverTests.cs @@ -0,0 +1,43 @@ +using System.Numerics; +using FluentAssertions; +using PerfectNumbers.Core.Cpu; +using Xunit; + +namespace PerfectNumbers.Core.Tests; + +public class BitTreeSolverTests +{ + [Fact] + public void BitTree_FindsKnownDivisor() + { + // 2^11-1 = 2047 = 23 * 89, so k=1 yields divisor 23. + var tester = new MersenneNumberDivisorByDivisorCpuTesterWithForBitTreeDivisorSetForCpuOrder + { + DivisorLimit = new BigInteger(1000), + MinK = BigInteger.One, + }; + + bool isPrime = tester.IsPrime(11, out bool exhausted, out BigInteger divisor); + + isPrime.Should().BeFalse(); + exhausted.Should().BeTrue(); + divisor.Should().Be(23); + } + + [Fact] + public void BitTree_ReportsPrimeWhenNoDivisorFound() + { + // 2^5-1 = 31 is prime. + var tester = new MersenneNumberDivisorByDivisorCpuTesterWithForBitTreeDivisorSetForCpuOrder + { + DivisorLimit = new BigInteger(10_000), + MinK = BigInteger.One, + }; + + bool isPrime = tester.IsPrime(5, out bool exhausted, out BigInteger divisor); + + isPrime.Should().BeTrue(); + exhausted.Should().BeTrue(); + divisor.Should().Be(BigInteger.Zero); + } +} diff --git a/PerfectNumbers.Core/Cpu/MersenneCombinedDivisorScanner.cs b/PerfectNumbers.Core/Cpu/MersenneCombinedDivisorScanner.cs new file mode 100644 index 0000000..8af9db0 --- /dev/null +++ b/PerfectNumbers.Core/Cpu/MersenneCombinedDivisorScanner.cs @@ -0,0 +1,1987 @@ +using System; +using System.Buffers; +using System.Collections.Generic; +using System.Numerics; +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core.Cpu; + +/// +/// Deterministic streaming (column-wise) reconstruction scan for M_p = 2^p - 1, +/// exploring divisors in the family q = (2p)*k + 1 without materializing full q during search. +/// +/// IMPORTANT: +/// - This scanner is exact for the explored (q,k) search tree up to maxQBits. +/// - However, exploring "all possible" k for very large maxQBits is combinatorially explosive. +/// The intent is to use strong pruning to rule out branches early. +/// - For practical use, keep maxQBits in a regime you can actually explore (hundreds or a few thousands), +/// and rely on pruning/caches. +/// +public static class MersenneCombinedDivisorScannerFullStreamExtended +{ + // ===== DEBUG SETTINGS ===== + private static readonly bool DebugEnabled = true; + private static readonly BigInteger DebugTargetQ = + BigInteger.Parse("1209708008767"); + private const int DebugMaxBlocks = 20; // wystarczy (q ma ~11 bloków) + + private static bool DebugPrefixMatches(BigInteger qNow, int producedNibbles) + { + int producedBits = producedNibbles * 4; + // target ma ~41 bitów; powyżej nie ma sensu maskować dalej + if (producedBits >= 128) producedBits = 128; // bezpiecznie, debug tylko + + BigInteger mask = (BigInteger.One << producedBits) - BigInteger.One; + return (qNow & mask) == (DebugTargetQ & mask); + } + + private const int MaxColumnChoices = 2; + + + public enum ScanResult + { + RuledOut = 0, + Candidate = 1, + FoundDivisor = MaxColumnChoices + } + + // Hard safety ceiling for in-memory bitsets/windows in this implementation. + // You can raise it, but exploration cost grows extremely fast with maxQBits. + public const int MaxSupportedQBits = 8192; + + // Block-0 filter: q is odd and q ≡ 1 or 7 (mod 8). + // NOTE: Do NOT apply any “not divisible by 5” filter based on low bits; it is not valid in base-2 local form. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool IsAllowedQ0Nibble(byte q0) + { + // odd + if ((q0 & 0x1) == 0) return false; + // q mod 8 in {1,7} + int mod8 = q0 & 0x7; + return mod8 == 0x1 || mod8 == 0x7; + } + + /// + /// Deterministic streaming reconstruction scan for M_p = 2^p - 1. + /// + /// The scan explores q in the family q = (2p)*k + 1, where k is an unbounded nonnegative integer, + /// but truncates q to maxQBits. The scan maintains a rolling window of a-bits of size maxQBits, + /// sufficient to evaluate the convolution exactly for all columns j up to p+maxQBits. + /// + /// Return values: + /// - FoundDivisor: found a q within maxQBits such that the reconstructed a satisfies a*q = 2^p - 1 (in all columns). + /// - RuledOut: no such q exists within the explored tree and limits. + /// - Candidate: scan did not find a divisor within the explored q-bit budget (maxQBits), + /// or maxQBits exceeds MaxSupportedQBits. + /// + /// Notes: + /// - This is an exact column-wise decision procedure for each explored branch, but exploration is exponential in maxQBits. + /// - Optional ModPow verification is performed only when FoundDivisor is reached. + /// + public static ScanResult TryFindDivisorByStreamingScan( + ulong prime, + int maxQBits, + bool verifyWithModPow, + out BigInteger foundQ) + { + foundQ = BigInteger.Zero; + + if (prime == 0) return ScanResult.RuledOut; + ArgumentOutOfRangeException.ThrowIfLessThan(maxQBits, MaxColumnChoices); + if (maxQBits > MaxSupportedQBits) return ScanResult.Candidate; + + int pBits = checked((int)prime); + int maxBlocks = (maxQBits + 3) / 4; + + var qGen = new QNibbleGenerator16(prime); + + // Layered (block-ordered) exploration: always exhaust smaller block indices first. + // This is important for correctness of conservative block-level pruning and for predictable progress. + var openByBlock = new Dictionary>(capacity: 4096); + var termByBlock = new Dictionary>(capacity: 4096); + var activeBlocks = new SortedSet(); + + // Per-block conservative global pruning stats (computed over the set of *actually reachable* reduced states at that block). + var blockPrune = new Dictionary(capacity: 1024); + + // Per-state allowed-a cache (must remain state-sensitive; see AllowedMaskKey). + var allowedMaskCache = new Dictionary(capacity: 1 << 16); + + void EnqueueState(State s) + { + var dict = s.KTerminated ? termByBlock : openByBlock; + if (!dict.TryGetValue(s.Block, out var q)) + { + q = new Queue(); + dict[s.Block] = q; + } + q.Enqueue(s); + activeBlocks.Add(s.Block); + } + + bool TryDequeueState(out State s) + { + s = default; + if (activeBlocks.Count == 0) return false; + int b = activeBlocks.Min; + + if (termByBlock.TryGetValue(b, out var tq) && tq.Count > 0) + { + s = tq.Dequeue(); + } + else if (openByBlock.TryGetValue(b, out var oq) && oq.Count > 0) + { + s = oq.Dequeue(); + } + else + { + // Should not happen, but keep structure consistent. + activeBlocks.Remove(b); + return false; + } + + // Cleanup active block if both queues are empty. + bool termEmpty = !termByBlock.TryGetValue(b, out var tq2) || tq2.Count == 0; + bool openEmpty = !openByBlock.TryGetValue(b, out var oq2) || oq2.Count == 0; + if (termEmpty && openEmpty) activeBlocks.Remove(b); + + return true; + } + + var init = State.CreateInitial(QNibbleGenerator16.InitialState, maxQBits); + EnqueueState(init); + + Span counts = stackalloc int[16]; + Span qNib = stackalloc byte[16]; + Span nextState = stackalloc QGenState16[16]; + Span aOrder = stackalloc byte[16]; // reusable buffer for nibble ordering (avoid stackalloc in loops) + int lastLoggedBlock = -1; + + State st; + while (TryDequeueState(out st)) + { + + // Intentionally no internal max-nodes limit. + + if (DebugEnabled && st.Block != lastLoggedBlock) + { + lastLoggedBlock = st.Block; + Console.WriteLine($"[PROGRESS] p={prime}, reached block={st.Block}"); + } + + // If we have finished producing all q nibbles (to maxQBits), finish scanning remaining columns. + if (st.Block >= maxBlocks) + { + if (FinishScanToEnd(pBits, maxQBits, ref st) && st.Carry == 0) + { + // Found a fully consistent reconstruction for this branch. + BigInteger qVal = st.QBits.ToBigInteger(); + if (qVal > BigInteger.One) + { + foundQ = qVal; + Console.WriteLine($"[stream-scan] Found divisor candidate q(bits<={maxQBits}) for p={prime}"); + if (verifyWithModPow) + { + bool ok = BigInteger.ModPow(MaxColumnChoices, prime, foundQ) == BigInteger.One; + Console.WriteLine($"[stream-scan] ModPow verification: {(ok ? "OK" : "FAILED")}"); + } + + return ScanResult.FoundDivisor; + } + } + + continue; + } + + int startBit = st.Block << 2; + + // Enumerate possible k-nibbles for this q-nibble; group by qNibble to avoid recomputation + // of allowed a choices (which depend on qNibble and carry). + + // Reset tables for this state. + qNib.Clear(); + nextState.Clear(); + counts.Clear(); + + // byte forcedKn; + // int ni = st.QState.NibbleIndex; + + // if (ni == 0) forcedKn = 0xF; + // else if (ni == 1) forcedKn = 0x1; + // else if (ni == 2) forcedKn = 0x1; + // else if (ni == 3) forcedKn = 0x1; + // else forcedKn = 0x0; + + // for (byte kn = 0; kn <= 15; kn++) + // { + // if (kn != forcedKn) continue; + + // for (int i = 0; i < 16; i++) qNib[i] = 0xFF; + // for (int i = 0; i < 16; i++) nextState[i] = default; // opcjonalnie + // qGen.ComputeNext(st.QState, forcedKn, out byte qn, out QGenState16 ns); + // qNib[forcedKn] = qn; + // nextState[forcedKn] = ns; + // counts[qn] = 1; + // } + + // Sentinel init: mark unused kn entries as invalid (0xFF) to avoid accidental matches when qn==0. + for (int i = 0; i < 16; i++) qNib[i] = 0xFF; + // nextState default is fine; we will only read it when qNib[kn] == qn (so it was written). + for (int i = 0; i < 16; i++) nextState[i] = default; + + // Jeśli k już zakończone, jedyny dozwolony kn=0 + if (st.KTerminated) + { + byte kn = 0; + qGen.ComputeNext(st.QState, kn, out byte qn, out QGenState16 ns); + + // Aktualizuj stan generatora + st.QState = ns; + + // Flush okna k: po terminacji trzeba wykonać dokładnie KWinLen kroków kn=0, + // żeby wszystkie historyczne nibble k przestały wpływać na q. + if (st.KFlushRemaining > 0) + st.KFlushRemaining--; + + // Dopisz q nibble do bitsetu q + ApplyQNibble(ref st, startBit, qn, maxQBits); + + // After k termination we keep streaming q with kn=0. + // We can only declare "end of q" once the k-window is flushed AND we observe a stable run of zero q-nibbles. + // This avoids premature freezing when a later carry/window interaction could still produce a non-zero nibble. + if (qn == 0) st.PostTermZeroQNibbles++; + else st.PostTermZeroQNibbles = 0; + + int requiredZeroRun = qGen.KWinLen + 1; + if (st.KFlushRemaining == 0 && st.QState.Carry16 == 0 && st.PostTermZeroQNibbles >= requiredZeroRun) + { + Console.WriteLine($"[QFROZEN] p={prime} at block={st.Block} nib={st.QState.NibbleIndex} carry16={st.QState.Carry16} zeroRun={st.PostTermZeroQNibbles}"); + + // Jump to final phase (FinishScanToEnd) without materializing full q. + var done = st.CloneShallow(); + done.Block = maxBlocks; + EnqueueState(done); + continue; + } + + qNib[kn] = qn; + nextState[kn] = ns; + counts[qn] = 1; + } + else + { + // Enumerate all possible k-nibbles at this position (optimized: compute base once per state). + qGen.ComputeAllNext(st.QState, qNib, nextState, counts); + + // Dodatkowa gałąź: “zakończ k TERAZ” (czyli od tego miejsca kn=0 zawsze). + // Realizujemy to poprzez dopuszczenie kn=0, ale z flagą KTerminated w dziecku. + // Nie musisz robić dodatkowego liczenia counts; obsłużymy to przy tworzeniu child. + } + + + // Conservative per-block pruning: + // Track which q-nibbles are impossible for *all reachable reduced states* at this block. + // This is safe because we explore blocks in increasing order; once we move past a block, no new states for it appear. + ulong qPrefixSig = st.QBits.Signature64UpTo(startBit - 1); + ulong aSigNow = st.AWindow.Signature64(); + var reducedKey = new ReducedStateKey(pBits, maxQBits, st.Block, st.Carry, qPrefixSig, aSigNow); + + if (!blockPrune.TryGetValue(st.Block, out var pruneStats)) + { + pruneStats = new BlockPruneStats(); + blockPrune[st.Block] = pruneStats; + } + pruneStats.RegisterState(reducedKey); + for (int qnVal = 0; qnVal < 16; qnVal++) + { + if (counts[qnVal] == 0) continue; + byte qn = (byte)qnVal; + + if (st.Block == 0 && !IsAllowedQ0Nibble(qn)) + continue; + + // Global (block-level) forbidden q-nibble mask, conservatively derived from all reachable reduced states. + if ((pruneStats.GlobalForbiddenMask & (1 << qn)) != 0) + continue; + + // Early contradiction pruning (interval-based carry propagation), inspired by BitContradictionSolver. + // Important: bits of 'a' outside the window are treated as Unknown (never as 0), so this pruning cannot + // introduce false negatives. + if (QuickRejectQNibbleByBounds(pBits, maxQBits, in st, startBit, qn)) + { + pruneStats.MarkRejected(reducedKey, qn); + continue; + } + + // Additional early pruning: derive a conservative candidate mask for the 4-bit a-nibble at this block + // by checking (using the same bounds logic) whether each a_j in this block can be 0/1. + // If a bit is forced (only one value feasible), we restrict the candidate a-nibble mask. + // We also derive a conservative *preference* mask/value for branch ordering (no pruning). + // Bits outside the a-window are treated as Unknown, so this cannot introduce false negatives. + var boundsInfo = ComputeCandidateANibbleMaskAndPreferencesByBounds(pBits, maxQBits, in st, startBit, qn); + ushort preMask = boundsInfo.CandidateMask; + if (preMask == 0) + { + pruneStats.MarkRejected(reducedKey, qn); + continue; + } + + // Precompute which a-nibbles are compatible with this q-nibble and current carry, + // using LUT for delta contribution plus exact base via the rolling a-window and current q-bitset. + int sigEndBit = Math.Min(startBit + 3, pBits + maxQBits); + ulong qSig = st.QBits.Signature64WithNibbleUpTo(sigEndBit, startBit, qn, maxQBits); + ulong aSig = st.AWindow.Signature64(); + var key = new AllowedMaskKey(pBits, maxQBits, st.Block, st.Carry, qn, qSig, aSig, preMask); + + if (!allowedMaskCache.TryGetValue(key, out ushort allowedAMask)) + { + allowedAMask = ComputeAllowedANibbleMaskFiltered(pBits, maxQBits, ref st, startBit, qn, preMask); + allowedMaskCache[key] = allowedAMask; + } + + if (allowedAMask == 0 && DebugEnabled && st.QState.NibbleIndex <= DebugMaxBlocks) + { + var tmp = st.CloneShallow(); + ApplyQNibble(ref tmp, startBit, qn, maxQBits); + BigInteger qTmp = tmp.QBits.ToBigInteger(); + + int produced = st.QState.NibbleIndex + 1; + if (DebugPrefixMatches(qTmp, produced)) + Console.WriteLine($"[DEBUG CUT-A-PREFIX] block={produced} qLow={qTmp} qNib=0x{qn:X}"); + } + + if (allowedAMask == 0) + { + pruneStats.MarkRejected(reducedKey, qn); + continue; + } + + // Expand all k-nibbles mapping to this q-nibble. + // We will apply branch ordering to a-nibbles using bounds-derived preferences (no pruning). + byte prefMask = boundsInfo.PrefMask; + byte prefValue = boundsInfo.PrefValue; + for (int kn = 0; kn < 16; kn++) + { + if (qNib[kn] != qn) continue; + if (st.Block == 0 && kn == 0) continue; + + // if (DebugEnabled && kn == 0 && qNib[kn] == qn) + // { + // Console.WriteLine($"[KN0] p={prime} block={st.Block} qn=0x{qn:X} canTerminateHere={!st.KTerminated}"); + // } + + var child = st.CloneShallow(); // share arrays via copy-on-write where safe (implemented as cloning bitsets) + child.QState = nextState[kn]; + + // Jeśli nie jesteśmy w trybie terminated i wybraliśmy kn==0, + // to rozgałęź na dwa warianty: + // - zwykły: k nadal może mieć wyższe nible, + // - terminated: uznajemy, że k kończy się tutaj. + bool canTerminateHere = !st.KTerminated && kn == 0 && st.Block > 0; + child.KTerminated = st.KTerminated; // default + + // wariant 1: zwykły (nie terminujemy) + ApplyQNibble(ref child, startBit, qn, maxQBits); + + // if (DebugEnabled && child.QState.NibbleIndex <= DebugMaxBlocks) + // { + // BigInteger qNow = child.QBits.ToBigInteger(); + // if (DebugPrefixMatches(qNow, child.QState.NibbleIndex)) + // { + // Console.WriteLine($"[DEBUG HIT-Q-ONLY-PREFIX] block={child.QState.NibbleIndex} qLow={qNow}"); + // } + // } + + // ===== wariant 1: normalnie kontynuujemy (KTerminated=false) ===== + int aCount = BuildOrderedNibbleList(allowedAMask, prefMask, prefValue, aOrder); + for (int ai = 0; ai < aCount; ai++) + { + int aNib = aOrder[ai]; + + var grand = child.CloneShallow(); + ApplyANibble(ref grand, startBit, (byte)aNib); + + foreach (var branched in StepColumnsForBlockMulti(pBits, maxQBits, grand, startBit)) + { + var pushState = branched; + pushState.Block = pushState.QState.NibbleIndex; + EnqueueState(pushState); + } + } + + // wariant 2: terminujemy k (tylko jeśli kn==0 i jeszcze nie terminated) + if (canTerminateHere) + { + var termChild = st.CloneShallow(); + termChild.QState = nextState[kn]; + termChild.KFlushRemaining = (byte)qGen.KWinLen; + termChild.PostTermZeroQNibbles = 0; + termChild.KTerminated = true; + + // Console.WriteLine($"[KTERM] p={prime} creating terminated branch at block={st.Block}, nib={st.QState.NibbleIndex}"); + + ApplyQNibble(ref termChild, startBit, qn, maxQBits); + + // kontynuujemy stan z terminacją k + int aCount2 = BuildOrderedNibbleList(allowedAMask, prefMask, prefValue, aOrder); + for (int ai2 = 0; ai2 < aCount2; ai2++) + { + int aNib = aOrder[ai2]; + + var grand = termChild.CloneShallow(); + ApplyANibble(ref grand, startBit, (byte)aNib); + + foreach (var branched in StepColumnsForBlockMulti(pBits, maxQBits, grand, startBit)) + { + var pushState = branched; + pushState.Block = pushState.QState.NibbleIndex; + EnqueueState(pushState); + } + } + } + } + } + } + + return ScanResult.RuledOut; + } + + + private static ushort ComputeAllowedANibbleMask(int pBits, int maxQBits, ref State st, int startBit, byte qNibble) + { + return ComputeAllowedANibbleMaskFiltered(pBits, maxQBits, ref st, startBit, qNibble, 0xFFFF); + } + + private static ushort ComputeAllowedANibbleMaskFiltered(int pBits, int maxQBits, ref State st, int startBit, byte qNibble, ushort candidateMask) + { + // Correctness-first implementation: + // - never mutate the caller state + // - pre-apply qNibble once to improve performance (it is the same for all aNib candidates) + + ushort mask = 0; + + var baseState = st.CloneShallow(); + ApplyQNibble(ref baseState, startBit, qNibble, maxQBits); + + for (int aNib = 0; aNib < 16; aNib++) + { + if ((candidateMask & (1 << aNib)) == 0) continue; + // In block0, enforce a0 bit0 = 1 (since a0=1). + if (startBit == 0 && (aNib & 1) == 0) continue; + + var tmp = baseState.CloneShallow(); + ApplyANibble(ref tmp, startBit, (byte)aNib); + + bool ok = false; + foreach (var _ in StepColumnsForBlockMulti(pBits, maxQBits, tmp, startBit)) + { + ok = true; + break; + } + if (ok) mask |= (ushort)(1 << aNib); + } + + return mask; + } + + // ------------------------------------------------------------------------ + // Early a-nibble candidate mask derivation (forced bits) via bounds + // ------------------------------------------------------------------------ + // This uses the same conservative parity/bounds logic as QuickRejectQNibbleByBounds, but instead of only + // rejecting qNibble, it attempts to determine whether each a_j in this 4-column block can be 0/1. + // If exactly one of {0,1} is feasible for a_j (under the conservative bounds model), that bit is treated as forced + // and used to restrict the candidate a-nibble space. + // + // IMPORTANT: This MUST be conservative. It may fail to detect some forced bits (false negatives), + // but it must never rule out a real solution (no false positives on forcing). Our feasibility checks use + // an over-approximate interval [minSum,maxSum] and therefore only declare "impossible" when parity cannot + // be achieved anywhere in that interval. + private readonly struct BoundsANibbleInfo + { + public readonly ushort CandidateMask; + public readonly byte PrefMask; // bits within the 4-bit nibble for which we have a preference + public readonly byte PrefValue; // preferred bit values for those positions + + public BoundsANibbleInfo(ushort candidateMask, byte prefMask, byte prefValue) + { + CandidateMask = candidateMask; + PrefMask = prefMask; + PrefValue = prefValue; + } + } + + private static BoundsANibbleInfo ComputeCandidateANibbleMaskAndPreferencesByBounds(int pBits, int maxQBits, in State st, int startBit, byte qNibble) + { + // Start with all 16 nibbles allowed. + ushort mask = 0xFFFF; + byte prefMask = 0; + byte prefValue = 0; + + long carryMin = st.Carry; + long carryMax = st.Carry; + + for (int j = startBit; j <= startBit + 3; j++) + { + if (j > pBits + maxQBits) break; + int requiredBit = (j < pBits) ? 1 : 0; + + // Compute forced/possible contributions from offsets t>0 (excluding a_j itself). + long forcedOther = 0; + long possibleOther = 0; + + int maxT = Math.Min(j, maxQBits - 1); + for (int t = 1; t <= maxT; t++) + { + if (GetQBitWithNibble(in st, startBit, qNibble, t, maxQBits) == 0) continue; + int aIndex = j - t; + if (aIndex < 0) continue; + + if (st.AWindow.TryGetKnownBit(aIndex, out int abit)) + { + if (abit != 0) + { + forcedOther++; + possibleOther++; + } + } + else + { + possibleOther++; + } + } + + bool ajKnown = st.AWindow.TryGetKnownBit(j, out int ajBit); + bool can0, can1; + if (ajKnown) + { + can0 = ajBit == 0 && IsParityReachableForCarryRange(carryMin, carryMax, forcedOther, possibleOther, 0, requiredBit); + can1 = ajBit == 1 && IsParityReachableForCarryRange(carryMin, carryMax, forcedOther, possibleOther, 1, requiredBit); + } + else + { + can0 = IsParityReachableForCarryRange(carryMin, carryMax, forcedOther, possibleOther, 0, requiredBit); + can1 = IsParityReachableForCarryRange(carryMin, carryMax, forcedOther, possibleOther, 1, requiredBit); + } + + if (!can0 && !can1) + return new BoundsANibbleInfo(0, 0, 0); + + // If only one value is feasible (under our conservative model), restrict the nibble mask. + int bitInNibble = j - startBit; + if (bitInNibble >= 0 && bitInNibble < 4) + { + if (!can0 && can1) + { + mask &= (ushort)ForceNibbleBit(bitInNibble, 1); + prefMask |= (byte)(1 << bitInNibble); + prefValue |= (byte)(1 << bitInNibble); + } + else if (can0 && !can1) + { + mask &= (ushort)ForceNibbleBit(bitInNibble, 0); + prefMask |= (byte)(1 << bitInNibble); + prefValue &= (byte)~(1 << bitInNibble); + } + else + { + // Both possible: derive a *preference* for branch ordering. + // Prefer the value that yields a narrower next carry range under the same conservative model. + long n0Min, n0Max, n1Min, n1Max; + bool ok0 = TryPropagateCarryBounds(carryMin, carryMax, forcedOther + 0, possibleOther + 0, requiredBit, out n0Min, out n0Max); + bool ok1 = TryPropagateCarryBounds(carryMin, carryMax, forcedOther + 1, possibleOther + 1, requiredBit, out n1Min, out n1Max); + // ok0/ok1 should match can0/can1, but keep defensive. + if (ok0 && ok1) + { + long w0 = n0Max - n0Min; + long w1 = n1Max - n1Min; + int prefer = (w1 < w0) ? 1 : (w0 < w1 ? 0 : 0); + prefMask |= (byte)(1 << bitInNibble); + if (prefer == 1) prefValue |= (byte)(1 << bitInNibble); + else prefValue &= (byte)~(1 << bitInNibble); + } + } + } + + // Update carry range for next column using *unknown* a_j if not known and not forced. + // For safety, we propagate using the full [min,max] interval that covers all feasible choices of a_j. + long forced = forcedOther; + long possible = possibleOther; + if (ajKnown) + { + if (ajBit != 0) { forced++; possible++; } + } + else + { + // if both 0 and 1 are feasible, treat as unknown -> possible++ + // if forced, we can reflect it in forced/possible. + if (!can0 && can1) { forced++; possible++; } + else if (can0 && !can1) { /* forced 0 -> nothing */ } + else { possible++; } + } + + if (!TryPropagateCarryBounds(carryMin, carryMax, forced, possible, requiredBit, out long nextMin, out long nextMax)) + return new BoundsANibbleInfo(0, 0, 0); + carryMin = nextMin; + carryMax = nextMax; + } + + return new BoundsANibbleInfo(mask, prefMask, prefValue); + } + + private static ushort ComputeCandidateANibbleMaskByBounds(int pBits, int maxQBits, in State st, int startBit, byte qNibble) + => ComputeCandidateANibbleMaskAndPreferencesByBounds(pBits, maxQBits, in st, startBit, qNibble).CandidateMask; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool IsParityReachableForCarryRange(long carryMin, long carryMax, long forcedOther, long possibleOther, int aJ, int requiredBit) + { + int parity = requiredBit & 1; + long minSum = carryMin + forcedOther + aJ; + long maxSum = carryMax + possibleOther + aJ; + long minAligned = AlignUpToParity(minSum, parity); + long maxAligned = AlignDownToParity(maxSum, parity); + return minAligned <= maxAligned; + } + + // Returns a 16-bit mask allowing only those a-nibbles whose bit 'b' equals requiredValue. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int ForceNibbleBit(int b, int requiredValue) + { + int mask = 0; + for (int n = 0; n < 16; n++) + { + if (((n >> b) & 1) == requiredValue) + mask |= 1 << n; + } + return mask; + } + + // Build an ordered list of a-nibbles from a bitmask, using bounds-derived preferences for branch ordering. + // This never prunes: it only reorders exploration to reach contradictions / valid branches earlier. + // + // Scoring: prefer nibbles matching preferred bit-values on positions indicated by prefMask. + private static int BuildOrderedNibbleList(ushort allowedMask, byte prefMask, byte prefValue, Span outList) + { + int count = 0; + for (int n = 0; n < 16; n++) + { + if ((allowedMask & (1 << n)) == 0) continue; + outList[count++] = (byte)n; + } + + if (count <= 1 || prefMask == 0) + return count; + + // Small selection sort (count <= 16). + for (int i = 0; i < count - 1; i++) + { + int best = i; + int bestScore = ScoreNibble(outList[i], prefMask, prefValue); + for (int j = i + 1; j < count; j++) + { + int s = ScoreNibble(outList[j], prefMask, prefValue); + if (s > bestScore) + { + best = j; + bestScore = s; + } + } + if (best != i) + { + (byte a, byte b) = (outList[i], outList[best]); + outList[i] = b; + outList[best] = a; + } + } + + return count; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int ScoreNibble(byte nib, byte prefMask, byte prefValue) + { + uint matches = (uint)(~(nib ^ prefValue) & prefMask) & 0xFu; + return (int)BitOperations.PopCount(matches); + } + + // ------------------------------------------------------------------------ + // Early contradiction pruning (interval-based carry propagation) + // ------------------------------------------------------------------------ + // This is a conservative pre-filter inspired by BitContradictionSolver.TryPropagateCarry. + // We do NOT commit any new a-bits here. We only prove that a given (state, qNibble) cannot + // satisfy the required parity constraints in any of the next 4 columns, even if all currently + // unknown a-bits were chosen adversarially. + // + // Key correctness constraint (per user request): a-bits outside the window are treated as UNKNOWN, + // never as 0. Therefore this pruning cannot eliminate a valid solution branch. + private static bool QuickRejectQNibbleByBounds(int pBits, int maxQBits, in State st, int startBit, byte qNibble) + { + long carryMin = st.Carry; + long carryMax = st.Carry; + + for (int j = startBit; j <= startBit + 3; j++) + { + if (j > pBits + maxQBits) break; + int requiredBit = (j < pBits) ? 1 : 0; + + // Compute forced/possible ones contributed to column j excluding carry. + // We model: sum = carry + a_j + Σ_{t>0, q_t=1} a_{j-t} + long forced = 0; + long possible = 0; + + int maxT = Math.Min(j, maxQBits - 1); + for (int t = 1; t <= maxT; t++) + { + if (GetQBitWithNibble(in st, startBit, qNibble, t, maxQBits) == 0) continue; + int aIndex = j - t; + if (aIndex < 0) continue; + + if (st.AWindow.TryGetKnownBit(aIndex, out int abit)) + { + if (abit != 0) + { + forced++; + possible++; + } + } + else + { + // Unknown a-bit outside window or not yet assigned. + possible++; + } + } + + // Include a_j itself. + if (st.AWindow.TryGetKnownBit(j, out int aJKnown)) + { + if (aJKnown != 0) + { + forced++; + possible++; + } + } + else + { + possible++; + } + + if (!TryPropagateCarryBounds(carryMin, carryMax, forced, possible, requiredBit, out long nextMin, out long nextMax)) + return true; + + carryMin = nextMin; + carryMax = nextMax; + } + + return false; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetQBitWithNibble(in State st, int startBit, byte qNibble, int bitPos, int maxQBits) + { + if (bitPos == 0) return 1; // enforce q0=1 always + if ((uint)bitPos >= (uint)maxQBits) return 0; + + int baseBit = st.QBits.Get(bitPos); + if (bitPos >= startBit && bitPos < startBit + 4) + { + int b = bitPos - startBit; + int nibBit = (qNibble >> b) & 1; + return baseBit | nibBit; + } + return baseBit; + } + + private static bool TryPropagateCarryBounds(long carryMin, long carryMax, long forcedOnes, long possibleOnes, int requiredBit, + out long nextCarryMin, out long nextCarryMax) + { + // Equivalent in spirit to BitContradictionSolver.TryPropagateCarry, but we track a carry *range* + // and return the union of feasible next carry values. + int parity = requiredBit & 1; + + nextCarryMin = long.MaxValue; + nextCarryMax = long.MinValue; + + for (long carry = carryMin; carry <= carryMax; carry++) + { + long minSum = carry + forcedOnes; + long maxSum = carry + possibleOnes; + + long minAligned = AlignUpToParity(minSum, parity); + long maxAligned = AlignDownToParity(maxSum, parity); + if (minAligned > maxAligned) continue; + + long cMin = (minAligned - requiredBit) >> 1; + long cMax = (maxAligned - requiredBit) >> 1; + + if (cMin < nextCarryMin) nextCarryMin = cMin; + if (cMax > nextCarryMax) nextCarryMax = cMax; + } + + if (nextCarryMin == long.MaxValue) + { + nextCarryMin = 0; + nextCarryMax = 0; + return false; + } + + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static long AlignUpToParity(long value, int parity) + => ((value & 1) == parity) ? value : value + 1; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static long AlignDownToParity(long value, int parity) + => ((value & 1) == parity) ? value : value - 1; + + private static bool StepColumnsForBlock(int pBits, int maxQBits, ref State st, int startBit) + { + int endBit = Math.Min(startBit + 3, pBits + maxQBits); + while (st.BitIndex <= endBit) + { + if (!StepSingleColumn(pBits, maxQBits, ref st)) + return false; + st.BitIndex++; + } + return true; + } + + private static IEnumerable StepColumnsForBlockMulti(int pBits, int maxQBits, State seed, int startBit) + { + int endBit = Math.Min(startBit + 3, pBits + maxQBits); + + // Lokalny work-stack: eksplorujemy tylko 4 kolumny. + // Maksymalnie 2^4 = 16 stanów na seed (zwykle znacznie mniej). + var work = new FixedCapacityStack(capacity: 32); + work.Push(seed); + + int[] choices = new int[MaxColumnChoices]; + while (work.Count > 0) + { + var st = work.Pop(); + + // Jeśli przeszliśmy już cały blok, zwracamy stan jako wynik. + if (st.BitIndex > endBit) + { + yield return st; + continue; + } + + int j = st.BitIndex; + int requiredBit = (j < pBits) ? 1 : 0; + + int sumOther = st.Carry; + + int maxT = Math.Min(j, maxQBits - 1); + + // policz wkład z q_t=1 (t>0): add a_{j-t} + // (bez capture ref – używamy lokalnych) + var qBits = st.QBits; + var aWin = st.AWindow; + + foreach (int t in qBits.EnumerateSetBitsUpTo(maxT)) + { + if (t == 0) continue; + int aIndex = j - t; + if (aIndex < 0) continue; + if (aWin.GetBitRelative(aIndex) != 0) + sumOther++; + } + + // planned? + bool hasPlanned = aWin.TryGetPlanned(j, out int plannedAJ); + + // Policz, które aJ są dopuszczalne dla tej kolumny + int choiceCount = 0; + + for (int aJ = 0; aJ <= 1; aJ++) + { + if (hasPlanned && aJ != plannedAJ) continue; + + int sum = sumOther + aJ; + if ((sum & 1) != requiredBit) + continue; + + choices[choiceCount++] = aJ; + } + + if (choiceCount == 0) + { + // dead end for this local state + continue; + } + + // Dla każdej dopuszczalnej wartości aJ twórz stan następny. + // Pierwszy wariant modyfikujemy "in place", kolejne klonujemy. + ScheduleState(work, choices[0], j, requiredBit, sumOther, st); + if (choiceCount == 2) + { + ScheduleState(work, choices[1], j, requiredBit, sumOther, st.CloneShallow()); + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static void ScheduleState(FixedCapacityStack work, int aJ, int j, int requiredBit, int sumOther, State next) + { + int sum = sumOther + aJ; + int nextCarry = (requiredBit == 1) ? ((sum - 1) >> 1) : (sum >> 1); + + next.Carry = nextCarry; + next.AWindow.PushBit(j, aJ); + next.BitIndex = j + 1; + + work.Push(next); + } + + + private static bool FinishScanToEnd(int pBits, int maxQBits, ref State st) + { + int end = pBits + 1; + + while (st.BitIndex <= end) + { + int startBit = (st.BitIndex >> MaxColumnChoices) << MaxColumnChoices; + bool advanced = false; + + foreach (var next in StepColumnsForBlockMulti(pBits, maxQBits, st, startBit)) + { + st = next; + advanced = true; + break; // deterministycznie wybieramy jedną spójną ścieżkę + } + + if (!advanced) + return false; + } + + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool StepSingleColumn(int pBits, int maxQBits, ref State st) + { + int j = st.BitIndex; + int requiredBit = (j < pBits) ? 1 : 0; + + int sumOther = st.Carry; + + // contributions from q_t=1 for t>0: add a_{j-t} + // iterate over set bits of q (excluding t=0) + int maxT = Math.Min(j, maxQBits - 1); + + // NOTE: 'st' is a ref parameter; do not capture it in a lambda (CS1628). + // Copy required members to locals. + ABitsWindow aWindow = st.AWindow; + int jj = j; + int baseJ = aWindow.CurrentJ; // last committed column (typically j-1 while processing column j) + + foreach (int t in st.QBits.EnumerateSetBitsUpTo(maxT)) + { + if (t == 0) continue; + int aIndex = jj - t; + if (aIndex < 0) continue; + if (aWindow.GetBitRelative(aIndex) != 0) + sumOther++; + } + + // decide a_j by parity, because q0=1 is always set. + bool anyOk = false; + int newCarry = 0; + + bool hasPlanned = st.AWindow.TryGetPlanned(j, out int plannedAJ); + for (int aJ = 0; aJ <= 1; aJ++) + { + if (hasPlanned && aJ != plannedAJ) continue; + + int sum = sumOther + aJ; + if ((sum & 1) != requiredBit) + continue; + + newCarry = (requiredBit == 1) ? ((sum - 1) >> 1) : (sum >> 1); + st.AWindow.PushBit(j, aJ); + anyOk = true; + break; + } + + if (!anyOk) + { + if (DebugEnabled && st.Block <= DebugMaxBlocks) + { + BigInteger qPrefix = st.QBits.ToBigInteger(); + if (qPrefix == DebugTargetQ) + { + Console.WriteLine( + $"[DEBUG CUT-BIT] block={st.Block}, col={j}, sumOther={sumOther}, carryIn={st.Carry}, hasPlanned={hasPlanned}, plannedAJ={plannedAJ}" + ); + } + } + return false; + } + + st.Carry = newCarry; + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void ApplyQNibble(ref State st, int startBit, byte qNibble, int maxQBits) + { + for (int b = 0; b < 4; b++) + { + int bitPos = startBit + b; + if (bitPos >= maxQBits) continue; + int bit = (qNibble >> b) & 1; + if (bit != 0) + st.QBits.Set(bitPos); + } + st.QBits.Set(0); // enforce q0=1 always + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void ApplyANibble(ref State st, int startBit, byte aNibble) + { + // Apply 4 bits a_{startBit..startBit+3} into the rolling window. + // The window is addressed by absolute j; we push bits in order as columns advance, not at nibble application time. + // Here we store "scheduled" bits into a sparse map (only within window), so StepSingleColumn can read them. + for (int b = 0; b < 4; b++) + { + int bitPos = startBit + b; + int bit = (aNibble >> b) & 1; + st.AWindow.SetPlanned(bitPos, bit); + } + } + + + private readonly struct AllowedMaskKey( + int pBits, + int maxQBits, + int block, + int carry, + byte qNibble, + ulong qSig, + ulong aSig, + ushort candidateMask) : IEquatable + { + public readonly int PBits = pBits; + public readonly int MaxQBits = maxQBits; + public readonly int Block = block; + public readonly int Carry = carry; + public readonly byte QNibble = qNibble; + public readonly ulong QSig = qSig; + public readonly ulong ASig = aSig; + public readonly ushort CandidateMask = candidateMask; + + public bool Equals(AllowedMaskKey other) + => PBits == other.PBits + && MaxQBits == other.MaxQBits + && Block == other.Block + && Carry == other.Carry + && QNibble == other.QNibble + && QSig == other.QSig + && ASig == other.ASig + && CandidateMask == other.CandidateMask; + + public override bool Equals(object? obj) => obj is AllowedMaskKey k && Equals(k); + + public override int GetHashCode() + => HashCode.Combine(PBits, MaxQBits, Block, Carry, QNibble, QSig, ASig, CandidateMask); + } + + + + // Reduced state key for conservative per-block pruning. + // This key must be strong enough to avoid conflating states that can yield different allowed masks. + private readonly struct ReducedStateKey( + int pBits, + int maxQBits, + int block, + int carry, + ulong qPrefixSig, + ulong aSig) : IEquatable + { + public readonly int PBits = pBits; + public readonly int MaxQBits = maxQBits; + public readonly int Block = block; + public readonly int Carry = carry; + public readonly ulong QPrefixSig = qPrefixSig; + public readonly ulong ASig = aSig; + + public bool Equals(ReducedStateKey other) + => PBits == other.PBits + && MaxQBits == other.MaxQBits + && Block == other.Block + && Carry == other.Carry + && QPrefixSig == other.QPrefixSig + && ASig == other.ASig; + + public override bool Equals(object? obj) => obj is ReducedStateKey k && Equals(k); + + public override int GetHashCode() => HashCode.Combine(PBits, MaxQBits, Block, Carry, QPrefixSig, ASig); + } + + // Conservative per-block pruning statistics: + // For a given block, if a q-nibble is rejected (allowed mask == 0) for every reachable reduced state at that block, + // we can safely skip it for the remainder of the run. + private sealed class BlockPruneStats + { + private readonly Dictionary _stateRejectedMask = new(capacity: 1024); + + public int SeenStates { get; private set; } + public ushort GlobalForbiddenMask { get; private set; } + + // For each qNibble d, count how many reduced states have already rejected it. + private readonly int[] _rejectCounts = new int[16]; + + public void RegisterState(in ReducedStateKey key) + { + if (_stateRejectedMask.ContainsKey(key)) return; + + _stateRejectedMask[key] = 0; + SeenStates++; + + // New state potentially invalidates previously "global" forbids; recompute conservatively. + RecomputeGlobalMask(); + } + + public void MarkRejected(in ReducedStateKey key, byte qNibble) + { + if (!_stateRejectedMask.TryGetValue(key, out ushort mask)) + { + // Should not happen if RegisterState() is called, but keep it safe. + _stateRejectedMask[key] = 0; + SeenStates++; + mask = 0; + } + + ushort bit = (ushort)(1 << qNibble); + if ((mask & bit) != 0) return; // already counted for this reduced state + + mask |= bit; + _stateRejectedMask[key] = mask; + + _rejectCounts[qNibble]++; + + // If all seen states reject this qNibble, we can forbid it globally at this block. + if (_rejectCounts[qNibble] == SeenStates) + GlobalForbiddenMask |= bit; + } + + private void RecomputeGlobalMask() + { + ushort m = 0; + for (int d = 0; d < 16; d++) + if (_rejectCounts[d] == SeenStates && SeenStates > 0) + m |= (ushort)(1 << d); + GlobalForbiddenMask = m; + } + } + + private struct State + { + public ABitsWindow AWindow; // rolling a window + planned bits map + public int BitIndex; + public int Block; + public int Carry; + public byte KFlushRemaining; + public bool KTerminated; + public byte PostTermZeroQNibbles; + public DynamicBitset QBits; // q bits [0..maxQBits-1] + public QGenState16 QState; + + public static State CreateInitial(QGenState16 initQ, int maxQBits) + { + var s = new State + { + Block = initQ.NibbleIndex, + BitIndex = 0, + Carry = 0, + QState = initQ, + QBits = new DynamicBitset(maxQBits), + AWindow = new ABitsWindow(maxQBits), + KFlushRemaining = 0, + KTerminated = false, + PostTermZeroQNibbles = 0, + }; + + s.QBits.Set(0); + // a0=1 + s.AWindow.SetPlanned(0, 1); + return s; + } + + // Clone bitsets/windows to isolate branches. + // For performance, replace with copy-on-write/pooling. + + public State CloneShallow() => new() + + { + Block = Block, + BitIndex = BitIndex, + Carry = Carry, + QState = QState, + QBits = QBits.CloneShared(), + AWindow = AWindow.Clone(), + KFlushRemaining = KFlushRemaining, + KTerminated = KTerminated, + PostTermZeroQNibbles = PostTermZeroQNibbles, + }; + } + + // ----------------------------- Bit primitives ----------------------------- + + /// Dynamic bitset with fast enumeration of set bits up to a limit. + + public sealed class DynamicBitset + { + private sealed class BitsetData + { + public readonly int Bits; + public ulong[] Words; + + // Indices of words that are non-zero (monotone increasing in this scanner). + public int[] NzWords; + public int NzCount; + + public BitsetData(int bits) + { + Bits = bits; + Words = new ulong[(bits + 63) >> 6]; + NzWords = new int[16]; + NzCount = 0; + } + + public BitsetData(int bits, ulong[] words, int[] nzWords, int nzCount) + { + Bits = bits; + Words = words; + NzWords = nzWords; + NzCount = nzCount; + } + } + + private readonly BitsetData _data; + private bool _shared; // copy-on-write guard + + public DynamicBitset(int bits) + { + _data = new BitsetData(bits); + _shared = false; + } + + private DynamicBitset(BitsetData data, bool shared) + { + _data = data; + _shared = shared; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private void EnsureWritable() + { + if (!_shared) return; + + // Deep copy underlying arrays. + var w = new ulong[_data.Words.Length]; + Array.Copy(_data.Words, w, w.Length); + + var nz = new int[_data.NzWords.Length]; + Array.Copy(_data.NzWords, nz, nz.Length); + + var copied = new BitsetData(_data.Bits, w, nz, _data.NzCount); + + // Switch this instance to its private copy. + // Note: _data is readonly, so we rely on the fact that CloneShared creates a *new* DynamicBitset + // instance for children. This instance becomes non-shared by clearing _shared; its _data remains unique + // because no other instance references 'copied'. + // + // To preserve that with readonly _data, we keep _data as reference to mutable arrays inside BitsetData. + // Here we mutate the arrays in-place by replacing their references. + _data.Words = copied.Words; + _data.NzWords = copied.NzWords; + _data.NzCount = copied.NzCount; + + _shared = false; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private void AddNzWord(int wIndex) + { + int n = _data.NzCount; + if ((uint)n >= (uint)_data.NzWords.Length) + { + // grow + int newLen = _data.NzWords.Length * 2; + if (newLen < 16) newLen = 16; + var tmp = new int[newLen]; + Array.Copy(_data.NzWords, tmp, _data.NzWords.Length); + _data.NzWords = tmp; + } + _data.NzWords[n] = wIndex; + _data.NzCount = n + 1; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Set(int bit) + { + int bits = _data.Bits; + if ((uint)bit >= (uint)bits) return; + + int w = bit >> 6; + int b = bit & 63; + ulong mask = 1UL << b; + + ulong before = _data.Words[w]; + if ((before & mask) != 0) return; // already set; no mutation needed + + EnsureWritable(); + + before = _data.Words[w]; + ulong after = before | mask; + _data.Words[w] = after; + + if (before == 0 && after != 0) + AddNzWord(w); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public int Get(int bit) + { + int bits = _data.Bits; + if ((uint)bit >= (uint)bits) return 0; + ulong v = _data.Words[bit >> 6]; + return (int)((v >> (bit & 63)) & 1UL); + } + + /// Deep clone (copies backing arrays). + public DynamicBitset Clone() + { + var w = new ulong[_data.Words.Length]; + Array.Copy(_data.Words, w, w.Length); + + var nz = new int[_data.NzWords.Length]; + Array.Copy(_data.NzWords, nz, nz.Length); + + return new DynamicBitset(new BitsetData(_data.Bits, w, nz, _data.NzCount), shared: false); + } + + /// + /// Shared clone (copy-on-write). Children share the same arrays until one of them mutates. + /// + public DynamicBitset CloneShared() + { + // Mark both this and the clone as shared so either mutation triggers EnsureWritable. + _shared = true; + return new DynamicBitset(_data, shared: true); + } + + public ulong Signature64UpTo(int maxBitInclusive) + { + if (maxBitInclusive < 0) return 0UL; + int bits = _data.Bits; + int maxBit = Math.Min(maxBitInclusive, bits - 1); + int maxWord = maxBit >> 6; + int lastBitInWord = maxBit & 63; + + const ulong FnvOffset = 1469598103934665603UL; + const ulong FnvPrime = 1099511628211UL; + + ulong h = FnvOffset; + for (int i = 0; i <= maxWord; i++) + { + ulong v = _data.Words[i]; + if (i == maxWord) + { + ulong mask = lastBitInWord == 63 ? ulong.MaxValue : ((1UL << (lastBitInWord + 1)) - 1UL); + v &= mask; + } + h ^= v; + h *= FnvPrime; + } + return h; + } + + public ulong Signature64WithNibble(int startBit, byte qNibble, int maxQBits) + => Signature64WithNibbleUpTo(_data.Bits - 1, startBit, qNibble, maxQBits); + + public ulong Signature64WithNibbleUpTo(int maxBitInclusive, int startBit, byte qNibble, int maxQBits) + { + const ulong FnvOffset = 1469598103934665603UL; + const ulong FnvPrime = 1099511628211UL; + + int bits = _data.Bits; + int maxBit = Math.Min(Math.Max(maxBitInclusive, 0), bits - 1); + int maxWord = maxBit >> 6; + + int w0 = startBit >> 6; + int w1 = (startBit + 3) >> 6; + + ulong add0 = 0; + ulong add1 = 0; + + for (int b = 0; b < 4; b++) + { + int bitPos = startBit + b; + if (bitPos >= maxQBits || bitPos > maxBit) break; + + if (((qNibble >> b) & 1) == 0) continue; + + int wi = bitPos >> 6; + int bi = bitPos & 63; + if (wi == w0) add0 |= 1UL << bi; + else add1 |= 1UL << bi; + } + + // enforce q0=1 always + if (w0 == 0) add0 |= 1UL; + else if (w1 == 0) add1 |= 1UL; + + ulong h = FnvOffset; + for (int i = 0; i <= maxWord; i++) + { + ulong v = _data.Words[i]; + if (i == w0) v |= add0; + if (i == w1 && w1 != w0) v |= add1; + + if (i == maxWord) + { + int lastBit = maxBit & 63; + ulong mask = lastBit == 63 ? ulong.MaxValue : ((1UL << (lastBit + 1)) - 1UL); + v &= mask; + } + + h ^= v; + h *= FnvPrime; + } + return h; + } + + public BigInteger ToBigInteger() + { + // BigInteger expects little-endian byte array. + byte[] bytes = new byte[_data.Words.Length * 8 + 1]; + int idx = 0; + for (int i = 0; i < _data.Words.Length; i++) + { + ulong v = _data.Words[i]; + for (int b = 0; b < 8; b++) + { + bytes[idx++] = (byte)(v & 0xFF); + v >>= 8; + } + } + return new BigInteger(bytes); + } + + /// + /// Allocation-free enumeration of set bits up to a limit (inclusive), iterating only over non-zero words. + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public SetBitsUpToEnumerable EnumerateSetBitsUpTo(int maxBitInclusive) + => new(_data.Words, _data.NzWords, _data.NzCount, _data.Bits, maxBitInclusive); + + public readonly struct SetBitsUpToEnumerable + { + private readonly ulong[] _words; + private readonly int[] _nzWords; + private readonly int _nzCount; + private readonly int _bits; + private readonly int _maxBitInclusive; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public SetBitsUpToEnumerable(ulong[] words, int[] nzWords, int nzCount, int bits, int maxBitInclusive) + { + _words = words; + _nzWords = nzWords; + _nzCount = nzCount; + _bits = bits; + _maxBitInclusive = maxBitInclusive; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public Enumerator GetEnumerator() => new(_words, _nzWords, _nzCount, _bits, _maxBitInclusive); + + public struct Enumerator + { + private readonly ulong[] _words; + private readonly int[] _nzWords; + private readonly int _nzCount; + private readonly int _maxWord; + private readonly int _lastBitInWord; + + private int _i; // index in nzWords list + private int _w; // current word index + private ulong _word; + private bool _started; + + public int Current { get; private set; } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public Enumerator(ulong[] words, int[] nzWords, int nzCount, int bits, int maxBitInclusive) + { + _words = words; + _nzWords = nzWords; + _nzCount = nzCount; + + if (maxBitInclusive < 0 || bits <= 0) + { + _maxWord = -1; + _lastBitInWord = 0; + } + else + { + int maxBit = Math.Min(maxBitInclusive, bits - 1); + _maxWord = maxBit >> 6; + _lastBitInWord = maxBit & 63; + } + + _i = 0; + _w = 0; + _word = 0; + _started = false; + Current = 0; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool MoveNext() + { + if (_maxWord < 0) return false; + + if (!_started) + { + _started = true; + _i = 0; + if (!AdvanceToNextNonZeroWord()) + return false; + } + + while (true) + { + while (_word != 0) + { + int tz = BitOperations.TrailingZeroCount(_word); + Current = (_w << 6) + tz; + _word &= _word - 1; + return true; + } + + _i++; + if (!AdvanceToNextNonZeroWord()) + return false; + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private bool AdvanceToNextNonZeroWord() + { + while (_i < _nzCount) + { + int w = _nzWords[_i]; + if (w > _maxWord) return false; + _w = w; + _word = LoadWord(_w); + if (_word != 0) return true; + _i++; + } + return false; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private ulong LoadWord(int w) + { + ulong word = _words[w]; + if (w == _maxWord) + { + ulong mask = _lastBitInWord == 63 ? ulong.MaxValue : ((1UL << (_lastBitInWord + 1)) - 1UL); + word &= mask; + } + return word; + } + } + } + + public void ForEachSetBitUpTo(int maxBitInclusive, Action action) + { + if (maxBitInclusive < 0) return; + + int bits = _data.Bits; + int maxBit = Math.Min(maxBitInclusive, bits - 1); + int maxWord = maxBit >> 6; + int lastBitInWord = maxBit & 63; + + int nzCount = _data.NzCount; + int[] nzWords = _data.NzWords; + ulong[] words = _data.Words; + + for (int i = 0; i < nzCount; i++) + { + int w = nzWords[i]; + if (w > maxWord) break; + + ulong word = words[w]; + if (w == maxWord) + { + ulong mask = lastBitInWord == 63 ? ulong.MaxValue : ((1UL << (lastBitInWord + 1)) - 1UL); + word &= mask; + } + + while (word != 0) + { + int tz = BitOperations.TrailingZeroCount(word); + int bit = (w << 6) + tz; + action(bit); + word &= word - 1; + } + } + } + } + + + /// + /// Rolling window of a-bits of length maxQBits with support for "planned" bits. + /// We separate: + /// - planned bits: explicit assignments for specific a_i within the window + /// - pushed bits: bits that have been committed as BitIndex advances + /// + /// This allows nibble-level assignment while still streaming column-by-column. + /// + + public sealed class ABitsWindow + { + private readonly int _windowBits; + private readonly int _wordCount; + + // Committed bits are stored in a circular (bit-level) ring buffer. + // headBit points to the position holding the most recently committed bit (distance 0). + private readonly ulong[] _ring; + private int _headBit; + private int _currentJ; // last committed index j (monotone increasing) + + // Planned a_i values (sparse) for not-yet-committed indices. + private readonly Dictionary _planned; + + public int CurrentJ => _currentJ; + + public ABitsWindow(int windowBits) + { + _windowBits = windowBits; + _wordCount = (windowBits + 63) >> 6; + _ring = new ulong[_wordCount]; + _headBit = -1; // "empty" + _currentJ = -1; + _planned = new Dictionary(capacity: 1024); + } + + private ABitsWindow(int windowBits, ulong[] ring, int headBit, int currentJ, Dictionary planned) + { + _windowBits = windowBits; + _wordCount = (windowBits + 63) >> 6; + _ring = ring; + _headBit = headBit; + _currentJ = currentJ; + _planned = planned; + } + + public ABitsWindow Clone() + { + var ring2 = new ulong[_ring.Length]; + Array.Copy(_ring, ring2, ring2.Length); + var planned2 = new Dictionary(_planned); + return new ABitsWindow(_windowBits, ring2, _headBit, _currentJ, planned2); + } + + /// + /// A compact signature of the current a-window state (committed ring + planned bits). + /// Used to keep AllowedMaskCache semantically correct: allowed-a depends on previously committed a-bits. + /// + public ulong Signature64() + { + const ulong FnvOffset = 1469598103934665603UL; + const ulong FnvPrime = 1099511628211UL; + + ulong h = FnvOffset; + + h ^= (ulong)(uint)_windowBits; h *= FnvPrime; + h ^= (ulong)(uint)_headBit; h *= FnvPrime; + h ^= (ulong)(uint)_currentJ; h *= FnvPrime; + + for (int i = 0; i < _ring.Length; i++) + { + h ^= _ring[i]; + h *= FnvPrime; + } + + // Planned bits: order-independent combine + ulong ph = 0; + foreach (var kv in _planned) + { + // mix index and bit (bit is 0/1) + ulong v = ((ulong)(uint)kv.Key << 1) ^ (uint)kv.Value; + ph ^= v + 0x9e3779b97f4a7c15UL + (ph << 6) + (ph >> 2); + } + h ^= ph; + h *= FnvPrime; + + return h; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void SetPlanned(int aIndex, int bit) + { + if (bit != 0) _planned[aIndex] = 1; + else _planned.Remove(aIndex); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void PushBit(int j, int bit) + { + // Commit bit for column j (a_j). + _planned.Remove(j); + _currentJ = j; + + // Advance ring head and write the new bit. + _headBit = (_headBit + 1) % _windowBits; + + int word = _headBit >> 6; + int b = _headBit & 63; + + _ring[word] &= ~(1UL << b); + if (bit != 0) _ring[word] |= 1UL << b; + + // planned map should not grow without bound; remove entries older than window. + int minKeep = j - _windowBits; + if (_planned.Count > 0) + { + var toRemove = new List(); + foreach (var kv in _planned) + if (kv.Key <= minKeep) toRemove.Add(kv.Key); + foreach (int k in toRemove) _planned.Remove(k); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public int GetBitRelative(int aIndex) + { + // planned must have priority even if aIndex > _currentJ + if (_planned.TryGetValue(aIndex, out int plannedBit)) + return plannedBit; + + int d = _currentJ - aIndex; + if (d < 0 || d >= _windowBits) return 0; + if (_headBit < 0) return 0; // empty + + int pos = _headBit - d; + pos %= _windowBits; + if (pos < 0) pos += _windowBits; + + int word = pos >> 6; + int bit = pos & 63; + return (int)((_ring[word] >> bit) & 1UL); + } + + /// + /// Try to read an a-bit if it is currently KNOWN (planned or already committed within the window). + /// Returns false when the bit is outside the committed window and not planned yet. + /// + /// IMPORTANT: Callers that implement contradiction pruning must treat "unknown" as Unknown (may be 0 or 1), + /// never as 0, to avoid false negatives. + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool TryGetKnownBit(int aIndex, out int bit) + { + // Planned has priority and is always a known assignment. + if (_planned.TryGetValue(aIndex, out int plannedBit)) + { + bit = plannedBit; + return true; + } + + // Not yet committed. + if (aIndex > _currentJ) + { + bit = 0; + return false; + } + + int d = _currentJ - aIndex; + if (d < 0 || d >= _windowBits || _headBit < 0) + { + bit = 0; + return false; + } + + int pos = _headBit - d; + pos %= _windowBits; + if (pos < 0) pos += _windowBits; + + int word = pos >> 6; + int b = pos & 63; + bit = (int)((_ring[word] >> b) & 1UL); + return true; + } + + public bool TryGetPlanned(int aIndex, out int bit) + => _planned.TryGetValue(aIndex, out bit); + } + + // ---------------------------------------------------------------------------- + // q nibble generator for q = (2p)*k + 1 in base 16, without materializing k/q. + // ---------------------------------------------------------------------------- + + public readonly struct QGenState16 + { + public readonly int NibbleIndex; // <-- DODAJ + public readonly uint Carry16; + public readonly ulong KWindowPacked; + + public QGenState16(int nibbleIndex, uint carry16, ulong kWindowPacked) + { + NibbleIndex = nibbleIndex; + Carry16 = carry16; + KWindowPacked = kWindowPacked; + } + } + + public sealed class QNibbleGenerator16 + { + private readonly byte[] _mDigits; // nibbles of m=2p in base 16, LSB-first + private readonly int _mLen; + private readonly int _kWinLen; // mLen-1 + private readonly ulong _kWinMask; + + public QNibbleGenerator16(ulong p) + { + ulong m = checked(2UL * p); + _mDigits = ToBase16DigitsLSB(m); + _mLen = _mDigits.Length; + _kWinLen = Math.Max(0, _mLen - 1); + _kWinMask = _kWinLen == 0 ? 0UL : ((1UL << (4 * _kWinLen)) - 1UL); + } + + public static QGenState16 InitialState => new(0, 0, 0); + public int KWinLen => _kWinLen; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void ComputeNext(QGenState16 state, byte kNibble, out byte qNibble, out QGenState16 nextState) + { + ulong sum = state.Carry16; + int nibbleIndex = state.NibbleIndex; + + // +1 only contributes to the very first nibble, once + if (nibbleIndex == 0) + sum += 1UL; + + // j=0 term + sum += (ulong)_mDigits[0] * kNibble; + + int maxJ = Math.Min(nibbleIndex, _mLen - 1); + for (int j = 1; j <= maxJ; j++) + { + byte mj = _mDigits[j]; + byte kPrev = GetKFromWindow(state.KWindowPacked, j - 1, _kWinLen); + sum += (ulong)mj * kPrev; + } + + qNibble = (byte)(sum & 0xFUL); + uint nextCarry = (uint)(sum >> 4); + + ulong nextWin = state.KWindowPacked; + if (_kWinLen > 0) + nextWin = ((nextWin << 4) | kNibble) & _kWinMask; + + nextState = new QGenState16( + nibbleIndex + 1, + nextCarry, + nextWin + ); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void ComputeAllNext( + QGenState16 state, + Span qNibByKn, // len=16 + Span nextByKn, // len=16 + Span countsByQn // len=16 + ) + { + ulong baseSum = state.Carry16; + int nibbleIndex = state.NibbleIndex; + + // +1 only contributes to the very first nibble, once + if (nibbleIndex == 0) + baseSum += 1UL; + + // Σ_{j=1..maxJ} mj * kPrev (part independent of current kn) + int maxJ = Math.Min(nibbleIndex, _mLen - 1); + for (int j = 1; j <= maxJ; j++) + { + byte mj = _mDigits[j]; + byte kPrev = GetKFromWindow(state.KWindowPacked, j - 1, _kWinLen); + baseSum += (ulong)mj * kPrev; + } + + ulong m0 = _mDigits[0]; + + for (byte kn = 0; kn <= 15; kn++) + { + ulong sum = baseSum + m0 * kn; + + byte qn = (byte)(sum & 0xFUL); + uint nextCarry = (uint)(sum >> 4); + + ulong nextWin = state.KWindowPacked; + if (_kWinLen > 0) + nextWin = ((nextWin << 4) | kn) & _kWinMask; + + qNibByKn[kn] = qn; + nextByKn[kn] = new QGenState16(nibbleIndex + 1, nextCarry, nextWin); + countsByQn[qn]++; + } + } + + + private static byte[] ToBase16DigitsLSB(ulong x) + { + if (x == 0) return [0]; + Span tmp = stackalloc byte[16]; + int len = 0; + while (x != 0) + { + tmp[len++] = (byte)(x & 0xF); + x >>= 4; + } + var r = new byte[len]; + for (int i = 0; i < len; i++) r[i] = tmp[i]; + return r; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte GetKFromWindow(ulong packed, int idxFromNewest, int winLen) + { + // newest nibble is stored at LSB (shift 0) + // idxFromNewest=0 => most recently appended nibble + if ((uint)idxFromNewest >= (uint)winLen) return 0; + int shift = idxFromNewest * 4; + return (byte)((packed >> shift) & 0xF); + } + + public static void SelfTest(ulong p) + { + var gen = new QNibbleGenerator16(p); + var st = InitialState; + + // k = 4383 = 0x111F => nibble LSB-first: F,1,1,1,0,0,... + byte[] kNibs = { 0xF, 0x1, 0x1, 0x1, 0x0, 0x0, 0x0, 0x0 }; + + BigInteger q = BigInteger.Zero; + for (int i = 0; i < 12; i++) + { + byte kn = i < kNibs.Length ? kNibs[i] : (byte)0; + gen.ComputeNext(st, kn, out byte qn, out var ns); + q |= (new BigInteger(qn) << (4 * i)); + st = ns; + } + + Console.WriteLine($"[SELFTEST] p={p} q(hex)={q.ToString("X")}"); + Console.WriteLine($"[SELFTEST] p={p} q(dec)={q}"); + } + } +} diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs index 17d3f4c..1c41004 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs @@ -1,9 +1,15 @@ +// #define DivisorSet_BitContradiction +#define DivisorSet_BitTree + using System; using System.Numerics; +using System.IO; +using System.Text; using System.Runtime.CompilerServices; using System.Globalization; using System.Collections.Generic; using System.Collections.Concurrent; +using System.Buffers; using PerfectNumbers.Core.Gpu; using PerfectNumbers.Core.Gpu.Accelerators; using PerfectNumbers.Core.ByDivisor; @@ -20,6 +26,7 @@ public enum DivisorSet Additive, TopDown, BitContradiction, + BitTree, } [EnumDependentTemplate(typeof(DivisorSet))] @@ -48,11 +55,14 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN private BigInteger _cheapResumeK; private GpuUInt128WorkSet _divisorScanGpuWorkSet; -private static readonly ConcurrentDictionary _checkedPow2MinusOne = new(); -private ulong _currentPrime; + private static readonly ConcurrentDictionary _checkedPow2MinusOne = new(); + private ulong _currentPrime; #if DivisorSet_BitContradiction private static Persistence.BitContradictionStateRepository? _bitContradictionStateRepository; #endif +#if DivisorSet_BitTree + private static Persistence.BitContradictionStateRepository? _bitTreeStateRepository; +#endif #if DivisorSet_TopDown private BigInteger _topDownCursor; #endif @@ -335,6 +345,9 @@ private ulong _currentPrime; #endif #if DivisorSet_BitContradiction _bitContradictionStateRepository = EnvironmentConfiguration.BitContradictionStateRepository; +#endif +#if DivisorSet_BitTree + _bitTreeStateRepository = EnvironmentConfiguration.BitContradictionStateRepository; #endif } @@ -344,11 +357,11 @@ private ulong _currentPrime; _currentPrime = prime; BigInteger allowedMax = MersenneNumberDivisorByDivisorTester.ComputeAllowedMaxDivisorBig(prime, DivisorLimit); - if (TryFindPowerOfTwoMinusOneDivisor(prime, allowedMax, out divisor)) - { - divisorsExhausted = true; - return false; - } + // if (TryFindPowerOfTwoMinusOneDivisor(prime, allowedMax, out divisor)) + // { + // divisorsExhausted = true; + // return false; + // } // The CPU by-divisor run always hands us primes with enormous divisor limits, so the fallback below never executes. // if (allowedMax < 3UL) @@ -413,9 +426,9 @@ private ulong _currentPrime; out BigInteger foundDivisor) { // TODO: Is minK ever < 1? Can this be optimized - normalizedMinK is always >= 1 after this - BigInteger normalizedMinK = minK >= BigInteger.One ? minK : BigInteger.One, - step = ((BigInteger)prime) << 1, - firstDivisor = (step * normalizedMinK) + BigInteger.One; + BigInteger normalizedMinK = minK >= BigInteger.One ? minK : BigInteger.One; + BigInteger step = ((BigInteger)prime) << 1; + BigInteger firstDivisor = (step * normalizedMinK) + BigInteger.One; bool primeLastIsSeven = (prime & 3UL) == 3UL; ushort primeDecimalMask = DivisorGenerator.GetDecimalMask(primeLastIsSeven); @@ -450,6 +463,17 @@ private ulong _currentPrime; out foundDivisor); return compositeBitContradiction; #endif +#if DivisorSet_BitTree + bool compositeBitTree = CheckDivisorsBitTree( + prime, + primeDecimalMask, + allowedMax, + normalizedMinK, + step, + out processedAll, + out foundDivisor); + return compositeBitTree; +#endif #if DivisorSet_TopDown bool compositeTopDown = CheckDivisorsTopDown( prime, @@ -708,6 +732,941 @@ private ulong _currentPrime; return false; } #endif +#if DivisorSet_BitTree + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private bool CheckDivisorsBitTree( + ulong prime, + ushort primeDecimalMask, + BigInteger allowedMax, + BigInteger minK, + BigInteger step, + out bool processedAll, + out BigInteger foundDivisor) + { + processedAll = true; + foundDivisor = BigInteger.Zero; + + // Default: frontier/windowed scan over (q,a) space. + MersenneCombinedDivisorScannerFullStreamExtended.QNibbleGenerator16.SelfTest(138000001); + MersenneCombinedDivisorScannerFullStreamExtended.ScanResult scanResult = MersenneCombinedDivisorScannerFullStreamExtended.TryFindDivisorByStreamingScan(prime, 128, false, out BigInteger branchDivisor); + // if (TryExactBitTreeCheck(prime, allowedMax, minK, step, out BigInteger branchDivisor)) + if (scanResult == MersenneCombinedDivisorScannerFullStreamExtended.ScanResult.FoundDivisor) + { + foundDivisor = branchDivisor; + return true; + } + else if (scanResult == MersenneCombinedDivisorScannerFullStreamExtended.ScanResult.Candidate) + { + processedAll = false; + } + + return false; + } + + private bool TryBitTreeBranchScan(ulong prime, BigInteger allowedMax, out BigInteger divisor) + { + divisor = BigInteger.Zero; + + if (prime == 0UL) + { + return false; + } + + int pBits = checked((int)prime); + var stack = new Stack<(int Index, int Carry, BigInteger Q, BigInteger A)>(); + // q0 = 1, a0 = 1 (to make r0 = 1) + stack.Push((Index: 1, Carry: 0, Q: BigInteger.One, A: BigInteger.One)); + + while (stack.Count > 0) + { + var state = stack.Pop(); + int j = state.Index; + int carry = state.Carry; + BigInteger qVal = state.Q; + BigInteger aVal = state.A; + + if (j >= pBits) + { + if (carry == 0) + { + if (allowedMax.IsZero || qVal <= allowedMax) + { + divisor = qVal; + return true; + } + } + continue; + } + + for (int qBit = 0; qBit <= 1; qBit++) + { + BigInteger qCandidate = qVal; + if (qBit == 1) + { + qCandidate |= (BigInteger.One << j); + } + + // Sum over assigned bits up to column j (a_i and q_{j-i} known up to j). + int sum = carry; + for (int i = 0; i <= j; i++) + { + bool aBit = ((aVal >> i) & BigInteger.One) == BigInteger.One; + bool qBitVal = ((qCandidate >> (j - i)) & BigInteger.One) == BigInteger.One; + if (aBit && qBitVal) + { + sum += 1; + } + } + + // Required result bit = 1 for columns < p + if ((sum & 1) != 1) + { + continue; + } + + // Set a_j to satisfy parity (already implied by sum&1==1) + BigInteger aNext = aVal; + if (((aVal >> j) & BigInteger.One) == BigInteger.Zero) + { + // keep as 0; sum already odd + } + + int carryNext = (sum - 1) >> 1; + + if (!allowedMax.IsZero && qCandidate > allowedMax) + { + continue; + } + + stack.Push((j + 1, carryNext, qCandidate, aNext)); + } + } + + return false; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static bool TryBitTreeFrontierScan(ulong prime, in BigInteger allowedMax, out BigInteger divisor) + { + divisor = BigInteger.Zero; + + const int WindowSize = 200_000_000; + const int MemoryFrontierLimit = 1_000_000; + const int DiskBatchMax = 100_000; + + int pBits; + try + { + pBits = checked((int)prime); + } + catch (OverflowException) + { + Console.WriteLine($"[bittree-frontier] exponent too large to check p={prime}"); + return false; + } + + int maxQBits = allowedMax.IsZero ? pBits : GetBitLengthPortable(allowedMax); + if (maxQBits <= 0) + { + return false; + } + + int hardLimit = pBits + maxQBits + 2; + Span initialOnes = stackalloc int[1]; + initialOnes[0] = 0; + + var frontier = new Queue(capacity: MemoryFrontierLimit); + BigInteger initialQ = BigInteger.One; + BigInteger initialA = BigInteger.One; + frontier.Enqueue(new BitTreeFrontierState( + column: 1, + carry: 0, + baseIndex: 0, + qSuffix: initialQ, + aSuffix: initialA, + maxQIndex: 0, + lastNonZeroA: 0, + qOnes: initialOnes.ToArray())); + using BitTreeFrontierRepository frontierRepo = BitTreeFrontierRepository.Open(Path.Combine("Checks"), "bittree.frontier.faster"); + + Console.WriteLine($"[bittree-frontier] p={prime} active_paths=1 column=1"); + + long prunedPaths = 0; + long pruneLogThreshold = 1_000_000; + const int ColumnLogStep = 1_000_000; + int lastLoggedColumn = 1; + void TrackPrune(int column) + { + prunedPaths++; + if (prunedPaths >= pruneLogThreshold) + { + Console.WriteLine($"[bittree-frontier] p={prime} pruned_paths={prunedPaths} column={column}"); + pruneLogThreshold <<= 1; + } + if (column - lastLoggedColumn >= ColumnLogStep) + { + lastLoggedColumn = column; + Console.WriteLine($"[bittree-frontier] p={prime} exploring_column={column}"); + } + } + + int maxFrontier = 1; + int printCount = 0; + + static byte[] SerializeState(in BitTreeFrontierState state) + { + using var ms = new MemoryStream(); + using var bw = new BinaryWriter(ms, Encoding.UTF8, leaveOpen: true); + bw.Write(state.Column); + bw.Write(state.Carry); + bw.Write(state.BaseIndex); + bw.Write(state.MaxQIndex); + bw.Write(state.LastNonZeroA); + bw.Write(state.QOnes.Length); + for (int i = 0; i < state.QOnes.Length; i++) + { + bw.Write(state.QOnes[i]); + } + byte[] qBytes = state.QSuffix.ToByteArray(isUnsigned: true, isBigEndian: false); + bw.Write(qBytes.Length); + bw.Write(qBytes); + byte[] aBytes = state.ASuffix.ToByteArray(isUnsigned: true, isBigEndian: false); + bw.Write(aBytes.Length); + bw.Write(aBytes); + bw.Flush(); + return ms.ToArray(); + } + + static BitTreeFrontierState DeserializeState(BinaryReader br) + { + int column = br.ReadInt32(); + int carry = br.ReadInt32(); + int baseIndex = br.ReadInt32(); + int maxQIndex = br.ReadInt32(); + int lastNonZeroA = br.ReadInt32(); + int onesLength = br.ReadInt32(); + int[] qOnes = new int[onesLength]; + for (int i = 0; i < onesLength; i++) + { + qOnes[i] = br.ReadInt32(); + } + int qLen = br.ReadInt32(); + byte[] qBytes = br.ReadBytes(qLen); + int aLen = br.ReadInt32(); + byte[] aBytes = br.ReadBytes(aLen); + BigInteger qSuffix = new BigInteger(qBytes, isUnsigned: true, isBigEndian: false); + BigInteger aSuffix = new BigInteger(aBytes, isUnsigned: true, isBigEndian: false); + return new BitTreeFrontierState(column, carry, baseIndex, qSuffix, aSuffix, maxQIndex, lastNonZeroA, qOnes); + } + + static byte[] SerializeBatch(IEnumerable states) + { + using var ms = new MemoryStream(); + using var bw = new BinaryWriter(ms, Encoding.UTF8, leaveOpen: true); + foreach (var s in states) + { + byte[] payload = SerializeState(in s); + bw.Write(payload.Length); + bw.Write(payload); + } + bw.Flush(); + return ms.ToArray(); + } + + static List DeserializeBatch(byte[] batch) + { + var list = new List(); + using var ms = new MemoryStream(batch); + using var br = new BinaryReader(ms, Encoding.UTF8, leaveOpen: true); + while (ms.Position < ms.Length) + { + int len = br.ReadInt32(); + byte[] payload = br.ReadBytes(len); + using var itemStream = new MemoryStream(payload); + using var itemReader = new BinaryReader(itemStream, Encoding.UTF8, leaveOpen: true); + list.Add(DeserializeState(itemReader)); + } + return list; + } + + void SpillFrontierIfNeeded() + { + if (frontier.Count <= MemoryFrontierLimit) + { + return; + } + + int toSpill = frontier.Count - MemoryFrontierLimit; + if (toSpill <= 0 || toSpill < DiskBatchMax) + { + return; + } + + var batchStates = new List(toSpill); + for (int i = 0; i < toSpill && frontier.Count > 0; i++) + { + batchStates.Add(frontier.Dequeue()); + } + + byte[] batchBytes = SerializeBatch(batchStates); + frontierRepo.AppendBatch(batchBytes); + // Console.WriteLine($"[bittree-frontier] spilled {batchStates.Count} states to disk; in-memory={frontier.Count}"); + } + + bool ReloadFrontierFromRepo() + { + if (frontier.Count > 0) + { + return true; + } + + if (frontierRepo.TryDequeueBatch(1, out List batches) && batches.Count > 0) + { + foreach (byte[] b in batches) + { + foreach (var s in DeserializeBatch(b)) + { + frontier.Enqueue(s); + } + } + Console.WriteLine($"[bittree-frontier] reloaded {frontier.Count} states from disk."); + return frontier.Count > 0; + } + + return false; + } + + + while (frontier.Count > 0) + { + BitTreeFrontierState state = frontier.Dequeue(); + if (state.Column > hardLimit) + { + TrackPrune(state.Column); + continue; + } + + state = SlideWindow(state, WindowSize); + + for (int qBitChoice = 0; qBitChoice <= 1; qBitChoice++) + { + BigInteger qCandidate = state.QSuffix; + int maxQIndex = state.MaxQIndex; + int[] qOnes = state.QOnes; + + if (qBitChoice == 1) + { + int bitIndex = state.Column - state.BaseIndex; + qCandidate |= BigInteger.One << bitIndex; + maxQIndex = Math.Max(maxQIndex, state.Column); + + int[] extended = new int[qOnes.Length + 1]; + Array.Copy(qOnes, extended, qOnes.Length); + extended[qOnes.Length] = state.Column; + qOnes = extended; + } + + BigInteger absoluteQ = qCandidate << state.BaseIndex; + if (!allowedMax.IsZero && absoluteQ > allowedMax) + { + TrackPrune(state.Column); + continue; + } + + int qBitLength = maxQIndex + 1; + if (qBitLength > maxQBits) + { + TrackPrune(state.Column); + continue; + } + + int maxAIndexAllowed = pBits - qBitLength; + if (maxAIndexAllowed < 0) + { + TrackPrune(state.Column); + continue; + } + + int sum = state.Carry; + for (int i = 0; i < qOnes.Length; i++) + { + int offset = qOnes[i]; + if (offset == 0 || offset > state.Column) + { + continue; + } + + int aIndex = state.Column - offset; + int aBitIndex = aIndex - state.BaseIndex; + if (aBitIndex >= 0 && aBitIndex < WindowSize && ((state.ASuffix >> aBitIndex) & BigInteger.One) == BigInteger.One) + { + sum += 1; + } + } + + int requiredParity = state.Column < pBits ? 1 : 0; + int parity = sum & 1; + int aBit = parity == requiredParity ? 0 : 1; + + BigInteger aCandidate = state.ASuffix; + int aBitIndexCurrent = state.Column - state.BaseIndex; + if (aBit == 1) + { + if (aBitIndexCurrent >= WindowSize) + { + TrackPrune(state.Column); + continue; + } + aCandidate |= BigInteger.One << aBitIndexCurrent; + } + + int lastNonZeroA = aBit == 1 ? state.Column : state.LastNonZeroA; + if (lastNonZeroA > maxAIndexAllowed) + { + TrackPrune(state.Column); + continue; + } + + int columnTotal = sum + aBit; + int carry = columnTotal >> 1; + + bool doneWithOnes = lastNonZeroA < 0 || (state.Column + 1 - lastNonZeroA) > maxQIndex; + if (state.Column + 1 >= pBits && carry == 0 && doneWithOnes) + { + divisor = absoluteQ; + return true; + } + + if (state.Column + 1 > hardLimit) + { + TrackPrune(state.Column); + continue; + } + + int[] nextOnes = qOnes; + if (qBitChoice == 1) + { + nextOnes = new int[qOnes.Length + 1]; + Array.Copy(qOnes, nextOnes, qOnes.Length); + nextOnes[^1] = state.Column; + } + + frontier.Enqueue(new BitTreeFrontierState( + column: state.Column + 1, + carry: carry, + baseIndex: state.BaseIndex, + qSuffix: qCandidate, + aSuffix: aCandidate, + maxQIndex: maxQIndex, + lastNonZeroA: lastNonZeroA, + qOnes: nextOnes)); + + SpillFrontierIfNeeded(); + + if (frontier.Count > maxFrontier) + { + maxFrontier = frontier.Count; + printCount++; + if (printCount == 10_000_000) + { + printCount = 0; + Console.WriteLine($"[bittree-frontier] p={prime} active_paths={maxFrontier} column={state.Column + 1}"); + } + } + } + + if (frontier.Count == 0 && !ReloadFrontierFromRepo()) + { + break; + } + } + + return false; + } + + private readonly struct BitTreeFrontierState + { + public BitTreeFrontierState( + int column, + int carry, + int baseIndex, + BigInteger qSuffix, + BigInteger aSuffix, + int maxQIndex, + int lastNonZeroA, + int[] qOnes) + { + Column = column; + Carry = carry; + BaseIndex = baseIndex; + QSuffix = qSuffix; + ASuffix = aSuffix; + MaxQIndex = maxQIndex; + LastNonZeroA = lastNonZeroA; + QOnes = qOnes; + } + + public readonly int Column; + public readonly int Carry; + public readonly int BaseIndex; + public readonly BigInteger QSuffix; + public readonly BigInteger ASuffix; + public readonly int MaxQIndex; + public readonly int LastNonZeroA; + public readonly int[] QOnes; + } + + private readonly struct BitTreeBacktrackFrame + { + public BitTreeBacktrackFrame(BitTreeFrontierState state, int nextChoice) + { + State = state; + NextChoice = nextChoice; + } + + public BitTreeFrontierState State { get; } + public int NextChoice { get; } + + public BitTreeBacktrackFrame WithNextChoice(int nextChoice) + { + return new BitTreeBacktrackFrame(State, nextChoice); + } +} + + + private static BitTreeFrontierState SlideWindow(in BitTreeFrontierState state, int windowSize) + { + int span = state.Column - state.BaseIndex + 1; + if (span < windowSize) + { + return state; + } + + int shift = span - (windowSize - 1); + int newBase = state.BaseIndex + shift; + + BigInteger qSuffix = state.QSuffix >> shift; + BigInteger aSuffix = state.ASuffix >> shift; + + List filtered = null; + for (int i = 0; i < state.QOnes.Length; i++) + { + int pos = state.QOnes[i]; + if (pos >= newBase) + { + filtered ??= new List(state.QOnes.Length - i); + filtered.Add(pos); + } + } + + int[] qOnes = filtered == null ? Array.Empty() : filtered.ToArray(); + + int lastNonZeroA = state.LastNonZeroA; + if (lastNonZeroA < newBase) + { + lastNonZeroA = -1; + } + + return new BitTreeFrontierState( + column: state.Column, + carry: state.Carry, + baseIndex: newBase, + qSuffix: qSuffix, + aSuffix: aSuffix, + maxQIndex: state.MaxQIndex, + lastNonZeroA: lastNonZeroA, + qOnes: qOnes); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private bool TryExactBitTreeCheck( + ulong prime, + BigInteger allowedMax, + BigInteger currentK, + BigInteger step, + out BigInteger foundDivisor) + { + foundDivisor = BigInteger.Zero; + + bool hasLimit = !allowedMax.IsZero; + BigInteger maxK = hasLimit ? (allowedMax - BigInteger.One) / step : BigInteger.Zero; + ulong computedCycle; + + bool decided; + bool divides; + BigInteger divisor = (step * currentK) + BigInteger.One; + while (true) + { +// if (divisor <= ulong.MaxValue) +// { +// ulong divisor64 = (ulong)divisor; +// MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor64); +// // We're reusing divides variable to limit registry pressure as primeOrderFailed flag. +// #if DEVICE_GPU +// decided = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentGpu( +// divisor64, +// Accelerator, +// prime, +// divisorData, +// out computedCycle, +// out divides); +// #elif DEVICE_HYBRID +// decided = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentHybrid( +// divisor64, +// Accelerator, +// prime, +// divisorData, +// out computedCycle, +// out divides); +// #else +// decided = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentCpu( +// divisor64, +// prime, +// divisorData, +// out computedCycle, +// out divides); +// #endif + +// if (!decided) +// { +// #if DEVICE_GPU +// computedCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthGpu( +// divisor64, +// Accelerator, +// divisorData, +// skipPrimeOrderHeuristic: divides); +// #elif DEVICE_HYBRID +// computedCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthHybrid( +// divisor64, +// Accelerator, +// divisorData, +// skipPrimeOrderHeuristic: divides); +// #else +// computedCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthCpu( +// divisor64, +// divisorData, +// skipPrimeOrderHeuristic: divides); +// #endif +// } + +// if (computedCycle > 0UL) +// { +// if ((prime % computedCycle) != 0UL) +// { +// goto MOVE_NEXT; +// } +// } +// else +// { +// computedCycle = prime; +// } +// } +// else +// { + computedCycle = prime; + // } + + decided = TryExactBitTreeCheck(computedCycle, divisor, out divides); + if (!decided) + { + throw new InvalidOperationException($"BitTree check did not decide for p={prime}, divisor={divisor} (k={currentK})."); + } + + if (divides) + { + foundDivisor = divisor; + return true; + } + + RecordState(currentK); + + MOVE_NEXT: + if (hasLimit && currentK >= maxK) + { + break; + } + + currentK += BigInteger.One; + divisor += step; + } + + return false; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private bool TryExactBitTreeCheck(ulong exponent, BigInteger divisor, out bool divides) + { + int qBitLength = GetBitLengthPortable(divisor); + int aBitLength = GetBitLengthPortable(exponent); + + ArrayPool intPool = ThreadStaticPools.IntPool; + int pBits = checked((int)exponent); + + int[] oneOffsetsArr = intPool.Rent(qBitLength); + int oneCount = 0; + int i; + for (i = 0; i < qBitLength; i++) + { + if (((divisor >> i) & BigInteger.One) == BigInteger.One) + { + oneOffsetsArr[oneCount++] = i; + } + } + + // This condition will never trigger on EvenPerfectBitScanner execution paths. We always provide valid q != 0. + // if (oneCount == 0) + // { + // divides = false; + // intPool.Return(oneOffsetsArr); + // return true; + // } + + // i's in oneOffsetsArr are built in ascending order by i. We don't need to sort it. + // Array.Sort(oneOffsetsArr, 0, oneCount); + ReadOnlySpan oneOffsets = oneOffsetsArr.AsSpan(0, oneCount); + int maxOffset = oneOffsets[^1]; + + int windowSize = qBitLength; + byte[] aWindowArray = ArrayPool.Shared.Rent(windowSize); + Span aWindow = aWindowArray.AsSpan(0, windowSize); + aWindow.Clear(); + + long lastNonZeroA = -1; + int carry = 0; + int column = 0; + + while (true) + { + bool targetOne = column < pBits; + int sum = carry; + + for (i = 0; i < oneCount; i++) + { + int offset = oneOffsets[i]; + if (offset == 0) + { + continue; // a_bit for current column, handled below + } + + if (offset > column) + { + break; + } + + int aIndex = column - offset; + byte aBit = aWindow[aIndex % windowSize]; + sum += aBit; + } + + int requiredParity = targetOne ? 1 : 0; + int parity = sum & 1; + byte aCurrent = (byte)(parity == requiredParity ? 0 : 1); + int columnTotal = sum + aCurrent; + + if ((columnTotal & 1) != requiredParity) + { + ArrayPool.Shared.Return(aWindowArray); + + divides = false; + intPool.Return(oneOffsetsArr); + return true; + } + + aWindow[column % windowSize] = aCurrent; + if (aCurrent == 1) + { + lastNonZeroA = column; + } + + carry = columnTotal >> 1; + column++; + + bool doneWithOnes = lastNonZeroA < 0 || column - lastNonZeroA > maxOffset; + if (column >= pBits && carry == 0 && doneWithOnes) + { + break; + } + + if (column > pBits + maxOffset + 2 && carry == 0 && doneWithOnes) + { + break; + } + + if (column > pBits + maxOffset + aBitLength + 2) + { + ArrayPool.Shared.Return(aWindowArray); + + divides = false; + intPool.Return(oneOffsetsArr); + return true; + } + } + + ArrayPool.Shared.Return(aWindowArray); + + divides = carry == 0; + Console.WriteLine($"[bittree] Prime={exponent}, Divisor={divisor} Decided=True Divides={divides}"); + if (!divides) + { + string message = $"[bittree] p={exponent} ruled out divisor={divisor} (parity/carry exhaustion)"; + _bitTreeStateRepository?.Upsert(exponent, message); + } + + intPool.Return(oneOffsetsArr); + return true; + } + + // [MethodImpl(MethodImplOptions.AggressiveInlining)] + // private bool TryExactBitTreeCheck(ulong prime, BigInteger divisor, ulong effectiveExponent, out bool divides) + // { + // divides = false; + + // // These conditions below will never be met on EvenPerfectBitScanner execution paths. + // // if (prime == 0 || divisor <= BigInteger.One) + // // { + // // Console.WriteLine($"[bittree] invalid inputs {prime} / {divisor}"); + // // return false; + // // } + + // // if ((divisor & BigInteger.One) == BigInteger.Zero) + // // { + // // divides = false; + // // return true; + // // } + + // int qBitLength = GetBitLengthPortable(divisor); + // if (prime <= (ulong)qBitLength) + // { + // divides = false; + // return true; + // } + + // ulong exponent = effectiveExponent > 0UL ? effectiveExponent : prime; + // int pBits; + // try + // { + // pBits = checked((int)exponent); + // } + // catch (OverflowException) + // { + // Console.WriteLine($"[bittree] exponent too large to check p={exponent}"); + // return false; + // } + + // int[] oneOffsetsArr = new int[qBitLength]; + // int oneCount = 0; + // for (int i = 0; i < qBitLength; i++) + // { + // if (((divisor >> i) & BigInteger.One) == BigInteger.One) + // { + // oneOffsetsArr[oneCount++] = i; + // } + // } + + // if (oneCount == 0) + // { + // divides = false; + // return true; + // } + + // Array.Sort(oneOffsetsArr, 0, oneCount); + // ReadOnlySpan oneOffsets = oneOffsetsArr.AsSpan(0, oneCount); + // int maxOffset = oneOffsets[^1]; + + // int windowSize = qBitLength; + // bool rented = windowSize <= 1024; + // byte[] aWindowArray = rented ? ArrayPool.Shared.Rent(1024) : new byte[windowSize]; + // Span aWindow = aWindowArray.AsSpan(0, windowSize); + // aWindow.Clear(); + + // long lastNonZeroA = -1; + // int carry = 0; + // int column = 0; + + // while (true) + // { + // bool targetOne = column < pBits; + // int sum = carry; + + // for (int idx = 0; idx < oneCount; idx++) + // { + // int offset = oneOffsets[idx]; + // if (offset == 0) + // { + // continue; // a_bit for current column, handled below + // } + + // if (offset > column) + // { + // break; + // } + + // int aIndex = column - offset; + // byte aBit = aWindow[aIndex % windowSize]; + // sum += aBit; + // } + + // int requiredParity = targetOne ? 1 : 0; + // int parity = sum & 1; + // byte aCurrent = (byte)(parity == requiredParity ? 0 : 1); + // int columnTotal = sum + aCurrent; + + // if ((columnTotal & 1) != requiredParity) + // { + // if (rented) + // { + // ArrayPool.Shared.Return(aWindowArray); + // } + // divides = false; + // return true; + // } + + // aWindow[column % windowSize] = aCurrent; + // if (aCurrent == 1) + // { + // lastNonZeroA = column; + // } + + // carry = columnTotal >> 1; + // column++; + + // bool doneWithOnes = lastNonZeroA < 0 || column - lastNonZeroA > maxOffset; + // if (column >= pBits && carry == 0 && doneWithOnes) + // { + // break; + // } + + // if (column > pBits + maxOffset + 2 && carry == 0 && doneWithOnes) + // { + // break; + // } + + // if (column > pBits + maxOffset + windowSize + 2) + // { + // if (rented) + // { + // ArrayPool.Shared.Return(aWindowArray); + // } + // divides = false; + // return true; + // } + // } + + // if (rented) + // { + // ArrayPool.Shared.Return(aWindowArray); + // } + // divides = carry == 0; + // Console.WriteLine($"[bittree] Prime={prime}, Divisor={divisor} Decided=True Divides={divides}"); + // if (!divides) + // { + // string message = $"[bittree] p={prime} ruled out divisor={divisor} (parity/carry exhaustion)"; + // _bitTreeStateRepository?.Upsert(prime, message); + // } + + // return true; + // } +#endif #if !DivisorSet_Pow2Groups && (DivisorSet_Predictive || DivisorSet_Percentile || DivisorSet_Additive) [MethodImpl(MethodImplOptions.AggressiveOptimization)] @@ -1121,7 +2080,8 @@ private ulong _currentPrime; { RecordState(currentK); BigInteger powResult = BigInteger.ModPow(BigIntegerNumbers.Two, primeBig, divisor); - if (powResult.IsOne && !IsMersenneValue(prime, divisor)) + // if (powResult.IsOne && !IsMersenneValue(prime, divisor)) + if (powResult.IsOne) { foundDivisor = divisor; return true; @@ -1746,7 +2706,8 @@ private ulong _currentPrime; } RecordState(currentK); - if (divisorCycle == prime && !IsMersenneValue(prime, candidate)) + // if (divisorCycle == prime && !IsMersenneValue(prime, candidate)) + if (divisorCycle == prime) { // A cycle equal to the tested exponent (which is prime in this path) guarantees that the candidate divides // the corresponding Mersenne number because the order of 2 modulo the divisor is exactly p. @@ -2310,7 +3271,8 @@ private ulong _currentPrime; BigInteger divisorBig = (BigInteger)divisor; RecordState((BigInteger)k); BigInteger powResult = BigInteger.ModPow(2, prime, divisorBig); - if (powResult.IsOne && !IsMersenneValue(prime, divisorBig)) + // if (powResult.IsOne && !IsMersenneValue(prime, divisorBig)) + if (powResult.IsOne) { foundDivisor = divisorBig; return true; @@ -2432,7 +3394,7 @@ private ulong _currentPrime; bool passesSmallModuli = remainder10 == 1 || remainder10 == 7; if (passesSmallModuli && ((primeDecimalMask >> remainder10) & 1) != 0) { - if (rem3.NextIsNotDivisible(divisor) && + if (rem3.NextIsNotDivisible(divisor) && rem7.NextIsNotDivisible(divisor) && rem11.NextIsNotDivisible(divisor) && rem13.NextIsNotDivisible(divisor) && @@ -2444,7 +3406,8 @@ private ulong _currentPrime; { RecordState(currentK); BigInteger powResult = BigInteger.ModPow(BigIntegerNumbers.Two, primeBig, divisor); - if (powResult.IsOne && !IsMersenneValue(prime, divisor)) + // if (powResult.IsOne && !IsMersenneValue(prime, divisor)) + if (powResult.IsOne) { foundDivisor = divisor; processedAll = true; @@ -2866,10 +3829,18 @@ private ulong _currentPrime; { if (k > _lastSavedK) { - _kRepository.Upsert(_currentPrime, k); - _lastSavedK = k; + int next = _stateCounter + 1; + if (next >= PerfectNumberConstants.ByDivisorStateSaveInterval) + { + _kRepository.Upsert(_currentPrime, k); + _lastSavedK = k; + _stateCounter = 0; + } + else + { + _stateCounter = next; + } } - _stateCounter = 0; } } @@ -3328,6 +4299,8 @@ private ulong _currentPrime; #endif #if DivisorSet_BitContradiction + private static readonly BigInteger Five = (BigInteger)5; + [MethodImpl(MethodImplOptions.AggressiveOptimization)] private bool CheckDivisorsBitContradiction( ulong prime, @@ -3348,7 +4321,88 @@ private ulong _currentPrime; while (true) { BigInteger divisor = (step * currentK) + BigInteger.One; - bool decided = TryExactBitContradictionCheck(prime, divisor, out bool divides); + if ( + divisor <= ulong.MaxValue && + !PrimeTesterByLastDigit.IsPrimeCpu((ulong)divisor) + ) + { + goto MOVE_NEXT; + } + + ulong effectiveExponent = 0UL; +// if (divisor <= ulong.MaxValue) +// { +// ulong divisor64 = (ulong)divisor; +// MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor64); +// #if DEVICE_GPU +// bool computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentGpu( +// divisor64, +// Accelerator, +// prime, +// divisorData, +// out ulong computedCycle, +// out bool primeOrderFailed); +// #elif DEVICE_HYBRID +// bool computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentHybrid( +// divisor64, +// Accelerator, +// prime, +// divisorData, +// out ulong computedCycle, +// out bool primeOrderFailed); +// #else +// bool computed = MersenneNumberDivisorByDivisorTester.TryCalculateCycleLengthForExponentCpu( +// divisor64, +// prime, +// divisorData, +// out ulong computedCycle, +// out bool primeOrderFailed); +// #endif +// ulong divisorCycle; +// if (!computed || computedCycle == 0UL) +// { +// #if DEVICE_GPU +// divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthGpu( +// divisor64, +// Accelerator, +// divisorData, +// skipPrimeOrderHeuristic: primeOrderFailed); +// #elif DEVICE_HYBRID +// divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthHybrid( +// divisor64, +// Accelerator, +// divisorData, +// skipPrimeOrderHeuristic: primeOrderFailed); +// #else +// divisorCycle = MersenneNumberDivisorByDivisorTester.CalculateCycleLengthCpu( +// divisor64, +// divisorData, +// skipPrimeOrderHeuristic: primeOrderFailed); +// #endif +// } +// else +// { +// divisorCycle = computedCycle; +// } + +// if (divisorCycle > 0UL && (prime % divisorCycle) != 0UL) +// { +// if (hasLimit && currentK >= maxK) +// { +// break; +// } + +// currentK += BigInteger.One; +// continue; +// } + +// if (divisorCycle > 0UL) +// { +// effectiveExponent = divisorCycle; +// } +// } + + bool decided = TryExactBitContradictionCheck(prime, divisor, effectiveExponent, out bool divides); if (!decided) { throw new InvalidOperationException($"BitContradiction check did not decide for p={prime}, divisor={divisor} (k={currentK})."); @@ -3360,6 +4414,8 @@ private ulong _currentPrime; return true; } + RecordState(currentK); + MOVE_NEXT: if (hasLimit && currentK >= maxK) { break; @@ -3424,7 +4480,7 @@ private ulong _currentPrime; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool TryExactBitContradictionCheck(ulong prime, BigInteger divisor, out bool divides) + private static bool TryExactBitContradictionCheck(ulong prime, BigInteger divisor, ulong effectiveExponent, out bool divides) { divides = false; @@ -3472,7 +4528,8 @@ private ulong _currentPrime; } ReadOnlySpan oneOffsetsSlice = oneOffsets[..oneCount]; - bool decided = BitContradictionSolver.TryCheckDivisibilityFromOneOffsets(oneOffsetsSlice, prime, out divides, out var reason); + ulong exponent = effectiveExponent > 0UL ? effectiveExponent : prime; + bool decided = BitContradictionSolver.TryCheckDivisibilityFromOneOffsets(oneOffsetsSlice, exponent, out divides, out var reason); if (decided && divides) { // Verify to avoid false positives on structurally admissible but non-dividing q. diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithTemplate.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithTemplate.cs index d6ce59e..655c50e 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithTemplate.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorPrimeScanSessionWithTemplate.cs @@ -27,6 +27,10 @@ using TesterGpu = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTeste using TesterCpu = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitContradictionDivisorSetForCpuOrder; using TesterHybrid = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitContradictionDivisorSetForHybridOrder; using TesterGpu = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitContradictionDivisorSetForGpuOrder; +#elif DivisorSet_BitTree +using TesterCpu = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitTreeDivisorSetForCpuOrder; +using TesterHybrid = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitTreeDivisorSetForHybridOrder; +using TesterGpu = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitTreeDivisorSetForGpuOrder; #else using TesterCpu = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForSequentialDivisorSetForCpuOrder; using TesterHybrid = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForSequentialDivisorSetForHybridOrder; @@ -91,6 +95,8 @@ internal struct MersenneNumberDivisorByDivisorPrimeScanSessionWithTemplate extension = ".topdown"; #elif DivisorSet_BitContradiction extension = ".bitcontradiction"; +#elif DivisorSet_BitTree + extension = ".bittree"; #endif string stateFile = prime.ToString(CultureInfo.InvariantCulture) + extension; diff --git a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs index 8374c15..683b2e5 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs @@ -149,7 +149,7 @@ public sealed class PrimeOrderCalculatorAccelerator private sealed class CalculatorKernels(Accelerator accelerator) { - public readonly Action, ArrayView1D, int, OrderKernelBuffers> CalculateOrderKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, OrderKernelBuffers>(PrimeOrderKernels.CalculateOrderKernel)).CreateLauncherDelegate, ArrayView1D, int, OrderKernelBuffers>>(); + // public readonly Action, ArrayView1D, int, OrderKernelBuffers> CalculateOrderKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, OrderKernelBuffers>(PrimeOrderKernels.CalculateOrderKernel)).CreateLauncherDelegate, ArrayView1D, int, OrderKernelBuffers>>(); public readonly Kernel CheckFactorsKernel = KernelUtil.GetKernel(accelerator.LoadStreamKernel, Stride1D.Dense>, ulong, ulong, ulong, ulong, ulong, ArrayView1D>(PrimeOrderKernels.CheckFactorsKernel)); @@ -177,7 +177,7 @@ public sealed class PrimeOrderCalculatorAccelerator public readonly Action, int, ulong, ArrayView> SpecialMaxKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, int, ulong, ArrayView>(PrimeOrderGpuHeuristics.EvaluateSpecialMaxCandidatesKernel)).CreateLauncherDelegate, int, ulong, ArrayView>>(); } - public readonly Action, ArrayView1D, int, OrderKernelBuffers> CalculateOrderKernelLauncher; + // public readonly Action, ArrayView1D, int, OrderKernelBuffers> CalculateOrderKernelLauncher; public readonly Kernel CheckFactorsKernel; public readonly Action, ulong, ulong, ulong, ulong, ulong, ArrayView1D> ConvertToStandardKernelLauncher; public readonly Action, ulong, ulong, ArrayView, ArrayView> HeuristicCombinedTrialDivisionKernelLauncher; @@ -421,9 +421,9 @@ public sealed class PrimeOrderCalculatorAccelerator DevicePrimesPow2LastThree = lastDigitSharedTables.DevicePrimesPow2LastThree.View; DevicePrimesPow2LastNine = lastDigitSharedTables.DevicePrimesPow2LastNine.View; - DivisorTables = HeuristicCombinedGpuTables - .GetStaticTables(acceleratorIndex) - .CreateViews(); + // DivisorTables = HeuristicCombinedGpuTables + // .GetStaticTables(acceleratorIndex) + // .CreateViews(); var smallPrimeFactorTables = SmallPrimeFactorGpuTables.GetStaticTables(acceleratorIndex); @@ -436,7 +436,7 @@ public sealed class PrimeOrderCalculatorAccelerator var kernels = _kernels[acceleratorIndex]; - CalculateOrderKernelLauncher = kernels.CalculateOrderKernelLauncher; + // CalculateOrderKernelLauncher = kernels.CalculateOrderKernelLauncher; CheckFactorsKernel = kernels.CheckFactorsKernel; ConvertToStandardKernelLauncher = kernels.ConvertToStandardKernelLauncher; HeuristicCombinedTrialDivisionKernelLauncher = kernels.HeuristicCombinedTrialDivisionKernelLauncher; diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 5294afc..cd8509e 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -272,27 +272,27 @@ internal static partial class PrimeOrderGpuHeuristics var kernelConfig = new CalculateOrderKernelConfig(previousValue, hasPrevious, limit, config.MaxPowChecks, config.StrictMode); ref var buffers = ref gpu.CalculateOrderKernelBuffers; - var kernelLauncher = gpu.CalculateOrderKernelLauncher; + // var kernelLauncher = gpu.CalculateOrderKernelLauncher; AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex); // TODO: Remove the cleaning after the order kernel is modified to always set the result. resultBufferView.MemSetToZero(stream); statusBufferView.MemSetToZero(stream); - kernelLauncher( - stream, - 1, - prime, - kernelConfig, - divisorData.Modulus, - divisorData.NPrime, - divisorData.MontgomeryOne, - divisorData.MontgomeryTwo, - divisorData.MontgomeryTwoSquared, - smallPrimesView, - smallSquaresView, - (int)smallPrimesView.Length, - buffers); + // kernelLauncher( + // stream, + // 1, + // prime, + // kernelConfig, + // divisorData.Modulus, + // divisorData.NPrime, + // divisorData.MontgomeryOne, + // divisorData.MontgomeryTwo, + // divisorData.MontgomeryTwoSquared, + // smallPrimesView, + // smallSquaresView, + // (int)smallPrimesView.Length, + // buffers); byte status = 0; diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index fed08ef..692eaa4 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -21,6 +21,9 @@ using MersenneNumberDivisorByDivisorCpuTesterWithAdditiveGpuOrder = PerfectNumbe using MersenneNumberDivisorByDivisorCpuTesterWithBitContradictionCpuOrder = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitContradictionDivisorSetForCpuOrder; using MersenneNumberDivisorByDivisorCpuTesterWithBitContradictionHybridOrder = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitContradictionDivisorSetForHybridOrder; using MersenneNumberDivisorByDivisorCpuTesterWithBitContradictionGpuOrder = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitContradictionDivisorSetForGpuOrder; +using MersenneNumberDivisorByDivisorCpuTesterWithBitTreeCpuOrder = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitTreeDivisorSetForCpuOrder; +using MersenneNumberDivisorByDivisorCpuTesterWithBitTreeHybridOrder = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitTreeDivisorSetForHybridOrder; +using MersenneNumberDivisorByDivisorCpuTesterWithBitTreeGpuOrder = PerfectNumbers.Core.Cpu.MersenneNumberDivisorByDivisorCpuTesterWithForBitTreeDivisorSetForGpuOrder; namespace PerfectNumbers.Core; @@ -88,6 +91,37 @@ public static class MersenneNumberDivisorByDivisorTester Action pr) => RunBitContradictionCpuOrderFiltered(filteredPrimes, prototype, workerCount, partitionSize, mc, cc, pr)); + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static void Run( + List candidates, + ref MersenneNumberDivisorByDivisorCpuTesterWithBitTreeCpuOrder tester, + Dictionary? previousResults, + ulong startPrime, + Action markComposite, + Action clearComposite, + Action printResult, + int threadCount, + int primesPerTask) + => RunStructTester( + candidates, + ref tester, + previousResults, + startPrime, + markComposite, + clearComposite, + printResult, + threadCount, + primesPerTask, + static ( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithBitTreeCpuOrder prototype, + int workerCount, + int partitionSize, + Action mc, + Action cc, + Action pr) => + RunBitTreeCpuOrderFiltered(filteredPrimes, prototype, workerCount, partitionSize, mc, cc, pr)); + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] public static void Run( List candidates, @@ -243,6 +277,37 @@ public static class MersenneNumberDivisorByDivisorTester Action pr) => RunBitContradictionHybridOrderFiltered(filteredPrimes, prototype, workerCount, partitionSize, mc, cc, pr)); + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static void Run( + List candidates, + ref MersenneNumberDivisorByDivisorCpuTesterWithBitTreeHybridOrder tester, + Dictionary? previousResults, + ulong startPrime, + Action markComposite, + Action clearComposite, + Action printResult, + int threadCount, + int primesPerTask) + => RunStructTester( + candidates, + ref tester, + previousResults, + startPrime, + markComposite, + clearComposite, + printResult, + threadCount, + primesPerTask, + static ( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithBitTreeHybridOrder prototype, + int workerCount, + int partitionSize, + Action mc, + Action cc, + Action pr) => + RunBitTreeHybridOrderFiltered(filteredPrimes, prototype, workerCount, partitionSize, mc, cc, pr)); + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] public static void Run( List candidates, @@ -398,6 +463,37 @@ public static class MersenneNumberDivisorByDivisorTester Action pr) => RunBitContradictionGpuOrderFiltered(filteredPrimes, prototype, workerCount, partitionSize, mc, cc, pr)); + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static void Run( + List candidates, + ref MersenneNumberDivisorByDivisorCpuTesterWithBitTreeGpuOrder tester, + Dictionary? previousResults, + ulong startPrime, + Action markComposite, + Action clearComposite, + Action printResult, + int threadCount, + int primesPerTask) + => RunStructTester( + candidates, + ref tester, + previousResults, + startPrime, + markComposite, + clearComposite, + printResult, + threadCount, + primesPerTask, + static ( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithBitTreeGpuOrder prototype, + int workerCount, + int partitionSize, + Action mc, + Action cc, + Action pr) => + RunBitTreeGpuOrderFiltered(filteredPrimes, prototype, workerCount, partitionSize, mc, cc, pr)); + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] public static void Run( List candidates, @@ -2575,6 +2671,150 @@ public static class MersenneNumberDivisorByDivisorTester })); } + private static void RunBitTreeCpuOrderFiltered( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithBitTreeCpuOrder prototype, + int workerCount, + int partitionSize, + Action markComposite, + Action clearComposite, + Action printResult) + { + MersenneNumberDivisorByDivisorCpuTesterWithBitTreeCpuOrder prototypeCopy = prototype; + int totalCount = filteredPrimes.Count; + if (partitionSize > totalCount) + { + partitionSize = totalCount; + } + + if (workerCount <= 1 || totalCount <= partitionSize) + { + var session = new Cpu.MersenneNumberDivisorByDivisorPrimeScanSessionWithForBitTreeDivisorSetForCpuOrder(prototypeCopy, markComposite, clearComposite, printResult); + for (int i = 0; i < totalCount; i++) + { + session.ProcessPrime(filteredPrimes[i]); + } + + return; + } + + var options = new ParallelOptions + { + MaxDegreeOfParallelism = workerCount, + }; + + Parallel.For( + 0, + (int)Math.Ceiling((double)totalCount / partitionSize), + options, + new Action(partitionIndex => + { + int startIndex = partitionIndex * partitionSize; + int count = Math.Min(partitionSize, totalCount - startIndex); + var session = new Cpu.MersenneNumberDivisorByDivisorPrimeScanSessionWithForBitTreeDivisorSetForCpuOrder(prototypeCopy, markComposite, clearComposite, printResult); + for (int i = 0; i < count; i++) + { + session.ProcessPrime(filteredPrimes[startIndex + i]); + } + })); + } + + private static void RunBitTreeHybridOrderFiltered( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithBitTreeHybridOrder prototype, + int workerCount, + int partitionSize, + Action markComposite, + Action clearComposite, + Action printResult) + { + MersenneNumberDivisorByDivisorCpuTesterWithBitTreeHybridOrder prototypeCopy = prototype; + int totalCount = filteredPrimes.Count; + if (partitionSize > totalCount) + { + partitionSize = totalCount; + } + + if (workerCount <= 1 || totalCount <= partitionSize) + { + var session = new Cpu.MersenneNumberDivisorByDivisorPrimeScanSessionWithForBitTreeDivisorSetForHybridOrder(prototypeCopy, markComposite, clearComposite, printResult); + for (int i = 0; i < totalCount; i++) + { + session.ProcessPrime(filteredPrimes[i]); + } + + return; + } + + var options = new ParallelOptions + { + MaxDegreeOfParallelism = workerCount, + }; + + Parallel.For( + 0, + (int)Math.Ceiling((double)totalCount / partitionSize), + options, + new Action(partitionIndex => + { + int startIndex = partitionIndex * partitionSize; + int count = Math.Min(partitionSize, totalCount - startIndex); + var session = new Cpu.MersenneNumberDivisorByDivisorPrimeScanSessionWithForBitTreeDivisorSetForHybridOrder(prototypeCopy, markComposite, clearComposite, printResult); + for (int i = 0; i < count; i++) + { + session.ProcessPrime(filteredPrimes[startIndex + i]); + } + })); + } + + private static void RunBitTreeGpuOrderFiltered( + List filteredPrimes, + in MersenneNumberDivisorByDivisorCpuTesterWithBitTreeGpuOrder prototype, + int workerCount, + int partitionSize, + Action markComposite, + Action clearComposite, + Action printResult) + { + MersenneNumberDivisorByDivisorCpuTesterWithBitTreeGpuOrder prototypeCopy = prototype; + int totalCount = filteredPrimes.Count; + if (partitionSize > totalCount) + { + partitionSize = totalCount; + } + + if (workerCount <= 1 || totalCount <= partitionSize) + { + var session = new Cpu.MersenneNumberDivisorByDivisorPrimeScanSessionWithForBitTreeDivisorSetForGpuOrder(prototypeCopy, markComposite, clearComposite, printResult); + for (int i = 0; i < totalCount; i++) + { + session.ProcessPrime(filteredPrimes[i]); + } + + return; + } + + var options = new ParallelOptions + { + MaxDegreeOfParallelism = workerCount, + }; + + Parallel.For( + 0, + (int)Math.Ceiling((double)totalCount / partitionSize), + options, + new Action(partitionIndex => + { + int startIndex = partitionIndex * partitionSize; + int count = Math.Min(partitionSize, totalCount - startIndex); + var session = new Cpu.MersenneNumberDivisorByDivisorPrimeScanSessionWithForBitTreeDivisorSetForGpuOrder(prototypeCopy, markComposite, clearComposite, printResult); + for (int i = 0; i < count; i++) + { + session.ProcessPrime(filteredPrimes[startIndex + i]); + } + })); + } + private static void RunBitContradictionHybridOrderFiltered( List filteredPrimes, in MersenneNumberDivisorByDivisorCpuTesterWithBitContradictionHybridOrder prototype, diff --git a/PerfectNumbers.Core/Persistence/BitTreeFrontierRepository.cs b/PerfectNumbers.Core/Persistence/BitTreeFrontierRepository.cs new file mode 100644 index 0000000..c89cb13 --- /dev/null +++ b/PerfectNumbers.Core/Persistence/BitTreeFrontierRepository.cs @@ -0,0 +1,129 @@ +using System; +using System.Collections.Generic; +using System.Globalization; +using System.IO; +using System.Numerics; +using System.Text; +using FASTER.core; + +namespace PerfectNumbers.Core.Persistence; + +/// +/// Stores BitTree frontier batches using Microsoft.FASTER.Core. +/// Each batch is a byte[] payload containing a sequence of serialized states. +/// Maintains a simple head/tail cursor stored inside the KV store. +/// +public sealed class BitTreeFrontierRepository : IDisposable +{ + private readonly FasterKV _store; + private readonly IDevice _logDevice; + private readonly IDevice _objectLogDevice; + private long _headIndex; + private long _tailIndex; + private const string HeadKey = "__head__"; + private const string TailKey = "__tail__"; + + public string DatabasePath { get; } + + public static BitTreeFrontierRepository Open(string directory, string repositoryFileName = "bittree.frontier.faster") + { + Directory.CreateDirectory(directory); + string dbPath = Path.Combine(directory, repositoryFileName); + return new BitTreeFrontierRepository(dbPath); + } + + private BitTreeFrontierRepository(string path) + { + DatabasePath = path; + _logDevice = Devices.CreateLogDevice(path + ".log"); + _objectLogDevice = Devices.CreateLogDevice(path + ".obj.log"); + _store = new FasterKV( + 1L << 20, + new LogSettings + { + LogDevice = _logDevice, + ObjectLogDevice = _objectLogDevice, + }); + + LoadCursors(); + } + + public void AppendBatch(ReadOnlySpan batch) + { + using var session = _store.For(new SimpleFunctions()).NewSession>(); + string key = BuildChunkKey(_tailIndex); + session.Upsert(key, batch.ToArray()); + _tailIndex++; + PersistCursors(session); + session.CompletePending(true); + _store.Log.FlushAndEvict(true); + } + + public bool TryDequeueBatch(int maxCount, out List batches) + { + batches = new List(); + using var session = _store.For(new SimpleFunctions()).NewSession>(); + + while (_headIndex < _tailIndex && batches.Count < maxCount) + { + string key = BuildChunkKey(_headIndex); + if (session.Read(key, out byte[]? stored).Found && stored is not null) + { + batches.Add(stored); + session.Delete(key); + } + + _headIndex++; + } + + + if (batches.Count > 0) + { + PersistCursors(session); + session.CompletePending(true); + _store.Log.FlushAndEvict(true); + return true; + } + + session.CompletePending(true); + _store.Log.FlushAndEvict(true); + return false; + } + + public void Dispose() + { + _store.Dispose(); + _logDevice.Dispose(); + _objectLogDevice.Dispose(); + } + + private void LoadCursors() + { + using var session = _store.For(new SimpleFunctions()).NewSession>(); + _headIndex = ReadCursor(session, HeadKey); + _tailIndex = ReadCursor(session, TailKey); + session.CompletePending(true); + } + + private static long ReadCursor(ClientSession> session, string key) + { + if (session.Read(key, out byte[]? stored).Found && stored is not null) + { + string text = Encoding.UTF8.GetString(stored); + if (long.TryParse(text, NumberStyles.Integer, CultureInfo.InvariantCulture, out long value)) + { + return value; + } + } + + return 0; + } + + private void PersistCursors(ClientSession> session) + { + session.Upsert(HeadKey, Encoding.UTF8.GetBytes(_headIndex.ToString(CultureInfo.InvariantCulture))); + session.Upsert(TailKey, Encoding.UTF8.GetBytes(_tailIndex.ToString(CultureInfo.InvariantCulture))); + } + + private static string BuildChunkKey(long index) => $"chunk-{index}"; +} diff --git a/docs/BitTree.md b/docs/BitTree.md new file mode 100644 index 0000000..584f88f --- /dev/null +++ b/docs/BitTree.md @@ -0,0 +1,52 @@ +# Bit-Tree Divisor Scan + +Plan and notes for the `DivisorSet.BitTree` / `--bydivisor-k-increment=bittree` mode. + +## Goal +- Divisor search via joint reconstruction of unknown bits in both the divisor \(q\) (only LSB=1 is known initially) and the cofactor \(a\), matching column-wise carries in \(q \cdot a = 2^p - 1\). +- Explore the full bit tree (all admissible \(q,a\) combinations) with early pruning from parity reachability, carry budgets, and length limits; no recursion (explicit stack/queue). +- Report primes \(p\) whose entire search space exhausts with no successful branch first. + +## Operational assumptions +- \(p\) is odd; we target prime-like divisors with \(q \equiv 1 \pmod{2p}\) but start from bit patterns (q LSB = 1) and let the tree explore. +- Work in joint (a,q) space, not tied to a single \(k\). All admissible \(q\) consistent with \(2pk+1\) should be reachable through the bit tree. +- Iterative DFS/BFS (explicit stack/queue), no recursion. Process columns streaming (no full 138M+ bit buffers). + +## Implementation steps +- [x] Add `BitTree` mode to CLI (`--bydivisor-k-increment`) and state repositories/files; ensure distinct state extensions. +- [x] Wire mode through CPU/Hybrid/GPU paths with `.bittree` state files. +- [x] Implement core bit-tree solver (frontier/BFS variant): + - Joint state: column index, carry, partial `q` (LSB=1), partial `a`, bounds on `a` length (≤ p - bitlen(q)+1), checks against `allowedMax`. + - Branching: choose next `q` bit (0/1), derive `a` bit from parity+carry; prune on length overflow, parity unreachable, carry cannot collapse, or `q` beyond allowedMax. + - Termination: success iff at/beyond column `p` with carry=0 and no further contributions from `(a,q)`; otherwise contradiction. +- [x] Integrate BitTree solver into `CheckDivisorsBitTree`; legacy per-`k` loop removed for this mode. +- [ ] Logging: emit Decided/Divides and contradiction reasons; throw if solver cannot decide. +- [ ] Tests: + - Unit tests on small \(p\) for success/failure paths. + - Integration test selecting `DivisorSet.BitTree` to ensure composites are detected and primes pass. +- [ ] Docs/help: mention the new mode in help text and README references. + +## Variants considered + +### Variant A (implemented baseline) +- Explicit stack-based DFS over `(a,q)` bits with reachability/pruning on parity and carry budgets. +- Processes columns sequentially without recursion; suited to small/medium `p` in current form. + +### Variant B (frontier/windowed exploration — in progress) +- Implemented: frontier/BFS over branches with sliding window (suffix) storage for `(a,q)` and per-branch carry, bounded by `allowedMax` and `p`; window slides forward to cap memory and drop obsolete history. +- Planned enhancement: tune window size/adaptive policy and add optional frontier compression (baseline/delta) for `p ≥ 138M` and many threads. + +### Variant C (frontier with shared baseline/deltas — planned) +- For each processed column maintain a shared baseline state (e.g., minimal carry/bit sums across the frontier). +- Each branch stores only its delta relative to the baseline (carry delta and the minimal differing `(a,q)` suffix), shrinking per-branch memory. +- When compressing the frontier after a window, pick the baseline that minimises aggregate delta size (e.g., minimal sum of carries/bit-counts across surviving branches). +- Objective: reduce memory further for highly parallel runs (many threads) while preserving full branch coverage. + +## Current status +- Variant A solver is integrated and invoked from `DivisorSet.BitTree` (branch scan path); legacy per-`k` loop is disabled for this mode. +- Variant B is not implemented yet; design notes above. + +## Next steps +- Implement frontier/windowed solver (Variant B) and add mode toggle if needed. +- Prototype baseline/delta compression for the frontier (Variant C) and measure memory savings under parallel load. +- Expand tests to cover both variants on small `p` and ensure integration keeps reporting decisions for all candidates.