BitContradiction GPU implementation fixes

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<ulong, Stride1D.Dense> debugBuffer = gpu.Accelerator.Allocate1D<ulong>(maxThreadCount * debugWordsPerSlot);
+		ArrayView1D<ulong, Stride1D.Dense> 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<ulong> debugQWords = stackalloc ulong[BitContradictionKernels.MaxQWordCount];
+					// Span<ulong> debugQWords = stackalloc ulong[BitContradictionKernels.MaxQWordCount];
-					if (!BitContradictionKernels.TryBuildQWordsFromBatchIndex(prime, batchIndexWords.AsSpan(i * wordCount, wordCount), wordCount, foundIndex, debugQWords, out int debugQBitLen))
+					// 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}");
+					// 	throw new InvalidOperationException($"BitContradiction GPU qWords overflow for p={prime}, k={candidateK}, batch={batchIndex}, index={foundIndex}");
-					}
+					// }
-					BigInteger debugQ = ReadBigIntegerFromWords(debugQWords);
+					// BigInteger debugQ = ReadBigIntegerFromWords(debugQWords);
-					if (debugQ != candidateQ)
+					// 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}");
+					// 	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
+					// #if DivisorSet_BitContradiction
-					int[]? qBits = _qBits;
+					// int[]? qBits = _qBits;
-					if (qBits == null || qBits.Length < (int)prime)
+					// if (qBits == null || qBits.Length < (int)prime)
-					{
+					// {
-						qBits = new int[(int)prime];
+					// 	qBits = new int[(int)prime];
-						_qBits = qBits;
+					// 	_qBits = qBits;
-					}
+					// }
-					bool cpuDecided = TryExactBitContradictionCheck(prime, candidateQ, candidateK, out bool cpuDivides);
+					// bool cpuDecided = TryExactBitContradictionCheck(prime, candidateQ, candidateK, out bool cpuDivides);
-					if (!cpuDecided)
+					// if (!cpuDecided)
-					{
+					// {
-						throw new InvalidOperationException($"BitContradiction CPU undecided for p={prime}, q={candidateQ}, k={candidateK}, batch={batchIndex}");
+					// 	throw new InvalidOperationException($"BitContradiction CPU undecided for p={prime}, q={candidateQ}, k={candidateK}, batch={batchIndex}");
-					}
+					// }
-					if (!cpuDivides)
+					// if (!cpuDivides)
-					{
+					// {
-						throw new InvalidOperationException($"BitContradiction CPU/GPU mismatch for p={prime}, q={candidateQ}, k={candidateK}, batch={batchIndex}");
+					// 	throw new InvalidOperationException($"BitContradiction CPU/GPU mismatch for p={prime}, q={candidateQ}, k={candidateK}, batch={batchIndex}");
-					}
+					// }
-					#endif
+					// #endif
-					bool dividesExact = BigInteger.ModPow(2, prime, candidateQ) == BigInteger.One;
+					// bool dividesExact = BigInteger.ModPow(2, prime, candidateQ) == BigInteger.One;
-					if (!dividesExact)
+					// if (!dividesExact)
-					{
+					// {
-						throw new InvalidOperationException($"BitContradiction CPU/GPU false positive for p={prime}, q={candidateQ}, k={candidateK}, batch={batchIndex}");
+					// 	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;
 	}
 }
+
+
+
+

PerfectNumbers.Core/Gpu/GpuKernelPool.cs

@@ -123,11 +123,11 @@ public class GpuKernelPool
 		{
 			InitOnce(ref kernels.BitContradiction, () =>
 			{
-				var loaded = accelerator.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<int, Stride1D.Dense>>(BitContradictionKernels.BitContradictionKernelScan);
+				var loaded = accelerator.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(BitContradictionKernels.BitContradictionKernelScan);
 
 				var kernel = KernelUtil.GetKernel(loaded);
 
-				return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<int, Stride1D.Dense>>>();
+				return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>>();
 			});
 		}
 

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
+    private struct CarryRange(long min, long max)
-    {
+	{
-        public long Min;
+        public long Min = min;
-        public long Max;
+        public long Max = max;
 
-        public CarryRange(long min, long max)
+		public static CarryRange Zero => new(0, 0);
-        {
-            Min = min;
-            Max = max;
-        }
-
-        public static CarryRange Zero => new(0, 0);
     }
 
     public static void BitContradictionKernelScan(
@@ -44,7 +39,8 @@ internal static class BitContradictionKernels
         ulong exponent,
         ArrayView1D<ulong, Stride1D.Dense> batchIndexWords,
         int countQ,
-        ArrayView1D<int, Stride1D.Dense> foundOut)
+        ArrayView1D<int, Stride1D.Dense> foundOut,
+        ArrayView1D<ulong, Stride1D.Dense> 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(qMaskWords, MaxQWordCount);
-            ClearWords(aOneWin, qWordCount);
+            ClearWords(aOneWin, MaxQWordCount);
-            ClearWords(aZeroWin, qWordCount);
+            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;
             }
         }
     }
+
+	/// <summary>
+	/// This method is used for testing purposes, only
+	/// </summary>
     internal static bool TryBuildQWordsFromBatchIndex(
         ulong exponent,
         ReadOnlySpan<ulong> 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<ulong> aZeroWin,
         ref int windowColumn,
         ArrayView<int> 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;
     }
 }
+
+