Merge pull request #72 from mlemanczyk/codex/implement-todos-in-mersennenumberdivisorbydivisorgputester

Fix residue tester filters and arithmetic utilities

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -26,6 +26,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
 
         private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<MontgomeryDivisorData>, ArrayView<byte>>> _kernelCache = new();
         private readonly ConcurrentDictionary<Accelerator, Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<byte>>> _kernelByPrimeCache = new();
+        private readonly ConcurrentBag<DivisorScanSession> _sessionPool = new();
 
 	private Action<Index1D, ulong, ArrayView<MontgomeryDivisorData>, ArrayView<byte>> GetKernel(Accelerator accelerator) =>
 									_kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<MontgomeryDivisorData>, ArrayView<byte>>(CheckKernel));
@@ -55,29 +56,77 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
 	}
 
 	public bool IsPrime(ulong prime, out bool divisorsExhausted)
+	{
+                ulong allowedMax;
+                bool useCycles;
+                int batchCapacity;
+
+                lock (_sync)
                 {
                         if (!_isConfigured)
                         {
                                 throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
                         }
 
-		ulong allowedMax = ComputeAllowedMaxDivisor(prime);
+                        allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit);
+                        useCycles = _useDivisorCycles;
+                        batchCapacity = _gpuBatchSize;
+                }
+
                 if (allowedMax < 3UL)
                 {
                         divisorsExhausted = true;
                         return true;
                 }
 
-		bool composite, coveredRange;
+                bool composite;
+                bool coveredRange;
+                ulong processedCount;
+                ulong lastProcessed;
 
-		// TODO: Let's rework this to check if the number is prime and check divisors outside of the lock. Lock and update the state afterwards.
+                var gpuLease = GpuContextPool.RentPreferred(preferCpu: false);
+                var accelerator = gpuLease.Accelerator;
+                var kernel = GetKernel(accelerator);
+
+                var divisorsBuffer = accelerator.Allocate1D<MontgomeryDivisorData>(batchCapacity);
+                var hitsBuffer = accelerator.Allocate1D<byte>(batchCapacity);
+                ulong[] divisors = ArrayPool<ulong>.Shared.Rent(batchCapacity);
+                byte[] hits = ArrayPool<byte>.Shared.Rent(batchCapacity);
+                MontgomeryDivisorData[] divisorData = ArrayPool<MontgomeryDivisorData>.Shared.Rent(batchCapacity);
+
+                try
+                {
+                        composite = CheckDivisors(
+                                prime,
+                                allowedMax,
+                                useCycles,
+                                accelerator,
+                                kernel,
+                                divisorsBuffer,
+                                hitsBuffer,
+                                divisors,
+                                hits,
+                                divisorData,
+                                out lastProcessed,
+                                out coveredRange,
+                                out processedCount
+                        );
+                }
+                finally
+                {
+                        ArrayPool<byte>.Shared.Return(hits, clearArray: true);
+                        ArrayPool<ulong>.Shared.Return(divisors, clearArray: true);
+                        ArrayPool<MontgomeryDivisorData>.Shared.Return(divisorData, clearArray: true);
+                        hitsBuffer.Dispose();
+                        divisorsBuffer.Dispose();
+                        gpuLease.Dispose();
+                }
+
+                if (processedCount > 0UL)
+                {
                         lock (_sync)
                         {
-			composite = CheckDivisors(prime, allowedMax, out ulong lastProcessed, out coveredRange);
+                                UpdateStatusUnsafe(lastProcessed, processedCount);
-
-			if (lastProcessed != 0UL)
-			{
-				ReportStatus(lastProcessed);
                         }
                 }
 
@@ -91,11 +140,23 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
 		return true;
 	}
 
