using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Algorithms; using ILGPU.Runtime; namespace PerfectNumbers.Core.Gpu; internal enum HeuristicGpuDivisorTableKind : byte { GroupA = 0, GroupBEnding1 = 1, GroupBEnding7 = 7, GroupBEnding9 = 9, Combined = 255, } internal static class PrimeTesterKernels { // GPU kernel: small-prime sieve only. Returns 1 if passes sieve (probable prime), 0 otherwise. public static void SmallPrimeSieveKernel( Index1D index, ArrayView numbers, ArrayView smallPrimes, ArrayView smallPrimesLastOne, ArrayView smallPrimesLastSeven, ArrayView smallPrimesLastThree, ArrayView smallPrimesLastNine, ArrayView smallPrimesPow2, ArrayView smallPrimesPow2LastOne, ArrayView smallPrimesPow2LastSeven, ArrayView smallPrimesPow2LastThree, ArrayView smallPrimesPow2LastNine, ArrayView results) { ulong n = numbers[index]; byte result = 1; ArrayView primes = smallPrimes; ArrayView primeSquares = smallPrimesPow2; ulong lastDigit = n % 10UL; switch (lastDigit) { case 1UL: primes = smallPrimesLastOne; primeSquares = smallPrimesPow2LastOne; break; case 3UL: primes = smallPrimesLastThree; primeSquares = smallPrimesPow2LastThree; break; case 7UL: primes = smallPrimesLastSeven; primeSquares = smallPrimesPow2LastSeven; break; case 9UL: primes = smallPrimesLastNine; primeSquares = smallPrimesPow2LastNine; break; } int length = (int)primes.Length; for (int i = 0; i < length; i++) { ulong primeSquare = primeSquares[i]; if (primeSquare > n) { break; } ulong prime = primes[i]; if (n % prime == 0UL) { result = 0; break; } } results[index] = result; } public static void HeuristicTrialDivisionKernel( Index1D index, ArrayView resultFlag, ulong n, ulong maxDivisorSquare, HeuristicGpuDivisorTableKind tableKind, HeuristicGpuDivisorTables tables) { byte nMod10 = (byte)(n % 10UL); ArrayView1D divisors = tables.SelectDivisors(tableKind, nMod10); ArrayView1D divisorSquares = tables.SelectDivisorSquares(tableKind, nMod10); int divisorLength = (int)divisorSquares.Length; int threadIndex = index; bool skipThread = threadIndex >= divisorLength; ulong divisorSquare = skipThread ? 0UL : divisorSquares[threadIndex]; skipThread = skipThread | (divisorSquare > maxDivisorSquare); ulong divisor = skipThread ? 1UL : divisors[threadIndex]; ulong remainder = skipThread ? 1UL : n % divisor; bool hasFactor = !skipThread & (remainder == 0UL); if (hasFactor) { resultFlag[0] = 1; } } public static void SharesFactorKernel(Index1D index, ArrayView numbers, ArrayView results) { ulong n = numbers[index]; ulong m = 63UL - (ulong)XMath.LeadingZeroCount(n); ulong gcd = BinaryGcdGpu(n, m); results[index] = gcd == 1UL ? (byte)0 : (byte)1; } private static ulong BinaryGcdGpu(ulong u, ulong v) { // TODO: Replace this inline GPU binary GCD with the kernel extracted from // GpuUInt128BinaryGcdBenchmarks via GpuKernelPool so device callers reuse the // fully unrolled ladder instead of this branchy fallback. bool uIsZero = u == 0UL; bool vIsZero = v == 0UL; bool eitherZero = uIsZero | vIsZero; ulong zeroResult = uIsZero ? v : u; ulong combined = u | v | 1UL; int shift = XMath.TrailingZeroCount(combined); ulong normalizedU = u >> XMath.TrailingZeroCount(u | 1UL); ulong normalizedV = v >> XMath.TrailingZeroCount(v | 1UL); ulong currentU = normalizedU; ulong currentV = normalizedV; bool loopCondition = currentV != 0UL; while (loopCondition) { currentV >>= XMath.TrailingZeroCount(currentV); ulong minValue = XMath.Min(currentU, currentV); ulong maxValue = XMath.Max(currentU, currentV); currentU = minValue; currentV = maxValue - minValue; loopCondition = currentV != 0UL; } ulong gcd = currentU << shift; return eitherZero ? zeroResult : gcd; } }