GPU calculations disabled again

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -546,7 +546,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         }
 
         ulong exponent = exponents[index];
-        hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
+        hits[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
     }
 
     public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -765,7 +765,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
         ulong exponent = exponents[index];
 
-        results[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true);
+        results[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: true);
     }
 
     public ulong DivisorLimit

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

@@ -809,7 +809,7 @@ internal static partial class PrimeOrderGpuHeuristics
             }
 
             ulong reduced = phi / factor;
-            if (reduced.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL)
+            if (ULongExtensions.Pow2MontgomeryModWindowedKernel(reduced, divisor, keepMontgomery: false) == 1UL)
             {
                 return false;
             }
@@ -1593,7 +1593,7 @@ internal static partial class PrimeOrderGpuHeuristics
 
     private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
     {
-        return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL;
+        return ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL;
     }
 
     private static ulong CalculateByDoublingKernel(ulong prime)
@@ -1793,7 +1793,7 @@ internal static partial class PrimeOrderGpuHeuristics
             return 0UL;
         }
 
-        return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false);
+        return ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisorData, keepMontgomery: false);
     }
 
     private static Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
@@ -1809,7 +1809,7 @@ internal static partial class PrimeOrderGpuHeuristics
     private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
     {
         ulong exponent = exponents[index];
-        remainders[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false);
+        remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false);
     }
 
     private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)

PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs

@@ -317,7 +317,7 @@ internal static partial class PrimeOrderCalculator
 				bool gpuStackRemainders = false;
 				GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable;
 