-	private bool CheckDivisors(ulong prime, ulong allowedMax, out ulong lastProcessed, out bool coveredRange)
+        private bool CheckDivisors(
+                ulong prime,
+                ulong allowedMax,
+                bool useCycles,
+                Accelerator accelerator,
+                Action<Index1D, ulong, ArrayView<MontgomeryDivisorData>, ArrayView<byte>> kernel,
+                MemoryBuffer1D<MontgomeryDivisorData, Stride1D.Dense> divisorsBuffer,
+                MemoryBuffer1D<byte, Stride1D.Dense> hitsBuffer,
+                ulong[] divisors,
+                byte[] hits,
+                MontgomeryDivisorData[] divisorData,
+                out ulong lastProcessed,
+                out bool coveredRange,
+                out ulong processedCount)
         {
-		// TODO: We should be able to calculate allowedMax outside in the parent method, unless it depends from some changing value.
-		allowedMax = Math.Min(allowedMax, _divisorLimit);
                 lastProcessed = 0UL;
+                processedCount = 0UL;
 
                 if (allowedMax < 3UL)
                 {
@@ -103,24 +164,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
                         return false;
                 }
 
-		// TODO: Rent these in the parent method before the loops - batch size doesn't change
+                int batchCapacity = (int)divisorsBuffer.Length;
-		var gpuLease = GpuContextPool.RentPreferred(preferCpu: false);
-		var accelerator = gpuLease.Accelerator;
-		var kernel = GetKernel(accelerator);
-
-		int batchCapacity = _gpuBatchSize;
-		// TODO: Rent buffers in the parent method before the loops - batch size doesn't change
-		var divisorsBuffer = accelerator.Allocate1D<MontgomeryDivisorData>(batchCapacity);
-		var hitsBuffer = accelerator.Allocate1D<byte>(batchCapacity);
-		ulong[] divisors = ArrayPool<ulong>.Shared.Rent(batchCapacity);
-		byte[] hits = ArrayPool<byte>.Shared.Rent(batchCapacity);
-		MontgomeryDivisorData[] divisorData = ArrayPool<MontgomeryDivisorData>.Shared.Rent(batchCapacity);
 
                 bool composite = false;
                 bool processedAll = false;
                 ulong cycle, divisor = 3UL, next;
                 int batchSize, i;
-		bool reachedEndInBatch, useCycles = _useDivisorCycles;
+                bool reachedEndInBatch;
 
                 Span<MontgomeryDivisorData> divisorDataSpan;
                 Span<ulong> divisorSpan;
@@ -136,24 +186,38 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
 
                         while (batchSize < batchCapacity && divisor <= allowedMax)
                         {
+                                ulong currentDivisor = divisor;
+                                ulong nextDivisor = currentDivisor + 2UL;
+                                processedCount++;
+
                                 if (useCycles)
                                 {
-					cycle = divisorCycles!.GetCycle(divisor);
+                                        cycle = divisorCycles!.GetCycle(currentDivisor);
                                         if (cycle == 0UL || prime % cycle != 0UL)
                                         {
-						ReportStatus(divisor);
+                                                lastProcessed = currentDivisor;
-					}
+
-					else
+                                                if (nextDivisor <= currentDivisor)
                                                 {
-						divisors[batchSize++] = divisor;
+                                                        reachedEndInBatch = true;
-					}
+                                                        break;
-				}
-				else
-				{
-					divisors[batchSize++] = divisor;
                                                 }
 
-				next = divisor + 2UL;
+                                                if (nextDivisor > allowedMax)
+                                                {
+                                                        divisor = nextDivisor;
+                                                        reachedEndInBatch = true;
+                                                        break;
+                                                }
+
+                                                divisor = nextDivisor;
+                                                continue;
+                                        }
+                                }
+
+                                divisors[batchSize++] = currentDivisor;
+                                lastProcessed = currentDivisor;
+                                next = nextDivisor;
 
                                 // This will be only true when we exceed ulong.MaxValue
                                 if (next <= divisor)
@@ -174,6 +238,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
 
                         if (batchSize == 0)
                         {
+                                if (reachedEndInBatch)
+                                {
+                                        processedAll = true;
+                                }
+
                                 break;
                         }
 
@@ -199,11 +268,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
                                 if (hitsSpan[i] != 0)
                                 {
                                         composite = true;
+                                        lastProcessed = divisorSpan[i];
                                         break;
                                 }
                         }
 
