diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks.cs index f26a60c..de961cf 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks.cs @@ -1,6 +1,7 @@ using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Jobs; using PerfectNumbers.Core; +using PerfectNumbers.Core.Gpu.Kernels; namespace EvenPerfectBitScanner.Benchmarks; diff --git a/PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs b/PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs index 89c6e38..f3ce19f 100644 --- a/PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/ULongExtensionsGpuTests.cs @@ -1,5 +1,6 @@ using FluentAssertions; using Xunit; +using PerfectNumbers.Core.Gpu.Kernels; namespace PerfectNumbers.Core.Tests.Gpu; diff --git a/PerfectNumbers.Core/Gpu/Accelerators/HeuristicCombinedGpuViews.cs b/PerfectNumbers.Core/Gpu/Accelerators/HeuristicCombinedGpuViews.cs index 65d70fd..a526214 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/HeuristicCombinedGpuViews.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/HeuristicCombinedGpuViews.cs @@ -14,32 +14,52 @@ public readonly struct HeuristicCombinedGpuViews( ArrayView1D combinedDivisorSquaresEnding7, ArrayView1D combinedDivisorSquaresEnding9) { - public readonly ArrayView1D CombinedDivisorsEnding1 = combinedDivisorsEnding1; - public readonly ArrayView1D CombinedDivisorsEnding3 = combinedDivisorsEnding3; - public readonly ArrayView1D CombinedDivisorsEnding7 = combinedDivisorsEnding7; - public readonly ArrayView1D CombinedDivisorsEnding9 = combinedDivisorsEnding9; - public readonly ArrayView1D CombinedDivisorSquaresEnding1 = combinedDivisorSquaresEnding1; - public readonly ArrayView1D CombinedDivisorSquaresEnding3 = combinedDivisorSquaresEnding3; - public readonly ArrayView1D CombinedDivisorSquaresEnding7 = combinedDivisorSquaresEnding7; - public readonly ArrayView1D CombinedDivisorSquaresEnding9 = combinedDivisorSquaresEnding9; + // public readonly ArrayView1D CombinedDivisorsEnding1 = combinedDivisorsEnding1; + // public readonly ArrayView1D CombinedDivisorsEnding3 = combinedDivisorsEnding3; + // public readonly ArrayView1D CombinedDivisorsEnding7 = combinedDivisorsEnding7; + // public readonly ArrayView1D CombinedDivisorsEnding9 = combinedDivisorsEnding9; + public readonly (ArrayView1D, ArrayView1D) CombinedDivisorsAndSquaresEnding1 = (combinedDivisorsEnding1, combinedDivisorSquaresEnding1); + + public readonly (ArrayView1D, ArrayView1D) CombinedDivisorsAndSquaresEnding3 = (combinedDivisorsEnding3, combinedDivisorSquaresEnding3); + + public readonly (ArrayView1D, ArrayView1D) CombinedDivisorsAndSquaresEnding7 = (combinedDivisorsEnding7, combinedDivisorSquaresEnding7); + + public readonly (ArrayView1D, ArrayView1D) CombinedDivisorsAndSquaresEnding9 = (combinedDivisorsEnding9, combinedDivisorSquaresEnding9); + + // public readonly ArrayView1D CombinedDivisorSquaresEnding1 = combinedDivisorSquaresEnding1; + // public readonly ArrayView1D CombinedDivisorSquaresEnding3 = combinedDivisorSquaresEnding3; + // public readonly ArrayView1D CombinedDivisorSquaresEnding7 = combinedDivisorSquaresEnding7; + // public readonly ArrayView1D CombinedDivisorSquaresEnding9 = combinedDivisorSquaresEnding9; + + private static readonly (ArrayView1D, ArrayView1D)EmptyCombinedDivisorsAndSquares = (ArrayView1D.Empty, ArrayView1D.Empty); + + // [MethodImpl(MethodImplOptions.AggressiveInlining)] + // public ArrayView1D SelectDivisors(byte ending) => ending switch + // { + // 1 => CombinedDivisorsEnding1, + // 3 => CombinedDivisorsEnding3, + // 7 => CombinedDivisorsEnding7, + // 9 => CombinedDivisorsEnding9, + // _ => ArrayView1D.Empty, + // }; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public ArrayView1D SelectDivisors(byte ending) => ending switch + public (ArrayView1D, ArrayView1D) SelectDivisorsAndSquares(byte ending) => ending switch { - 1 => CombinedDivisorsEnding1, - 3 => CombinedDivisorsEnding3, - 7 => CombinedDivisorsEnding7, - 9 => CombinedDivisorsEnding9, - _ => ArrayView1D.Empty, + 1 => CombinedDivisorsAndSquaresEnding1, + 3 => CombinedDivisorsAndSquaresEnding3, + 7 => CombinedDivisorsAndSquaresEnding7, + 9 => CombinedDivisorsAndSquaresEnding9, + _ => EmptyCombinedDivisorsAndSquares, }; - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public ArrayView1D SelectDivisorSquares(byte ending) => ending switch - { - 1 => CombinedDivisorSquaresEnding1, - 3 => CombinedDivisorSquaresEnding3, - 7 => CombinedDivisorSquaresEnding7, - 9 => CombinedDivisorSquaresEnding9, - _ => ArrayView1D.Empty, - }; + // [MethodImpl(MethodImplOptions.AggressiveInlining)] + // public ArrayView1D SelectDivisorSquares(byte ending) => ending switch + // { + // 1 => CombinedDivisorSquaresEnding1, + // 3 => CombinedDivisorSquaresEnding3, + // 7 => CombinedDivisorSquaresEnding7, + // 9 => CombinedDivisorSquaresEnding9, + // _ => ArrayView1D.Empty, + // }; } diff --git a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs index 71e8a71..1ff5958 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs @@ -78,7 +78,7 @@ public sealed class PrimeOrderCalculatorAccelerator // GpuStaticTableInitializer.EnsureStaticTables(accelerator, kernels, stream); stream.Synchronize(); stream.Dispose(); - accelerator.Synchronize(); + GC.Collect(GC.MaxGeneration, GCCollectionMode.Aggressive | GCCollectionMode.Forced, blocking: true, compacting: true); } } diff --git a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs index 0f74f40..2cc7926 100644 --- a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs +++ b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs @@ -1,4 +1,5 @@ using System.Runtime.CompilerServices; +using PerfectNumbers.Core.Gpu.Kernels; namespace PerfectNumbers.Core.Gpu; @@ -92,7 +93,7 @@ internal struct ExponentRemainderStepperGpu [MethodImpl(MethodImplOptions.AggressiveInlining)] private void InitializeStateGpu(ulong exponent) { - ulong residue = ULongExtensions.Pow2ModWindowedGpuKernel(exponent, _modulus); + ulong residue = exponent.Pow2ModWindowedGpuKernel(_modulus); _currentResidue = new GpuUInt128(residue); PreviousExponent = exponent; } @@ -101,7 +102,7 @@ internal struct ExponentRemainderStepperGpu private void AdvanceStateGpu(ulong exponent) { ulong delta = exponent - PreviousExponent; - ulong multiplier = ULongExtensions.Pow2ModWindowedGpuKernel(delta, _modulus); + ulong multiplier = delta.Pow2ModWindowedGpuKernel(_modulus); _currentResidue.MulMod(multiplier, _modulusWide); PreviousExponent = exponent; } diff --git a/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs index 9e114a1..3b09b99 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/Pow2MontgomeryKernels.cs @@ -1,10 +1,11 @@ using System.Runtime.CompilerServices; using ILGPU; +using ILGPU.Algorithms; using ILGPU.Runtime; namespace PerfectNumbers.Core.Gpu.Kernels; -internal static class Pow2MontgomeryKernels +public static class Pow2MontgomeryKernels { [MethodImpl(MethodImplOptions.AggressiveInlining)] private static GpuUInt128 AdvancePolynomialGpu(GpuUInt128 x, ulong c, ulong modulus) @@ -14,6 +15,355 @@ internal static class Pow2MontgomeryKernels return x; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ComputeMontgomeryOddPowerGpu(ulong exponent, in MontgomeryDivisorDataGpu divisor, ulong modulus, ulong nPrime) + { + ulong baseValue = divisor.MontgomeryTwo; + ulong power = divisor.MontgomeryOne; + ulong remaining = exponent; + + while (remaining != 0UL) + { + ulong bit = remaining & 1UL; + ulong multiplied = MontgomeryMultiplyGpu(power, baseValue, modulus, nPrime); + ulong mask = (ulong)-(long)bit; + power = (power & ~mask) | (multiplied & mask); + + remaining >>= 1; + if (remaining == 0UL) + { + break; + } + + baseValue = MontgomeryMultiplyGpu(baseValue, baseValue, modulus, nPrime); + } + + return power; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetPortableBitLengthGpu(ulong value) => 64 - XMath.LeadingZeroCount(value); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetWindowSizeGpu(int n) + { + int w = PerfectNumberConstants.Pow2WindowSize; + + w = (n <= 671) ? 7 : w; + w = (n <= 239) ? 6 : w; + w = (n <= 79) ? 5 : w; + w = (n <= 23) ? 4 : w; + + int m = (n >= 1) ? n : 1; + w = (n <= 6) ? m : w; + + return w; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetNextSetBitIndexGpu(ulong exponent, int startIndex) + { + ulong guard = (ulong)(((long)startIndex - 64) >> 63); + int shift = startIndex & 63; + ulong mask = (~0UL << shift) & guard; + ulong masked = exponent & mask; + + return XMath.TrailingZeroCount(masked); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void InitializeStandardOddPowers(ulong modulus, Span oddPowers) + { + // oddPowers spans at least one entry for production workloads; keep the guard commented out to avoid redundant checks in the hot loop. + // if (oddPowers.IsEmpty) + // { + // return; + // } + + // EvenPerfectBitScanner feeds odd prime moduli (≥ 3) here, so the base value stays within range without a modulo reduction. + ulong baseValue = 2UL; + oddPowers[0] = baseValue; + if (oddPowers.Length == 1) + { + return; + } + + ulong square = baseValue.MulModGpu(baseValue, modulus); + for (int i = 1; i < oddPowers.Length; i++) + { + ulong previous = oddPowers[i - 1]; + oddPowers[i] = previous.MulModGpu(square, modulus); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong MontgomeryMultiplyGpu(this ulong a, ulong b, ulong modulus, ulong nPrime) + { + // Mirror the struct-based MontMul64 helper measured as the second-fastest GPU-compatible + // option in GpuUInt128Montgomery64Benchmarks so accelerator kernels reuse the validated path. + MultiplyPartsGpu(a, b, out ulong productHigh, out ulong productLow); + MultiplyPartsGpu(productLow, nPrime, out _, out ulong mLow); + MultiplyPartsGpu(mLow, modulus, out ulong mTimesModulusHigh, out ulong mTimesModulusLow); + + ulong sumLow = unchecked(productLow + mTimesModulusLow); + ulong carry = sumLow < productLow ? 1UL : 0UL; + ulong sumHigh = unchecked(productHigh + mTimesModulusHigh + carry); + + if (sumHigh >= modulus) + { + sumHigh -= modulus; + } + + return sumHigh; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void MultiplyPartsGpu(ulong left, ulong right, out ulong high, out ulong low) + { + ulong leftLow = (uint)left; + ulong leftHigh = left >> 32; + ulong rightLow = (uint)right; + ulong rightHigh = right >> 32; + + ulong lowProduct = unchecked(leftLow * rightLow); + ulong cross1 = unchecked(leftHigh * rightLow); + ulong cross2 = unchecked(leftLow * rightHigh); + ulong highProduct = unchecked(leftHigh * rightHigh); + + ulong carry = unchecked((lowProduct >> 32) + (uint)cross1 + (uint)cross2); + low = unchecked((lowProduct & 0xFFFFFFFFUL) | (carry << 32)); + high = unchecked(highProduct + (cross1 >> 32) + (cross2 >> 32) + (carry >> 32)); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong Pow2MontgomeryModWindowedGpuKeepMontgomery(in MontgomeryDivisorDataGpu divisor, ulong exponent) + { + return Pow2MontgomeryModWindowedGpuMontgomeryResult(divisor, exponent); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModWindowedGpuConvertToStandard(this MontgomeryDivisorDataGpu divisor, ulong exponent) + { + ulong modulus = divisor.Modulus; + ulong nPrime = divisor.NPrime; + ulong montgomeryResult = Pow2MontgomeryModWindowedGpuMontgomeryResult(divisor, exponent); + return MontgomeryMultiplyGpu(montgomeryResult, 1UL, modulus, nPrime); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong Pow2MontgomeryModWindowedGpuMontgomeryResult(in MontgomeryDivisorDataGpu divisor, ulong exponent) + { + ulong modulus = divisor.Modulus; + // Reuse the GPU-compatible MontgomeryMultiply helper highlighted in MontgomeryMultiplyBenchmarks so the accelerator path + // mirrors the fastest CPU reduction without relying on UInt128 intrinsics. + // We shouldn't hit it in production code. If we do, we're doing something wrong. + // if (exponent == 0UL) + // { + // return divisor.MontgomeryOne; + // } + + // We barely ever hit it in production code. It's e.g. 772 calls out of billions + // if (exponent <= Pow2WindowFallbackThreshold) + // { + // return Pow2MontgomeryModSingleBit(exponent, divisor); + // } + + int bitLength = GetPortableBitLengthGpu(exponent); + int windowSize = GetWindowSizeGpu(bitLength); + ulong result = divisor.MontgomeryOne; + ulong nPrime = divisor.NPrime; + + int index = bitLength - 1; + while (index >= 0) + { + ulong currentBit = (exponent >> index) & 1UL; + ulong squared = MontgomeryMultiplyGpu(result, result, modulus, nPrime); + bool processWindow = currentBit != 0UL; + + result = processWindow ? result : squared; + index -= (int)(currentBit ^ 1UL); + + int windowStartCandidate = index - windowSize + 1; + // if (windowStartCandidate < 0) + // { + // windowStartCandidate = 0; + // } + // Clamp the negative offset without branching so the GPU loop stays divergence-free. + // This still handles windows that would otherwise extend past the most significant bit. + int negativeMask = windowStartCandidate >> 31; + windowStartCandidate &= ~negativeMask; + + int windowStart = processWindow ? GetNextSetBitIndexGpu(exponent, windowStartCandidate) : windowStartCandidate; + int windowLength = processWindow ? index - windowStart + 1 : 0; + for (int square = 0; square < windowLength; square++) + { + result = MontgomeryMultiplyGpu(result, result, modulus, nPrime); + } + + result = processWindow + ? MontgomeryMultiplyGpu( + result, + ComputeMontgomeryOddPowerGpu( + (exponent >> windowStart) & ((1UL << windowLength) - 1UL), + divisor, + modulus, + nPrime), + modulus, + nPrime) + : result; + + index = processWindow ? windowStart - 1 : index; + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus) + { + // EvenPerfectBitScanner never routes moduli ≤ 1 through this path; keep the guard commented out so the hot loop stays branch-free. + // if (modulus <= 1UL) + // { + // return 0UL; + // } + + // EvenPerfectBitScanner always works with positive exponents here; leave the defensive check commented for future reference. + // if (exponent == 0UL) + // { + // return 1UL % modulus; + // } + + int bitLength = GetPortableBitLengthGpu(exponent); + int windowSize = GetWindowSizeGpu(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + // oddPowerCount is always ≥ 1 for the production workloads handled by EvenPerfectBitScanner. + // if (oddPowerCount <= 0) + // { + // return 1UL; + // } + + ulong[] oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount); + Span oddPowers = oddPowersArray.AsSpan(0, oddPowerCount); + + InitializeStandardOddPowers(modulus, oddPowers); + + ulong result = 1UL; + int index = bitLength - 1; + + while (index >= 0) + { + ulong currentBit = (exponent >> index) & 1UL; + ulong squared = result.MulModGpu(result, modulus); + result = currentBit == 0UL ? squared : result; + index = currentBit == 0UL ? index - 1 : index; + if (currentBit == 0UL) + { + continue; + } + + int windowStart = index - windowSize + 1; + // EvenPerfectBitScanner guarantees windowStart stays non-negative for production exponents, so keep the guard commented out to avoid extra branching. + // if (windowStart < 0) + // { + // windowStart = 0; + // } + + windowStart = GetNextSetBitIndexGpu(exponent, windowStart); + + int windowLength = index - windowStart + 1; + for (int square = 0; square < windowLength; square++) + { + result = result.MulModGpu(result, modulus); + } + + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + int tableIndex = (int)((windowValue - 1UL) >> 1); + ulong multiplier = oddPowers[tableIndex]; + result = result.MulModGpu(multiplier, modulus); + + index = windowStart - 1; + } + + ThreadStaticPools.UlongPool.Return(oddPowersArray); + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2ModWindowedGpuKernel(this ulong exponent, ulong modulus) + { + int bitLength = GetPortableBitLengthGpu(exponent); + int windowSize = GetWindowSizeGpu(bitLength); + ReadOnlyGpuUInt128 modulusWide = new(modulus); + GpuUInt128 result = new(1UL); + int index = bitLength - 1; + + while (index >= 0) + { + ulong currentBit = (exponent >> index) & 1UL; + if (currentBit == 0UL) + { + result.MulMod(result.AsReadOnly(), modulusWide); + index--; + continue; + } + + int windowStartCandidate = index - windowSize + 1; + int negativeMask = windowStartCandidate >> 31; + windowStartCandidate &= ~negativeMask; + int windowStart = GetNextSetBitIndexGpu(exponent, windowStartCandidate); + int windowLength = index - windowStart + 1; + + for (int square = 0; square < windowLength; square++) + { + result.MulMod(result.AsReadOnly(), modulusWide); + } + + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + ulong multiplier = ComputeWindowedOddPowerGpuKernel(windowValue, modulusWide); + result.MulMod(multiplier, modulusWide); + + index = windowStart - 1; + } + + return result.Low; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ComputeWindowedOddPowerGpuKernel(ulong windowValue, in ReadOnlyGpuUInt128 modulus) + { + GpuUInt128 result = new(2UL); + if (windowValue == 1UL) + { + // Single-bit windows show up whenever the bit scan finds an isolated exponent bit, + // which happens routinely on the GPU path, so keep the fast return enabled. + return 2UL; + } + + ulong remaining = (windowValue - 1UL) >> 1; + GpuUInt128 squareBase = new(2UL); + squareBase.MulMod(squareBase.AsReadOnly(), modulus); + + while (remaining != 0UL) + { + if ((remaining & 1UL) != 0UL) + { + result.MulMod(squareBase.AsReadOnly(), modulus); + } + + remaining >>= 1; + if (remaining != 0UL) + { + squareBase.MulMod(squareBase.AsReadOnly(), modulus); + } + } + + return result.Low; + } + public static void Pow2MontgomeryKernelKeepMontgomery( Index1D index, ArrayView1D exponents, @@ -25,7 +375,7 @@ internal static class Pow2MontgomeryKernels ArrayView1D results) { MontgomeryDivisorDataGpu divisor = new(divisorModulus, divisorNPrime, divisorMontgomeryOne, divisorMontgomeryTwo, divisorMontgomeryTwoSquared); - results[0] = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(divisor, exponents[0]); + results[0] = Pow2MontgomeryModWindowedGpuKeepMontgomery(divisor, exponents[0]); } public static void Pow2MontgomeryKernelConvertToStandard( @@ -39,7 +389,7 @@ internal static class Pow2MontgomeryKernels ArrayView1D results) { MontgomeryDivisorDataGpu divisor = new(divisorModulus, divisorNPrime, divisorMontgomeryOne, divisorMontgomeryTwo, divisorMontgomeryTwoSquared); - results[0] = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponents[0]); + results[0] = Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponents[0]); } diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs index 45e69a5..19c3ba1 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderKernels.cs @@ -2,6 +2,7 @@ using System.Numerics; using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Runtime; +using PerfectNumbers.Core.Gpu.Kernels; namespace PerfectNumbers.Core.Gpu; @@ -38,7 +39,7 @@ internal static partial class PrimeOrderGpuHeuristics } // ulong candidate = order / primeFactor; - if (ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisorData, candidate) == 1UL) + if (divisorData.Pow2MontgomeryModWindowedGpuConvertToStandard(candidate) == 1UL) { order = candidate; continue; @@ -403,7 +404,7 @@ internal static partial class PrimeOrderGpuHeuristics } ulong reduced = phi / factor; - if (ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, reduced) == 1UL) + if (divisor.Pow2MontgomeryModWindowedGpuConvertToStandard(reduced) == 1UL) { return false; } @@ -1187,7 +1188,7 @@ internal static partial class PrimeOrderGpuHeuristics private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorDataGpu divisor) { - return ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent) == 1UL; + return divisor.Pow2MontgomeryModWindowedGpuConvertToStandard(exponent) == 1UL; } private static ulong CalculateByDoublingKernel(ulong prime) @@ -1215,7 +1216,7 @@ internal static partial class PrimeOrderGpuHeuristics { ulong exponent = exponents[index]; MontgomeryDivisorDataGpu divisor = new(divisorModulus, divisorNPrime, divisorMontgomeryOne, divisorMontgomeryTwo, divisorMontgomeryTwoSquared); - remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent); + remainders[index] = divisor.Pow2MontgomeryModWindowedGpuConvertToStandard(exponent); } /// This kernel always sets the result of the corresponding element. Callers don't need to clear the output buffers. diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs index 3c03c6c..0ede5dd 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs @@ -86,28 +86,17 @@ internal static class PrimeTesterKernels ulong maxDivisorSquare, HeuristicCombinedGpuViews tables) { - // byte nMod10 = (byte)(n % 10UL); - ArrayView1D divisors = tables.SelectDivisors(nMod10); - ArrayView1D divisorSquares = tables.SelectDivisorSquares(nMod10); + var (divisors, divisorSquares) = tables.SelectDivisorsAndSquares(nMod10); - int divisorLength = (int)divisorSquares.Length; int threadIndex = index; - if (threadIndex >= divisorLength) + if (divisorSquares[threadIndex] > maxDivisorSquare) { return; } - ulong divisorSquare = divisorSquares[threadIndex]; - if (divisorSquare > maxDivisorSquare) + if (n % divisors[threadIndex] == 0UL) { - return; - } - - ulong divisor = divisors[threadIndex]; - divisor = n % divisor; - if (divisor == 0UL) - { - resultFlag[0] = 1; + Atomic.Exchange(ref resultFlag[0], 1); } } diff --git a/PerfectNumbers.Core/Gpu/MontgomeryDivisorDataGpu.cs b/PerfectNumbers.Core/Gpu/MontgomeryDivisorDataGpu.cs index fcb4573..b384908 100644 --- a/PerfectNumbers.Core/Gpu/MontgomeryDivisorDataGpu.cs +++ b/PerfectNumbers.Core/Gpu/MontgomeryDivisorDataGpu.cs @@ -8,3 +8,4 @@ public readonly struct MontgomeryDivisorDataGpu(ulong modulus, ulong nPrime, ulo public readonly ulong MontgomeryTwo = montgomeryTwo; public readonly ulong MontgomeryTwoSquared = montgomeryTwoSquared; } + \ No newline at end of file diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index ddd453d..7f2d36f 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -5,6 +5,8 @@ using ILGPU; using ILGPU.Runtime; using ILGPU.Runtime.OpenCL; using PerfectNumbers.Core.Gpu.Accelerators; +using PerfectNumbers.Core.Gpu.Kernels; + namespace PerfectNumbers.Core.Gpu; internal enum GpuPow2ModStatus @@ -648,7 +650,7 @@ internal static partial class PrimeOrderGpuHeuristics return 0UL; } - return ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisorData, exponent); + return divisorData.Pow2MontgomeryModWindowedGpuConvertToStandard(exponent); } private static Kernel GetPow2ModKernel(int acceleratorIndex) diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index 1f6743f..2033992 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -13,6 +13,7 @@ public static class PerfectNumberConstants public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; public const int MaxOddPowersCount = 128; public const int PooledArrayThreshold = 64; + public const int Pow2WindowSize = 8; public static readonly int RollingAccelerators = 298; //SharedGpuContext.Device.MaxNumThreadsPerGroup; public const int ThreadsByAccelerator = 1; // TODO: Load these limits from the benchmark-driven configuration so CPU and GPU scans stay aligned with the optimal diff --git a/PerfectNumbers.Core/ULongExtensions.Cpu.cs b/PerfectNumbers.Core/ULongExtensions.Cpu.cs index f0e70c5..0e14ffe 100644 --- a/PerfectNumbers.Core/ULongExtensions.Cpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Cpu.cs @@ -261,7 +261,7 @@ public static partial class ULongExtensions return 7; } - return Pow2WindowSize; + return PerfectNumberConstants.Pow2WindowSize; } } diff --git a/PerfectNumbers.Core/ULongExtensions.Gpu.cs b/PerfectNumbers.Core/ULongExtensions.Gpu.cs index ba698cf..b509f1f 100644 --- a/PerfectNumbers.Core/ULongExtensions.Gpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Gpu.cs @@ -43,274 +43,6 @@ public static partial class ULongExtensions return state.MulMod(right, modulus); } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong MontgomeryMultiplyGpu(this ulong a, ulong b, ulong modulus, ulong nPrime) - { - // Mirror the struct-based MontMul64 helper measured as the second-fastest GPU-compatible - // option in GpuUInt128Montgomery64Benchmarks so accelerator kernels reuse the validated path. - MultiplyPartsGpu(a, b, out ulong productHigh, out ulong productLow); - MultiplyPartsGpu(productLow, nPrime, out _, out ulong mLow); - MultiplyPartsGpu(mLow, modulus, out ulong mTimesModulusHigh, out ulong mTimesModulusLow); - - ulong sumLow = unchecked(productLow + mTimesModulusLow); - ulong carry = sumLow < productLow ? 1UL : 0UL; - ulong sumHigh = unchecked(productHigh + mTimesModulusHigh + carry); - - if (sumHigh >= modulus) - { - sumHigh -= modulus; - } - - return sumHigh; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void MultiplyPartsGpu(ulong left, ulong right, out ulong high, out ulong low) - { - ulong leftLow = (uint)left; - ulong leftHigh = left >> 32; - ulong rightLow = (uint)right; - ulong rightHigh = right >> 32; - - ulong lowProduct = unchecked(leftLow * rightLow); - ulong cross1 = unchecked(leftHigh * rightLow); - ulong cross2 = unchecked(leftLow * rightHigh); - ulong highProduct = unchecked(leftHigh * rightHigh); - - ulong carry = unchecked((lowProduct >> 32) + (uint)cross1 + (uint)cross2); - low = unchecked((lowProduct & 0xFFFFFFFFUL) | (carry << 32)); - high = unchecked(highProduct + (cross1 >> 32) + (cross2 >> 32) + (carry >> 32)); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus) - { - // EvenPerfectBitScanner never routes moduli ≤ 1 through this path; keep the guard commented out so the hot loop stays branch-free. - // if (modulus <= 1UL) - // { - // return 0UL; - // } - - // EvenPerfectBitScanner always works with positive exponents here; leave the defensive check commented for future reference. - // if (exponent == 0UL) - // { - // return 1UL % modulus; - // } - - int bitLength = GetPortableBitLengthGpu(exponent); - int windowSize = GetWindowSizeGpu(bitLength); - int oddPowerCount = 1 << (windowSize - 1); - - // oddPowerCount is always ≥ 1 for the production workloads handled by EvenPerfectBitScanner. - // if (oddPowerCount <= 0) - // { - // return 1UL; - // } - - ulong[] oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount); - Span oddPowers = oddPowersArray.AsSpan(0, oddPowerCount); - - InitializeStandardOddPowers(modulus, oddPowers); - - ulong result = 1UL; - int index = bitLength - 1; - - while (index >= 0) - { - ulong currentBit = (exponent >> index) & 1UL; - ulong squared = result.MulModGpu(result, modulus); - result = currentBit == 0UL ? squared : result; - index = currentBit == 0UL ? index - 1 : index; - if (currentBit == 0UL) - { - continue; - } - - int windowStart = index - windowSize + 1; - // EvenPerfectBitScanner guarantees windowStart stays non-negative for production exponents, so keep the guard commented out to avoid extra branching. - // if (windowStart < 0) - // { - // windowStart = 0; - // } - - windowStart = GetNextSetBitIndexGpu(exponent, windowStart); - - int windowLength = index - windowStart + 1; - for (int square = 0; square < windowLength; square++) - { - result = result.MulModGpu(result, modulus); - } - - ulong mask = (1UL << windowLength) - 1UL; - ulong windowValue = (exponent >> windowStart) & mask; - int tableIndex = (int)((windowValue - 1UL) >> 1); - ulong multiplier = oddPowers[tableIndex]; - result = result.MulModGpu(multiplier, modulus); - - index = windowStart - 1; - } - - ThreadStaticPools.UlongPool.Return(oddPowersArray); - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong Pow2ModWindowedGpuKernel(ulong exponent, ulong modulus) - { - int bitLength = GetPortableBitLengthGpu(exponent); - int windowSize = GetWindowSizeGpu(bitLength); - ReadOnlyGpuUInt128 modulusWide = new(modulus); - GpuUInt128 result = new(1UL); - int index = bitLength - 1; - - while (index >= 0) - { - ulong currentBit = (exponent >> index) & 1UL; - if (currentBit == 0UL) - { - result.MulMod(result.AsReadOnly(), modulusWide); - index--; - continue; - } - - int windowStartCandidate = index - windowSize + 1; - int negativeMask = windowStartCandidate >> 31; - windowStartCandidate &= ~negativeMask; - int windowStart = GetNextSetBitIndexGpu(exponent, windowStartCandidate); - int windowLength = index - windowStart + 1; - - for (int square = 0; square < windowLength; square++) - { - result.MulMod(result.AsReadOnly(), modulusWide); - } - - ulong mask = (1UL << windowLength) - 1UL; - ulong windowValue = (exponent >> windowStart) & mask; - ulong multiplier = ComputeWindowedOddPowerGpuKernel(windowValue, modulusWide); - result.MulMod(multiplier, modulusWide); - - index = windowStart - 1; - } - - return result.Low; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong ComputeWindowedOddPowerGpuKernel(ulong windowValue, in ReadOnlyGpuUInt128 modulus) - { - GpuUInt128 result = new(2UL); - if (windowValue == 1UL) - { - // Single-bit windows show up whenever the bit scan finds an isolated exponent bit, - // which happens routinely on the GPU path, so keep the fast return enabled. - return 2UL; - } - - ulong remaining = (windowValue - 1UL) >> 1; - GpuUInt128 squareBase = new(2UL); - squareBase.MulMod(squareBase.AsReadOnly(), modulus); - - while (remaining != 0UL) - { - if ((remaining & 1UL) != 0UL) - { - result.MulMod(squareBase.AsReadOnly(), modulus); - } - - remaining >>= 1; - if (remaining != 0UL) - { - squareBase.MulMod(squareBase.AsReadOnly(), modulus); - } - } - - return result.Low; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong Pow2MontgomeryModWindowedGpuKeepMontgomery(in MontgomeryDivisorDataGpu divisor, ulong exponent) - { - return Pow2MontgomeryModWindowedGpuMontgomeryResult(divisor, exponent); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong Pow2MontgomeryModWindowedGpuConvertToStandard(in MontgomeryDivisorDataGpu divisor, ulong exponent) - { - ulong modulus = divisor.Modulus; - ulong nPrime = divisor.NPrime; - ulong montgomeryResult = Pow2MontgomeryModWindowedGpuMontgomeryResult(divisor, exponent); - return MontgomeryMultiplyGpu(montgomeryResult, 1UL, modulus, nPrime); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong Pow2MontgomeryModWindowedGpuMontgomeryResult(in MontgomeryDivisorDataGpu divisor, ulong exponent) - { - ulong modulus = divisor.Modulus; - // Reuse the GPU-compatible MontgomeryMultiply helper highlighted in MontgomeryMultiplyBenchmarks so the accelerator path - // mirrors the fastest CPU reduction without relying on UInt128 intrinsics. - // We shouldn't hit it in production code. If we do, we're doing something wrong. - // if (exponent == 0UL) - // { - // return divisor.MontgomeryOne; - // } - - // We barely ever hit it in production code. It's e.g. 772 calls out of billions - // if (exponent <= Pow2WindowFallbackThreshold) - // { - // return Pow2MontgomeryModSingleBit(exponent, divisor); - // } - - int bitLength = GetPortableBitLengthGpu(exponent); - int windowSize = GetWindowSizeGpu(bitLength); - ulong result = divisor.MontgomeryOne; - ulong nPrime = divisor.NPrime; - - int index = bitLength - 1; - while (index >= 0) - { - ulong currentBit = (exponent >> index) & 1UL; - ulong squared = MontgomeryMultiplyGpu(result, result, modulus, nPrime); - bool processWindow = currentBit != 0UL; - - result = processWindow ? result : squared; - index -= (int)(currentBit ^ 1UL); - - int windowStartCandidate = index - windowSize + 1; - // if (windowStartCandidate < 0) - // { - // windowStartCandidate = 0; - // } - // Clamp the negative offset without branching so the GPU loop stays divergence-free. - // This still handles windows that would otherwise extend past the most significant bit. - int negativeMask = windowStartCandidate >> 31; - windowStartCandidate &= ~negativeMask; - - int windowStart = processWindow ? GetNextSetBitIndexGpu(exponent, windowStartCandidate) : windowStartCandidate; - int windowLength = processWindow ? index - windowStart + 1 : 0; - for (int square = 0; square < windowLength; square++) - { - result = MontgomeryMultiplyGpu(result, result, modulus, nPrime); - } - - result = processWindow - ? MontgomeryMultiplyGpu( - result, - ComputeMontgomeryOddPowerGpu( - (exponent >> windowStart) & ((1UL << windowLength) - 1UL), - divisor, - modulus, - nPrime), - modulus, - nPrime) - : result; - - index = processWindow ? windowStart - 1 : index; - } - - return result; - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2MontgomeryModWindowedKeepGpu(this ulong exponent, in MontgomeryDivisorData divisor) => Pow2MontgomeryGpuCalculator.CalculateKeep(exponent, divisor); @@ -323,87 +55,6 @@ public static partial class ULongExtensions public static ulong Pow2MontgomeryModWindowedConvertGpu(this ulong exponent, PrimeOrderCalculatorAccelerator gpu, in MontgomeryDivisorData divisor) => Pow2MontgomeryGpuCalculator.CalculateConvert(gpu, exponent, divisor); - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void InitializeStandardOddPowers(ulong modulus, Span oddPowers) - { - // oddPowers spans at least one entry for production workloads; keep the guard commented out to avoid redundant checks in the hot loop. - // if (oddPowers.IsEmpty) - // { - // return; - // } - - // EvenPerfectBitScanner feeds odd prime moduli (≥ 3) here, so the base value stays within range without a modulo reduction. - ulong baseValue = 2UL; - oddPowers[0] = baseValue; - if (oddPowers.Length == 1) - { - return; - } - - ulong square = baseValue.MulModGpu(baseValue, modulus); - for (int i = 1; i < oddPowers.Length; i++) - { - ulong previous = oddPowers[i - 1]; - oddPowers[i] = previous.MulModGpu(square, modulus); - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong ComputeMontgomeryOddPowerGpu(ulong exponent, in MontgomeryDivisorDataGpu divisor, ulong modulus, ulong nPrime) - { - ulong baseValue = divisor.MontgomeryTwo; - ulong power = divisor.MontgomeryOne; - ulong remaining = exponent; - - while (remaining != 0UL) - { - ulong bit = remaining & 1UL; - ulong multiplied = MontgomeryMultiplyGpu(power, baseValue, modulus, nPrime); - ulong mask = (ulong)-(long)bit; - power = (power & ~mask) | (multiplied & mask); - - remaining >>= 1; - if (remaining == 0UL) - { - break; - } - - baseValue = MontgomeryMultiplyGpu(baseValue, baseValue, modulus, nPrime); - } - - return power; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static int GetPortableBitLengthGpu(ulong value) => 64 - XMath.LeadingZeroCount(value); - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static int GetWindowSizeGpu(int n) - { - int w = Pow2WindowSize; - - w = (n <= 671) ? 7 : w; - w = (n <= 239) ? 6 : w; - w = (n <= 79) ? 5 : w; - w = (n <= 23) ? 4 : w; - - int m = (n >= 1) ? n : 1; - w = (n <= 6) ? m : w; - - return w; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static int GetNextSetBitIndexGpu(ulong exponent, int startIndex) - { - ulong guard = (ulong)(((long)startIndex - 64) >> 63); - int shift = startIndex & 63; - ulong mask = (~0UL << shift) & guard; - ulong masked = exponent & mask; - - return XMath.TrailingZeroCount(masked); - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) { diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 64e87f2..ce5bb40 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -5,11 +5,10 @@ namespace PerfectNumbers.Core; public static partial class ULongExtensions { - private const int Pow2WindowSize = 8; private const ulong Pow2WindowFallbackThreshold = 32UL; - private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1); + private const int Pow2WindowOddPowerCount = 1 << (PerfectNumberConstants.Pow2WindowSize - 1); - private static readonly UInt128 Pow2WindowedModulusThreshold = (UInt128)PerfectNumberConstants.MaxQForDivisorCycles; + private static readonly UInt128 Pow2WindowedModulusThreshold = (UInt128)PerfectNumberConstants.MaxQForDivisorCycles; private static readonly byte[] BitStatsZeroCount = new byte[256]; private static readonly byte[] BitStatsPrefixZero = new byte[256];