-				// if (allowGpuBatch && IsGpuPow2Allowed)
+				if (allowGpuBatch && IsGpuPow2Allowed)
 				{
 					if (batchSize <= StackGpuBatchSize)
 					{
@@ -638,21 +638,21 @@ internal static partial class PrimeOrderCalculator
 
 	private static bool Pow2EqualsOneCpu(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
 	{
-		// if (IsGpuPow2Allowed)
-		// {
-			ulong remainder = exponent.Pow2MontgomeryModWindowedGpu(divisorData, false);
-			// GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData);
-			// if (status == GpuPow2ModStatus.Success)
-			// {
+		if (IsGpuPow2Allowed)
+		{
+			// ulong remainder = exponent.Pow2MontgomeryModWindowedGpu(divisorData, false);
+			GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData);
+			if (status == GpuPow2ModStatus.Success)
+			{
 				return remainder == 1UL;
-			// }
-		// }
+			}
+		}
 
-		// return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
+		return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
 	}
 
 	// [ThreadStatic]
-	private static readonly Dictionary<ulong, bool> _primalityCache = new(16_777_216);
+	private static readonly Dictionary<ulong, bool> _primalityCache = new(134_217_728);
 	private static readonly ReaderWriterLockSlim _primalityCacheLock = new(LockRecursionPolicy.NoRecursion);
 
 	// [ThreadStatic]

PerfectNumbers.Core/ULongExtensions.cs

@@ -1,6 +1,9 @@
+using System.Collections.Concurrent;
 using System.Numerics;
 using System.Runtime.CompilerServices;
+using ILGPU;
 using ILGPU.Algorithms;
+using ILGPU.Runtime;
 using PerfectNumbers.Core.Gpu;
 
 namespace PerfectNumbers.Core;
@@ -238,20 +241,20 @@ public static class ULongExtensions
 		// faster, while GPU code can keep using GpuUInt128.MulHigh.
 		ulong xLow = (uint)x;
 		ulong xHigh = x >> 32;
-        ulong yLow = (uint)y;
+		ulong result = (uint)y;
 		ulong yHigh = y >> 32;
 
 		ulong w1 = xLow * yHigh;
-        ulong w2 = xHigh * yLow;
-        ulong w3 = xLow * yLow;
+		ulong w2 = xHigh * result;
+		xLow *= result;
 
 		// Matching the layout used in GpuUInt128.MulHigh: introducing the
 		// intermediate result looks like one extra store, but it lets RyuJIT keep
 		// the accumulated high word entirely in registers. Without this explicit
 		// local the JIT spills the partial sum, which is where the performance
 		// regression in the benchmarks came from.
-        ulong result = (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32);
-        result += ((w3 >> 32) + (uint)w1 + (uint)w2) >> 32;
+		result = (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32);
+		result += ((xLow >> 32) + (uint)w1 + (uint)w2) >> 32;
 		return result;
 	}
 
@@ -313,9 +316,9 @@ public static class ULongExtensions
 	{
 		ulong tLow = unchecked(a * b);
 		ulong m = unchecked(tLow * nPrime);
-        ulong mTimesModulusLow = unchecked(m * modulus);
+		ulong result = unchecked(m * modulus);
 
-        ulong result = unchecked(a.MulHigh(b) + m.MulHigh(modulus) + (unchecked(tLow + mTimesModulusLow) < tLow ? 1UL : 0UL));
+		result = unchecked(a.MulHigh(b) + m.MulHigh(modulus) + (unchecked(tLow + result) < tLow ? 1UL : 0UL));
 		if (result >= modulus)
 		{
 			result -= modulus;
@@ -409,13 +412,33 @@ public static class ULongExtensions
 			return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery);
 		}
 
+		if (Pow2MontgomeryGpuExecutor.TryExecute(exponent, divisor, keepMontgomery, out ulong gpuResult))
+		{
+			return gpuResult;
+		}
+
+		return Pow2MontgomeryModWindowedCpu(exponent, divisor, keepMontgomery);
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	internal static ulong Pow2MontgomeryModWindowedKernel(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
+	{
+		ulong modulus = divisor.Modulus;
+		// This will never happen in production code
+		// if (exponent == 0UL)
+		// {
+		//     return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus;
+		// }
+
+		if (exponent <= Pow2WindowFallbackThreshold)
+		{
+			return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery);
+		}
+
 		int bitLength = GetPortableBitLength(exponent);
 		int windowSize = GetWindowSize(bitLength);
-        int oddPowerCount = 1 << (windowSize - 1);
-
 		ulong result = divisor.MontgomeryOne;
 		ulong nPrime = divisor.NPrime;
-        ulong[] oddPowers = InitializeMontgomeryOddPowersGpu(divisor, modulus, nPrime, oddPowerCount);
 
 		int index = bitLength - 1;
 		while (index >= 0)
@@ -428,37 +451,122 @@ public static class ULongExtensions
 			}
 
 			int windowStart = index - windowSize + 1;
-            if (windowStart < 0)
-            {
-                windowStart = 0;
-            }
+			// This can never happen in production code
+			// if (windowStart < 0)
+			// {
+			//     windowStart = 0;
+			// }
 
 			while (((exponent >> windowStart) & 1UL) == 0UL)
 			{
 				windowStart++;
 			}
 
-            int windowLength = index - windowStart + 1;
-            for (int square = 0; square < windowLength; square++)
+			bitLength = index - windowStart + 1;
+			for (int square = 0; square < bitLength; square++)
 			{
 				result = result.MontgomeryMultiply(result, modulus, nPrime);
 			}
 
-            ulong mask = (1UL << windowLength) - 1UL;
-            ulong windowValue = (exponent >> windowStart) & mask;
-            int tableIndex = (int)((windowValue - 1UL) >> 1);
-            ulong multiplier = oddPowers[tableIndex];
-            result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
+			ulong windowValue = (1UL << bitLength) - 1UL;
+			windowValue = (exponent >> windowStart) & windowValue;
+			windowValue = ComputeMontgomeryOddPower(windowValue, divisor, modulus, nPrime);
+			result = result.MontgomeryMultiply(windowValue, modulus, nPrime);
 
 			index = windowStart - 1;
 		}
 
-        if (keepMontgomery)
+		if (!keepMontgomery)
 		{
+			result = result.MontgomeryMultiply(1UL, modulus, nPrime);
+		}
+
 		return result;
 	}
 