-			ReportStatus(divisorSpan[i < batchSize ? i : i - 1]);
+                        if (!composite)
+                        {
+                                lastProcessed = divisorSpan[batchSize - 1];
+                        }
 			if (composite)
 			{
 				break;
@@ -218,25 +291,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
 
                 coveredRange = composite || processedAll || divisor > allowedMax;
 
-		ArrayPool<byte>.Shared.Return(hits, clearArray: true);
-		ArrayPool<ulong>.Shared.Return(divisors, clearArray: true);
-		ArrayPool<MontgomeryDivisorData>.Shared.Return(divisorData, clearArray: true);
-		hitsBuffer.Dispose();
-		divisorsBuffer.Dispose();
-		gpuLease.Dispose();
-
                 return composite;
         }
 
-	private void ReportStatus(ulong divisor)
-	{
-		if (++_lastStatusDivisor == PerfectNumberConstants.ConsoleInterval)
-		{
-			Console.WriteLine($"...processed by-divisor candidate = {divisor}");
-			_lastStatusDivisor = 0;
-		}
-	}
-
 	private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
 	{
 		if (maxPrime <= 1UL)
@@ -252,7 +309,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
 		return (1UL << (int)(maxPrime - 1UL)) - 1UL;
 	}
 
-	private ulong ComputeAllowedMaxDivisor(ulong prime)
+        private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit)
         {
                 if (prime <= 1UL)
                 {
@@ -261,10 +318,34 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
 
                 if (prime - 1UL >= 64UL)
                 {
-			return _divisorLimit;
+                        return divisorLimit;
                 }
 
-		return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, _divisorLimit);
+                return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
+        }
+
+        private void UpdateStatusUnsafe(ulong lastProcessed, ulong processedCount)
+        {
+                if (processedCount == 0UL)
+                {
+                        return;
+                }
+
+                ulong interval = (ulong)PerfectNumberConstants.ConsoleInterval;
+                if (interval == 0UL)
+                {
+                        _lastStatusDivisor = 0UL;
+                        return;
+                }
+
+                ulong total = _lastStatusDivisor + processedCount;
+                ulong messages = total / interval;
+                _lastStatusDivisor = total % interval;
+
+                for (ulong i = 0; i < messages; i++)
+                {
+                        Console.WriteLine($"...processed by-divisor candidate = {lastProcessed}");
+                }
         }
 
 	private static void CheckKernel(Index1D index, ulong prime, ArrayView<MontgomeryDivisorData> divisors, ArrayView<byte> hits)
@@ -280,7 +361,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
 		hits[index] = Pow2Mod(prime, divisor) == 1UL ? (byte)1 : (byte)0;
 	}
 
-	// TODO: Let's create a pool for these and let's rent them and return, when not needed
         public sealed class DivisorScanSession : IDisposable
         {
                 private readonly MersenneNumberDivisorByDivisorGpuTester _owner;
@@ -307,6 +387,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
                         _indexHost = ArrayPool<int>.Shared.Rent(owner._gpuBatchSize);
                 }
 
+                internal void Reset()
+                {
+                        _disposed = false;
+                }
+
                 public void CheckDivisor(ulong divisor, ReadOnlySpan<ulong> primes, Span<byte> hits)
                 {
                         if (_disposed)
@@ -411,15 +496,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
                         }
 
                         _disposed = true;
-			ArrayPool<ulong>.Shared.Return(_primesHost, clearArray: true);
+                        _owner.ReturnSession(this);
-			ArrayPool<byte>.Shared.Return(_hitsHost, clearArray: true);
-			ArrayPool<int>.Shared.Return(_indexHost, clearArray: true);
-			_primesBuffer.Dispose();
-			_hitsBuffer.Dispose();
-			_lease.Dispose();
                 }
         }
 
+        internal void ReturnSession(DivisorScanSession session)
+        {
+                _sessionPool.Add(session);
+        }
+
 	private static ulong Pow2Mod(ulong exponent, in MontgomeryDivisorData divisor)
 	{
 		ulong modulus = divisor.Modulus;
@@ -544,7 +629,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
                                 throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
                         }
 
