diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 229f5dd..97cfb4e 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -515,7 +515,9 @@ internal static class Program { _byDivisorTester = mersenneOnGpu ? new MersenneNumberDivisorByDivisorGpuTester() - : new MersenneNumberDivisorByDivisorCpuTester(); + : new MersenneNumberDivisorByDivisorCpuTester( + useGpuPow2Mod: !GpuContextPool.ForceCpu, + useGpuCycleComputation: orderOnGpu); _byDivisorTester.BatchSize = scanBatchSize; } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index c744bf6..24ab14a 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -3,6 +3,7 @@ using System.Collections.Concurrent; using System.Collections.Generic; using System.Numerics; using PerfectNumbers.Core; +using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core.Cpu; @@ -10,11 +11,18 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv { private readonly object _sync = new(); private readonly ConcurrentBag _sessionPool = new(); + private readonly bool _useGpuPow2Mod; + private readonly bool _useGpuCycleComputation; private ulong _divisorLimit; private ulong _lastStatusDivisor; private bool _isConfigured; private int _batchSize = 1_024; + public MersenneNumberDivisorByDivisorCpuTester(bool useGpuPow2Mod = false, bool useGpuCycleComputation = false) + { + _useGpuPow2Mod = useGpuPow2Mod; + _useGpuCycleComputation = useGpuCycleComputation; + } public int BatchSize { @@ -190,6 +198,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } Dictionary? factorCache = null; + Dictionary? cycleRemainderCache = null; + CycleHint cycleHint = default; DivisorCycleCache cycleCache = DivisorCycleCache.Shared; byte step10 = (byte)(step % 10UL); @@ -223,26 +233,14 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv if (admissible && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0) { MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(candidate); - ulong divisorCycle; - - if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles) + if (TryResolveDivisorCycle(candidate, prime, divisorData, cycleCache, ref factorCache, ref cycleHint, out ulong divisorCycle) && divisorCycle != 0UL) { - divisorCycle = cycleCache.GetCycleLength(candidate); - } - else - { - factorCache ??= new Dictionary(8); - if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL) + if (CheckDivisor(prime, candidate, divisorCycle, divisorData, ref cycleRemainderCache, ref cycleHint) != 0) { - divisorCycle = cycleCache.GetCycleLength(candidate); + processedAll = true; + return true; } } - - if (divisorCycle != 0UL && CheckDivisor(prime, divisorCycle, divisorData) != 0) - { - processedAll = true; - return true; - } } divisor += step; @@ -258,12 +256,152 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } - private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) + private byte CheckDivisor( + ulong prime, + ulong divisor, + ulong divisorCycle, + in MontgomeryDivisorData divisorData, + ref Dictionary? cycleRemainders, + ref CycleHint cycleHint) { - ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData); + if (divisorCycle == 0UL) + { + throw new InvalidOperationException($"Cycle length resolution returned zero for divisor {divisor}."); + } + + Dictionary? remainderCache = cycleRemainders; + if (_useGpuPow2Mod && remainderCache is null) + { + remainderCache = new Dictionary(8); + cycleRemainders = remainderCache; + } + + ulong reducedExponent; + if (remainderCache is not null && remainderCache.TryGetValue(divisorCycle, out ulong cachedRemainder)) + { + reducedExponent = cachedRemainder; + } + else if (cycleHint.HasCycle && cycleHint.LastCycleLength == divisorCycle) + { + reducedExponent = cycleHint.LastReducedExponent; + } + else + { + reducedExponent = prime % divisorCycle; + if (remainderCache is not null) + { + remainderCache[divisorCycle] = reducedExponent; + } + } + + cycleHint.LastCycleLength = divisorCycle; + cycleHint.LastReducedExponent = reducedExponent; + cycleHint.HasCycle = true; + + if (_useGpuPow2Mod && remainderCache is not null) + { + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(reducedExponent, divisor, out ulong remainder, divisorData); + if (status == GpuPow2ModStatus.Success) + { + return remainder == 1UL ? (byte)1 : (byte)0; + } + + PrimeOrderGpuHeuristics.ReportPow2Failure( + "by-divisor residue scan", + divisor, + status); + } + + ulong residue = reducedExponent.Pow2MontgomeryModFromCycleRemainder(divisorData); return residue == 1UL ? (byte)1 : (byte)0; } + private bool TryResolveDivisorCycle( + ulong divisor, + ulong prime, + in MontgomeryDivisorData divisorData, + DivisorCycleCache cycleCache, + ref Dictionary? factorCache, + ref CycleHint cycleHint, + out ulong cycleLength) + { + if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) + { + cycleLength = cycleCache.GetCycleLength(divisor); + if (cycleLength == 0UL) + { + throw new InvalidOperationException($"Missing cached cycle length for divisor {divisor}."); + } + + cycleHint.LastCycleLength = cycleLength; + cycleHint.HasCycle = true; + return true; + } + + if (_useGpuCycleComputation && TryCalculateCycleLengthGpu(divisor, prime, divisorData, ref cycleHint, out cycleLength)) + { + return cycleLength != 0UL; + } + + factorCache ??= new Dictionary(8); + if (MersenneDivisorCycles.TryCalculateCycleLengthForExponent( + divisor, + prime, + divisorData, + factorCache, + out cycleLength, + allowGpuPow2: _useGpuPow2Mod) && cycleLength != 0UL) + { + cycleHint.LastCycleLength = cycleLength; + cycleHint.HasCycle = true; + return true; + } + + cycleLength = cycleCache.GetCycleLength(divisor); + if (cycleLength == 0UL) + { + return false; + } + + cycleHint.LastCycleLength = cycleLength; + cycleHint.HasCycle = true; + return true; + } + + private bool TryCalculateCycleLengthGpu( + ulong divisor, + ulong prime, + in MontgomeryDivisorData divisorData, + ref CycleHint cycleHint, + out ulong cycleLength) + { + PrimeOrderCalculator.PrimeOrderSearchConfig config = PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault; + ulong? previousOrder = cycleHint.HasCycle ? cycleHint.LastCycleLength : null; + + try + { + if (!PrimeOrderGpuHeuristics.TryCalculateOrder(divisor, previousOrder, config, divisorData, out PrimeOrderCalculator.PrimeOrderResult result)) + { + throw new InvalidOperationException("GPU requested fallback."); + } + + cycleLength = result.Order; + if (cycleLength == 0UL) + { + return false; + } + + cycleHint.LastCycleLength = cycleLength; + cycleHint.HasCycle = true; + return true; + } + catch (InvalidOperationException ex) + { + Console.Error.WriteLine($"GPU cycle length computation failed for divisor {divisor}: {ex.Message}"); + throw; + } + } + private static byte AddMod(byte value, byte delta, byte modulus) { int sum = value + delta; @@ -310,6 +448,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return (1UL << (int)(maxPrime - 1UL)) - 1UL; } + private struct CycleHint + { + public ulong LastCycleLength; + public ulong LastReducedExponent; + public bool HasCycle; + } + private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit) { if (prime <= 1UL) diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index fa4251e..bb09d15 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -4,10 +4,12 @@ using System.Collections.Concurrent; using ILGPU; using ILGPU.Runtime; using System.Numerics; +using System.Globalization; using PerfectNumbers.Core; using System.Runtime.InteropServices; using System.Runtime.CompilerServices; using ILGPU.Runtime.OpenCL; + namespace PerfectNumbers.Core.Gpu; internal enum GpuPow2ModStatus @@ -27,6 +29,25 @@ internal static partial class PrimeOrderGpuHeuristics private static readonly ConcurrentDictionary OrderKernelCache = new(); private static readonly ConcurrentDictionary SmallPrimeDeviceCaches = new(); + private static void ReportPow2FailureCore(string context, string modulusDescription, GpuPow2ModStatus status) + { + string message = $"GPU pow2mod failed with status {status} for divisor {modulusDescription} during {context}."; + Console.Error.WriteLine(message); + throw new InvalidOperationException(message); + } + + internal static void ReportPow2Failure(string context, ulong modulus, GpuPow2ModStatus status) + { + string modulusDescription = modulus.ToString(CultureInfo.InvariantCulture); + ReportPow2FailureCore(context, modulusDescription, status); + } + + internal static void ReportPow2Failure(string context, UInt128 modulus, GpuPow2ModStatus status) + { + string modulusDescription = modulus.ToString(CultureInfo.InvariantCulture); + ReportPow2FailureCore(context, modulusDescription, status); + } + public readonly struct OrderKernelConfig(ulong previousOrder, byte hasPreviousOrder, uint smallFactorLimit, int maxPowChecks, int mode) { public readonly ulong PreviousOrder = previousOrder; diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index ee77260..3d8fe7f 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -265,11 +265,12 @@ public class MersenneDivisorCycles } public static bool TryCalculateCycleLengthForExponent( - ulong divisor, + ulong divisor, ulong exponent, in MontgomeryDivisorData divisorData, Dictionary? factorCache, - out ulong cycleLength) + out ulong cycleLength, + bool allowGpuPow2 = false) { cycleLength = 0UL; @@ -290,7 +291,8 @@ public class MersenneDivisorCycles previousOrder: null, divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, - PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); + PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu, + allowGpuPow2); if (orderResult.Order != 0UL) { cycleLength = orderResult.Order; @@ -327,7 +329,7 @@ public class MersenneDivisorCycles return false; } - cycleLength = ReduceOrder(divisorData, phi, factorCounts); + cycleLength = ReduceOrder(divisorData, phi, factorCounts, allowGpuPow2); return true; } @@ -471,7 +473,7 @@ public class MersenneDivisorCycles } } - private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary factorCounts) + private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary factorCounts, bool allowGpuPow2) { if (factorCounts.Count == 0) { @@ -505,7 +507,7 @@ public class MersenneDivisorCycles } ulong candidate = order / prime; - if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL) + if (Pow2EqualsOne(candidate, divisorData, allowGpuPow2)) { order = candidate; continue; @@ -523,6 +525,26 @@ public class MersenneDivisorCycles } } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool Pow2EqualsOne(ulong exponent, in MontgomeryDivisorData divisorData, bool allowGpuPow2) + { + if (allowGpuPow2) + { + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, divisorData.Modulus, out ulong remainder, divisorData); + if (status == GpuPow2ModStatus.Success) + { + return remainder == 1UL; + } + + PrimeOrderGpuHeuristics.ReportPow2Failure( + "divisor cycle reduction", + divisorData.Modulus, + status); + } + + return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL; + } + private static ulong PollardRho64(ulong n) { if ((n & 1UL) == 0UL) @@ -924,7 +946,7 @@ public class MersenneDivisorCycles return CalculateCycleLengthFallback(divisor); } - internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength) + internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength, bool allowGpuPow2 = false) { if ((divisor & (divisor - 1UL)) == 0UL) { @@ -945,7 +967,8 @@ public class MersenneDivisorCycles previousOrder: null, divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, - PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); + PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu, + allowGpuPow2); if (orderResult.Order != 0UL) { cycleLength = orderResult.Order; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs index 9e1fd27..594f2fa 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs @@ -31,7 +31,8 @@ internal static partial class PrimeOrderCalculator PrimeOrderSearchConfig config, PrimeOrderHeuristicDevice device) { - using var scope = UsePow2Mode(device); + bool enableGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu; + using var scope = UsePow2Mode(device, enableGpuPow2); MontgomeryDivisorData divisorData; if (prime <= ulong.MaxValue) { @@ -51,7 +52,7 @@ internal static partial class PrimeOrderCalculator ulong prime64 = (ulong)prime; divisorData = MontgomeryDivisorData.FromModulus(prime64); - PrimeOrderResult result = Calculate(prime64, previous, divisorData, config, device); + PrimeOrderResult result = Calculate(prime64, previous, divisorData, config, device, enableGpuPow2); return new PrimeOrderResultWide(result.Status, result.Order); } else @@ -1030,6 +1031,11 @@ internal static partial class PrimeOrderCalculator { return wideRemainder; } + + PrimeOrderGpuHeuristics.ReportPow2Failure( + "wide prime order pow2 check", + modulus, + wideStatus); } return exponent.Pow2MontgomeryModWindowed(modulus); diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index f079c42..d898425 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -90,13 +90,13 @@ internal static partial class PrimeOrderCalculator private readonly bool _previousAllow; private readonly PrimeOrderHeuristicDevice _previousDevice; - public Pow2ModeScope(PrimeOrderHeuristicDevice device) + public Pow2ModeScope(PrimeOrderHeuristicDevice device, bool allowGpuPow2) { _previousInitialized = s_pow2ModeInitialized; _previousAllow = s_allowGpuPow2; _previousDevice = s_deviceMode; s_pow2ModeInitialized = true; - s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu; + s_allowGpuPow2 = allowGpuPow2; s_deviceMode = device; } @@ -124,7 +124,7 @@ internal static partial class PrimeOrderCalculator } } - private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device) => new(device); + private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device, bool allowGpuPow2) => new(device, allowGpuPow2); private static DebugLoggingScope UseDebugLogging(bool enabled) => new(enabled); @@ -154,9 +154,11 @@ internal static partial class PrimeOrderCalculator ulong? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config, - PrimeOrderHeuristicDevice device) + PrimeOrderHeuristicDevice device, + bool allowGpuPow2 = false) { - using var scope = UsePow2Mode(device); + bool enableGpuPow2 = allowGpuPow2 || device == PrimeOrderHeuristicDevice.Gpu; + using var scope = UsePow2Mode(device, enableGpuPow2); return CalculateInternal(prime, previousOrder, divisorData, config); } @@ -443,7 +445,6 @@ internal static partial class PrimeOrderCalculator int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks; int powUsed = 0; int candidateCount = candidates.Count; - bool allowGpuBatch = true; Span candidateSpan = CollectionsMarshal.AsSpan(candidates); DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)"); @@ -467,9 +468,11 @@ internal static partial class PrimeOrderCalculator bool gpuSuccess = false; bool gpuStackRemainders = false; GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable; + bool attemptedGpu = false; - if (allowGpuBatch && IsGpuPow2Allowed) + if (IsGpuPow2Allowed) { + attemptedGpu = true; if (batchSize <= StackGpuBatchSize) { Span localRemainders = stackGpuRemainders.Slice(0, batchSize); @@ -496,11 +499,14 @@ internal static partial class PrimeOrderCalculator gpuPool = null; } } + } - if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable)) - { - allowGpuBatch = false; - } + if (attemptedGpu && !gpuSuccess) + { + PrimeOrderGpuHeuristics.ReportPow2Failure( + "prime order candidate scan", + prime, + status); } for (int i = 0; i < batchSize && powUsed < powBudget; i++) @@ -780,6 +786,11 @@ internal static partial class PrimeOrderCalculator { return remainder == 1UL; } + + PrimeOrderGpuHeuristics.ReportPow2Failure( + "prime order pow2 check", + prime, + status); } return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;