diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index ddd2498..4e38fe4 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -112,18 +112,25 @@ private static ThreadLocal? MersenneTesters; } Console.WriteLine("Initializing"); - int gpuPrimeThreads = Math.Max(1, _cliArguments.GpuPrimeThreads); + int gpuPrimeThreads = Math.Max(1, _cliArguments.GpuPrimeThreads * PerfectNumberConstants.ThreadsByAccelerator); int threadCount = Math.Max(1, _cliArguments.ThreadCount); EnvironmentConfiguration.MinK = _cliArguments.MinK; EnvironmentConfiguration.GpuBatchSize = Math.Max(1, _cliArguments.ScanBatchSize); EnvironmentConfiguration.GpuRatio = _cliArguments.GpuRatio; EnvironmentConfiguration.MersenneDevice = _cliArguments.MersenneDevice; - EnvironmentConfiguration.OrderDevice = _cliArguments.OrderDevice; + ComputationDevice orderDevice = _cliArguments.OrderDevice; + if (_cliArguments.MersenneDevice == ComputationDevice.Gpu && + _cliArguments.ByDivisorKIncrementMode == ByDivisorKIncrementMode.BitContradiction) + { + orderDevice = ComputationDevice.Gpu; + } + + EnvironmentConfiguration.OrderDevice = orderDevice; EnvironmentConfiguration.UsePow2GroupDivisors = _cliArguments.ByDivisorKIncrementMode == ByDivisorKIncrementMode.Pow2Groups; EnvironmentConfiguration.ByDivisorSpecialRange = _cliArguments.ByDivisorSpecialRange; - bool useGpu = _cliArguments.OrderDevice == ComputationDevice.Hybrid || _cliArguments.OrderDevice == ComputationDevice.Gpu || + bool useGpu = orderDevice == ComputationDevice.Hybrid || orderDevice == ComputationDevice.Gpu || _cliArguments.MersenneDevice == ComputationDevice.Hybrid || _cliArguments.MersenneDevice == ComputationDevice.Gpu; EnvironmentConfiguration.RollingAccelerators = useGpu @@ -132,6 +139,10 @@ private static ThreadLocal? MersenneTesters; EnvironmentConfiguration.Initialize(); + PrimeOrderCalculatorAccelerator.SkipOrderTablesAndKernels = _cliArguments.MersenneDevice == ComputationDevice.Gpu && + _cliArguments.UseByDivisor && + _cliArguments.ByDivisorKIncrementMode == ByDivisorKIncrementMode.BitContradiction; + // NttGpuMath.GpuTransformBackend = _cliArguments.NttBackend; // NttGpuMath.ReductionMode = _cliArguments.ModReductionMode; _runPrimesOnCpu = _cliArguments.ForcePrimeKernelsOnCpu; @@ -148,7 +159,16 @@ private static ThreadLocal? MersenneTesters; if (useGpu) { Console.WriteLine("Warming up GPU kernels"); - PrimeTester.WarmUpGpuKernels(EnvironmentConfiguration.RollingAccelerators); + if (_cliArguments.UseByDivisor && + _cliArguments.ByDivisorKIncrementMode == ByDivisorKIncrementMode.BitContradiction && + _cliArguments.MersenneDevice == ComputationDevice.Gpu) + { + PrimeTester.WarmUpGpuKernelsForBitContradiction(EnvironmentConfiguration.RollingAccelerators); + } + else + { + PrimeTester.WarmUpGpuKernels(EnvironmentConfiguration.RollingAccelerators); + } } // PrimeTester.WarmUpGpuKernels(gpuPrimeThreads >> 4); @@ -328,11 +348,6 @@ private static ThreadLocal? MersenneTesters; } else if (useByDivisor) { - if (mersenneOnGpu) - { - throw new NotSupportedException("--mersenne-device=gpu is not supported for --mersenne=bydivisor on this execution path."); - } - ByDivisorKIncrementMode kIncrementMode = _cliArguments.ByDivisorKIncrementMode; if (EnvironmentConfiguration.UseCpuOrder) { diff --git a/PerfectNumbers.Core.Tests/BitContradictionIntegrationTests.cs b/PerfectNumbers.Core.Tests/BitContradictionIntegrationTests.cs index b3364fb..83a0bae 100644 --- a/PerfectNumbers.Core.Tests/BitContradictionIntegrationTests.cs +++ b/PerfectNumbers.Core.Tests/BitContradictionIntegrationTests.cs @@ -5,6 +5,7 @@ using Xunit; namespace PerfectNumbers.Core.Tests; +[Trait("Category", "Fast")] public class BitContradictionIntegrationTests { // 2^11 - 1 = 2047 = 23 * 89 @@ -39,11 +40,21 @@ public class BitContradictionIntegrationTests divides.Should().BeFalse(); } + [Fact] + public void BitContradictionAcceptsKnownLargeDivisor() + { + BigInteger divisor = BigInteger.Parse("1209708008767"); + ulong p = 138000001; + bool decided = Invoke(divisor, p, out bool divides); + decided.Should().BeTrue(); + divides.Should().BeTrue(); + } + [Fact] public void BitContradictionRejectsNonDivisorWithDecision() { - // 2^11 - 1 is not divisible by 25 - BigInteger divisor = 25; + // 2^11 - 1 is not divisible by 45 (k=2, q=22*2+1). + BigInteger divisor = 45; ulong p = 11; bool decided = Invoke(divisor, p, out bool divides); decided.Should().BeTrue("small p should be decided exactly"); diff --git a/PerfectNumbers.Core.Tests/BitContradictionKernelsTests.cs b/PerfectNumbers.Core.Tests/BitContradictionKernelsTests.cs new file mode 100644 index 0000000..554b8e7 --- /dev/null +++ b/PerfectNumbers.Core.Tests/BitContradictionKernelsTests.cs @@ -0,0 +1,80 @@ +using System; +using System.Numerics; +using FluentAssertions; +using PerfectNumbers.Core.Gpu.Kernels; +using Xunit; + +namespace PerfectNumbers.Core.Tests; + +public class BitContradictionKernelsTests +{ + [Theory] + [Trait("Category", "Fast")] + [InlineData(31UL, 0UL, 0)] + [InlineData(31UL, 0UL, 5)] + [InlineData(61UL, 7UL, 3)] + [InlineData(36721UL, 1UL, 0)] + public void TryBuildQWordsFromBatchIndex_matches_BigInteger(ulong exponent, ulong batchIndex, int batchOffset) + { + int wordCount = BitContradictionKernels.MaxQWordCount; + ulong[] batchIndexWords = new ulong[wordCount]; + batchIndexWords[0] = batchIndex; + + ulong[] qWords = new ulong[wordCount]; + bool ok = BitContradictionKernels.TryBuildQWordsFromBatchIndex( + exponent, + batchIndexWords, + wordCount, + batchOffset, + qWords, + out int qBitLen); + + ok.Should().BeTrue(); + + BigInteger expectedK = ((BigInteger)batchIndex * BitContradictionKernels.BatchCount) + BigInteger.One + batchOffset; + BigInteger expectedQ = (expectedK * (exponent << 1)) + BigInteger.One; + BigInteger actualQ = ReadBigIntegerFromWords(qWords); + + actualQ.Should().Be(expectedQ); + qBitLen.Should().Be(GetBitLength(expectedQ)); + } + + private static BigInteger ReadBigIntegerFromWords(ReadOnlySpan words) + { + byte[] bytes = new byte[words.Length * sizeof(ulong)]; + for (int w = 0; w < words.Length; w++) + { + ulong word = words[w]; + int byteIndex = w << 3; + bytes[byteIndex] = (byte)word; + bytes[byteIndex + 1] = (byte)(word >> 8); + bytes[byteIndex + 2] = (byte)(word >> 16); + bytes[byteIndex + 3] = (byte)(word >> 24); + bytes[byteIndex + 4] = (byte)(word >> 32); + bytes[byteIndex + 5] = (byte)(word >> 40); + bytes[byteIndex + 6] = (byte)(word >> 48); + bytes[byteIndex + 7] = (byte)(word >> 56); + } + + return new BigInteger(bytes, isUnsigned: true, isBigEndian: false); + } + + private static int GetBitLength(BigInteger value) + { + if (value.IsZero) + { + return 0; + } + + byte[] bytes = value.ToByteArray(isUnsigned: true, isBigEndian: false); + int last = bytes[^1]; + int bits = (bytes.Length - 1) * 8; + while (last != 0) + { + bits++; + last >>= 1; + } + + return bits; + } +} diff --git a/PerfectNumbers.Core.Tests/BitContradictionSolverTests.cs b/PerfectNumbers.Core.Tests/BitContradictionSolverTests.cs index 7f4a309..cfb365f 100644 --- a/PerfectNumbers.Core.Tests/BitContradictionSolverTests.cs +++ b/PerfectNumbers.Core.Tests/BitContradictionSolverTests.cs @@ -16,7 +16,7 @@ public class BitContradictionSolverTests ColumnBounds bounds = ComputeColumnBounds(a, q, 1); bounds.ForcedOnes.Should().Be(1); - bounds.PossibleOnes.Should().Be(1); + bounds.PossibleOnes.Should().Be(2); } [Fact] diff --git a/PerfectNumbers.Core/ByDivisor/BitContradictionSolverWithAMultiplier.cs b/PerfectNumbers.Core/ByDivisor/BitContradictionSolverWithAMultiplier.cs index fccd16d..4d3f7d8 100644 --- a/PerfectNumbers.Core/ByDivisor/BitContradictionSolverWithAMultiplier.cs +++ b/PerfectNumbers.Core/ByDivisor/BitContradictionSolverWithAMultiplier.cs @@ -1,5 +1,6 @@ using System.Runtime.CompilerServices; using System.Numerics; +using PerfectNumbers.Core; namespace PerfectNumbers.Core.ByDivisor; @@ -119,6 +120,7 @@ internal sealed class PrivateWorkSet /// interval-based carry propagation described in the BitContradictions plan document. /// The solver itself will build on top of these primitives. /// +[DeviceDependentTemplate(typeof(ComputationDevice))] internal static class BitContradictionSolverWithAMultiplier { private const int PingPongBlockColumns = 4096; @@ -1510,6 +1512,11 @@ internal static class BitContradictionSolverWithAMultiplier out bool divides) #endif { + if (!_initialized) + { + Initialize(); + } + PrivateWorkSet workSet = _privateWorkSet; divides = false; @@ -2070,6 +2077,22 @@ internal static class BitContradictionSolverWithAMultiplier // } divides = true; + if (p <= 63 && divides) + { + ulong qLow = BuildQLow64(qOneOffsets); + if ((qLow & 1UL) == 0UL) + { + divides = false; + } + else + { + BigInteger qExact = new BigInteger(qLow); + if (BigInteger.ModPow(2, (int)p, qExact) != BigInteger.One) + { + divides = false; + } + } + } #if DETAILED_LOG reason = ContradictionReason.None; PrintStats(); diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs index 59fd888..1706ce2 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs @@ -11,7 +11,9 @@ using System.Globalization; using System.Collections.Generic; using System.Collections.Concurrent; using System.Buffers; +using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; +using PerfectNumbers.Core.Gpu.Kernels; using PerfectNumbers.Core.Gpu.Accelerators; using PerfectNumbers.Core.ByDivisor; using PerfectNumbers.Core.Persistence; @@ -58,6 +60,226 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN private GpuUInt128WorkSet _divisorScanGpuWorkSet; private static readonly ConcurrentDictionary _checkedPow2MinusOne = new(); private ulong _currentPrime; + +#if DEVICE_GPU + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private bool CheckDivisorsBitContradictionGpu( + ulong prime, + ushort primeDecimalMask, + BigInteger allowedMax, + BigInteger minK, + BigInteger step, + out bool processedAll, + out BigInteger foundDivisor) + { + processedAll = false; + foundDivisor = BigInteger.Zero; + + BigInteger currentK = minK; + bool hasLimit = !allowedMax.IsZero; + BigInteger maxK = hasLimit ? (allowedMax - BigInteger.One) / step : BigInteger.Zero; + + BigInteger maxQ = (BigInteger.One << BitContradictionKernels.MaxQBitLength) - BigInteger.One; + BigInteger maxKByBits = (maxQ - BigInteger.One) / step; + if (!hasLimit || maxKByBits < maxK) + { + maxK = maxKByBits; + } + + if (maxK < currentK) + { + processedAll = hasLimit; + return false; + } + + PrimeOrderCalculatorAccelerator gpu = Accelerator; + int maxThreadCount = PerfectNumberConstants.BitContradictionGpuThreads; + int batchSize = PerfectNumberConstants.BitContradictionGpuBatchCount; + BigInteger batchSizeBig = batchSize; + int wordCount = BitContradictionKernels.MaxQWordCount; + int inputCapacity = maxThreadCount * wordCount; + gpu.EnsureCapacity(PerfectNumberConstants.DefaultFactorsBuffer, inputCapacity); + gpu.EnsurePartialFactorCapacity(inputCapacity); + + var stream = AcceleratorStreamPool.Rent(gpu.AcceleratorIndex); + var kernels = GpuKernelPool.GetOrAddKernels(gpu.AcceleratorIndex, stream, KernelType.BitContradictionScan); + var kernel = kernels.BitContradiction!; + + int[] found = new int[maxThreadCount]; + ulong[] batchIndexWords = new ulong[inputCapacity]; + + BigInteger currentBatch = (currentK - BigInteger.One) / batchSizeBig; + currentK = (currentBatch * batchSizeBig) + BigInteger.One; + BigInteger maxBatch = (maxK - BigInteger.One) / batchSizeBig; + + while (true) + { + + if (hasLimit && currentBatch > maxBatch) + { + processedAll = true; + break; + } + + BigInteger remainingBatches = maxBatch - currentBatch + BigInteger.One; + if (remainingBatches <= BigInteger.Zero) + { + processedAll = hasLimit; + break; + } + + int threadCount = maxThreadCount; + if (remainingBatches < threadCount) + { + threadCount = (int)remainingBatches; + } + + if (currentBatch > ulong.MaxValue) + { + throw new InvalidOperationException($"BitContradiction GPU batch index overflow for p={prime}, batch={currentBatch}"); + } + + BigInteger maxBatchIndex = currentBatch + (threadCount - 1); + if (maxBatchIndex > ulong.MaxValue) + { + throw new InvalidOperationException($"BitContradiction GPU batch index overflow for p={prime}, batch={currentBatch}, threads={threadCount}"); + } + + int countQ = batchSize; + + for (int t = 0; t < threadCount; t++) + { + BigInteger batchIndex = currentBatch + t; + WriteBigIntegerToWords(batchIndex, batchIndexWords.AsSpan(t * wordCount, wordCount)); + } + + gpu.InputView.CopyFromCPU(stream, batchIndexWords.AsSpan(0, threadCount * wordCount)); + kernel(stream, threadCount, prime, gpu.InputView, countQ, gpu.OutputIntView); + stream.Synchronize(); + gpu.OutputIntView.CopyToCPU(stream, found.AsSpan(0, threadCount)); + + for (int i = 0; i < threadCount; i++) + { + int foundIndex = found[i]; + if (foundIndex >= 0) + { + BigInteger batchIndex = currentBatch + i; + BigInteger candidateK = (batchIndex * batchSizeBig) + BigInteger.One + foundIndex; + if (hasLimit && candidateK > maxK) + { + processedAll = true; + AcceleratorStreamPool.Return(gpu.AcceleratorIndex, stream); + return false; + } + + BigInteger candidateQ = (step * candidateK) + BigInteger.One; + Span debugQWords = stackalloc ulong[BitContradictionKernels.MaxQWordCount]; + if (!BitContradictionKernels.TryBuildQWordsFromBatchIndex(prime, batchIndexWords.AsSpan(i * wordCount, wordCount), wordCount, foundIndex, debugQWords, out int debugQBitLen)) + { + throw new InvalidOperationException($"BitContradiction GPU qWords overflow for p={prime}, k={candidateK}, batch={batchIndex}, index={foundIndex}"); + } + BigInteger debugQ = ReadBigIntegerFromWords(debugQWords); + if (debugQ != candidateQ) + { + throw new InvalidOperationException($"BitContradiction GPU qWords mismatch for p={prime}, k={candidateK}, batch={batchIndex}, index={foundIndex}, cpuQ={candidateQ}, gpuQ={debugQ}, gpuQBits={debugQBitLen}"); + } + #if DivisorSet_BitContradiction + int[]? qBits = _qBits; + if (qBits == null || qBits.Length < (int)prime) + { + qBits = new int[(int)prime]; + _qBits = qBits; + } + bool cpuDecided = TryExactBitContradictionCheck(prime, candidateQ, candidateK, out bool cpuDivides); + if (!cpuDecided) + { + throw new InvalidOperationException($"BitContradiction CPU undecided for p={prime}, q={candidateQ}, k={candidateK}, batch={batchIndex}"); + } + if (!cpuDivides) + { + throw new InvalidOperationException($"BitContradiction CPU/GPU mismatch for p={prime}, q={candidateQ}, k={candidateK}, batch={batchIndex}"); + } + #endif + bool dividesExact = BigInteger.ModPow(2, prime, candidateQ) == BigInteger.One; + if (!dividesExact) + { + throw new InvalidOperationException($"BitContradiction CPU/GPU false positive for p={prime}, q={candidateQ}, k={candidateK}, batch={batchIndex}"); + } + foundDivisor = candidateQ; + AcceleratorStreamPool.Return(gpu.AcceleratorIndex, stream); + return true; + } + } + + ContinueOuter: + RecordState(currentK); + if (hasLimit && currentK > maxK) + { + processedAll = true; + break; + } + + currentBatch += threadCount; + currentK = (currentBatch * batchSizeBig) + BigInteger.One; + } + + AcceleratorStreamPool.Return(gpu.AcceleratorIndex, stream); + return false; + } + + private static void WriteBigIntegerToWords(BigInteger value, Span words) + { + Span bytes = stackalloc byte[BitContradictionKernels.MaxQWordCount * 8]; + if (!value.TryWriteBytes(bytes, out int written, isUnsigned: true, isBigEndian: false)) + { + written = 0; + } + + int wordCount = words.Length; + for (int w = 0; w < wordCount; w++) + { + ulong word = 0UL; + int byteIndex = w << 3; + int limit = written - byteIndex; + if (limit > 0) + { + if (limit > 8) + { + limit = 8; + } + + for (int b = 0; b < limit; b++) + { + word |= (ulong)bytes[byteIndex + b] << (b << 3); + } + } + + words[w] = word; + } + } + + private static BigInteger ReadBigIntegerFromWords(ReadOnlySpan words) + { + Span bytes = stackalloc byte[BitContradictionKernels.MaxQWordCount * 8]; + int wordCount = words.Length; + for (int w = 0; w < wordCount; w++) + { + ulong word = words[w]; + int byteIndex = w << 3; + bytes[byteIndex] = (byte)word; + bytes[byteIndex + 1] = (byte)(word >> 8); + bytes[byteIndex + 2] = (byte)(word >> 16); + bytes[byteIndex + 3] = (byte)(word >> 24); + bytes[byteIndex + 4] = (byte)(word >> 32); + bytes[byteIndex + 5] = (byte)(word >> 40); + bytes[byteIndex + 6] = (byte)(word >> 48); + bytes[byteIndex + 7] = (byte)(word >> 56); + } + + return new BigInteger(bytes, isUnsigned: true, isBigEndian: false); + } +#endif + [ThreadStatic] private static Dictionary? _pow2DivisorPowModCache; #if DETAILED_LOG @@ -578,7 +800,6 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN out processedAll, out foundDivisor); } -#endif // return CheckDivisorsLarge( // prime, @@ -596,6 +817,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN step, out processedAll, out foundDivisor); +#endif } #if DivisorSet_Pow2Groups @@ -1421,6 +1643,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN return true; } + ContinueOuter: RecordState(currentK); MOVE_NEXT: @@ -1732,6 +1955,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN { processedAll = false; foundDivisor = BigInteger.Zero; + bool coversRange = true; if (minK > ulong.MaxValue) { @@ -2130,7 +2354,8 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN { if (((primeDecimalMask >> remainder10) & 1) != 0) { - RecordState(currentK); + ContinueOuter: + RecordState(currentK); BigInteger powResult = BigInteger.ModPow(BigIntegerNumbers.Two, primeBig, divisor); // if (powResult.IsOne && !IsMersenneValue(prime, divisor)) if (powResult.IsOne) @@ -2757,7 +2982,8 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN #endif } - RecordState(currentK); + ContinueOuter: + RecordState(currentK); // if (divisorCycle == prime && !IsMersenneValue(prime, candidate)) if (divisorCycle == prime) { @@ -2877,7 +3103,8 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN if (residueStepper.IsAdmissible()) { // Stopwatch sw = Stopwatch.StartNew(); - RecordState(currentK); + ContinueOuter: + RecordState(currentK); ulong powMod; bool isPow2 = IsPowerOfTwo(currentK); #if DETAILED_LOG @@ -2978,7 +3205,8 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN byte remainder8 = rem8.ComputeNext(divisor); if (remainder8 == 1 || remainder8 == 7) { - RecordState(currentK); + ContinueOuter: + RecordState(currentK); BigInteger powResult = BigInteger.ModPow(BigIntegerNumbers.Two, primeBig, divisor); if (powResult.IsOne) { @@ -3677,7 +3905,8 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN byte remainder8 = rem8.ComputeNext(divisor); if (remainder8 == 1 || remainder8 == 7) { - RecordState(currentK); + ContinueOuter: + RecordState(currentK); BigInteger powResult = BigInteger.ModPow(BigIntegerNumbers.Two, primeBig, divisor); // if (powResult.IsOne && !IsMersenneValue(prime, divisor)) if (powResult.IsOne) @@ -4031,7 +4260,8 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN divisorCycle = computedCycle; } - RecordState(currentK); + ContinueOuter: + RecordState(currentK); if (divisorCycle == prime && !IsMersenneValue(prime, divisor)) { foundDivisor = divisor; @@ -4585,6 +4815,10 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN out bool processedAll, out BigInteger foundDivisor) { +#if DEVICE_GPU + return CheckDivisorsBitContradictionGpu(prime, primeDecimalMask, allowedMax, minK, step, out processedAll, out foundDivisor); +#else + processedAll = true; foundDivisor = BigInteger.Zero; @@ -4643,6 +4877,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN } MOVE_NEXT: + ContinueOuter: RecordState(currentK); if (hasLimit && currentK >= maxK) { @@ -4707,6 +4942,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN // step, // out processedAll, // out foundDivisor); +#endif } [ThreadStatic] diff --git a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs index 683b2e5..44952cf 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs @@ -56,7 +56,11 @@ public sealed class PrimeOrderCalculatorAccelerator internal static void WarmUp() { Accelerator[] accelerators = AcceleratorPool.Shared.Accelerators; - int acceleratorCount = accelerators.Length; + int acceleratorCount = EnvironmentConfiguration.RollingAccelerators; + if (acceleratorCount > accelerators.Length) + { + acceleratorCount = accelerators.Length; + } for (var i = 0; i < acceleratorCount; i++) { Console.WriteLine($"Preparing accelerator {i}..."); @@ -64,6 +68,13 @@ public sealed class PrimeOrderCalculatorAccelerator AcceleratorStreamPool.WarmUp(i); // Don't take this from the pool as quick uploads of data to the accelerator consumes much of GPU's memory and throws. AcceleratorStream stream = accelerator.CreateStream(); + if (SkipOrderTablesAndKernels) + { + stream.Synchronize(); + stream.Dispose(); + continue; + } + // TODO: Review which tables are needed when final execution path is defined. LastDigitGpuTables.WarmUp(i, stream); _kernels[i] = new CalculatorKernels(accelerator); @@ -99,6 +110,9 @@ public sealed class PrimeOrderCalculatorAccelerator private static readonly Accelerator[] _accelerators = AcceleratorPool.Shared.Accelerators; + public static bool SkipOrderTablesAndKernels; + + public readonly Accelerator Accelerator; public readonly int AcceleratorIndex; @@ -411,6 +425,28 @@ public sealed class PrimeOrderCalculatorAccelerator CalculateOrderWideExponentBufferView = CalculateOrderWideExponentBuffer.View; CalculateOrderWideRemainderBufferView = CalculateOrderWideRemainderBuffer.View; + if (SkipOrderTablesAndKernels) + { + DevicePrimesLastOne = default; + DevicePrimesLastSeven = default; + DevicePrimesLastThree = default; + DevicePrimesLastNine = default; + DevicePrimesPow2LastOne = default; + DevicePrimesPow2LastSeven = default; + DevicePrimesPow2LastThree = default; + DevicePrimesPow2LastNine = default; + SmallPrimeFactorPrimes = default; + SmallPrimeFactorSquares = default; + CheckFactorsKernel = default; + ConvertToStandardKernelLauncher = default; + HeuristicCombinedTrialDivisionKernelLauncher = default; + Pow2ModKernelLauncher = default; + Pow2ModWideKernelLauncher = default; + SmallPrimeFactorKernelLauncher = default; + SpecialMaxKernelLauncher = default; + return; + } + LastDigitGpuTables lastDigitSharedTables = LastDigitGpuTables.GetStaticTables(acceleratorIndex); DevicePrimesLastOne = lastDigitSharedTables.DevicePrimesLastOne.View; DevicePrimesLastSeven = lastDigitSharedTables.DevicePrimesLastSeven.View; diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index 5e4c243..52a6aa0 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -2,6 +2,7 @@ using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Runtime; using PerfectNumbers.Core.Gpu.Accelerators; +using PerfectNumbers.Core.Gpu.Kernels; namespace PerfectNumbers.Core.Gpu; @@ -26,7 +27,7 @@ public class GpuKernelPool internal static KernelContainer GetOrAddKernels(int acceleratorIndex, AcceleratorStream stream, KernelType kernelType) { - var pool = _pool ??= []; + var pool = _pool ??= new KernelContainer[_accelerators.Length]; if (pool[acceleratorIndex] is not { } kernels) { kernels = new(); @@ -58,6 +59,11 @@ public class GpuKernelPool InitPow2ModOrderKernelScan(accelerator, kernels); } + if (kernelType.HasFlag(KernelType.BitContradictionScan)) + { + InitBitContradictionKernel(accelerator, kernels); + } + PreloadStaticTables(accelerator, kernels, stream, kernelType); pool[acceleratorIndex] = kernels; @@ -113,6 +119,18 @@ public class GpuKernelPool }); } + static void InitBitContradictionKernel(Accelerator accelerator, KernelContainer kernels) + { + InitOnce(ref kernels.BitContradiction, () => + { + var loaded = accelerator.LoadAutoGroupedStreamKernel, int, ArrayView1D>(BitContradictionKernels.BitContradictionKernelScan); + + var kernel = KernelUtil.GetKernel(loaded); + + return kernel.CreateLauncherDelegate, int, ArrayView1D>>(); + }); + } + static void InitPow2ModOrderKernelScan(Accelerator accelerator, KernelContainer kernels) { InitOnce(ref kernels.Pow2ModOrder, () => diff --git a/PerfectNumbers.Core/Gpu/KernelContainer.cs b/PerfectNumbers.Core/Gpu/KernelContainer.cs index 65c4594..f6e06d2 100644 --- a/PerfectNumbers.Core/Gpu/KernelContainer.cs +++ b/PerfectNumbers.Core/Gpu/KernelContainer.cs @@ -22,6 +22,9 @@ public sealed class KernelContainer public Action, ArrayView1D>? Pow2ModOrder; + public Action, int, + ArrayView1D>? BitContradiction; + public Kernel? SmallPrimeFactor; public Action, int, MontgomeryDivisorData, ArrayView1D, ArrayView1D>? SpecialMax; diff --git a/PerfectNumbers.Core/Gpu/KernelType.cs b/PerfectNumbers.Core/Gpu/KernelType.cs index 654808d..fb88a3c 100644 --- a/PerfectNumbers.Core/Gpu/KernelType.cs +++ b/PerfectNumbers.Core/Gpu/KernelType.cs @@ -10,5 +10,6 @@ public enum KernelType OrderKernelScan, IncrementalKernelScan, IncrementalOrderKernelScan, - Pow2ModOrderKernelScan + Pow2ModOrderKernelScan, + BitContradictionScan } diff --git a/PerfectNumbers.Core/Gpu/Kernels/BitContradictionKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/BitContradictionKernels.cs new file mode 100644 index 0000000..a55a6f0 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Kernels/BitContradictionKernels.cs @@ -0,0 +1,1480 @@ + +using ILGPU; +using ILGPU.Algorithms; +using ILGPU.Runtime; +using PerfectNumbers.Core; + +namespace PerfectNumbers.Core.Gpu.Kernels; + +internal static class BitContradictionKernels +{ + public const int BatchCount = PerfectNumberConstants.BitContradictionGpuBatchCount; + public const int ForcedALowBits = 1024; + public const int ForcedALowWords = 16; + public const int MaxQBitLength = 4096; + public const int MaxQWordCount = MaxQBitLength / 64; + public const int MaxQOffsets = 4096; + public const int MaxAWordCount = MaxQBitLength / 64; + + private const int InverseIterations = 10; + private const int HighBitCarryPrefilterColumns = 16; + private const int TopDownBorrowColumns = 4096; + private const int TopDownBorrowMinUnknown = 12; + private const int TailCarryBatchColumns = 16384; + private const int DeltaColumnsAtOnce = 8; + private const int DeltaLength = 1 << DeltaColumnsAtOnce; + private const int IntervalBufferLength = 64; + + private struct CarryRange + { + public long Min; + public long Max; + + public CarryRange(long min, long max) + { + Min = min; + Max = max; + } + + public static CarryRange Zero => new(0, 0); + } + + public static void BitContradictionKernelScan( + Index1D index, + ulong exponent, + ArrayView1D batchIndexWords, + int countQ, + ArrayView1D foundOut) + { + int slot = index; + foundOut[slot] = -1; + + if (countQ <= 0) + { + return; + } + + int localCountQ = countQ; + ulong twoP = exponent << 1; + int startOffset = slot * MaxQWordCount; + + var baseKWords = LocalMemory.Allocate(MaxQWordCount); + var kWords = LocalMemory.Allocate(MaxQWordCount); + var qWords = LocalMemory.Allocate(MaxQWordCount); + var qOffsets = LocalMemory.Allocate(MaxQOffsets); + var qPrefix = LocalMemory.Allocate(ForcedALowWords); + var reversePrefix = LocalMemory.Allocate(ForcedALowWords); + var prefixMask = LocalMemory.Allocate(ForcedALowWords); + var invWords = LocalMemory.Allocate(ForcedALowWords); + var tmpWords = LocalMemory.Allocate(ForcedALowWords); + var tmpWords2 = LocalMemory.Allocate(ForcedALowWords); + var aLowWords = LocalMemory.Allocate(ForcedALowWords); + var delta8 = LocalMemory.Allocate(DeltaLength); + + var knownOne0 = LocalMemory.Allocate(MaxAWordCount); + var knownOne1 = LocalMemory.Allocate(MaxAWordCount); + var knownZero0 = LocalMemory.Allocate(MaxAWordCount); + var knownZero1 = LocalMemory.Allocate(MaxAWordCount); + var qMaskWords = LocalMemory.Allocate(MaxQWordCount); + var aOneWin = LocalMemory.Allocate(MaxQWordCount); + var aZeroWin = LocalMemory.Allocate(MaxQWordCount); + + var lo0 = LocalMemory.Allocate(IntervalBufferLength); + var hi0 = LocalMemory.Allocate(IntervalBufferLength); + var lo1 = LocalMemory.Allocate(IntervalBufferLength); + var hi1 = LocalMemory.Allocate(IntervalBufferLength); + + for (int i = 0; i < MaxQWordCount; i++) + { + baseKWords[i] = batchIndexWords[startOffset + i]; + } + + if (MultiplyWordsByUInt64(baseKWords, kWords, MaxQWordCount, (ulong)BatchCount) != 0UL) + { + return; + } + + if (AddUInt64ToWords(kWords, MaxQWordCount, 1UL) != 0UL) + { + return; + } + + CopyWords(kWords, baseKWords, MaxQWordCount); + + for (int batch = 0; batch < localCountQ; batch++) + { + CopyWords(baseKWords, kWords, MaxQWordCount); + if (batch != 0) + { + if (AddUInt64ToWords(kWords, MaxQWordCount, (ulong)batch) != 0UL) + { + break; + } + } + + if (MultiplyWordsByUInt64(kWords, qWords, MaxQWordCount, twoP) != 0UL) + { + continue; + } + + if (AddUInt64ToWords(qWords, MaxQWordCount, 1UL) != 0UL) + { + continue; + } + + int qBitLen = GetBitLength(qWords, MaxQWordCount); + if (qBitLen <= 0 || qBitLen > MaxQBitLength) + { + continue; + } + + if ((ulong)qBitLen >= exponent) + { + continue; + } + + int pLong = (int)exponent; + int maxAllowedA = pLong - qBitLen; + if (maxAllowedA < 0) + { + continue; + } + + int qWordCount = (qBitLen + 63) >> 6; + int offsetCount; + if (!BuildQOneOffsetsWords(qWords, qWordCount, qOffsets, out offsetCount)) + { + continue; + } + + if (offsetCount <= 0) + { + continue; + } + + int windowSize = qBitLen < ForcedALowBits ? ForcedALowBits : qBitLen; + windowSize = (windowSize + 63) & ~63; + int aWordCount = (windowSize + 63) >> 6; + + ClearWords(knownOne0, aWordCount); + ClearWords(knownOne1, aWordCount); + ClearWords(knownZero0, aWordCount); + ClearWords(knownZero1, aWordCount); + ClearWords(qMaskWords, qWordCount); + ClearWords(aOneWin, qWordCount); + ClearWords(aZeroWin, qWordCount); + + ClearWords(qPrefix, ForcedALowWords); + CopyWords(qWords, qPrefix, ForcedALowWords); + CopyWords(qPrefix, invWords, ForcedALowWords); + ComputeInverseMod2k(qPrefix, invWords, tmpWords, tmpWords2, ForcedALowWords); + NegateMod2k(invWords, aLowWords, ForcedALowWords); + + for (int i = 0; i < qWordCount; i++) + { + qMaskWords[i] = qWords[i]; + } + + ulong lastWordMask = (qBitLen & 63) == 0 ? ulong.MaxValue : ((1UL << (qBitLen & 63)) - 1UL); + if (qWordCount > 0) + { + qMaskWords[qWordCount - 1] &= lastWordMask; + } + + int maxFixed = maxAllowedA < (ForcedALowBits - 1) ? maxAllowedA : (ForcedALowBits - 1); + if (maxFixed >= 0) + { + int fullWords = maxFixed >> 6; + int lastBit = maxFixed & 63; + for (int w = 0; w < fullWords; w++) + { + ulong wordValue = aLowWords[w]; + knownOne0[w] = wordValue; + knownZero0[w] = ~wordValue; + } + + if (fullWords < aWordCount) + { + ulong wordValue = aLowWords[fullWords]; + ulong mask = lastBit == 63 ? ulong.MaxValue : ((1UL << (lastBit + 1)) - 1UL); + knownOne0[fullWords] = wordValue & mask; + knownZero0[fullWords] = (~wordValue) & mask; + } + } + + int maxKnownA = maxFixed >= 0 ? maxFixed : -1; + + int maxFixedPrefilter = 512 - 1; + if (maxFixedPrefilter > maxKnownA) + { + maxFixedPrefilter = maxKnownA; + } + + if (maxFixedPrefilter > maxAllowedA) + { + maxFixedPrefilter = maxAllowedA; + } + + if (maxFixedPrefilter >= 8) + { + ReverseBits1024(qPrefix, reversePrefix); + + bool prefixOk = true; + CarryRange carryLow = CarryRange.Zero; + for (int column = 0; column <= maxFixedPrefilter; column++) + { + int forced = 0; + if (column <= 1023) + { + int shift = 1023 - column; + int wordShift = shift >> 6; + int bitShift = shift & 63; + for (int word = 0; word < ForcedALowWords; word++) + { + int src = word + wordShift; + ulong value = 0UL; + if (src < ForcedALowWords) + { + value = reversePrefix[src] >> bitShift; + if (bitShift != 0 && src + 1 < ForcedALowWords) + { + value |= reversePrefix[src + 1] << (64 - bitShift); + } + } + + prefixMask[word] = value; + } + + for (int word = 0; word < ForcedALowWords; word++) + { + forced += XMath.PopCount(knownOne0[word] & prefixMask[word]); + } + } + else + { + for (int i = 0; i < offsetCount; i++) + { + int t = qOffsets[i]; + if (t > column) + { + break; + } + + int aIndex = column - t; + int word = aIndex >> 6; + ulong mask = 1UL << (aIndex & 63); + if ((knownOne0[word] & mask) != 0UL) + { + forced++; + } + } + } + + if (!TryPropagateCarry(ref carryLow, forced, forced, 1)) + { + prefixOk = false; + break; + } + } + + if (!prefixOk) + { + continue; + } + } + + int segmentBase = 0; + int windowColumn = 0; + + int initialState = GetAKnownStateRowAware(0, segmentBase, windowSize, knownOne0, knownZero0, knownOne1, knownZero1); + if (initialState == 1) + { + aOneWin[0] |= 1UL; + } + else if (initialState == 2) + { + aZeroWin[0] |= 1UL; + } + + if (!TryRunHighBitAndBorrowPrefiltersCombinedGpu(qOffsets, offsetCount, pLong, maxAllowedA, lo0, hi0, lo1, hi1)) + { + continue; + } + + BuildStableUnknownDelta8(offsetCount, delta8); + + CarryRange carry = CarryRange.Zero; + bool ok = TryProcessBottomUpBlockRowAwareGpu( + qOffsets, + offsetCount, + qMaskWords, + qWordCount, + lastWordMask, + maxAllowedA, + 0, + pLong, + ref carry, + ref maxKnownA, + ref segmentBase, + knownOne0, + knownOne1, + knownZero0, + knownZero1, + windowSize, + aWordCount, + aOneWin, + aZeroWin, + ref windowColumn, + delta8, + offsetCount); + + if (ok) + { + foundOut[slot] = batch; + return; + } + } + } + internal static bool TryBuildQWordsFromBatchIndex( + ulong exponent, + ReadOnlySpan batchIndexWords, + int wordCount, + int batchOffset, + Span qWords, + out int qBitLen) + { + Span baseKWords = stackalloc ulong[wordCount]; + Span kWords = stackalloc ulong[wordCount]; + + for (int i = 0; i < wordCount; i++) + { + baseKWords[i] = batchIndexWords[i]; + } + + if (MultiplyWordsByUInt64(baseKWords, kWords, wordCount, (ulong)BatchCount) != 0UL) + { + qBitLen = 0; + return false; + } + + if (AddUInt64ToWords(kWords, wordCount, 1UL) != 0UL) + { + qBitLen = 0; + return false; + } + + if (batchOffset != 0) + { + if (AddUInt64ToWords(kWords, wordCount, (ulong)batchOffset) != 0UL) + { + qBitLen = 0; + return false; + } + } + + if (MultiplyWordsByUInt64(kWords, qWords, wordCount, exponent << 1) != 0UL) + { + qBitLen = 0; + return false; + } + + if (AddUInt64ToWords(qWords, wordCount, 1UL) != 0UL) + { + qBitLen = 0; + return false; + } + + qBitLen = GetBitLength(qWords, wordCount); + return true; + } + + private static ulong AddUInt64ToWords(Span words, int wordCount, ulong value) + { + ulong carry = value; + for (int i = 0; i < wordCount; i++) + { + if (carry == 0UL) + { + return 0UL; + } + + ulong sum = words[i] + carry; + words[i] = sum; + carry = sum < carry ? 1UL : 0UL; + } + + return carry; + } + + private static ulong MultiplyWordsByUInt64(ReadOnlySpan source, Span destination, int wordCount, ulong factor) + { + ulong carry = 0UL; + ulong carryExtra = 0UL; + for (int i = 0; i < wordCount; i++) + { + ulong high; + ulong low; + MultiplyPartsGpu(source[i], factor, out high, out low); + + ulong sum = low + carry; + ulong carryOut = sum < low ? 1UL : 0UL; + if (carryExtra != 0UL) + { + sum += 1UL; + if (sum == 0UL) + { + carryOut++; + } + } + + destination[i] = sum; + ulong nextCarry = high + carryOut; + carryExtra = nextCarry < high ? 1UL : 0UL; + carry = nextCarry; + } + + return carry + carryExtra; + } + + private static int GetBitLength(ReadOnlySpan words, int wordCount) + { + for (int i = wordCount - 1; i >= 0; i--) + { + ulong word = words[i]; + if (word != 0UL) + { + int leading = XMath.LeadingZeroCount(word); + return (i << 6) + (64 - leading); + } + } + + return 0; + } + + private static ulong AddUInt64ToWords(ArrayView words, int wordCount, ulong value) + { + ulong carry = value; + for (int i = 0; i < wordCount; i++) + { + if (carry == 0UL) + { + return 0UL; + } + + ulong sum = words[i] + carry; + words[i] = sum; + carry = sum < carry ? 1UL : 0UL; + } + + return carry; + } + private static ulong MultiplyWordsByUInt64(ArrayView source, ArrayView destination, int wordCount, ulong factor) + { + ulong carry = 0UL; + ulong carryExtra = 0UL; + for (int i = 0; i < wordCount; i++) + { + ulong high; + ulong low; + MultiplyPartsGpu(source[i], factor, out high, out low); + + ulong sum = low + carry; + ulong carryOut = sum < low ? 1UL : 0UL; + if (carryExtra != 0UL) + { + sum += 1UL; + if (sum == 0UL) + { + carryOut++; + } + } + + destination[i] = sum; + ulong nextCarry = high + carryOut; + carryExtra = nextCarry < high ? 1UL : 0UL; + carry = nextCarry; + } + + return carry + carryExtra; + } + + private static ulong ReverseBits64(ulong value) + { + value = ((value >> 1) & 0x5555555555555555UL) | ((value & 0x5555555555555555UL) << 1); + value = ((value >> 2) & 0x3333333333333333UL) | ((value & 0x3333333333333333UL) << 2); + value = ((value >> 4) & 0x0F0F0F0F0F0F0F0FUL) | ((value & 0x0F0F0F0F0F0F0F0FUL) << 4); + value = ((value >> 8) & 0x00FF00FF00FF00FFUL) | ((value & 0x00FF00FF00FF00FFUL) << 8); + value = ((value >> 16) & 0x0000FFFF0000FFFFUL) | ((value & 0x0000FFFF0000FFFFUL) << 16); + value = (value >> 32) | (value << 32); + return value; + } + + private static void ReverseBits1024(ArrayView prefix, ArrayView reversed) + { + reversed[0] = ReverseBits64(prefix[15]); + reversed[1] = ReverseBits64(prefix[14]); + reversed[2] = ReverseBits64(prefix[13]); + reversed[3] = ReverseBits64(prefix[12]); + reversed[4] = ReverseBits64(prefix[11]); + reversed[5] = ReverseBits64(prefix[10]); + reversed[6] = ReverseBits64(prefix[9]); + reversed[7] = ReverseBits64(prefix[8]); + reversed[8] = ReverseBits64(prefix[7]); + reversed[9] = ReverseBits64(prefix[6]); + reversed[10] = ReverseBits64(prefix[5]); + reversed[11] = ReverseBits64(prefix[4]); + reversed[12] = ReverseBits64(prefix[3]); + reversed[13] = ReverseBits64(prefix[2]); + reversed[14] = ReverseBits64(prefix[1]); + reversed[15] = ReverseBits64(prefix[0]); + } + + private static void MultiplyPartsGpu(ulong left, ulong right, out ulong high, out ulong low) + { + ulong leftLow = (uint)left; + ulong leftHigh = left >> 32; + ulong rightLow = (uint)right; + ulong rightHigh = right >> 32; + + ulong lowProduct = unchecked(leftLow * rightLow); + ulong cross1 = unchecked(leftHigh * rightLow); + ulong cross2 = unchecked(leftLow * rightHigh); + ulong highProduct = unchecked(leftHigh * rightHigh); + + ulong carry = unchecked((lowProduct >> 32) + (uint)cross1 + (uint)cross2); + low = unchecked((lowProduct & 0xFFFFFFFFUL) | (carry << 32)); + high = unchecked(highProduct + (cross1 >> 32) + (cross2 >> 32) + (carry >> 32)); + } + + private static int GetBitLength(ArrayView words, int wordCount) + { + for (int i = wordCount - 1; i >= 0; i--) + { + ulong word = words[i]; + if (word != 0UL) + { + int leading = XMath.LeadingZeroCount(word); + return (i << 6) + (64 - leading); + } + } + + return 0; + } + + private static bool BuildQOneOffsetsWords( + ArrayView words, + int wordCount, + ArrayView offsets, + out int offsetCount) + { + offsetCount = 0; + for (int wordIndex = 0; wordIndex < wordCount; wordIndex++) + { + ulong word = words[wordIndex]; + while (word != 0UL) + { + int bit = XMath.TrailingZeroCount(word); + if (offsetCount >= MaxQOffsets) + { + return false; + } + + offsets[offsetCount] = (wordIndex << 6) + bit; + offsetCount++; + word &= word - 1UL; + } + } + + return true; + } + + private static void FindBounds(ArrayView values, int valueCount, int lowerTarget, int upperTarget, out int lower, out int upper) + { + int lo = 0; + int hi = valueCount; + int mid; + while (lo < hi) + { + mid = lo + ((hi - lo) >> 1); + if (values[mid] < lowerTarget) + { + lo = mid + 1; + } + else + { + hi = mid; + } + } + + lower = lo; + hi = valueCount; + while (lo < hi) + { + mid = lo + ((hi - lo) >> 1); + if (values[mid] <= upperTarget) + { + lo = mid + 1; + } + else + { + hi = mid; + } + } + + upper = lo; + } + private static bool TryRunHighBitAndBorrowPrefiltersCombinedGpu( + ArrayView qOffsets, + int qOffsetCount, + int pLong, + int maxAllowedA, + ArrayView lo0, + ArrayView hi0, + ArrayView lo1, + ArrayView hi1) + { + int column = pLong - 1; + int columnPlusOne = column + 1; + + int highBitColumns = HighBitCarryPrefilterColumns; + if (highBitColumns > columnPlusOne) + { + highBitColumns = columnPlusOne; + } + + int borrowColumns = 0; + bool runBorrowPrefilter = qOffsetCount >= TopDownBorrowMinUnknown; + if (runBorrowPrefilter) + { + borrowColumns = TopDownBorrowColumns; + if (borrowColumns > columnPlusOne) + { + borrowColumns = columnPlusOne; + } + } + + int sequentialColumns = highBitColumns > borrowColumns ? highBitColumns : borrowColumns; + + var outLo = lo0; + var outHi = hi0; + var inLo = lo1; + var inHi = hi1; + + int outCount = 1; + outLo[0] = 0; + outHi[0] = 0; + + int inCount = 0; + bool flip = false; + CarryRange borrow = CarryRange.Zero; + + int step = 0; + int upperT; + int lowerT = column - maxAllowedA; + if (lowerT < 0) + { + lowerT = 0; + } + + FindBounds(qOffsets, qOffsetCount, lowerT, column, out lowerT, out upperT); + int n = upperT - lowerT; + if (n < 0) + { + return false; + } + + if (runBorrowPrefilter && n < TopDownBorrowMinUnknown) + { + runBorrowPrefilter = false; + if (highBitColumns < sequentialColumns) + { + sequentialColumns = highBitColumns; + } + } + + while (true) + { + if (step < highBitColumns) + { + inCount = 0; + for (int i = 0; i < outCount; i++) + { + long rMin = outLo[i]; + long rMax = outHi[i]; + + if (rMax < 0) continue; + if (rMin < 0) rMin = 0; + if (rMax < rMin) continue; + + const long TwicePlusOneOverflowLimit = (long.MaxValue - 1) >> 1; + if (rMin > TwicePlusOneOverflowLimit || rMax > TwicePlusOneOverflowLimit) + { + return true; + } + + rMin = (rMin << 1) + 1; + rMax = (rMax << 1) + 1; + + if (n <= 0) + { + AddInterval(ref inCount, inLo, inHi, rMin, rMax); + if (inCount == 0) return true; + continue; + } + + if (n == 1) + { + AddInterval(ref inCount, inLo, inHi, rMin - 1, rMax); + if (inCount == 0) return true; + continue; + } + + if (rMax <= n) + { + AddInterval(ref inCount, inLo, inHi, 0, rMax); + if (inCount == 0) return true; + } + else if (rMin > n) + { + AddInterval(ref inCount, inLo, inHi, rMin - n, rMax); + if (inCount == 0) return true; + } + else + { + AddInterval(ref inCount, inLo, inHi, 0, rMax); + if (inCount == 0) return true; + } + } + + if (inCount <= 0) + { + return false; + } + + outCount = inCount; + flip = !flip; + if (!flip) + { + outLo = lo0; outHi = hi0; + inLo = lo1; inHi = hi1; + } + else + { + outLo = lo1; outHi = hi1; + inLo = lo0; inHi = hi0; + } + } + + if (runBorrowPrefilter && step < borrowColumns) + { + if (!TryPropagateCarry(ref borrow, 0, n, 1)) + { + return false; + } + } + + step++; + column--; + if (step >= sequentialColumns || column < 0) + { + break; + } + + lowerT = column - maxAllowedA; + if (lowerT < 0) + { + lowerT = 0; + } + + FindBounds(qOffsets, qOffsetCount, lowerT, column, out lowerT, out upperT); + n = upperT - lowerT; + if (n < 0) + { + return false; + } + } + + return true; + } + + private static void AddInterval(ref int count, ArrayView loArr, ArrayView hiArr, long lo, long hi) + { + if (hi < lo) + { + return; + } + + if (lo < 0) lo = 0; + if (hi < 0) return; + + if (count == 0) + { + loArr[0] = lo; + hiArr[0] = hi; + count = 1; + return; + } + + long lastHi = hiArr[count - 1]; + if (lo <= lastHi + 1) + { + if (hi > lastHi) + { + hiArr[count - 1] = hi; + } + + return; + } + + if (count >= IntervalBufferLength) + { + count = 0; + return; + } + + loArr[count] = lo; + hiArr[count] = hi; + count++; + } + private static void CopyWords(ArrayView source, ArrayView destination, int wordCount) + { + for (int i = 0; i < wordCount; i++) + { + destination[i] = source[i]; + } + } + + private static void ClearWords(ArrayView words, int wordCount) + { + for (int i = 0; i < wordCount; i++) + { + words[i] = 0UL; + } + } + + private static void MulMod2k(ArrayView left, ArrayView right, ArrayView result, int wordCount) + { + ClearWords(result, wordCount); + + for (int i = 0; i < wordCount; i++) + { + ulong carry = 0UL; + ulong carryExtra = 0UL; + ulong leftValue = left[i]; + int maxJ = wordCount - i; + for (int j = 0; j < maxJ; j++) + { + int target = i + j; + ulong high; + ulong low; + MultiplyPartsGpu(leftValue, right[j], out high, out low); + + ulong sum = result[target] + low; + ulong carryOut = sum < low ? 1UL : 0UL; + sum += carry; + if (sum < carry) + { + carryOut++; + } + if (carryExtra != 0UL) + { + sum += 1UL; + if (sum == 0UL) + { + carryOut++; + } + } + + result[target] = sum; + ulong nextCarry = high + carryOut; + carryExtra = nextCarry < high ? 1UL : 0UL; + carry = nextCarry; + } + } + } + + private static void ComputeInverseMod2k( + ArrayView modulus, + ArrayView inverse, + ArrayView tmp, + ArrayView tmp2, + int wordCount) + { + CopyWords(modulus, inverse, wordCount); + + for (int i = 0; i < InverseIterations; i++) + { + MulMod2k(modulus, inverse, tmp, wordCount); + SubtractFromTwo(tmp, tmp2, wordCount); + MulMod2k(inverse, tmp2, tmp, wordCount); + CopyWords(tmp, inverse, wordCount); + } + } + + private static void SubtractFromTwo(ArrayView value, ArrayView result, int wordCount) + { + ulong borrow = 0UL; + for (int i = 0; i < wordCount; i++) + { + ulong minuend = i == 0 ? 2UL : 0UL; + ulong sub = value[i]; + ulong res; + if (borrow == 0UL) + { + res = minuend - sub; + borrow = minuend < sub ? 1UL : 0UL; + } + else + { + res = minuend - sub - 1UL; + borrow = minuend <= sub ? 1UL : 0UL; + } + + result[i] = res; + } + } + + private static void NegateMod2k(ArrayView value, ArrayView result, int wordCount) + { + ulong carry = 1UL; + for (int i = 0; i < wordCount; i++) + { + ulong inverted = ~value[i]; + ulong sum = inverted + carry; + result[i] = sum; + carry = sum < inverted ? 1UL : 0UL; + } + } + private static void BuildStableUnknownDelta8(int unknown, ArrayView delta) + { + for (int low = 0; low < DeltaLength; low++) + { + long value = low; + for (int i = 0; i < DeltaColumnsAtOnce; i++) + { + int sMax = (((unknown ^ (int)value) & 1) != 0) ? unknown : (unknown - 1); + value = (value + sMax - 1) >> 1; + } + + delta[low] = (int)((value << DeltaColumnsAtOnce) - low); + } + } + + private static bool TryAdvanceZeroColumns(ref CarryRange carry, int columnCount) + { + const int Step = 30; + long carryMin = carry.Min; + long carryMax = carry.Max; + int add; + while (columnCount >= Step) + { + add = (1 << Step) - 1; + carryMin = (carryMin + add) >> Step; + carryMax >>= Step; + if (carryMin > carryMax) + { + return false; + } + + columnCount -= Step; + } + + if (columnCount > 0) + { + add = (1 << columnCount) - 1; + carryMin = (carryMin + add) >> columnCount; + carryMax >>= columnCount; + if (carryMin > carryMax) + { + return false; + } + } + + carry = new CarryRange(carryMin, carryMax); + return true; + } + + private static bool TryAdvanceStableUnknownGpu(ref CarryRange carry, int unknown, int columnCount, ArrayView delta) + { + if (columnCount <= 0) + { + return true; + } + + if (unknown <= 0) + { + return TryAdvanceZeroColumns(ref carry, columnCount); + } + + long carryMin = columnCount >= 63 ? 0 : (carry.Min >> columnCount); + long carryMax = carry.Max; + int remaining = columnCount; + + while (remaining > 0) + { + int sMax = (((unknown ^ (int)carryMax) & 1) != 0) ? unknown : unknown - 1; + carryMax = (carryMax + sMax - 1) >> 1; + remaining--; + } + + if (carryMax < carryMin) + { + return false; + } + + carry = new CarryRange(carryMin, carryMax); + return true; + } + + private static bool TryPropagateCarry(ref CarryRange carry, int min, int max, int requiredParity) + { + long minSum = carry.Min + min; + long maxSum = carry.Max + max; + + long alignedMin = AlignUpToParity(minSum, requiredParity); + long alignedMax = AlignDownToParity(maxSum, requiredParity); + if (alignedMin > alignedMax) + { + return false; + } + + carry.Min = (alignedMin - requiredParity) >> 1; + carry.Max = (alignedMax - requiredParity) >> 1; + return true; + } + + private static long AlignUpToParity(long value, int parity) => (value & 1L) == parity ? value : value + 1L; + + private static long AlignDownToParity(long value, int parity) => (value & 1L) == parity ? value : value - 1L; + private static void GetQPrefixMask(int startColumn, int qBitLen, int qWordCount, out int words, out ulong lastMask) + { + if (startColumn >= qBitLen - 1) + { + words = qWordCount; + lastMask = ulong.MaxValue; + return; + } + + if (startColumn < 0) + { + words = 0; + lastMask = 0; + return; + } + + int wordIndex = startColumn >> 6; + words = wordIndex + 1; + + int bit = startColumn & 63; + lastMask = bit == 63 ? ulong.MaxValue : ((1UL << (bit + 1)) - 1UL); + } + + private static int GetAKnownStateRowAware( + int aIndex, + int segmentBase, + int windowSize, + ArrayView knownOne0, + ArrayView knownZero0, + ArrayView knownOne1, + ArrayView knownZero1) + { + if (aIndex < 0) + { + return 2; + } + + int rel; + int word; + ulong mask; + if (aIndex >= segmentBase && aIndex < segmentBase + windowSize) + { + rel = aIndex - segmentBase; + word = rel >> 6; + mask = 1UL << (rel & 63); + if ((knownOne0[word] & mask) != 0UL) return 1; + if ((knownZero0[word] & mask) != 0UL) return 2; + return 0; + } + + int segment1Start = segmentBase - windowSize; + if (aIndex >= segment1Start && aIndex < segmentBase) + { + rel = aIndex - segment1Start; + word = rel >> 6; + mask = 1UL << (rel & 63); + if ((knownOne1[word] & mask) != 0UL) return 1; + if ((knownZero1[word] & mask) != 0UL) return 2; + return 0; + } + + return 0; + } + + private static void SetAKnownRowAware( + int aIndex, + bool isOne, + int segmentBase, + int windowSize, + ArrayView knownOne0, + ArrayView knownZero0, + ArrayView knownOne1, + ArrayView knownZero1) + { + if (aIndex < 0) + { + return; + } + + int rel; + int word; + ulong mask; + if (aIndex >= segmentBase && aIndex < segmentBase + windowSize) + { + rel = aIndex - segmentBase; + word = rel >> 6; + mask = 1UL << (rel & 63); + if (isOne) + { + knownOne0[word] |= mask; + knownZero0[word] &= ~mask; + } + else + { + knownZero0[word] |= mask; + knownOne0[word] &= ~mask; + } + + return; + } + + int segment1Start = segmentBase - windowSize; + if (aIndex >= segment1Start && aIndex < segmentBase) + { + rel = aIndex - segment1Start; + word = rel >> 6; + mask = 1UL << (rel & 63); + if (isOne) + { + knownOne1[word] |= mask; + knownZero1[word] &= ~mask; + } + else + { + knownZero1[word] |= mask; + knownOne1[word] &= ~mask; + } + } + } + + private static void ShiftLeft1InPlace(ArrayView wordsOne, ArrayView wordsZero, int wordCount, ulong lastWordMask) + { + ulong carryOne = 0UL; + ulong carryZero = 0UL; + + for (int i = 0; i < wordCount; i++) + { + ulong word = wordsOne[i]; + ulong newCarry = word >> 63; + wordsOne[i] = (word << 1) | carryOne; + carryOne = newCarry; + + word = wordsZero[i]; + newCarry = word >> 63; + wordsZero[i] = (word << 1) | carryZero; + carryZero = newCarry; + } + + if (wordCount > 0) + { + int last = wordCount - 1; + wordsOne[last] &= lastWordMask; + wordsZero[last] &= lastWordMask; + } + } + + private static void ComputeForcedRemainingAndCleanPrefix( + ArrayView qMask, + ArrayView oneWin, + ArrayView zeroWin, + int wordCount, + ulong lastMask, + out int forced, + out int remaining) + { + forced = 0; + remaining = 0; + + int lastIndex = wordCount - 1; + for (int i = 0; i < wordCount; i++) + { + ulong q = qMask[i]; + if (i == lastIndex) + { + q &= lastMask; + } + + ulong one = oneWin[i] & q; + ulong known = (oneWin[i] | zeroWin[i]); + ulong unknown = q & ~known; + + forced += XMath.PopCount(one); + remaining += XMath.PopCount(unknown); + } + } + + private static bool TryFindSingleUnknownInQPrefix( + ArrayView qMask, + ArrayView oneWin, + ArrayView zeroWin, + int wordCount, + ulong lastMask, + out int offset) + { + offset = -1; + int found = 0; + int lastIndex = wordCount - 1; + for (int w = 0; w < wordCount; w++) + { + ulong q = qMask[w]; + if (w == lastIndex) + { + q &= lastMask; + } + + ulong known = oneWin[w] | zeroWin[w]; + ulong unknown = q & ~known; + if (unknown == 0UL) + { + continue; + } + + found += XMath.PopCount(unknown); + if (found > 1) + { + return false; + } + + int bit = XMath.TrailingZeroCount(unknown); + offset = (w << 6) + bit; + } + + return found == 1; + } + private static bool TryProcessBottomUpBlockRowAwareGpu( + ArrayView qOffsets, + int qOffsetCount, + ArrayView qMaskWords, + int qWordCount, + ulong lastWordMask, + int maxAllowedA, + int startColumn, + int columnCount, + ref CarryRange carry, + ref int maxKnownA, + ref int segmentBase, + ArrayView knownOne0, + ArrayView knownOne1, + ArrayView knownZero0, + ArrayView knownZero1, + int windowSize, + int wordCount, + ArrayView aOneWin, + ArrayView aZeroWin, + ref int windowColumn, + ArrayView stableUnknownDelta, + int stableUnknown) + { + bool needRebuildWindow = false; + int endColumn = startColumn + columnCount; + int endColumnMinusOne = endColumn - 1; + int maxOffset = qOffsets[qOffsetCount - 1]; + int maxOffsetPlusOne = maxOffset + 1; + + if (windowColumn != startColumn) + { + return false; + } + + int qWordCountMinusOne = qWordCount - 1; + + while (startColumn < endColumn) + { + int bound = (startColumn / windowSize) * windowSize; + if (bound != segmentBase) + { + if (bound == segmentBase + windowSize) + { + var temp = knownOne0; + knownOne0 = knownOne1; + knownOne1 = temp; + temp = knownZero0; + knownZero0 = knownZero1; + knownZero1 = temp; + ClearWords(knownOne0, wordCount); + ClearWords(knownZero0, wordCount); + } + else + { + ClearWords(knownOne0, wordCount); + ClearWords(knownOne1, wordCount); + ClearWords(knownZero0, wordCount); + ClearWords(knownZero1, wordCount); + needRebuildWindow = true; + } + + segmentBase = bound; + if (needRebuildWindow) + { + ClearWords(aOneWin, qWordCount); + ClearWords(aZeroWin, qWordCount); + + for (int step = 0; step <= maxOffset; step++) + { + int aIndex = startColumn - step; + int state = aIndex < 0 || aIndex > maxAllowedA + ? 2 + : GetAKnownStateRowAware(aIndex, segmentBase, windowSize, knownOne0, knownZero0, knownOne1, knownZero1); + + if (state == 1) + { + aOneWin[step >> 6] |= 1UL << (step & 63); + } + else if (state == 2) + { + aZeroWin[step >> 6] |= 1UL << (step & 63); + } + } + + if (qWordCount > 0) + { + aOneWin[qWordCountMinusOne] &= lastWordMask; + aZeroWin[qWordCountMinusOne] &= lastWordMask; + } + + needRebuildWindow = false; + } + } + + int qWordsEff = qWordCount; + ulong qLastEff = lastWordMask; + if (startColumn < maxOffset) + { + GetQPrefixMask(startColumn, maxOffsetPlusOne, qWordCount, out qWordsEff, out qLastEff); + } + + if (startColumn >= maxOffset && startColumn <= maxAllowedA) + { + int remaining = endColumnMinusOne; + if (remaining > maxAllowedA) + { + remaining = maxAllowedA; + } + + remaining = remaining - startColumn + 1; + if (remaining > 0) + { + while (remaining >= TailCarryBatchColumns) + { + if (!TryAdvanceStableUnknownGpu(ref carry, stableUnknown, TailCarryBatchColumns, stableUnknownDelta)) + { + return false; + } + + startColumn += TailCarryBatchColumns; + windowColumn = startColumn; + remaining -= TailCarryBatchColumns; + } + + if (remaining > 0) + { + if (!TryAdvanceStableUnknownGpu(ref carry, stableUnknown, remaining, stableUnknownDelta)) + { + return false; + } + + startColumn += remaining; + windowColumn = startColumn; + } + + ClearWords(aOneWin, qWordCount); + ClearWords(aZeroWin, qWordCount); + continue; + } + } + + ComputeForcedRemainingAndCleanPrefix(qMaskWords, aOneWin, aZeroWin, qWordsEff, qLastEff, out int forced, out int remainingCols); + + if (remainingCols == 1 && ((carry.Min ^ carry.Max) & 1L) == 0) + { + int step; + if (TryFindSingleUnknownInQPrefix(qMaskWords, aOneWin, aZeroWin, qWordsEff, qLastEff, out step)) + { + int aIndex = startColumn - step; + const int requiredBit = 1; + int parity = (int)(carry.Min & 1L); + int needOne = (requiredBit ^ parity ^ (forced & 1)) & 1; + + int state = GetAKnownStateRowAware(aIndex, segmentBase, windowSize, knownOne0, knownZero0, knownOne1, knownZero1); + if (needOne != 0) + { + if (state == 2) + { + return false; + } + + SetAKnownRowAware(aIndex, true, segmentBase, windowSize, knownOne0, knownZero0, knownOne1, knownZero1); + int w = step >> 6; + int b = step & 63; + aOneWin[w] |= 1UL << b; + forced++; + } + else + { + if (state == 1) + { + return false; + } + + SetAKnownRowAware(aIndex, false, segmentBase, windowSize, knownOne0, knownZero0, knownOne1, knownZero1); + int w = step >> 6; + int b = step & 63; + aZeroWin[w] |= 1UL << b; + } + + remainingCols = 0; + if (aIndex > maxKnownA) + { + maxKnownA = aIndex; + } + } + } + + int max = forced + remainingCols; + if (!TryPropagateCarry(ref carry, forced, max, 1)) + { + return false; + } + + int nextA = startColumn + 1; + ShiftLeft1InPlace(aOneWin, aZeroWin, qWordCount, lastWordMask); + + int nextState = nextA > maxAllowedA + ? 2 + : GetAKnownStateRowAware(nextA, segmentBase, windowSize, knownOne0, knownZero0, knownOne1, knownZero1); + + if (nextState == 1) + { + aOneWin[0] |= 1UL; + } + else if (nextState == 2) + { + aZeroWin[0] |= 1UL; + } + + startColumn++; + windowColumn = startColumn; + } + + return true; + } +} diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index e2320d9..40c50ff 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -20,12 +20,14 @@ public static class PerfectNumberConstants public const int PooledArrayThreshold = 64; public const int Pow2WindowSize = 8; public const int RollingAccelerators = 298; - public const int ThreadsByAccelerator = 1; + public const int ThreadsByAccelerator = 128; public const int ByDivisorStateSaveInterval = 1_000; public const string ByDivisorStateDirectory = "Checks"; public const int PollardRhoMilliseconds = 100; public const int ByDivisorMersenneTailCount = 8; public const int ByDivisorMersenneBitLimit = 4096; + public const int BitContradictionGpuBatchCount = 16_777_216; + public const int BitContradictionGpuThreads = 4096; // TODO: Load these limits from the benchmark-driven configuration so CPU and GPU scans stay aligned with the optimal // divisor-cycle datasets we generate offline. // TODO: Promote these magic numbers into a runtime profile derived from EvenPerfectBitScanner.Benchmarks so we can retune diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index c1f44ea..31872f0 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -1,6 +1,7 @@ using System.Numerics; using System.Runtime.CompilerServices; using ILGPU.Runtime; +using PerfectNumbers.Core.Gpu; using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core; @@ -157,6 +158,45 @@ public sealed class PrimeTester private static readonly object GpuWarmUpLock = new(); private static int WarmedGpuLeaseCount; + private static int WarmedBitContradictionGpuLeaseCount; + + public static void WarmUpGpuKernelsForBitContradiction(int threadCount) + { + int target = threadCount; + if (target == 0) + { + target = threadCount; + } + + lock (GpuWarmUpLock) + { + PrimeOrderCalculatorAccelerator.SkipOrderTablesAndKernels = true; + + if (target <= WarmedBitContradictionGpuLeaseCount) + { + return; + } + + Accelerator[] accelerators = AcceleratorPool.Shared.Accelerators; + int acceleratorCount = EnvironmentConfiguration.RollingAccelerators; + if (acceleratorCount > accelerators.Length) + { + acceleratorCount = accelerators.Length; + } + for (var i = 0; i < acceleratorCount; i++) + { + Console.WriteLine($"Preparing accelerator {i}..."); + AcceleratorStreamPool.WarmUp(i); + Accelerator accelerator = accelerators[i]; + AcceleratorStream stream = accelerator.CreateStream(); + _ = GpuKernelPool.GetOrAddKernels(i, stream, KernelType.BitContradictionScan); + stream.Synchronize(); + stream.Dispose(); + } + + WarmedBitContradictionGpuLeaseCount = target; + } + } public static void WarmUpGpuKernels(int threadCount) { @@ -173,6 +213,7 @@ public sealed class PrimeTester return; } + PrimeOrderCalculatorAccelerator.SkipOrderTablesAndKernels = false; PrimeOrderCalculatorAccelerator.WarmUp(); WarmedGpuLeaseCount = target; }