-			return ComputeAllowedMaxDivisor(prime);
+                        return ComputeAllowedMaxDivisor(prime, _divisorLimit);
                 }
         }
 
@@ -557,6 +642,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
                                 throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
                         }
 
+                        if (_sessionPool.TryTake(out DivisorScanSession? session))
+                        {
+                                session.Reset();
+                                return session;
+                        }
+
                         return new DivisorScanSession(this);
                 }
         }

PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs

@@ -22,12 +22,11 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
 		ResiduePrimeViews primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
 		int batchSize = GpuConstants.ScanBatchSize;
                 UInt128 kStart = 1UL;
+                UInt128 limit = maxK + UInt128.One;
 		byte last = lastIsSeven ? (byte)1 : (byte)0;
 		var kernel = gpuLease.Pow2ModKernel;
-		ulong step3 = exponent % 3UL,
+                twoP.Mod10_8_5_3(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
-			  step5 = exponent % 5UL,
+                step10 = step10.Mod10();
-			  step8 = exponent % 8UL,
-			  step10 = exponent.Mod10();
 
 		GpuUInt128 twoPGpu = (GpuUInt128)twoP;
 
@@ -46,9 +45,9 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
 
 		try
 		{
-			while (kStart < maxK && Volatile.Read(ref isPrime))
+                        while (kStart < limit && Volatile.Read(ref isPrime))
                         {
-				UInt128 remaining = maxK - kStart;
+                                UInt128 remaining = limit - kStart;
 				int currentSize = remaining > batchSize128 ? batchSize : (int)remaining;
 				Span<ulong> orders = orderArray.AsSpan(0, currentSize);
 				ref ulong ordersRef = ref MemoryMarshal.GetReference(orders);
@@ -84,7 +83,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
 
 				UInt128 processed = (UInt128)currentSize;
                                 kStart += processed;
-				if (kStart < maxK)
+                                if (kStart < limit)
                                 {
                                         if (processed == batchSize128)
                                         {

PerfectNumbers.Core/UInt128Extensions.cs

@@ -148,10 +148,16 @@ public static class UInt128Extensions
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static ulong Mod7(this UInt128 value)
         {
-		ulong remainder = (((ulong)value) % 7UL) + (((ulong)(value >> 64)) % 7UL) << 1;
+                ulong low = (ulong)value % 7UL;
-		while (remainder >= 7UL)
+                ulong high = (ulong)(value >> 64) % 7UL;
+                ulong remainder = low + (high * 2UL);
+                if (remainder >= 7UL)
                 {
                         remainder -= 7UL;
+                        if (remainder >= 7UL)
+                        {
+                                remainder -= 7UL;
+                        }
                 }
 
                 return remainder;
@@ -164,8 +170,9 @@ public static class UInt128Extensions
         public static ulong Mod10(this UInt128 value128)
         {
                 // Split and fold under mod 10: 2^64 ≡ 6 (mod 10)
-		ulong value = (ulong)value128;
+                ulong low = (ulong)value128;
-		return ((value % 10UL) + ((value >> 64) % 10UL) * 6UL) % 10UL;
+                ulong high = (ulong)(value128 >> 64);
+                return ((low % 10UL) + ((high % 10UL) * 6UL)) % 10UL;
         }
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]

PerfectNumbers.Core/ULongExtensions.cs

@@ -48,13 +48,14 @@ public static class ULongExtensions
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static ulong FastDiv64(this ulong value, ulong divisor, ulong mul)
         {
-		mul = value.MulHigh(mul);
+                ulong quotient = (ulong)(((UInt128)value * mul) >> 64);
-		if ((UInt128)value - mul.Mul64(divisor) >= divisor)
+                UInt128 remainder = (UInt128)value - ((UInt128)quotient * divisor);
+                if (remainder >= divisor)
                 {
-			mul++;
+                        quotient++;
                 }
 
-		return mul;
+                return quotient;
         }
 
 	public const ulong WordBitMask = 0xFFFFUL;