From e30c83584bab4bf3fa7ec781cdd5d40350dfecf0 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Sun, 25 Jan 2026 19:04:56 +0100 Subject: [PATCH] BitContradiction GPU implementation fixes --- ...erDivisorByDivisorCpuTesterWithTemplate.cs | 115 +++++++++++----- PerfectNumbers.Core/Gpu/GpuKernelPool.cs | 4 +- .../Gpu/Kernels/BitContradictionKernels.cs | 129 +++++++++++++++--- 3 files changed, 195 insertions(+), 53 deletions(-) diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs index 1706ce2..6be2ec4 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs @@ -11,6 +11,7 @@ using System.Globalization; using System.Collections.Generic; using System.Collections.Concurrent; using System.Buffers; +using ILGPU; using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; using PerfectNumbers.Core.Gpu.Kernels; @@ -107,6 +108,9 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN int[] found = new int[maxThreadCount]; ulong[] batchIndexWords = new ulong[inputCapacity]; + int debugWordsPerSlot = BitContradictionKernels.DebugWordCountPerSlot; + using MemoryBuffer1D debugBuffer = gpu.Accelerator.Allocate1D(maxThreadCount * debugWordsPerSlot); + ArrayView1D debugView = debugBuffer.View; BigInteger currentBatch = (currentK - BigInteger.One) / batchSizeBig; currentK = (currentBatch * batchSizeBig) + BigInteger.One; @@ -154,7 +158,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN } gpu.InputView.CopyFromCPU(stream, batchIndexWords.AsSpan(0, threadCount * wordCount)); - kernel(stream, threadCount, prime, gpu.InputView, countQ, gpu.OutputIntView); + kernel(stream, threadCount, prime, gpu.InputView, countQ, gpu.OutputIntView, debugView); stream.Synchronize(); gpu.OutputIntView.CopyToCPU(stream, found.AsSpan(0, threadCount)); @@ -173,38 +177,79 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN } 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}"); - } + // 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}"); + if (prime == 138000001UL && candidateQ == new BigInteger(276000003)) + { + ulong[] debugWords = new ulong[debugWordsPerSlot]; + debugView.SubView(i * debugWordsPerSlot, debugWordsPerSlot).CopyToCPU(stream, debugWords); + BigInteger gpuQ = ReadBigIntegerFromWords(debugWords.AsSpan(0, 16)); + BigInteger gpuInv = ReadBigIntegerFromWords(debugWords.AsSpan(16, 16)); + BigInteger gpuALow = ReadBigIntegerFromWords(debugWords.AsSpan(32, 16)); + ulong gpuQBits = debugWords[48]; + ulong gpuQWordCount = debugWords[49]; + ulong gpuQMask0 = debugWords[50]; + ulong gpuQMask1 = debugWords[51]; + ulong gpuKnownOne0 = debugWords[52]; + ulong gpuKnownZero0 = debugWords[53]; + ulong gpuAOneWin0 = debugWords[54]; + ulong gpuAZeroWin0 = debugWords[55]; + long gpuCarryMin = unchecked((long)debugWords[56]); + long gpuCarryMax = unchecked((long)debugWords[57]); + long gpuMaxAllowedA = unchecked((long)debugWords[58]); + long gpuMaxFixed = unchecked((long)debugWords[59]); + long gpuWindowSize = unchecked((long)debugWords[60]); + long gpuAWordCount = unchecked((long)debugWords[61]); + long gpuMaxKnownA = unchecked((long)debugWords[62]); + long gpuMaxFixedPrefilter = unchecked((long)debugWords[63]); + ulong preKnownOne0 = debugWords[64]; + ulong preKnownZero0 = debugWords[65]; + ulong preAOneWin0 = debugWords[66]; + ulong preAZeroWin0 = debugWords[67]; + long preCarryMin = unchecked((long)debugWords[68]); + long preCarryMax = unchecked((long)debugWords[69]); + long preSegmentBase = unchecked((long)debugWords[70]); + long preWindowColumn = unchecked((long)debugWords[71]); + long gpuOffsetCount = unchecked((long)debugWords[72]); + long gpuMaxOffset = unchecked((long)debugWords[73]); + long gpuFirstForced = unchecked((long)debugWords[74]); + long gpuFirstRemaining = unchecked((long)debugWords[75]); + long gpuFirstQWordsEff = unchecked((long)debugWords[76]); + ulong gpuFirstQLastEff = debugWords[77]; + long gpuFirstCarryMin = unchecked((long)debugWords[78]); + long gpuFirstCarryMax = unchecked((long)debugWords[79]); + throw new InvalidOperationException($"BitContradiction GPU debug for p={prime}, q={candidateQ}, k={candidateK}, qBits={gpuQBits}, qWordCount={gpuQWordCount}, qMask0=0x{gpuQMask0:X16}, qMask1=0x{gpuQMask1:X16}, knownOne0=0x{gpuKnownOne0:X16}, knownZero0=0x{gpuKnownZero0:X16}, aOneWin0=0x{gpuAOneWin0:X16}, aZeroWin0=0x{gpuAZeroWin0:X16}, carryMin={gpuCarryMin}, carryMax={gpuCarryMax}, maxAllowedA={gpuMaxAllowedA}, maxFixed={gpuMaxFixed}, windowSize={gpuWindowSize}, aWordCount={gpuAWordCount}, maxKnownA={gpuMaxKnownA}, maxFixedPrefilter={gpuMaxFixedPrefilter}, preKnownOne0=0x{preKnownOne0:X16}, preKnownZero0=0x{preKnownZero0:X16}, preAOneWin0=0x{preAOneWin0:X16}, preAZeroWin0=0x{preAZeroWin0:X16}, preCarryMin={preCarryMin}, preCarryMax={preCarryMax}, preSegmentBase={preSegmentBase}, preWindowColumn={preWindowColumn}, offsetCount={gpuOffsetCount}, maxOffset={gpuMaxOffset}, firstForced={gpuFirstForced}, firstRemaining={gpuFirstRemaining}, firstQWordsEff={gpuFirstQWordsEff}, firstQLastEff=0x{gpuFirstQLastEff:X16}, firstCarryMin={gpuFirstCarryMin}, firstCarryMax={gpuFirstCarryMax}, gpuQ={gpuQ}, gpuInv={gpuInv}, gpuALow={gpuALow}"); + } + foundDivisor = candidateQ; AcceleratorStreamPool.Return(gpu.AcceleratorIndex, stream); return true; @@ -5201,3 +5246,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN return false; } } + + + + diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index 52a6aa0..d6aec2b 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -123,11 +123,11 @@ public class GpuKernelPool { InitOnce(ref kernels.BitContradiction, () => { - var loaded = accelerator.LoadAutoGroupedStreamKernel, int, ArrayView1D>(BitContradictionKernels.BitContradictionKernelScan); + var loaded = accelerator.LoadAutoGroupedStreamKernel, int, ArrayView1D, ArrayView1D>(BitContradictionKernels.BitContradictionKernelScan); var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, int, ArrayView1D>>(); + return kernel.CreateLauncherDelegate, int, ArrayView1D, ArrayView1D>>(); }); } diff --git a/PerfectNumbers.Core/Gpu/Kernels/BitContradictionKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/BitContradictionKernels.cs index a55a6f0..799babb 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/BitContradictionKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/BitContradictionKernels.cs @@ -21,22 +21,17 @@ internal static class BitContradictionKernels private const int TopDownBorrowColumns = 4096; private const int TopDownBorrowMinUnknown = 12; private const int TailCarryBatchColumns = 16384; - private const int DeltaColumnsAtOnce = 8; + private const int DeltaColumnsAtOnce = 16; private const int DeltaLength = 1 << DeltaColumnsAtOnce; + public const int DebugWordCountPerSlot = 80; private const int IntervalBufferLength = 64; - private struct CarryRange - { - public long Min; - public long Max; + private struct CarryRange(long min, long max) + { + public long Min = min; + public long Max = max; - public CarryRange(long min, long max) - { - Min = min; - Max = max; - } - - public static CarryRange Zero => new(0, 0); + public static CarryRange Zero => new(0, 0); } public static void BitContradictionKernelScan( @@ -44,7 +39,8 @@ internal static class BitContradictionKernels ulong exponent, ArrayView1D batchIndexWords, int countQ, - ArrayView1D foundOut) + ArrayView1D foundOut, + ArrayView1D debugOut) { int slot = index; foundOut[slot] = -1; @@ -142,6 +138,7 @@ internal static class BitContradictionKernels int qWordCount = (qBitLen + 63) >> 6; int offsetCount; + int maxOffsetValue = 0; if (!BuildQOneOffsetsWords(qWords, qWordCount, qOffsets, out offsetCount)) { continue; @@ -149,6 +146,7 @@ internal static class BitContradictionKernels if (offsetCount <= 0) { + maxOffsetValue = qOffsets[offsetCount - 1]; continue; } @@ -160,9 +158,9 @@ internal static class BitContradictionKernels ClearWords(knownOne1, aWordCount); ClearWords(knownZero0, aWordCount); ClearWords(knownZero1, aWordCount); - ClearWords(qMaskWords, qWordCount); - ClearWords(aOneWin, qWordCount); - ClearWords(aZeroWin, qWordCount); + ClearWords(qMaskWords, MaxQWordCount); + ClearWords(aOneWin, MaxQWordCount); + ClearWords(aZeroWin, MaxQWordCount); ClearWords(qPrefix, ForcedALowWords); CopyWords(qWords, qPrefix, ForcedALowWords); @@ -203,6 +201,12 @@ internal static class BitContradictionKernels } int maxKnownA = maxFixed >= 0 ? maxFixed : -1; + int firstForced = -1; + int firstRemaining = -1; + int firstQWordsEff = 0; + ulong firstQLastEff = 0; + long firstCarryMin = 0; + long firstCarryMax = 0; int maxFixedPrefilter = 512 - 1; if (maxFixedPrefilter > maxKnownA) @@ -304,6 +308,14 @@ internal static class BitContradictionKernels BuildStableUnknownDelta8(offsetCount, delta8); CarryRange carry = CarryRange.Zero; + ulong preKnownOne0 = knownOne0[0]; + ulong preKnownZero0 = knownZero0[0]; + ulong preAOneWin0 = aOneWin[0]; + ulong preAZeroWin0 = aZeroWin[0]; + long preCarryMin = 0L; + long preCarryMax = 0L; + int preSegmentBase = segmentBase; + int preWindowColumn = windowColumn; bool ok = TryProcessBottomUpBlockRowAwareGpu( qOffsets, offsetCount, @@ -326,15 +338,71 @@ internal static class BitContradictionKernels aZeroWin, ref windowColumn, delta8, - offsetCount); + offsetCount, + ref firstForced, + ref firstRemaining, + ref firstQWordsEff, + ref firstQLastEff, + ref firstCarryMin, + ref firstCarryMax); if (ok) { + int debugBase = slot * DebugWordCountPerSlot; + for (int w = 0; w < ForcedALowWords; w++) + { + debugOut[debugBase + w] = qWords[w]; + } + for (int w = 0; w < ForcedALowWords; w++) + { + debugOut[debugBase + ForcedALowWords + w] = invWords[w]; + } + for (int w = 0; w < ForcedALowWords; w++) + { + debugOut[debugBase + (ForcedALowWords * 2) + w] = aLowWords[w]; + } + debugOut[debugBase + (ForcedALowWords * 3)] = (ulong)qBitLen; + int extraBase = debugBase + (ForcedALowWords * 3) + 1; + debugOut[extraBase] = (ulong)qWordCount; + debugOut[extraBase + 1] = qMaskWords[0]; + debugOut[extraBase + 2] = qMaskWords[1]; + debugOut[extraBase + 3] = knownOne0[0]; + debugOut[extraBase + 4] = knownZero0[0]; + debugOut[extraBase + 5] = aOneWin[0]; + debugOut[extraBase + 6] = aZeroWin[0]; + debugOut[extraBase + 7] = (ulong)carry.Min; + debugOut[extraBase + 8] = (ulong)carry.Max; + debugOut[extraBase + 9] = (ulong)maxAllowedA; + debugOut[extraBase + 10] = (ulong)maxFixed; + debugOut[extraBase + 11] = (ulong)windowSize; + debugOut[extraBase + 12] = (ulong)aWordCount; + debugOut[extraBase + 13] = (ulong)maxKnownA; + debugOut[extraBase + 14] = (ulong)maxFixedPrefilter; + debugOut[extraBase + 15] = preKnownOne0; + debugOut[extraBase + 16] = preKnownZero0; + debugOut[extraBase + 17] = preAOneWin0; + debugOut[extraBase + 18] = preAZeroWin0; + debugOut[extraBase + 19] = (ulong)preCarryMin; + debugOut[extraBase + 20] = (ulong)preCarryMax; + debugOut[extraBase + 21] = (ulong)preSegmentBase; + debugOut[extraBase + 22] = (ulong)preWindowColumn; + debugOut[extraBase + 23] = (ulong)offsetCount; + debugOut[extraBase + 24] = (ulong)maxOffsetValue; + debugOut[extraBase + 25] = (ulong)firstForced; + debugOut[extraBase + 26] = (ulong)firstRemaining; + debugOut[extraBase + 27] = (ulong)firstQWordsEff; + debugOut[extraBase + 28] = firstQLastEff; + debugOut[extraBase + 29] = (ulong)firstCarryMin; + debugOut[extraBase + 30] = (ulong)firstCarryMax; foundOut[slot] = batch; return; } } } + + /// + /// This method is used for testing purposes, only + /// internal static bool TryBuildQWordsFromBatchIndex( ulong exponent, ReadOnlySpan batchIndexWords, @@ -1010,6 +1078,14 @@ internal static class BitContradictionKernels long carryMax = carry.Max; int remaining = columnCount; + const long CarryDecimalMask = (1L << DeltaColumnsAtOnce) - 1L; + while (remaining >= DeltaColumnsAtOnce) + { + int carryMask = (int)(carryMax & CarryDecimalMask); + carryMax = (carryMax + delta[carryMask]) >> DeltaColumnsAtOnce; + remaining -= DeltaColumnsAtOnce; + } + while (remaining > 0) { int sMax = (((unknown ^ (int)carryMax) & 1) != 0) ? unknown : unknown - 1; @@ -1282,7 +1358,13 @@ internal static class BitContradictionKernels ArrayView aZeroWin, ref int windowColumn, ArrayView stableUnknownDelta, - int stableUnknown) + int stableUnknown, + ref int firstForced, + ref int firstRemaining, + ref int firstQWordsEff, + ref ulong firstQLastEff, + ref long firstCarryMin, + ref long firstCarryMax) { bool needRebuildWindow = false; int endColumn = startColumn + columnCount; @@ -1451,6 +1533,15 @@ internal static class BitContradictionKernels int max = forced + remainingCols; if (!TryPropagateCarry(ref carry, forced, max, 1)) + if (startColumn == 0 && firstForced < 0) + { + firstForced = forced; + firstRemaining = remainingCols; + firstQWordsEff = qWordsEff; + firstQLastEff = qLastEff; + firstCarryMin = carry.Min; + firstCarryMax = carry.Max; + } { return false; } @@ -1478,3 +1569,5 @@ internal static class BitContradictionKernels return true; } } + +