-        return result.MontgomeryMultiply(1UL, modulus, nPrime);
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private static ulong ComputeMontgomeryOddPower(ulong exponent, in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime)
+	{
+		ulong baseValue = divisor.MontgomeryTwo;
+		ulong power = divisor.MontgomeryOne;
+		ulong remaining = exponent;
+
+		while (remaining != 0UL)
+		{
+			if ((remaining & 1UL) != 0UL)
+			{
+				power = power.MontgomeryMultiply(baseValue, modulus, nPrime);
+			}
+
+			remaining >>= 1;
+			if (remaining == 0UL)
+			{
+				break;
+			}
+
+			baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime);
+		}
+
+		return power;
+	}
+
+	private static class Pow2MontgomeryGpuExecutor
+	{
+		private static readonly ConcurrentDictionary<Accelerator, Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, byte, ArrayView1D<ulong, Stride1D.Dense>>> KernelCache = new();
+
+		public static bool TryExecute(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery, out ulong result)
+		{
+			result = 0UL;
+
+			GpuKernelLease lease;
+			try
+			{
+				lease = GpuKernelPool.GetKernel(useGpuOrder: true);
+			}
+			catch (Exception)
+			{
+				return false;
+			}
+
+			var execution = lease.EnterExecutionScope();
+			Accelerator accelerator = lease.Accelerator;
+			if (accelerator.AcceleratorType == AcceleratorType.CPU)
+			{
+				return false;
+			}
+
+			var kernel = KernelCache.GetOrAdd(accelerator, static accel =>
+			{
+				var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, byte, ArrayView1D<ulong, Stride1D.Dense>>(Pow2MontgomeryKernel);
+				var launcher = KernelUtil.GetKernel(loaded);
+				return launcher.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, byte, ArrayView1D<ulong, Stride1D.Dense>>>();
+			});
+
+			var exponentBuffer = accelerator.Allocate1D<ulong>(1);
+			var resultBuffer = accelerator.Allocate1D<ulong>(1);
+
+			exponentBuffer.View.CopyFromCPU(ref exponent, 1);
+			// We could clear the buffer, but it'll be overwritten anyway
+			// resultBuffer.MemSetToZero();
+
+			byte keepFlag = keepMontgomery ? (byte)1 : (byte)0;
+			AcceleratorStream stream = lease.Stream;
+			kernel(stream, 1, exponentBuffer.View, divisor, keepFlag, resultBuffer.View);
+			stream.Synchronize();
+
+			resultBuffer.View.CopyToCPU(ref result, 1);
+			resultBuffer.Dispose();
+			exponentBuffer.Dispose();
+			execution.Dispose();
+			lease.Dispose();
+
+			return true;
+		}
+
+		private static void Pow2MontgomeryKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, byte keepMontgomery, ArrayView1D<ulong, Stride1D.Dense> results)
+		{
+			results[index] = Pow2MontgomeryModWindowedKernel(exponents[index], divisor, keepMontgomery != 0);
+		}
 	}
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -492,43 +600,25 @@ public static class ULongExtensions
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	private static int GetPortableBitLength(ulong value)
 	{
-        if (value == 0UL)
-        {
-            return 0;
-        }
+		// This will never happen in production code
+		// if (value == 0UL)
+		// {
+		//     return 0;
+		// }
 
 		return 64 - XMath.LeadingZeroCount(value);
 	}
 
-    private static int GetWindowSize(int bitLength)
+	private static int GetWindowSize(int bitLength) => bitLength switch
 	{
-        if (bitLength <= 6)
-        {
-            return Math.Max(bitLength, 1);
-        }
-
-        if (bitLength <= 23)
-        {
-            return 4;
-        }
-
-        if (bitLength <= 79)
-        {
-            return 5;
-        }
-
-        if (bitLength <= 239)
-        {
-            return 6;
-        }
-
-        if (bitLength <= 671)
-        {
-            return 7;
-        }
-
-        return Pow2WindowSize;
-    }
+		<= 6 => bitLength,// bitLength is always >= 1 in production code
+						  // return Math.Max(bitLength, 1);
+		<= 23 => 4,
+		<= 79 => 5,
+		<= 239 => 6,
+		<= 671 => 7,
+		_ => Pow2WindowSize,
+	};
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span<ulong> oddPowers)
@@ -549,26 +639,6 @@ public static class ULongExtensions
 		}
 	}
 
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong[] InitializeMontgomeryOddPowersGpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
-    {
-        ulong[] oddPowers = new ulong[PerfectNumberConstants.MaxOddPowersCount];
-        oddPowers[0] = divisor.MontgomeryTwo;
-        if (oddPowerCount == 1)
-        {
-            return oddPowers;
-        }
-
-        ulong square = divisor.MontgomeryTwoSquared;
-        for (int i = 1; i < oddPowerCount; i++)
-        {
-            ulong previous = oddPowers[i - 1];
-            oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime);
-        }
-
-        return oddPowers;
-    }
-
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
 	{