Various micro- and medium- optimizations

This commit is contained in:
Marek Lemańczyk committed 2025-10-07 22:14:34 +02:00
1 parent 2eb622030e
commit cd1debe9d6
25 files changed
+4403 -3210

No files matched your search

@@ -39,7 +39,7 @@ public class Pow2MontgomeryModBenchmarks
ulong smallModulus = NextSmallOddModulus();
_smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus);
_smallExponents[i] = NextSmallExponent();
_smallCycles[i] = MersenneDivisorCycles.CalculateCycleLength(smallModulus);
_smallCycles[i] = MersenneDivisorCycles.CalculateCycleLength(smallModulus, MontgomeryDivisorData.FromModulus(smallModulus));
(ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle();
_largeDivisors[i] = CreateMontgomeryDivisorData(largeModulus);
@@ -261,9 +261,10 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
private ulong CalculateCycleLengthWithHeuristics(ulong modulus, ulong fallbackOrder = 0UL)
{
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus);
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus);
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
}
// When the below is uncommented, the benchmark never completes due to the primality test cost on very large random set.
@@ -281,6 +282,7 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
modulus,
_previousPrimeOrder,
divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
if (orderResult.Order != 0UL)
@@ -292,7 +294,7 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
_previousPrimeOrder = null;
Console.WriteLine($"Heuristic failed for {modulus}, falling back to full cycle calculation");
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus);
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
}
private UInt128 CalculateCycleLengthWithHeuristics(UInt128 modulus)
File diff suppressed because it is too large. Load diff
+9 -7
View File
@@ -1,5 +1,6 @@
using BenchmarkDotNet.Configs;
using BenchmarkDotNet.Running;
using PerfectNumbers.Core.Gpu;
namespace EvenPerfectBitScanner.Benchmarks;
@@ -7,15 +8,16 @@ public static class Program
{
public static void Main(string[] args)
{
var config = DefaultConfig.Instance.WithOptions(ConfigOptions.JoinSummary);
BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly).Run(args, config);
// var config = DefaultConfig.Instance.WithOptions(ConfigOptions.JoinSummary | ConfigOptions.DisableOptimizationsValidator);
// BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly).Run(args, config);
// Pow2MontgomeryModCycleComputationBenchmarks benchmarks = new()
// {
// Scale = Pow2MontgomeryModCycleComputationBenchmarks.InputScale.Large
// };
Pow2MontgomeryModCycleComputationBenchmarks benchmarks = new()
{
Scale = Pow2MontgomeryModCycleComputationBenchmarks.InputScale.Large
};
// benchmarks.Setup();
GpuPrimeWorkLimiter.SetLimit(1);
benchmarks.Setup();
// Console.WriteLine("Press ENTER to run");
// Console.ReadLine();
@@ -25,8 +25,8 @@ public class PrimeOrderGpuHeuristicsTests
try
{
PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0;
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow);
remainder.Should().Be(0UL);
@@ -40,7 +40,9 @@ public class PrimeOrderGpuHeuristicsTests
[Fact]
public void TryPow2Mod_returns_success_for_supported_prime()
{
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, 193UL, out ulong remainder);
const ulong prime = 193UL;
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder, divisorData);
status.Should().Be(GpuPow2ModStatus.Success);
remainder.Should().Be(2UL);
@@ -54,7 +56,8 @@ public class PrimeOrderGpuHeuristicsTests
PrimeOrderGpuHeuristics.OverrideCapabilitiesForTesting(capability);
try
{
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow);
remainder.Should().Be(0UL);
@@ -73,12 +76,14 @@ public class PrimeOrderGpuHeuristicsTests
PrimeOrderGpuHeuristics.OverrideCapabilitiesForTesting(capability);
try
{
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(16UL, 193UL, out ulong remainder);
const ulong prime = 193UL;
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(16UL, prime, out ulong remainder, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow);
remainder.Should().Be(0UL);
GpuPow2ModStatus followUp = PrimeOrderGpuHeuristics.TryPow2Mod(4UL, 193UL, out ulong followUpRemainder);
GpuPow2ModStatus followUp = PrimeOrderGpuHeuristics.TryPow2Mod(4UL, prime, out ulong followUpRemainder, divisorData);
followUp.Should().Be(GpuPow2ModStatus.Success);
followUpRemainder.Should().Be(16UL);
}
@@ -98,7 +103,8 @@ public class PrimeOrderGpuHeuristicsTests
ulong[] remainders = new ulong[2];
ulong[] exponents = { 1UL, 2UL };
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow);
remainders[0].Should().Be(0UL);
@@ -115,7 +121,9 @@ public class PrimeOrderGpuHeuristicsTests
{
ulong[] remainders = { 123UL, 456UL };
ulong[] exponents = { 1UL, 2UL };
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, 257UL, remainders);
const ulong prime = 257UL;
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData);
status.Should().Be(GpuPow2ModStatus.Success);
remainders[0].Should().Be(2UL);
@@ -133,7 +141,8 @@ public class PrimeOrderGpuHeuristicsTests
ulong[] exponents = { 1UL, 2UL };
ulong[] remainders = { 5UL, 6UL };
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow);
remainders[0].Should().Be(0UL);
@@ -156,7 +165,9 @@ public class PrimeOrderGpuHeuristicsTests
ulong[] exponents = { 2UL, 16UL };
ulong[] remainders = { 9UL, 11UL };
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, 193UL, remainders);
const ulong prime = 193UL;
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData);
status.Should().Be(GpuPow2ModStatus.Overflow);
remainders[0].Should().Be(0UL);
@@ -322,16 +333,22 @@ public class PrimeOrderGpuHeuristicsTests
}
}
[Fact]
public void TryPow2ModBatch_throws_when_remainder_span_is_too_small()
{
ulong[] remainders = new ulong[1];
ulong[] exponents = { 1UL, 2UL };
// This will never happen in production code because the caller ensures the remainder span is large enough.
// [Fact]
// public void TryPow2ModBatch_throws_when_remainder_span_is_too_small()
// {
// ulong[] remainders = new ulong[1];
// ulong[] exponents = { 1UL, 2UL };
Action act = () => PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, 5UL, remainders);
// Action act = () =>
// {
// const ulong prime = 5UL;
// MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
// PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData);
// };
act.Should().Throw<ArgumentException>();
}
// act.Should().Throw<ArgumentException>();
// }
[Fact]
public void Calculate_falls_back_to_cpu_when_gpu_overflow_is_marked()
@@ -341,7 +358,8 @@ public class PrimeOrderGpuHeuristicsTests
{
PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0;
PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(prime, null, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(prime, null, divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
result.Order.Should().Be(3UL);
@@ -38,7 +38,8 @@ public class MersenneDivisorCyclesGenerateGpuTests
foreach (var pair in pairs)
{
expected.Should().Contain(pair.divisor);
pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLength(pair.divisor));
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(pair.divisor);
pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLength(pair.divisor, divisorData));
}
}
finally
@@ -86,7 +87,8 @@ public class MersenneDivisorCyclesGenerateGpuTests
foreach (var pair in pairs)
{
uniqueDivisors.Add(pair.divisor).Should().BeTrue();
pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLength(pair.divisor));
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(pair.divisor);
pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLength(pair.divisor, divisorData));
}
uniqueDivisors.Should().BeEquivalentTo(expected);
@@ -18,7 +18,8 @@ public class MersenneDivisorCyclesTests
[MemberData(nameof(Divisors))]
public void CalculateCycleLength_returns_expected_value(ulong divisor, ulong expected)
{
MersenneDivisorCycles.CalculateCycleLength(divisor).Should().Be(expected);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData).Should().Be(expected);
MersenneDivisorCycles.CalculateCycleLengthGpu(divisor).Should().Be(expected);
}
}
@@ -35,13 +35,13 @@ public class MersenneNumberDivisorCpuTesterTests
ulong[] primes = { 5UL, 7UL, 11UL, 13UL };
byte[] hits = new byte[primes.Length];
ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL);
session.CheckDivisor(23UL, cycle23, primes, hits);
ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL));
session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
Array.Fill(hits, (byte)0);
ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL);
session.CheckDivisor(31UL, cycle31, primes, hits);
ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL));
session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits);
hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 });
}
}
@@ -146,13 +146,13 @@ public class MersenneNumberDivisorGpuTesterTests
ulong[] primes = { 31UL, 37UL, 41UL, 43UL };
byte[] hits = new byte[primes.Length];
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL);
session.CheckDivisor(223UL, cycle223, primes, hits);
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0 });
Array.Fill(hits, (byte)0);
ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLength(13367UL);
session.CheckDivisor(13367UL, cycle13367, primes, hits);
ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLength(13367UL, MontgomeryDivisorData.FromModulus(13367UL));
session.CheckDivisor(13367UL, MontgomeryDivisorData.FromModulus(13367UL), cycle13367, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
}
@@ -171,8 +171,8 @@ public class MersenneNumberDivisorGpuTesterTests
ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
byte[] hits = new byte[primes.Length];
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL);
session.CheckDivisor(223UL, cycle223, primes, hits);
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
}
@@ -192,8 +192,8 @@ public class MersenneNumberDivisorGpuTesterTests
ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
byte[] hits = new byte[primes.Length];
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL);
session.CheckDivisor(223UL, cycle223, primes, hits);
ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
}
@@ -48,6 +48,7 @@ public class PrimeOrderCalculatorTests
PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
prime,
previousOrder: null,
MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
@@ -64,11 +65,13 @@ public class PrimeOrderCalculatorTests
PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
prime,
previousOrder,
MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
prime,
previousOrder: null,
MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault);
heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
@@ -85,6 +88,7 @@ public class PrimeOrderCalculatorTests
PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
prime,
previousOrder: null,
MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
@@ -104,6 +108,7 @@ public class PrimeOrderCalculatorTests
PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
prime: 13UL,
previousOrder: null,
MontgomeryDivisorData.FromModulus(13UL),
config);
result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
@@ -120,9 +125,11 @@ public class PrimeOrderCalculatorTests
maxPowChecks: 1,
mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL);
PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
prime: 239UL,
previousOrder: null,
divisorData,
heuristicConfig);
heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
@@ -137,6 +144,7 @@ public class PrimeOrderCalculatorTests
PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
prime: 239UL,
previousOrder: null,
divisorData,
strictConfig);
strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
@@ -194,7 +194,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
processedCount++;
lastProcessed = candidate;
if (MersenneDivisorCycles.CycleEqualsExponentForMersenneCandidate(candidate, prime))
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate);
if (MersenneDivisorCycles.CycleEqualsExponentForMersenneCandidate(candidate, divisorData, prime))
{
processedAll = true;
return true;
@@ -272,7 +273,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
_disposed = false;
}
public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits)
public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits)
{
if (_disposed)
{
@@ -288,7 +289,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
for (int i = 0; i < length; i++)
{
ulong prime = primes[i];
bool divides = MersenneDivisorCycles.CycleEqualsExponentForMersenneCandidate(divisor, prime);
bool divides = MersenneDivisorCycles.CycleEqualsExponentForMersenneCandidate(divisor, divisorData, prime);
hits[i] = divides ? (byte)1 : (byte)0;
}
}
+4 -2
View File
@@ -109,14 +109,16 @@ public sealed class DivisorCycleCache
ArrayPool<int>.Shared.Return(rentedMissing, clearArray: false);
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void ComputeCyclesCpu(ReadOnlySpan<ulong> divisors, Span<ulong> cycles, ReadOnlySpan<int> indices)
{
for (int i = 0; i < indices.Length; i++)
{
int targetIndex = indices[i];
cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLength(divisors[targetIndex]);
ulong divisor = divisors[targetIndex];
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData);
}
}
+257 -254
View File
@@ -25,16 +25,16 @@ public readonly struct ResidueAutomatonArgs
public sealed class KernelContainer
{
// Serializes first-time initialization of kernels/buffers per accelerator.
public Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>? Order;
public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>? Incremental;
// Serializes first-time initialization of kernels/buffers per accelerator.
public Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>? Order;
public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>? Incremental;
public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>,
ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>,
ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>? Pow2Mod;
public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>? IncrementalOrder;
public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>? IncrementalOrder;
public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>? Pow2ModOrder;
@@ -45,22 +45,22 @@ public sealed class KernelContainer
public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimesLastSeven;
public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimesPow2LastSeven;
public static T InitOnce<T>(ref T? slot, Func<T> factory) where T : class
public static T InitOnce<T>(ref T? slot, Func<T> factory) where T : class
{
var current = Volatile.Read(ref slot);
if (current is not null)
{
var current = Volatile.Read(ref slot);
if (current is not null)
{
return current;
}
current = Volatile.Read(ref slot);
if (current is null)
{
current = factory();
Volatile.Write(ref slot, current);
}
return current;
return current;
}
current = Volatile.Read(ref slot);
if (current is null)
{
current = factory();
Volatile.Write(ref slot, current);
}
return current;
}
}
public readonly struct ResiduePrimeViews(
@@ -77,122 +77,120 @@ public readonly struct ResiduePrimeViews(
public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels)
{
private bool disposedValue;
private IDisposable? _limiter = limiter;
private GpuContextLease? _gpu = gpu;
private KernelContainer _kernels = kernels;
private AcceleratorStream? _stream;
private object? _executionLock = gpu.ExecutionLock;
private bool disposedValue;
private IDisposable? _limiter = limiter;
private GpuContextLease? _gpu = gpu;
private KernelContainer _kernels = kernels;
private AcceleratorStream? _stream;
private object? _executionLock = gpu.ExecutionLock;
public readonly Accelerator Accelerator => _gpu?.Accelerator ?? throw new NullReferenceException("GPU context is already released");
public readonly Accelerator Accelerator => _gpu?.Accelerator ?? throw new NullReferenceException("GPU context is already released");
public AcceleratorStream Stream
public AcceleratorStream Stream
{
get
{
get
if (_stream is { } stream)
{
return stream;
}
// It feels that we can combine these into 1 return statement, but that increases the max stack frame size to 3 from 2
stream = Accelerator.CreateStream();
_stream = stream;
return stream;
}
}
public readonly ExecutionScope EnterExecutionScope() => new ExecutionScope(_executionLock);
public readonly Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>> OrderKernel
{
get
{
var accel = Accelerator; // avoid capturing 'this' in lambda
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Order, () =>
{
if (_stream is { } stream)
{
return stream;
}
stream = Accelerator.CreateStream();
_stream = stream;
return stream;
}
}
public readonly ExecutionScope EnterExecutionScope()
{
return new ExecutionScope(_executionLock);
}
public readonly Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>> OrderKernel
{
get
{
var accel = Accelerator; // avoid capturing 'this' in lambda
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Order, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>(OrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>>();
});
}
}
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>(OrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>>();
});
}
}
}
public readonly Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>> Pow2ModKernel
{
get
{
get
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () =>
{
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
}
public readonly Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>> IncrementalKernel
{
get
{
get
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Incremental, () =>
{
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Incremental, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
}
public readonly Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>> IncrementalOrderKernel
{
get
{
get
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.IncrementalOrder, () =>
{
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.IncrementalOrder, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalOrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalOrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
}
public readonly Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>> Pow2ModOrderKernel
{
get
{
get
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Pow2ModOrder, () =>
{
var accel = Accelerator;
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Pow2ModOrder, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModOrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModOrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
}
}
}
// TODO: Plumb the small-cycles device buffer into all kernels that can benefit
// (some already accept it). Consider a compact type (byte/ushort) for memory footprint.
// TODO: Plumb the small-cycles device buffer into all kernels that can benefit
// (some already accept it). Consider a compact type (byte/ushort) for memory footprint.
private static void IncrementalKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong divMul,
ulong q0m10, ulong q0m8, ulong q0m3, ulong q0m5, ArrayView<ulong> orders,
@@ -206,21 +204,21 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
// residue automaton
ulong step10 = (exponent.Mod10() << 1).Mod10();
ulong r10 = q0m10 + (step10 * idx).Mod10();
r10 -= (r10 >= 10UL) ? 10UL : 0UL;
bool shouldCheck = r10 != 5UL; // skip q ≡ 5 (mod 10); other residues handled by r8/r3/r5
if (shouldCheck)
{
r10 -= (r10 >= 10UL) ? 10UL : 0UL;
bool shouldCheck = r10 != 5UL; // skip q ≡ 5 (mod 10); other residues handled by r8/r3/r5
if (shouldCheck)
{
ulong step8 = ((exponent & 7UL) << 1) & 7UL;
ulong r8 = (q0m8 + ((step8 * idx) & 7UL)) & 7UL;
if (r8 != 1UL && r8 != 7UL)
{
shouldCheck = false;
}
if (r8 != 1UL && r8 != 7UL)
{
shouldCheck = false;
}
else
{
// TODO: Swap these `%` filters to the shared Mod helpers so the residue automaton matches the benchmarked
// bitmask-based implementation for GPU workloads.
ulong step3 = ((exponent % 3UL) << 1) % 3UL;
// TODO: Swap these `%` filters to the shared Mod helpers so the residue automaton matches the benchmarked
// bitmask-based implementation for GPU workloads.
ulong step3 = ((exponent % 3UL) << 1) % 3UL;
ulong r3 = q0m3 + (step3 * idxMod3); r3 -= (r3 >= 3UL) ? 3UL : 0UL;
if (r3 == 0UL)
{
@@ -236,19 +234,20 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
}
}
}
}
}
if (!shouldCheck)
{
orders[index] = 0UL;
return;
}
if (!shouldCheck)
{
orders[index] = 0UL;
return;
}
// TODO: Is this expected that the value of kStart is modified?
GpuUInt128 k = kStart + (GpuUInt128)idx;
GpuUInt128 q = twoP;
q.Mul64(k);
q.Add(GpuUInt128.One);
// Small-cycles in-kernel early rejection from device table
if (q.High == 0UL && q.Low < (ulong)smallCycles.Length)
{
@@ -260,7 +259,9 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
}
}
// TODO: q will be lowered by GpuUInt128.One here. Is this expected? Maybe that's why we have issues with the incremental kernel? Review other places where we use GpuUInt128 operators.
GpuUInt128 phi = q - GpuUInt128.One;
// TODO: Do we have the check here because we don't support high values or this is expected result?
if (phi.High != 0UL)
{
orders[index] = 0UL;
@@ -332,29 +333,29 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
orders[index] = exponent;
}
private void Dispose(bool disposing)
private void Dispose(bool disposing)
{
if (!disposedValue)
{
if (!disposedValue)
{
if (disposing)
{
_stream?.Dispose();
_stream = null;
if (disposing)
{
_stream?.Dispose();
_stream = null;
_gpu?.Dispose();
_gpu = null;
_gpu?.Dispose();
_gpu = null;
_executionLock = null;
_executionLock = null;
_limiter?.Dispose();
_limiter = null;
}
_limiter?.Dispose();
_limiter = null;
}
// free unmanaged resources (unmanaged objects) and override finalizer
// set large fields to null
disposedValue = true;
}
}
// free unmanaged resources (unmanaged objects) and override finalizer
// set large fields to null
disposedValue = true;
}
}
private static void Pow2ModKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _,
ResidueAutomatonArgs ra, ArrayView<ulong> orders,
@@ -404,13 +405,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
return;
}
GpuUInt128 k = kStart + (GpuUInt128)idx;
// TODO: kStart is modified after this. Is this expected?
kStart.Add(idx);
GpuUInt128 q = twoP;
q.Mul64(k);
q.Mul64(kStart);
q.Add(GpuUInt128.One);
if (q.High == 0UL && q.Low < (ulong)smallCycles.Length)
{
ulong cycle = smallCycles[(int)q.Low];
// cycle should always be initialized if we're within array limit in production code
if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL)
{
orders[index] = 0UL;
@@ -517,40 +520,40 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
ulong idxMod3 = idx % 3UL;
ulong idxMod5 = idx % 5UL;
ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL;
bool shouldCheck = r10 != 5UL;
if (shouldCheck)
{
ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL;
if (r8 != 1UL && r8 != 7UL)
{
shouldCheck = false;
}
else
{
ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3);
if (r3 >= 6UL) r3 -= 6UL;
if (r3 >= 3UL) r3 -= 3UL;
if (r3 == 0UL)
{
shouldCheck = false;
}
else
{
ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5);
if (r5 >= 15UL) r5 -= 15UL;
if (r5 >= 10UL) r5 -= 10UL;
if (r5 >= 5UL) r5 -= 5UL;
if (r5 == 0UL)
{
shouldCheck = false;
}
}
}
}
if (!shouldCheck)
{
return;
}
bool shouldCheck = r10 != 5UL;
if (shouldCheck)
{
ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL;
if (r8 != 1UL && r8 != 7UL)
{
shouldCheck = false;
}
else
{
ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3);
if (r3 >= 6UL) r3 -= 6UL;
if (r3 >= 3UL) r3 -= 3UL;
if (r3 == 0UL)
{
shouldCheck = false;
}
else
{
ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5);
if (r5 >= 15UL) r5 -= 15UL;
if (r5 >= 10UL) r5 -= 10UL;
if (r5 >= 5UL) r5 -= 5UL;
if (r5 == 0UL)
{
shouldCheck = false;
}
}
}
}
if (!shouldCheck)
{
return;
}
GpuUInt128 k = kStart + (GpuUInt128)idx;
GpuUInt128 q = twoP;
@@ -603,40 +606,40 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
ulong idxMod3 = idx % 3UL;
ulong idxMod5 = idx % 5UL;
ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL;
bool shouldCheck = r10 != 5UL;
if (shouldCheck)
{
ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL;
if (r8 != 1UL && r8 != 7UL)
{
shouldCheck = false;
}
else
{
ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3);
if (r3 >= 6UL) r3 -= 6UL;
if (r3 >= 3UL) r3 -= 3UL;
if (r3 == 0UL)
{
shouldCheck = false;
}
else
{
ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5);
if (r5 >= 15UL) r5 -= 15UL;
if (r5 >= 10UL) r5 -= 10UL;
if (r5 >= 5UL) r5 -= 5UL;
if (r5 == 0UL)
{
shouldCheck = false;
}
}
}
}
if (!shouldCheck)
{
return;
}
bool shouldCheck = r10 != 5UL;
if (shouldCheck)
{
ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL;
if (r8 != 1UL && r8 != 7UL)
{
shouldCheck = false;
}
else
{
ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3);
if (r3 >= 6UL) r3 -= 6UL;
if (r3 >= 3UL) r3 -= 3UL;
if (r3 == 0UL)
{
shouldCheck = false;
}
else
{
ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5);
if (r5 >= 15UL) r5 -= 15UL;
if (r5 >= 10UL) r5 -= 10UL;
if (r5 >= 5UL) r5 -= 5UL;
if (r5 == 0UL)
{
shouldCheck = false;
}
}
}
}
if (!shouldCheck)
{
return;
}
GpuUInt128 k = kStart + (GpuUInt128)idx;
GpuUInt128 q = twoP;
@@ -661,35 +664,35 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
Atomic.Or(ref found[0], 1);
}
public void Dispose()
{
// Do not change this code. Put cleanup code in 'Dispose(bool disposing)' method
Dispose(disposing: true);
// Enabled this when the struct is converted to a class
// GC.SuppressFinalize(this);
}
public readonly struct ExecutionScope : IDisposable
{
private readonly object? _lock;
public ExecutionScope(object? sync)
{
_lock = sync;
if (_lock is not null)
{
Monitor.Enter(_lock);
}
}
public void Dispose()
{
// Do not change this code. Put cleanup code in 'Dispose(bool disposing)' method
Dispose(disposing: true);
// Enabled this when the struct is converted to a class
// GC.SuppressFinalize(this);
}
public readonly struct ExecutionScope : IDisposable
{
private readonly object? _lock;
public ExecutionScope(object? sync)
{
_lock = sync;
if (_lock is not null)
{
Monitor.Enter(_lock);
}
}
public void Dispose()
{
if (_lock is not null)
{
Monitor.Exit(_lock);
}
}
if (_lock is not null)
{
Monitor.Exit(_lock);
}
}
}
}
public class GpuKernelPool
@@ -771,19 +774,19 @@ public class GpuKernelPool
}
}
public static GpuKernelLease GetKernel(bool useGpuOrder)
{
var limiter = GpuPrimeWorkLimiter.Acquire();
var gpu = RentPreferred(preferCpu: !useGpuOrder);
var accelerator = gpu.Accelerator;
var kernels = GetKernels(accelerator);
return new(limiter, gpu, kernels);
}
public static GpuKernelLease GetKernel(bool useGpuOrder)
{
var limiter = GpuPrimeWorkLimiter.Acquire();
var gpu = RentPreferred(preferCpu: !useGpuOrder);
var accelerator = gpu.Accelerator;
var kernels = GetKernels(accelerator);
return new(limiter, gpu, kernels);
}
/// <summary>
/// Runs a GPU action with an acquired accelerator and stream.
/// </summary>
/// <param name="action">Action to run with (Accelerator, Stream).</param>
/// <summary>
/// Runs a GPU action with an acquired accelerator and stream.
/// </summary>
/// <param name="action">Action to run with (Accelerator, Stream).</param>
public static void Run(Action<Accelerator, AcceleratorStream> action)
{
var lease = GetKernel(useGpuOrder: true);
+65 -737
View File
@@ -8,7 +8,6 @@ namespace PerfectNumbers.Core.Gpu;
public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
{
public ulong High;
public ulong Low;
private const int NativeModuloChunkBits = 8;
@@ -21,690 +20,8 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1);
private const int Pow2WindowMaxExponent = (1 << Pow2WindowSizeBits) - 1;
private struct Pow2OddPowerTable
{
public GpuUInt128 Element0;
public GpuUInt128 Element1;
public GpuUInt128 Element2;
public GpuUInt128 Element3;
public GpuUInt128 Element4;
public GpuUInt128 Element5;
public GpuUInt128 Element6;
public GpuUInt128 Element7;
public GpuUInt128 Element8;
public GpuUInt128 Element9;
public GpuUInt128 Element10;
public GpuUInt128 Element11;
public GpuUInt128 Element12;
public GpuUInt128 Element13;
public GpuUInt128 Element14;
public GpuUInt128 Element15;
public GpuUInt128 Element16;
public GpuUInt128 Element17;
public GpuUInt128 Element18;
public GpuUInt128 Element19;
public GpuUInt128 Element20;
public GpuUInt128 Element21;
public GpuUInt128 Element22;
public GpuUInt128 Element23;
public GpuUInt128 Element24;
public GpuUInt128 Element25;
public GpuUInt128 Element26;
public GpuUInt128 Element27;
public GpuUInt128 Element28;
public GpuUInt128 Element29;
public GpuUInt128 Element30;
public GpuUInt128 Element31;
public GpuUInt128 Element32;
public GpuUInt128 Element33;
public GpuUInt128 Element34;
public GpuUInt128 Element35;
public GpuUInt128 Element36;
public GpuUInt128 Element37;
public GpuUInt128 Element38;
public GpuUInt128 Element39;
public GpuUInt128 Element40;
public GpuUInt128 Element41;
public GpuUInt128 Element42;
public GpuUInt128 Element43;
public GpuUInt128 Element44;
public GpuUInt128 Element45;
public GpuUInt128 Element46;
public GpuUInt128 Element47;
public GpuUInt128 Element48;
public GpuUInt128 Element49;
public GpuUInt128 Element50;
public GpuUInt128 Element51;
public GpuUInt128 Element52;
public GpuUInt128 Element53;
public GpuUInt128 Element54;
public GpuUInt128 Element55;
public GpuUInt128 Element56;
public GpuUInt128 Element57;
public GpuUInt128 Element58;
public GpuUInt128 Element59;
public GpuUInt128 Element60;
public GpuUInt128 Element61;
public GpuUInt128 Element62;
public GpuUInt128 Element63;
public GpuUInt128 Element64;
public GpuUInt128 Element65;
public GpuUInt128 Element66;
public GpuUInt128 Element67;
public GpuUInt128 Element68;
public GpuUInt128 Element69;
public GpuUInt128 Element70;
public GpuUInt128 Element71;
public GpuUInt128 Element72;
public GpuUInt128 Element73;
public GpuUInt128 Element74;
public GpuUInt128 Element75;
public GpuUInt128 Element76;
public GpuUInt128 Element77;
public GpuUInt128 Element78;
public GpuUInt128 Element79;
public GpuUInt128 Element80;
public GpuUInt128 Element81;
public GpuUInt128 Element82;
public GpuUInt128 Element83;
public GpuUInt128 Element84;
public GpuUInt128 Element85;
public GpuUInt128 Element86;
public GpuUInt128 Element87;
public GpuUInt128 Element88;
public GpuUInt128 Element89;
public GpuUInt128 Element90;
public GpuUInt128 Element91;
public GpuUInt128 Element92;
public GpuUInt128 Element93;
public GpuUInt128 Element94;
public GpuUInt128 Element95;
public GpuUInt128 Element96;
public GpuUInt128 Element97;
public GpuUInt128 Element98;
public GpuUInt128 Element99;
public GpuUInt128 Element100;
public GpuUInt128 Element101;
public GpuUInt128 Element102;
public GpuUInt128 Element103;
public GpuUInt128 Element104;
public GpuUInt128 Element105;
public GpuUInt128 Element106;
public GpuUInt128 Element107;
public GpuUInt128 Element108;
public GpuUInt128 Element109;
public GpuUInt128 Element110;
public GpuUInt128 Element111;
public GpuUInt128 Element112;
public GpuUInt128 Element113;
public GpuUInt128 Element114;
public GpuUInt128 Element115;
public GpuUInt128 Element116;
public GpuUInt128 Element117;
public GpuUInt128 Element118;
public GpuUInt128 Element119;
public GpuUInt128 Element120;
public GpuUInt128 Element121;
public GpuUInt128 Element122;
public GpuUInt128 Element123;
public GpuUInt128 Element124;
public GpuUInt128 Element125;
public GpuUInt128 Element126;
public GpuUInt128 Element127;
public GpuUInt128 this[int index]
{
readonly get
{
return GetElement(index);
}
set
{
SetElement(index, value);
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private readonly GpuUInt128 GetElement(int index) => index switch
{
0 => Element0,
1 => Element1,
2 => Element2,
3 => Element3,
4 => Element4,
5 => Element5,
6 => Element6,
7 => Element7,
8 => Element8,
9 => Element9,
10 => Element10,
11 => Element11,
12 => Element12,
13 => Element13,
14 => Element14,
15 => Element15,
16 => Element16,
17 => Element17,
18 => Element18,
19 => Element19,
20 => Element20,
21 => Element21,
22 => Element22,
23 => Element23,
24 => Element24,
25 => Element25,
26 => Element26,
27 => Element27,
28 => Element28,
29 => Element29,
30 => Element30,
31 => Element31,
32 => Element32,
33 => Element33,
34 => Element34,
35 => Element35,
36 => Element36,
37 => Element37,
38 => Element38,
39 => Element39,
40 => Element40,
41 => Element41,
42 => Element42,
43 => Element43,
44 => Element44,
45 => Element45,
46 => Element46,
47 => Element47,
48 => Element48,
49 => Element49,
50 => Element50,
51 => Element51,
52 => Element52,
53 => Element53,
54 => Element54,
55 => Element55,
56 => Element56,
57 => Element57,
58 => Element58,
59 => Element59,
60 => Element60,
61 => Element61,
62 => Element62,
63 => Element63,
64 => Element64,
65 => Element65,
66 => Element66,
67 => Element67,
68 => Element68,
69 => Element69,
70 => Element70,
71 => Element71,
72 => Element72,
73 => Element73,
74 => Element74,
75 => Element75,
76 => Element76,
77 => Element77,
78 => Element78,
79 => Element79,
80 => Element80,
81 => Element81,
82 => Element82,
83 => Element83,
84 => Element84,
85 => Element85,
86 => Element86,
87 => Element87,
88 => Element88,
89 => Element89,
90 => Element90,
91 => Element91,
92 => Element92,
93 => Element93,
94 => Element94,
95 => Element95,
96 => Element96,
97 => Element97,
98 => Element98,
99 => Element99,
100 => Element100,
101 => Element101,
102 => Element102,
103 => Element103,
104 => Element104,
105 => Element105,
106 => Element106,
107 => Element107,
108 => Element108,
109 => Element109,
110 => Element110,
111 => Element111,
112 => Element112,
113 => Element113,
114 => Element114,
115 => Element115,
116 => Element116,
117 => Element117,
118 => Element118,
119 => Element119,
120 => Element120,
121 => Element121,
122 => Element122,
123 => Element123,
124 => Element124,
125 => Element125,
126 => Element126,
127 => Element127,
_ => Element0,
};
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void SetElement(int index, GpuUInt128 value)
{
switch (index)
{
case 0:
Element0 = value;
return;
case 1:
Element1 = value;
return;
case 2:
Element2 = value;
return;
case 3:
Element3 = value;
return;
case 4:
Element4 = value;
return;
case 5:
Element5 = value;
return;
case 6:
Element6 = value;
return;
case 7:
Element7 = value;
return;
case 8:
Element8 = value;
return;
case 9:
Element9 = value;
return;
case 10:
Element10 = value;
return;
case 11:
Element11 = value;
return;
case 12:
Element12 = value;
return;
case 13:
Element13 = value;
return;
case 14:
Element14 = value;
return;
case 15:
Element15 = value;
return;
case 16:
Element16 = value;
return;
case 17:
Element17 = value;
return;
case 18:
Element18 = value;
return;
case 19:
Element19 = value;
return;
case 20:
Element20 = value;
return;
case 21:
Element21 = value;
return;
case 22:
Element22 = value;
return;
case 23:
Element23 = value;
return;
case 24:
Element24 = value;
return;
case 25:
Element25 = value;
return;
case 26:
Element26 = value;
return;
case 27:
Element27 = value;
return;
case 28:
Element28 = value;
return;
case 29:
Element29 = value;
return;
case 30:
Element30 = value;
return;
case 31:
Element31 = value;
return;
case 32:
Element32 = value;
return;
case 33:
Element33 = value;
return;
case 34:
Element34 = value;
return;
case 35:
Element35 = value;
return;
case 36:
Element36 = value;
return;
case 37:
Element37 = value;
return;
case 38:
Element38 = value;
return;
case 39:
Element39 = value;
return;
case 40:
Element40 = value;
return;
case 41:
Element41 = value;
return;
case 42:
Element42 = value;
return;
case 43:
Element43 = value;
return;
case 44:
Element44 = value;
return;
case 45:
Element45 = value;
return;
case 46:
Element46 = value;
return;
case 47:
Element47 = value;
return;
case 48:
Element48 = value;
return;
case 49:
Element49 = value;
return;
case 50:
Element50 = value;
return;
case 51:
Element51 = value;
return;
case 52:
Element52 = value;
return;
case 53:
Element53 = value;
return;
case 54:
Element54 = value;
return;
case 55:
Element55 = value;
return;
case 56:
Element56 = value;
return;
case 57:
Element57 = value;
return;
case 58:
Element58 = value;
return;
case 59:
Element59 = value;
return;
case 60:
Element60 = value;
return;
case 61:
Element61 = value;
return;
case 62:
Element62 = value;
return;
case 63:
Element63 = value;
return;
case 64:
Element64 = value;
return;
case 65:
Element65 = value;
return;
case 66:
Element66 = value;
return;
case 67:
Element67 = value;
return;
case 68:
Element68 = value;
return;
case 69:
Element69 = value;
return;
case 70:
Element70 = value;
return;
case 71:
Element71 = value;
return;
case 72:
Element72 = value;
return;
case 73:
Element73 = value;
return;
case 74:
Element74 = value;
return;
case 75:
Element75 = value;
return;
case 76:
Element76 = value;
return;
case 77:
Element77 = value;
return;
case 78:
Element78 = value;
return;
case 79:
Element79 = value;
return;
case 80:
Element80 = value;
return;
case 81:
Element81 = value;
return;
case 82:
Element82 = value;
return;
case 83:
Element83 = value;
return;
case 84:
Element84 = value;
return;
case 85:
Element85 = value;
return;
case 86:
Element86 = value;
return;
case 87:
Element87 = value;
return;
case 88:
Element88 = value;
return;
case 89:
Element89 = value;
return;
case 90:
Element90 = value;
return;
case 91:
Element91 = value;
return;
case 92:
Element92 = value;
return;
case 93:
Element93 = value;
return;
case 94:
Element94 = value;
return;
case 95:
Element95 = value;
return;
case 96:
Element96 = value;
return;
case 97:
Element97 = value;
return;
case 98:
Element98 = value;
return;
case 99:
Element99 = value;
return;
case 100:
Element100 = value;
return;
case 101:
Element101 = value;
return;
case 102:
Element102 = value;
return;
case 103:
Element103 = value;
return;
case 104:
Element104 = value;
return;
case 105:
Element105 = value;
return;
case 106:
Element106 = value;
return;
case 107:
Element107 = value;
return;
case 108:
Element108 = value;
return;
case 109:
Element109 = value;
return;
case 110:
Element110 = value;
return;
case 111:
Element111 = value;
return;
case 112:
Element112 = value;
return;
case 113:
Element113 = value;
return;
case 114:
Element114 = value;
return;
case 115:
Element115 = value;
return;
case 116:
Element116 = value;
return;
case 117:
Element117 = value;
return;
case 118:
Element118 = value;
return;
case 119:
Element119 = value;
return;
case 120:
Element120 = value;
return;
case 121:
Element121 = value;
return;
case 122:
Element122 = value;
return;
case 123:
Element123 = value;
return;
case 124:
Element124 = value;
return;
case 125:
Element125 = value;
return;
case 126:
Element126 = value;
return;
case 127:
Element127 = value;
return;
default:
// ILGPU kernels cannot throw exceptions, and callers guarantee the index range.
return;
}
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static int GetScalarBitLength(ulong value)
{
if (value == 0UL)
{
return 0;
}
return 64 - XMath.LeadingZeroCount(value);
}
private static int GetScalarBitLength(ulong value) => value == 0UL ? 0 : 64 - XMath.LeadingZeroCount(value);
public static readonly GpuUInt128 Zero = new();
public static readonly GpuUInt128 One = new(1UL);
@@ -743,23 +60,21 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public readonly int GetBitLength()
{
if (High != 0UL)
ulong temp = High;
if (temp != 0UL)
{
return 64 + (64 - BitOperations.LeadingZeroCount(High));
return 64 + (64 - BitOperations.LeadingZeroCount(temp));
}
ulong low = Low;
if (low == 0UL)
temp = Low;
if (temp == 0UL)
{
return 0;
}
return 64 - BitOperations.LeadingZeroCount(low);
return 64 - BitOperations.LeadingZeroCount(temp);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static int GetBitLength(GpuUInt128 value) => value.GetBitLength();
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static implicit operator GpuUInt128(UInt128 value) => new(value);
@@ -770,6 +85,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
public static implicit operator UInt128(GpuUInt128 value) =>
((UInt128)value.High << 64) | value.Low;
// TODO: The operators modifying the state of left is intuitive and can easily lead to bugs. Change this back to return a new instance, but replace where possible with method calls, if we want to modify the state of the instance.
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GpuUInt128 operator +(GpuUInt128 left, GpuUInt128 right)
{
@@ -809,16 +125,17 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
ulong high = value.High;
ulong low = value.Low;
int shiftedValue;
if (shift >= 64)
{
int longShift = shift - 64;
value.High = low << longShift;
shiftedValue = shift - 64;
value.High = low << shiftedValue;
value.Low = 0UL;
return value;
}
int inverseShift = 64 - shift;
value.High = (high << shift) | (low >> inverseShift);
shiftedValue = 64 - shift;
value.High = (high << shift) | (low >> shiftedValue);
value.Low = low << shift;
return value;
}
@@ -835,16 +152,17 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
ulong high = value.High;
ulong low = value.Low;
int shiftedValue;
if (shift >= 64)
{
int longShift = shift - 64;
value.Low = longShift == 0 ? high : high >> longShift;
shiftedValue = shift - 64;
value.Low = shiftedValue == 0 ? high : high >> shiftedValue;
value.High = 0UL;
return value;
}
int inverseShift = 64 - shift;
value.Low = (low >> shift) | (high << inverseShift);
shiftedValue = 64 - shift;
value.Low = (low >> shift) | (high << shiftedValue);
value.High = high >> shift;
return value;
}
@@ -913,17 +231,20 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
public void Add(GpuUInt128 other)
{
ulong originalLow = Low;
// We're reusing originalLow for carry to avoid an extra local.
ulong low = originalLow + other.Low;
ulong carry = low < originalLow ? 1UL : 0UL;
High = High + other.High + carry;
originalLow = low < originalLow ? 1UL : 0UL;
Low = low;
// We're reusing low for the final sum to avoid an extra local and to keep max stack frame at 3 instead of 4 with one statement.
low = High + other.High;
High = low + originalLow;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Add(ulong value)
{
ulong low = Low + value;
High += (low < Low ? 1UL : 0UL);
High += low < Low ? 1UL : 0UL;
Low = low;
}
@@ -962,10 +283,11 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GpuUInt128 Pow2Minus1Mod(ulong exponent, GpuUInt128 modulus)
{
if (modulus.IsZero)
{
return Zero;
}
// This should never happen in production code
// if (modulus.IsZero)
// {
// return Zero;
// }
GpuUInt128 pow = Pow2Mod(exponent, modulus);
if (pow.IsZero)
@@ -1034,30 +356,34 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong MulHigh(ulong x, ulong y)
{
// TODO: Review if we can reuse any other existing variables to reduce registry pressure.
ulong xLow = (uint)x;
ulong yLow = (uint)y;
ulong lowProduct = xLow * yLow;
ulong yHigh = y >> 32;
ulong cross = xLow * yHigh;
// We're reusing xLow for the cross product to avoid an extra local.
xLow *= yHigh;
ulong xHigh = x >> 32;
ulong result = xHigh * yHigh;
ulong temp = xHigh * yLow;
// We're reusing yHigh for the result to void an extra local.
yHigh *= xHigh;
// We're reusing yLow for the other cross product to avoid an extra local.
yLow *= xHigh;
// Keeping the wide partial sum in a dedicated local prevents the JIT from
// materialising it on the stack before the final carry propagation. The
// additional store looks redundant in C#, but it shortens the generated
// instruction sequence by avoiding an extra temporary and results in a
// measurable throughput win in the MulHigh benchmarks.
result += cross >> 32;
result += temp >> 32;
yHigh += (xLow >> 32);
// We've more instructions with this but the max method stack frame is 3 instead of 4.
yHigh += (yLow >> 32);
lowProduct >>= 32;
lowProduct += (uint)cross;
lowProduct += (uint)temp;
result += lowProduct >> 32;
return result;
lowProduct += (uint)xLow + (uint)yLow;
yHigh += lowProduct >> 32;
return yHigh;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -1099,23 +425,25 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GpuUInt128 Pow2Mod(ulong exponent, GpuUInt128 modulus)
{
if (modulus.IsZero || modulus == One)
{
return new();
}
// This should never happen in production code.
// if (modulus.IsZero || modulus == One)
// {
// return new();
// }
if (exponent == 0UL)
{
return new GpuUInt128(1UL);
}
// This should never happen in production code.
// if (exponent == 0UL)
// {
// return new GpuUInt128(1UL);
// }
Pow2OddPowerTable oddPowers = default;
GpuUInt128 power = new GpuUInt128(1UL);
GpuUInt128[] oddPowers = new GpuUInt128[PerfectNumberConstants.MaxOddPowersCount];
GpuUInt128 power = new(1UL);
int oddIndex = 0;
for (int bit = 1; bit <= Pow2WindowMaxExponent; bit++)
{
power += power;
power.Add(power);
if (power.CompareTo(modulus) >= 0)
{
power.Sub(modulus);
@@ -1128,7 +456,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
}
}
GpuUInt128 result = new GpuUInt128(1UL);
GpuUInt128 result = GpuUInt128.One;
int bitLength = GetScalarBitLength(exponent);
int index = bitLength - 1;
@@ -2144,14 +1472,14 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
}
}
// TODO: Check if the TODO below is still relevant.
// TODO: Check if the TODO below is still relevant.
// TODO(MOD-OPT): Montgomery/Barrett integration plan
// - Introduce caches of modulus-dependent constants:
// * Montgomery: n' (-(n^{-1}) mod 2^64 or 2^128), R2 = (R^2 mod n)
// * Barrett: mu = floor(2^k / n) for k ∈ {128, 192, 256}
// - Add fast-path for 64-bit NTT primes (modulus.High == 0UL) using pure 64-bit ops.
// - Expose helpers to retrieve/calc constants once per modulus and reuse in kernels.
// - Wire these into MulMod and SquareMod hot paths under feature toggles.
// - Ensure ILGPU compatibility (no BigInteger, no % inside kernels).
// TODO(MOD-OPT): Montgomery/Barrett integration plan
// - Introduce caches of modulus-dependent constants:
// * Montgomery: n' (-(n^{-1}) mod 2^64 or 2^128), R2 = (R^2 mod n)
// * Barrett: mu = floor(2^k / n) for k ∈ {128, 192, 256}
// - Add fast-path for 64-bit NTT primes (modulus.High == 0UL) using pure 64-bit ops.
// - Expose helpers to retrieve/calc constants once per modulus and reuse in kernels.
// - Wire these into MulMod and SquareMod hot paths under feature toggles.
// - Ensure ILGPU compatibility (no BigInteger, no % inside kernels).
@@ -375,7 +375,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
for (int i = 0; i < batchSize; i++)
{
divisorDataSpan[i] = MontgomeryDivisorDataCache.Get(divisorSpan[i]);
divisorDataSpan[i] = MontgomeryDivisorData.FromModulus(divisorSpan[i]);
}
divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize);
@@ -538,7 +538,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
_disposed = false;
}
public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits)
public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits)
{
if (_disposed)
{
@@ -551,7 +551,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
return;
}
MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor);
if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL)
{
hits.Clear();
@@ -0,0 +1,681 @@
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core.Gpu;
internal static partial class PrimeOrderGpuHeuristics
{
internal readonly struct Pow2OddPowerTableGpuUInt128
{
public Pow2OddPowerTableGpuUInt128(GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount)
{
Element0 = baseValue;
if (oddPowerCount == 1)
{
return;
}
// Reusing baseValue to hold base^2 for the shared odd-power ladder that follows.
baseValue.MulMod(baseValue, modulus);
// TODO: We can calculate baseValue % modulus before loop and use it to increase ladderEntry calculation speed - we'll reuse the base for incremental calculations.
GpuUInt128 current = baseValue;
// We're manually assigning each field to prevent the compiler to initialize each field twice due to auto-initialization. We're using the action to lower the code base size.
for (int i = 1; i < oddPowerCount; i++)
{
current.MulMod(baseValue, modulus);
switch (i)
{
case 1:
Element1 = current;
break;
case 2:
Element2 = current;
break;
case 3:
Element3 = current;
break;
case 4:
Element4 = current;
break;
case 5:
Element5 = current;
break;
case 6:
Element6 = current;
break;
case 7:
Element7 = current;
break;
case 8:
Element8 = current;
break;
case 9:
Element9 = current;
break;
case 10:
Element10 = current;
break;
case 11:
Element11 = current;
break;
case 12:
Element12 = current;
break;
case 13:
Element13 = current;
break;
case 14:
Element14 = current;
break;
case 15:
Element15 = current;
break;
case 16:
Element16 = current;
break;
case 17:
Element17 = current;
break;
case 18:
Element18 = current;
break;
case 19:
Element19 = current;
break;
case 20:
Element20 = current;
break;
case 21:
Element21 = current;
break;
case 22:
Element22 = current;
break;
case 23:
Element23 = current;
break;
case 24:
Element24 = current;
break;
case 25:
Element25 = current;
break;
case 26:
Element26 = current;
break;
case 27:
Element27 = current;
break;
case 28:
Element28 = current;
break;
case 29:
Element29 = current;
break;
case 30:
Element30 = current;
break;
case 31:
Element31 = current;
break;
case 32:
Element32 = current;
break;
case 33:
Element33 = current;
break;
case 34:
Element34 = current;
break;
case 35:
Element35 = current;
break;
case 36:
Element36 = current;
break;
case 37:
Element37 = current;
break;
case 38:
Element38 = current;
break;
case 39:
Element39 = current;
break;
case 40:
Element40 = current;
break;
case 41:
Element41 = current;
break;
case 42:
Element42 = current;
break;
case 43:
Element43 = current;
break;
case 44:
Element44 = current;
break;
case 45:
Element45 = current;
break;
case 46:
Element46 = current;
break;
case 47:
Element47 = current;
break;
case 48:
Element48 = current;
break;
case 49:
Element49 = current;
break;
case 50:
Element50 = current;
break;
case 51:
Element51 = current;
break;
case 52:
Element52 = current;
break;
case 53:
Element53 = current;
break;
case 54:
Element54 = current;
break;
case 55:
Element55 = current;
break;
case 56:
Element56 = current;
break;
case 57:
Element57 = current;
break;
case 58:
Element58 = current;
break;
case 59:
Element59 = current;
break;
case 60:
Element60 = current;
break;
case 61:
Element61 = current;
break;
case 62:
Element62 = current;
break;
case 63:
Element63 = current;
break;
case 64:
Element64 = current;
break;
case 65:
Element65 = current;
break;
case 66:
Element66 = current;
break;
case 67:
Element67 = current;
break;
case 68:
Element68 = current;
break;
case 69:
Element69 = current;
break;
case 70:
Element70 = current;
break;
case 71:
Element71 = current;
break;
case 72:
Element72 = current;
break;
case 73:
Element73 = current;
break;
case 74:
Element74 = current;
break;
case 75:
Element75 = current;
break;
case 76:
Element76 = current;
break;
case 77:
Element77 = current;
break;
case 78:
Element78 = current;
break;
case 79:
Element79 = current;
break;
case 80:
Element80 = current;
break;
case 81:
Element81 = current;
break;
case 82:
Element82 = current;
break;
case 83:
Element83 = current;
break;
case 84:
Element84 = current;
break;
case 85:
Element85 = current;
break;
case 86:
Element86 = current;
break;
case 87:
Element87 = current;
break;
case 88:
Element88 = current;
break;
case 89:
Element89 = current;
break;
case 90:
Element90 = current;
break;
case 91:
Element91 = current;
break;
case 92:
Element92 = current;
break;
case 93:
Element93 = current;
break;
case 94:
Element94 = current;
break;
case 95:
Element95 = current;
break;
case 96:
Element96 = current;
break;
case 97:
Element97 = current;
break;
case 98:
Element98 = current;
break;
case 99:
Element99 = current;
break;
case 100:
Element100 = current;
break;
case 101:
Element101 = current;
break;
case 102:
Element102 = current;
break;
case 103:
Element103 = current;
break;
case 104:
Element104 = current;
break;
case 105:
Element105 = current;
break;
case 106:
Element106 = current;
break;
case 107:
Element107 = current;
break;
case 108:
Element108 = current;
break;
case 109:
Element109 = current;
break;
case 110:
Element110 = current;
break;
case 111:
Element111 = current;
break;
case 112:
Element112 = current;
break;
case 113:
Element113 = current;
break;
case 114:
Element114 = current;
break;
case 115:
Element115 = current;
break;
case 116:
Element116 = current;
break;
case 117:
Element117 = current;
break;
case 118:
Element118 = current;
break;
case 119:
Element119 = current;
break;
case 120:
Element120 = current;
break;
case 121:
Element121 = current;
break;
case 122:
Element122 = current;
break;
case 123:
Element123 = current;
break;
case 124:
Element124 = current;
break;
case 125:
Element125 = current;
break;
case 126:
Element126 = current;
break;
case 127:
Element127 = current;
break;
}
}
}
public readonly GpuUInt128 Element0;
public readonly GpuUInt128 Element1;
public readonly GpuUInt128 Element2;
public readonly GpuUInt128 Element3;
public readonly GpuUInt128 Element4;
public readonly GpuUInt128 Element5;
public readonly GpuUInt128 Element6;
public readonly GpuUInt128 Element7;
public readonly GpuUInt128 Element8;
public readonly GpuUInt128 Element9;
public readonly GpuUInt128 Element10;
public readonly GpuUInt128 Element11;
public readonly GpuUInt128 Element12;
public readonly GpuUInt128 Element13;
public readonly GpuUInt128 Element14;
public readonly GpuUInt128 Element15;
public readonly GpuUInt128 Element16;
public readonly GpuUInt128 Element17;
public readonly GpuUInt128 Element18;
public readonly GpuUInt128 Element19;
public readonly GpuUInt128 Element20;
public readonly GpuUInt128 Element21;
public readonly GpuUInt128 Element22;
public readonly GpuUInt128 Element23;
public readonly GpuUInt128 Element24;
public readonly GpuUInt128 Element25;
public readonly GpuUInt128 Element26;
public readonly GpuUInt128 Element27;
public readonly GpuUInt128 Element28;
public readonly GpuUInt128 Element29;
public readonly GpuUInt128 Element30;
public readonly GpuUInt128 Element31;
public readonly GpuUInt128 Element32;
public readonly GpuUInt128 Element33;
public readonly GpuUInt128 Element34;
public readonly GpuUInt128 Element35;
public readonly GpuUInt128 Element36;
public readonly GpuUInt128 Element37;
public readonly GpuUInt128 Element38;
public readonly GpuUInt128 Element39;
public readonly GpuUInt128 Element40;
public readonly GpuUInt128 Element41;
public readonly GpuUInt128 Element42;
public readonly GpuUInt128 Element43;
public readonly GpuUInt128 Element44;
public readonly GpuUInt128 Element45;
public readonly GpuUInt128 Element46;
public readonly GpuUInt128 Element47;
public readonly GpuUInt128 Element48;
public readonly GpuUInt128 Element49;
public readonly GpuUInt128 Element50;
public readonly GpuUInt128 Element51;
public readonly GpuUInt128 Element52;
public readonly GpuUInt128 Element53;
public readonly GpuUInt128 Element54;
public readonly GpuUInt128 Element55;
public readonly GpuUInt128 Element56;
public readonly GpuUInt128 Element57;
public readonly GpuUInt128 Element58;
public readonly GpuUInt128 Element59;
public readonly GpuUInt128 Element60;
public readonly GpuUInt128 Element61;
public readonly GpuUInt128 Element62;
public readonly GpuUInt128 Element63;
public readonly GpuUInt128 Element64;
public readonly GpuUInt128 Element65;
public readonly GpuUInt128 Element66;
public readonly GpuUInt128 Element67;
public readonly GpuUInt128 Element68;
public readonly GpuUInt128 Element69;
public readonly GpuUInt128 Element70;
public readonly GpuUInt128 Element71;
public readonly GpuUInt128 Element72;
public readonly GpuUInt128 Element73;
public readonly GpuUInt128 Element74;
public readonly GpuUInt128 Element75;
public readonly GpuUInt128 Element76;
public readonly GpuUInt128 Element77;
public readonly GpuUInt128 Element78;
public readonly GpuUInt128 Element79;
public readonly GpuUInt128 Element80;
public readonly GpuUInt128 Element81;
public readonly GpuUInt128 Element82;
public readonly GpuUInt128 Element83;
public readonly GpuUInt128 Element84;
public readonly GpuUInt128 Element85;
public readonly GpuUInt128 Element86;
public readonly GpuUInt128 Element87;
public readonly GpuUInt128 Element88;
public readonly GpuUInt128 Element89;
public readonly GpuUInt128 Element90;
public readonly GpuUInt128 Element91;
public readonly GpuUInt128 Element92;
public readonly GpuUInt128 Element93;
public readonly GpuUInt128 Element94;
public readonly GpuUInt128 Element95;
public readonly GpuUInt128 Element96;
public readonly GpuUInt128 Element97;
public readonly GpuUInt128 Element98;
public readonly GpuUInt128 Element99;
public readonly GpuUInt128 Element100;
public readonly GpuUInt128 Element101;
public readonly GpuUInt128 Element102;
public readonly GpuUInt128 Element103;
public readonly GpuUInt128 Element104;
public readonly GpuUInt128 Element105;
public readonly GpuUInt128 Element106;
public readonly GpuUInt128 Element107;
public readonly GpuUInt128 Element108;
public readonly GpuUInt128 Element109;
public readonly GpuUInt128 Element110;
public readonly GpuUInt128 Element111;
public readonly GpuUInt128 Element112;
public readonly GpuUInt128 Element113;
public readonly GpuUInt128 Element114;
public readonly GpuUInt128 Element115;
public readonly GpuUInt128 Element116;
public readonly GpuUInt128 Element117;
public readonly GpuUInt128 Element118;
public readonly GpuUInt128 Element119;
public readonly GpuUInt128 Element120;
public readonly GpuUInt128 Element121;
public readonly GpuUInt128 Element122;
public readonly GpuUInt128 Element123;
public readonly GpuUInt128 Element124;
public readonly GpuUInt128 Element125;
public readonly GpuUInt128 Element126;
public readonly GpuUInt128 Element127;
public readonly GpuUInt128 this[int index]
{
[MethodImpl(MethodImplOptions.AggressiveInlining)]
get
{
return index switch
{
0 => Element0,
1 => Element1,
2 => Element2,
3 => Element3,
4 => Element4,
5 => Element5,
6 => Element6,
7 => Element7,
8 => Element8,
9 => Element9,
10 => Element10,
11 => Element11,
12 => Element12,
13 => Element13,
14 => Element14,
15 => Element15,
16 => Element16,
17 => Element17,
18 => Element18,
19 => Element19,
20 => Element20,
21 => Element21,
22 => Element22,
23 => Element23,
24 => Element24,
25 => Element25,
26 => Element26,
27 => Element27,
28 => Element28,
29 => Element29,
30 => Element30,
31 => Element31,
32 => Element32,
33 => Element33,
34 => Element34,
35 => Element35,
36 => Element36,
37 => Element37,
38 => Element38,
39 => Element39,
40 => Element40,
41 => Element41,
42 => Element42,
43 => Element43,
44 => Element44,
45 => Element45,
46 => Element46,
47 => Element47,
48 => Element48,
49 => Element49,
50 => Element50,
51 => Element51,
52 => Element52,
53 => Element53,
54 => Element54,
55 => Element55,
56 => Element56,
57 => Element57,
58 => Element58,
59 => Element59,
60 => Element60,
61 => Element61,
62 => Element62,
63 => Element63,
64 => Element64,
65 => Element65,
66 => Element66,
67 => Element67,
68 => Element68,
69 => Element69,
70 => Element70,
71 => Element71,
72 => Element72,
73 => Element73,
74 => Element74,
75 => Element75,
76 => Element76,
77 => Element77,
78 => Element78,
79 => Element79,
80 => Element80,
81 => Element81,
82 => Element82,
83 => Element83,
84 => Element84,
85 => Element85,
86 => Element86,
87 => Element87,
88 => Element88,
89 => Element89,
90 => Element90,
91 => Element91,
92 => Element92,
93 => Element93,
94 => Element94,
95 => Element95,
96 => Element96,
97 => Element97,
98 => Element98,
99 => Element99,
100 => Element100,
101 => Element101,
102 => Element102,
103 => Element103,
104 => Element104,
105 => Element105,
106 => Element106,
107 => Element107,
108 => Element108,
109 => Element109,
110 => Element110,
111 => Element111,
112 => Element112,
113 => Element113,
114 => Element114,
115 => Element115,
116 => Element116,
117 => Element117,
118 => Element118,
119 => Element119,
120 => Element120,
121 => Element121,
122 => Element122,
123 => Element123,
124 => Element124,
125 => Element125,
126 => Element126,
127 => Element127,
_ => Element0,
};
}
}
}
}
File diff suppressed because it is too large. Load diff
@@ -20,6 +20,6 @@ public interface IMersenneNumberDivisorByDivisorTester
public interface IDivisorScanSession : IDisposable
{
void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits);
void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits);
}
}
File diff suppressed because it is too large. Load diff
+24 -23
View File
@@ -14,29 +14,6 @@ public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong
public readonly ulong MontgomeryTwo = montgomeryTwo;
public readonly ulong MontgomeryTwoSquared = montgomeryTwoSquared;
}
internal static class MontgomeryDivisorDataCache
{
private static readonly ConcurrentDictionary<ulong, MontgomeryDivisorData> Cache = new();
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static MontgomeryDivisorData Get(ulong modulus) => Cache.GetOrAdd(modulus, static m => Create(m));
private static MontgomeryDivisorData Create(ulong modulus)
{
ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64x2, modulus);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return new MontgomeryDivisorData(
modulus,
nPrime,
montgomeryOne,
montgomeryTwo,
montgomeryTwoSquared);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong ComputeMontgomeryResidue(UInt128 value, ulong modulus) => (ulong)(value % modulus);
@@ -53,4 +30,28 @@ internal static class MontgomeryDivisorDataCache
inv *= unchecked(2UL - modulus * inv);
return unchecked(0UL - inv);
}
public static MontgomeryDivisorData FromModulus(ulong modulus)
{
ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64x2, modulus);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return new MontgomeryDivisorData(
modulus,
nPrime,
montgomeryOne,
montgomeryTwo,
montgomeryTwoSquared);
}
}
internal static class MontgomeryDivisorDataCache
{
private static readonly ConcurrentDictionary<ulong, MontgomeryDivisorData> Cache = new();
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static MontgomeryDivisorData Get(ulong modulus) => Cache.GetOrAdd(modulus, static m => MontgomeryDivisorData.FromModulus(m));
}
@@ -7,6 +7,7 @@ public static class PerfectNumberConstants
public const ulong ExtraDivisorCycleSearchLimit = 64UL;
public const int MaxQForDivisorCycles = 4_000_000;
public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000;
public const int MaxOddPowersCount = 128;
// TODO: Load these limits from the benchmark-driven configuration so CPU and GPU scans stay aligned with the optimal
// divisor-cycle datasets we generate offline.
// TODO: Promote these magic numbers into a runtime profile derived from EvenPerfectBitScanner.Benchmarks so we can retune
@@ -32,6 +32,7 @@ internal static partial class PrimeOrderCalculator
PrimeOrderHeuristicDevice device)
{
using var scope = UsePow2Mode(device);
MontgomeryDivisorData divisorData;
if (prime <= ulong.MaxValue)
{
ulong? previous = null;
@@ -48,14 +49,20 @@ internal static partial class PrimeOrderCalculator
}
}
PrimeOrderResult result = Calculate((ulong)prime, previous, config, device);
ulong prime64 = (ulong)prime;
divisorData = MontgomeryDivisorData.FromModulus(prime64);
PrimeOrderResult result = Calculate(prime64, previous, divisorData, config, device);
return new PrimeOrderResultWide(result.Status, result.Order);
}
else
{
divisorData = default;
}
return CalculateWideInternal(prime, previousOrder, config);
return CalculateWideInternal(prime, previousOrder, divisorData, config);
}
private static PrimeOrderResultWide CalculateWideInternal(UInt128 prime, UInt128? previousOrder, PrimeOrderSearchConfig config)
private static PrimeOrderResultWide CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
{
if (prime <= UInt128.One)
{
@@ -72,53 +79,54 @@ internal static partial class PrimeOrderCalculator
PartialFactorResult128 phiFactors = PartialFactorWide(phi, config);
if (phiFactors.Factors is null)
{
return FinishStrictlyWide(prime, config.Mode);
return FinishStrictlyWide(prime, divisorData, config.Mode);
}
return RunHeuristicPipelineWide(prime, previousOrder, config, phi, phiFactors);
return RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
}
private static PrimeOrderResultWide RunHeuristicPipelineWide(
UInt128 prime,
UInt128? previousOrder,
in MontgomeryDivisorData divisorData,
PrimeOrderSearchConfig config,
UInt128 phi,
PartialFactorResult128 phiFactors)
{
if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, phiFactors))
if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, divisorData, phiFactors))
{
return new PrimeOrderResultWide(PrimeOrderStatus.Found, phi);
}
UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi);
candidateOrder = ExponentLoweringWide(candidateOrder, prime, phiFactors);
UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi, divisorData);
candidateOrder = ExponentLoweringWide(candidateOrder, prime, divisorData, phiFactors);
if (TryConfirmOrderWide(prime, candidateOrder, config))
if (TryConfirmOrderWide(prime, candidateOrder, divisorData, config))
{
return new PrimeOrderResultWide(PrimeOrderStatus.Found, candidateOrder);
}
if (config.Mode == PrimeOrderMode.Strict)
{
return FinishStrictlyWide(prime, PrimeOrderMode.Strict);
return FinishStrictlyWide(prime, divisorData, PrimeOrderMode.Strict);
}
if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, config, out UInt128 order))
if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order))
{
return new PrimeOrderResultWide(PrimeOrderStatus.Found, order);
}
return FinishStrictlyWide(prime, config.Mode);
return FinishStrictlyWide(prime, divisorData, config.Mode);
}
private static PrimeOrderResultWide FinishStrictlyWide(UInt128 prime, PrimeOrderMode mode)
private static PrimeOrderResultWide FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData, PrimeOrderMode mode)
{
UInt128 strictOrder = CalculateByFactorizationWide(prime);
UInt128 strictOrder = CalculateByFactorizationWide(prime, divisorData);
PrimeOrderStatus status = mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved;
return new PrimeOrderResultWide(status, strictOrder);
}
private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, PartialFactorResult128 factors)
private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors)
{
ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors;
int length = factors.Count;
@@ -126,7 +134,7 @@ internal static partial class PrimeOrderCalculator
{
UInt128 factor = factorSpan[i].Value;
UInt128 reduced = phi / factor;
if (Pow2ModWide(reduced, prime) == UInt128.One)
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
{
return false;
}
@@ -135,14 +143,14 @@ internal static partial class PrimeOrderCalculator
return true;
}
private static UInt128 InitializeStartingOrderWide(UInt128 prime, UInt128 phi)
private static UInt128 InitializeStartingOrderWide(UInt128 prime, UInt128 phi, in MontgomeryDivisorData divisorData)
{
UInt128 order = phi;
UInt128 mod8 = prime & (UInt128)7UL;
if (mod8 == UInt128.One || mod8 == (UInt128)7UL)
{
UInt128 half = phi >> 1;
if (Pow2ModWide(half, prime) == UInt128.One)
if (Pow2ModWide(half, prime, divisorData) == UInt128.One)
{
order = half;
}
@@ -151,7 +159,7 @@ internal static partial class PrimeOrderCalculator
return order;
}
private static UInt128 ExponentLoweringWide(UInt128 order, UInt128 prime, PartialFactorResult128 factors)
private static UInt128 ExponentLoweringWide(UInt128 order, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors)
{
FactorEntry128[]? tempArray = null;
try
@@ -184,7 +192,7 @@ internal static partial class PrimeOrderCalculator
if ((order % primeFactor) == UInt128.Zero)
{
UInt128 reduced = order / primeFactor;
if (Pow2ModWide(reduced, prime) == UInt128.One)
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
{
order = reduced;
continue;
@@ -206,14 +214,14 @@ internal static partial class PrimeOrderCalculator
}
}
private static bool TryConfirmOrderWide(UInt128 prime, UInt128 order, PrimeOrderSearchConfig config)
private static bool TryConfirmOrderWide(UInt128 prime, UInt128 order, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
{
if (order == UInt128.Zero)
{
return false;
}
if (Pow2ModWide(order, prime) != UInt128.One)
if (Pow2ModWide(order, prime, divisorData) != UInt128.One)
{
return false;
}
@@ -253,7 +261,7 @@ internal static partial class PrimeOrderCalculator
}
reduced /= primeFactor;
if (Pow2ModWide(reduced, prime) == UInt128.One)
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
{
return false;
}
@@ -267,6 +275,7 @@ internal static partial class PrimeOrderCalculator
UInt128 prime,
UInt128 order,
UInt128? previousOrder,
in MontgomeryDivisorData divisorData,
PrimeOrderSearchConfig config,
out UInt128 result)
{
@@ -322,12 +331,12 @@ internal static partial class PrimeOrderCalculator
UInt128 candidate = candidates[i];
powUsed++;
if (Pow2ModWide(candidate, prime) != UInt128.One)
if (Pow2ModWide(candidate, prime, divisorData) != UInt128.One)
{
continue;
}
if (!TryConfirmCandidateWide(prime, candidate, config, ref powUsed, powBudget))
if (!TryConfirmCandidateWide(prime, candidate, config, ref powUsed, powBudget, divisorData))
{
continue;
}
@@ -511,7 +520,7 @@ internal static partial class PrimeOrderCalculator
}
}
private static bool TryConfirmCandidateWide(UInt128 prime, UInt128 candidate, PrimeOrderSearchConfig config, ref int powUsed, int powBudget)
private static bool TryConfirmCandidateWide(UInt128 prime, UInt128 candidate, PrimeOrderSearchConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData divisorData)
{
PartialFactorResult128 factorization = PartialFactorWide(candidate, config);
if (factorization.Factors is null)
@@ -554,7 +563,7 @@ internal static partial class PrimeOrderCalculator
}
powUsed++;
if (Pow2ModWide(reduced, prime) == UInt128.One)
if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One)
{
return false;
}
@@ -826,7 +835,7 @@ internal static partial class PrimeOrderCalculator
return 64 + BitOperations.TrailingZeroCount(high);
}
private static UInt128 CalculateByFactorizationWide(UInt128 prime)
private static UInt128 CalculateByFactorizationWide(in UInt128 prime, in MontgomeryDivisorData divisorData)
{
UInt128 phi = prime - UInt128.One;
Dictionary<UInt128, int> counts = new(capacity: 8);
@@ -850,7 +859,7 @@ internal static partial class PrimeOrderCalculator
if ((order % primeFactor) != UInt128.Zero)
{
UInt128 candidate = order / primeFactor;
if (Pow2ModWide(candidate, prime) == UInt128.One)
if (Pow2ModWide(candidate, prime, divisorData) == UInt128.One)
{
order = candidate;
}
@@ -996,7 +1005,7 @@ internal static partial class PrimeOrderCalculator
return true;
}
private static UInt128 Pow2ModWide(UInt128 exponent, UInt128 modulus)
private static UInt128 Pow2ModWide(UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData divisorData)
{
if (modulus == UInt128.One)
{
@@ -1009,7 +1018,7 @@ internal static partial class PrimeOrderCalculator
{
ulong prime64 = (ulong)modulus;
ulong exponent64 = (ulong)exponent;
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent64, prime64, out ulong remainder);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent64, prime64, out ulong remainder, divisorData);
if (status == GpuPow2ModStatus.Success)
{
return remainder;
+10 -10
View File
@@ -146,20 +146,21 @@ internal static partial class PrimeOrderCalculator
}
}
public static PrimeOrderResult Calculate(ulong prime, ulong? previousOrder, PrimeOrderSearchConfig config)
=> Calculate(prime, previousOrder, config, PrimeOrderHeuristicDevice.Gpu);
public static PrimeOrderResult Calculate(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
=> Calculate(prime, previousOrder, divisorData, config, PrimeOrderHeuristicDevice.Gpu);
public static PrimeOrderResult Calculate(
ulong prime,
ulong? previousOrder,
in MontgomeryDivisorData divisorData,
PrimeOrderSearchConfig config,
PrimeOrderHeuristicDevice device)
{
using var scope = UsePow2Mode(device);
return CalculateInternal(prime, previousOrder, config);
return CalculateInternal(prime, previousOrder, divisorData, config);
}
private static PrimeOrderResult CalculateInternal(ulong prime, ulong? previousOrder, PrimeOrderSearchConfig config)
private static PrimeOrderResult CalculateInternal(ulong prime, ulong? previousOrder, MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
{
if (prime <= 3UL)
{
@@ -167,7 +168,6 @@ internal static partial class PrimeOrderCalculator
}
ulong phi = prime - 1UL;
MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(prime);
if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out PrimeOrderResult gpuResult))
{
@@ -448,8 +448,8 @@ internal static partial class PrimeOrderCalculator
DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
int index = 0;
const int MaxGpuBatchSize = 256;
const int StackGpuBatchSize = 32;
const int MaxGpuBatchSize = 1024;
const int StackGpuBatchSize = 1024;
Span<ulong> stackGpuRemainders = stackalloc ulong[StackGpuBatchSize];
while (index < candidateCount && powUsed < powBudget)
{
@@ -473,7 +473,7 @@ internal static partial class PrimeOrderCalculator
if (batchSize <= StackGpuBatchSize)
{
Span<ulong> localRemainders = stackGpuRemainders.Slice(0, batchSize);
status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders);
status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData);
if (status == GpuPow2ModStatus.Success)
{
gpuSuccess = true;
@@ -484,7 +484,7 @@ internal static partial class PrimeOrderCalculator
{
gpuPool = ArrayPool<ulong>.Shared.Rent(batchSize);
Span<ulong> pooledRemainders = gpuPool.AsSpan(0, batchSize);
status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders);
status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData);
if (status == GpuPow2ModStatus.Success)
{
pooledGpuRemainders = pooledRemainders;
@@ -775,7 +775,7 @@ internal static partial class PrimeOrderCalculator
{
if (IsGpuPow2Allowed)
{
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder);
GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData);
if (status == GpuPow2ModStatus.Success)
{
return remainder == 1UL;
+4
View File
@@ -1,3 +1,5 @@
using PerfectNumbers.Core.Gpu;
namespace PerfectNumbers.Core
{
public static class UInt128Numbers
@@ -9,6 +11,8 @@ namespace PerfectNumbers.Core
public static readonly UInt128 OneShiftedLeft64 = UInt128.One << 64;
public static readonly UInt128 OneShiftedLeft64x2 = Two * OneShiftedLeft64;
public static readonly GpuUInt128 TwoGpu = new(2UL);
// TODO: Populate this table from the UInt128 residue benchmarks so GPU kernels can pull precomputed constants
// without re-deriving them at runtime.
// TODO: Extend the table with cached Mod3/Mod5 folding constants so both CPU and GPU residue helpers can reuse
+6 -678
View File
@@ -12,680 +12,6 @@ public static class ULongExtensions
private const ulong Pow2WindowFallbackThreshold = 32UL;
private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1);
private struct Pow2OddPowerTable
{
public ulong Element0;
public ulong Element1;
public ulong Element2;
public ulong Element3;
public ulong Element4;
public ulong Element5;
public ulong Element6;
public ulong Element7;
public ulong Element8;
public ulong Element9;
public ulong Element10;
public ulong Element11;
public ulong Element12;
public ulong Element13;
public ulong Element14;
public ulong Element15;
public ulong Element16;
public ulong Element17;
public ulong Element18;
public ulong Element19;
public ulong Element20;
public ulong Element21;
public ulong Element22;
public ulong Element23;
public ulong Element24;
public ulong Element25;
public ulong Element26;
public ulong Element27;
public ulong Element28;
public ulong Element29;
public ulong Element30;
public ulong Element31;
public ulong Element32;
public ulong Element33;
public ulong Element34;
public ulong Element35;
public ulong Element36;
public ulong Element37;
public ulong Element38;
public ulong Element39;
public ulong Element40;
public ulong Element41;
public ulong Element42;
public ulong Element43;
public ulong Element44;
public ulong Element45;
public ulong Element46;
public ulong Element47;
public ulong Element48;
public ulong Element49;
public ulong Element50;
public ulong Element51;
public ulong Element52;
public ulong Element53;
public ulong Element54;
public ulong Element55;
public ulong Element56;
public ulong Element57;
public ulong Element58;
public ulong Element59;
public ulong Element60;
public ulong Element61;
public ulong Element62;
public ulong Element63;
public ulong Element64;
public ulong Element65;
public ulong Element66;
public ulong Element67;
public ulong Element68;
public ulong Element69;
public ulong Element70;
public ulong Element71;
public ulong Element72;
public ulong Element73;
public ulong Element74;
public ulong Element75;
public ulong Element76;
public ulong Element77;
public ulong Element78;
public ulong Element79;
public ulong Element80;
public ulong Element81;
public ulong Element82;
public ulong Element83;
public ulong Element84;
public ulong Element85;
public ulong Element86;
public ulong Element87;
public ulong Element88;
public ulong Element89;
public ulong Element90;
public ulong Element91;
public ulong Element92;
public ulong Element93;
public ulong Element94;
public ulong Element95;
public ulong Element96;
public ulong Element97;
public ulong Element98;
public ulong Element99;
public ulong Element100;
public ulong Element101;
public ulong Element102;
public ulong Element103;
public ulong Element104;
public ulong Element105;
public ulong Element106;
public ulong Element107;
public ulong Element108;
public ulong Element109;
public ulong Element110;
public ulong Element111;
public ulong Element112;
public ulong Element113;
public ulong Element114;
public ulong Element115;
public ulong Element116;
public ulong Element117;
public ulong Element118;
public ulong Element119;
public ulong Element120;
public ulong Element121;
public ulong Element122;
public ulong Element123;
public ulong Element124;
public ulong Element125;
public ulong Element126;
public ulong Element127;
public ulong this[int index]
{
get
{
return GetElement(index);
}
set
{
SetElement(index, value);
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private ulong GetElement(int index) => index switch
{
0 => Element0,
1 => Element1,
2 => Element2,
3 => Element3,
4 => Element4,
5 => Element5,
6 => Element6,
7 => Element7,
8 => Element8,
9 => Element9,
10 => Element10,
11 => Element11,
12 => Element12,
13 => Element13,
14 => Element14,
15 => Element15,
16 => Element16,
17 => Element17,
18 => Element18,
19 => Element19,
20 => Element20,
21 => Element21,
22 => Element22,
23 => Element23,
24 => Element24,
25 => Element25,
26 => Element26,
27 => Element27,
28 => Element28,
29 => Element29,
30 => Element30,
31 => Element31,
32 => Element32,
33 => Element33,
34 => Element34,
35 => Element35,
36 => Element36,
37 => Element37,
38 => Element38,
39 => Element39,
40 => Element40,
41 => Element41,
42 => Element42,
43 => Element43,
44 => Element44,
45 => Element45,
46 => Element46,
47 => Element47,
48 => Element48,
49 => Element49,
50 => Element50,
51 => Element51,
52 => Element52,
53 => Element53,
54 => Element54,
55 => Element55,
56 => Element56,
57 => Element57,
58 => Element58,
59 => Element59,
60 => Element60,
61 => Element61,
62 => Element62,
63 => Element63,
64 => Element64,
65 => Element65,
66 => Element66,
67 => Element67,
68 => Element68,
69 => Element69,
70 => Element70,
71 => Element71,
72 => Element72,
73 => Element73,
74 => Element74,
75 => Element75,
76 => Element76,
77 => Element77,
78 => Element78,
79 => Element79,
80 => Element80,
81 => Element81,
82 => Element82,
83 => Element83,
84 => Element84,
85 => Element85,
86 => Element86,
87 => Element87,
88 => Element88,
89 => Element89,
90 => Element90,
91 => Element91,
92 => Element92,
93 => Element93,
94 => Element94,
95 => Element95,
96 => Element96,
97 => Element97,
98 => Element98,
99 => Element99,
100 => Element100,
101 => Element101,
102 => Element102,
103 => Element103,
104 => Element104,
105 => Element105,
106 => Element106,
107 => Element107,
108 => Element108,
109 => Element109,
110 => Element110,
111 => Element111,
112 => Element112,
113 => Element113,
114 => Element114,
115 => Element115,
116 => Element116,
117 => Element117,
118 => Element118,
119 => Element119,
120 => Element120,
121 => Element121,
122 => Element122,
123 => Element123,
124 => Element124,
125 => Element125,
126 => Element126,
127 => Element127,
_ => Element0,
};
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void SetElement(int index, ulong value)
{
switch (index)
{
case 0:
Element0 = value;
return;
case 1:
Element1 = value;
return;
case 2:
Element2 = value;
return;
case 3:
Element3 = value;
return;
case 4:
Element4 = value;
return;
case 5:
Element5 = value;
return;
case 6:
Element6 = value;
return;
case 7:
Element7 = value;
return;
case 8:
Element8 = value;
return;
case 9:
Element9 = value;
return;
case 10:
Element10 = value;
return;
case 11:
Element11 = value;
return;
case 12:
Element12 = value;
return;
case 13:
Element13 = value;
return;
case 14:
Element14 = value;
return;
case 15:
Element15 = value;
return;
case 16:
Element16 = value;
return;
case 17:
Element17 = value;
return;
case 18:
Element18 = value;
return;
case 19:
Element19 = value;
return;
case 20:
Element20 = value;
return;
case 21:
Element21 = value;
return;
case 22:
Element22 = value;
return;
case 23:
Element23 = value;
return;
case 24:
Element24 = value;
return;
case 25:
Element25 = value;
return;
case 26:
Element26 = value;
return;
case 27:
Element27 = value;
return;
case 28:
Element28 = value;
return;
case 29:
Element29 = value;
return;
case 30:
Element30 = value;
return;
case 31:
Element31 = value;
return;
case 32:
Element32 = value;
return;
case 33:
Element33 = value;
return;
case 34:
Element34 = value;
return;
case 35:
Element35 = value;
return;
case 36:
Element36 = value;
return;
case 37:
Element37 = value;
return;
case 38:
Element38 = value;
return;
case 39:
Element39 = value;
return;
case 40:
Element40 = value;
return;
case 41:
Element41 = value;
return;
case 42:
Element42 = value;
return;
case 43:
Element43 = value;
return;
case 44:
Element44 = value;
return;
case 45:
Element45 = value;
return;
case 46:
Element46 = value;
return;
case 47:
Element47 = value;
return;
case 48:
Element48 = value;
return;
case 49:
Element49 = value;
return;
case 50:
Element50 = value;
return;
case 51:
Element51 = value;
return;
case 52:
Element52 = value;
return;
case 53:
Element53 = value;
return;
case 54:
Element54 = value;
return;
case 55:
Element55 = value;
return;
case 56:
Element56 = value;
return;
case 57:
Element57 = value;
return;
case 58:
Element58 = value;
return;
case 59:
Element59 = value;
return;
case 60:
Element60 = value;
return;
case 61:
Element61 = value;
return;
case 62:
Element62 = value;
return;
case 63:
Element63 = value;
return;
case 64:
Element64 = value;
return;
case 65:
Element65 = value;
return;
case 66:
Element66 = value;
return;
case 67:
Element67 = value;
return;
case 68:
Element68 = value;
return;
case 69:
Element69 = value;
return;
case 70:
Element70 = value;
return;
case 71:
Element71 = value;
return;
case 72:
Element72 = value;
return;
case 73:
Element73 = value;
return;
case 74:
Element74 = value;
return;
case 75:
Element75 = value;
return;
case 76:
Element76 = value;
return;
case 77:
Element77 = value;
return;
case 78:
Element78 = value;
return;
case 79:
Element79 = value;
return;
case 80:
Element80 = value;
return;
case 81:
Element81 = value;
return;
case 82:
Element82 = value;
return;
case 83:
Element83 = value;
return;
case 84:
Element84 = value;
return;
case 85:
Element85 = value;
return;
case 86:
Element86 = value;
return;
case 87:
Element87 = value;
return;
case 88:
Element88 = value;
return;
case 89:
Element89 = value;
return;
case 90:
Element90 = value;
return;
case 91:
Element91 = value;
return;
case 92:
Element92 = value;
return;
case 93:
Element93 = value;
return;
case 94:
Element94 = value;
return;
case 95:
Element95 = value;
return;
case 96:
Element96 = value;
return;
case 97:
Element97 = value;
return;
case 98:
Element98 = value;
return;
case 99:
Element99 = value;
return;
case 100:
Element100 = value;
return;
case 101:
Element101 = value;
return;
case 102:
Element102 = value;
return;
case 103:
Element103 = value;
return;
case 104:
Element104 = value;
return;
case 105:
Element105 = value;
return;
case 106:
Element106 = value;
return;
case 107:
Element107 = value;
return;
case 108:
Element108 = value;
return;
case 109:
Element109 = value;
return;
case 110:
Element110 = value;
return;
case 111:
Element111 = value;
return;
case 112:
Element112 = value;
return;
case 113:
Element113 = value;
return;
case 114:
Element114 = value;
return;
case 115:
Element115 = value;
return;
case 116:
Element116 = value;
return;
case 117:
Element117 = value;
return;
case 118:
Element118 = value;
return;
case 119:
Element119 = value;
return;
case 120:
Element120 = value;
return;
case 121:
Element121 = value;
return;
case 122:
Element122 = value;
return;
case 123:
Element123 = value;
return;
case 124:
Element124 = value;
return;
case 125:
Element125 = value;
return;
case 126:
Element126 = value;
return;
case 127:
Element127 = value;
return;
default:
// ILGPU kernels cannot throw exceptions, and callers guarantee the index range.
return;
}
}
}
public static ulong CalculateOrder(this ulong q)
{
if (q <= 2UL)
@@ -1019,8 +345,7 @@ public static class ULongExtensions
ulong result = divisor.MontgomeryOne;
ulong nPrime = divisor.NPrime;
Pow2OddPowerTable oddPowers = default;
InitializeMontgomeryOddPowers(divisor, modulus, nPrime, ref oddPowers, oddPowerCount);
ulong[] oddPowers = InitializeMontgomeryOddPowers(divisor, modulus, nPrime, oddPowerCount);
int index = bitLength - 1;
while (index >= 0)
@@ -1136,12 +461,13 @@ public static class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void InitializeMontgomeryOddPowers(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, ref Pow2OddPowerTable oddPowers, int oddPowerCount)
private static ulong[] InitializeMontgomeryOddPowers(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
{
ulong[] oddPowers = new ulong[PerfectNumberConstants.MaxOddPowersCount];
oddPowers[0] = divisor.MontgomeryTwo;
if (oddPowerCount == 1)
{
return;
return oddPowers;
}
ulong square = divisor.MontgomeryTwoSquared;
@@ -1150,6 +476,8 @@ public static class ULongExtensions
ulong previous = oddPowers[i - 1];
oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime);
}
return oddPowers;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]