333 lines
11 KiB
C#
333 lines
11 KiB
C#
using System.Runtime.CompilerServices;
|
|
using PerfectNumbers.Core.Gpu;
|
|
|
|
namespace PerfectNumbers.Core.Cpu;
|
|
|
|
public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester
|
|
{
|
|
private readonly object _sync = new();
|
|
|
|
private ulong _divisorLimit;
|
|
private ulong _lastStatusDivisor;
|
|
private int _batchSize = 1_024;
|
|
|
|
[ThreadStatic]
|
|
private static GpuUInt128WorkSet _divisorScanGpuWorkSet;
|
|
|
|
private struct GpuUInt128WorkSet
|
|
{
|
|
public GpuUInt128 Step;
|
|
public GpuUInt128 Divisor;
|
|
public GpuUInt128 Limit;
|
|
}
|
|
|
|
public int BatchSize
|
|
{
|
|
get => _batchSize;
|
|
set => _batchSize = Math.Max(1, value);
|
|
}
|
|
|
|
public void ConfigureFromMaxPrime(ulong maxPrime)
|
|
{
|
|
_divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
|
|
_lastStatusDivisor = 0UL;
|
|
}
|
|
|
|
public ulong DivisorLimit
|
|
{
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
get
|
|
{
|
|
return _divisorLimit;
|
|
}
|
|
}
|
|
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
public ulong GetAllowedMaxDivisor(ulong prime) => ComputeAllowedMaxDivisor(prime, _divisorLimit);
|
|
|
|
public bool IsPrime(ulong prime, out bool divisorsExhausted)
|
|
{
|
|
ulong allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit);
|
|
|
|
// The CPU by-divisor run always hands us primes with enormous divisor limits, so the fallback below never executes.
|
|
// if (allowedMax < 3UL)
|
|
// {
|
|
// // EvenPerfectBitScanner routes primes below the small-divisor cutoff to the GPU path, so the CPU path still sees
|
|
// // trivial candidates during targeted tests. Short-circuit here to keep those runs aligned with the production flow.
|
|
// divisorsExhausted = true;
|
|
// return true;
|
|
// }
|
|
|
|
ulong processedCount;
|
|
bool processedAll;
|
|
|
|
bool composite = CheckDivisors(
|
|
prime,
|
|
allowedMax,
|
|
out processedAll,
|
|
out processedCount);
|
|
|
|
if (processedCount > 0UL)
|
|
{
|
|
lock (_sync)
|
|
{
|
|
UpdateStatusUnsafe(processedCount);
|
|
}
|
|
}
|
|
|
|
if (composite)
|
|
{
|
|
divisorsExhausted = true;
|
|
return false;
|
|
}
|
|
|
|
divisorsExhausted = processedAll || composite;
|
|
return true;
|
|
}
|
|
|
|
public void PrepareCandidates(in ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues)
|
|
{
|
|
ulong divisorLimit = _divisorLimit;
|
|
int length = primes.Length;
|
|
for (int index = 0; index < length; index++)
|
|
{
|
|
allowedMaxValues[index] = ComputeAllowedMaxDivisor(primes[index], divisorLimit);
|
|
}
|
|
}
|
|
|
|
public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession()
|
|
{
|
|
MersenneCpuDivisorScanSession? session = ThreadStaticPools.RentMersenneCpuDivisorSession();
|
|
if (session is not null)
|
|
{
|
|
return session;
|
|
}
|
|
|
|
return new MersenneCpuDivisorScanSession();
|
|
}
|
|
|
|
private bool CheckDivisors(
|
|
ulong prime,
|
|
ulong allowedMax,
|
|
out bool processedAll,
|
|
out ulong processedCount)
|
|
{
|
|
processedCount = 0UL;
|
|
processedAll = false;
|
|
|
|
// The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs.
|
|
// Keeping the guard commented out documents the reasoning for benchmarks and tests.
|
|
// if (allowedMax < 3UL)
|
|
// {
|
|
// return false;
|
|
// }
|
|
|
|
ref GpuUInt128WorkSet workSet = ref _divisorScanGpuWorkSet;
|
|
|
|
ref GpuUInt128 step = ref workSet.Step;
|
|
step.High = 0UL;
|
|
step.Low = prime;
|
|
step.ShiftLeft(1);
|
|
|
|
ref GpuUInt128 limit = ref workSet.Limit;
|
|
limit.High = 0UL;
|
|
limit.Low = allowedMax;
|
|
|
|
ref GpuUInt128 divisor = ref workSet.Divisor;
|
|
divisor.High = step.High;
|
|
divisor.Low = step.Low;
|
|
divisor.Add(1UL);
|
|
if (divisor.CompareTo(limit) > 0)
|
|
{
|
|
processedAll = true;
|
|
divisor.High = 0UL;
|
|
divisor.Low = 0UL;
|
|
step.High = 0UL;
|
|
step.Low = 0UL;
|
|
limit.High = 0UL;
|
|
limit.Low = 0UL;
|
|
return false;
|
|
}
|
|
|
|
FactorCacheLease factorCacheLease = new FactorCacheLease();
|
|
Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? factorCache = null;
|
|
|
|
DivisorCycleCache cycleCache = DivisorCycleCache.Shared;
|
|
|
|
ulong stepHigh = step.High;
|
|
ulong stepLow = step.Low;
|
|
byte step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL);
|
|
byte step8 = (byte)(stepLow % 8UL);
|
|
byte step5 = (byte)(((stepHigh % 5UL) + (stepLow % 5UL)) % 5UL);
|
|
byte step3 = (byte)(((stepHigh % 3UL) + (stepLow % 3UL)) % 3UL);
|
|
byte step7 = (byte)((((stepHigh % 7UL) * 2UL) + (stepLow % 7UL)) % 7UL);
|
|
byte step11 = (byte)((((stepHigh % 11UL) * 5UL) + (stepLow % 11UL)) % 11UL);
|
|
|
|
ulong divisorHigh = divisor.High;
|
|
ulong divisorLow = divisor.Low;
|
|
byte remainder10 = (byte)((((divisorHigh % 10UL) * 6UL) + (divisorLow % 10UL)) % 10UL);
|
|
byte remainder8 = (byte)(divisorLow % 8UL);
|
|
byte remainder5 = (byte)(((divisorHigh % 5UL) + (divisorLow % 5UL)) % 5UL);
|
|
byte remainder3 = (byte)(((divisorHigh % 3UL) + (divisorLow % 3UL)) % 3UL);
|
|
byte remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL);
|
|
byte remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL);
|
|
|
|
// Keep the divisibility filters aligned with the divisor-cycle generator so the
|
|
// CPU path never requests cycles that were skipped during cache creation.
|
|
bool lastIsSeven = (prime & 3UL) == 3UL;
|
|
|
|
while (divisor.CompareTo(limit) <= 0)
|
|
{
|
|
ulong candidate = divisor.Low;
|
|
processedCount++;
|
|
|
|
bool admissible = lastIsSeven
|
|
? (remainder10 == 3 || remainder10 == 7 || remainder10 == 9)
|
|
: (remainder10 == 1 || remainder10 == 3 || remainder10 == 9);
|
|
|
|
if (admissible && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0)
|
|
{
|
|
MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(candidate);
|
|
ulong divisorCycle;
|
|
// Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs.
|
|
// if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles)
|
|
// {
|
|
// divisorCycle = cycleCache.GetCycleLength(candidate);
|
|
// }
|
|
// else
|
|
{
|
|
if (factorCache is null)
|
|
{
|
|
factorCacheLease.EnsureInitialized(ref factorCache);
|
|
}
|
|
|
|
if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(
|
|
candidate,
|
|
prime,
|
|
divisorData,
|
|
factorCache,
|
|
out ulong computedCycle,
|
|
out bool primeOrderFailed) || computedCycle == 0UL)
|
|
{
|
|
divisorCycle = cycleCache.GetCycleLength(candidate, skipPrimeOrderHeuristic: primeOrderFailed);
|
|
}
|
|
else
|
|
{
|
|
divisorCycle = computedCycle;
|
|
}
|
|
}
|
|
|
|
if (divisorCycle == prime)
|
|
{
|
|
// A cycle equal to the tested exponent (which is prime in this path) guarantees that the candidate divides
|
|
// the corresponding Mersenne number because the order of 2 modulo the divisor is exactly p.
|
|
factorCacheLease.Dispose();
|
|
factorCache = null;
|
|
divisor.High = 0UL;
|
|
divisor.Low = 0UL;
|
|
step.High = 0UL;
|
|
step.Low = 0UL;
|
|
limit.High = 0UL;
|
|
limit.Low = 0UL;
|
|
processedAll = true;
|
|
return true;
|
|
}
|
|
|
|
if (divisorCycle == 0UL)
|
|
{
|
|
Console.WriteLine($"Divisor cycle was not calculated for {prime}");
|
|
}
|
|
}
|
|
|
|
divisor.Add(step);
|
|
remainder10 = AddMod(remainder10, step10, (byte)10);
|
|
remainder8 = AddMod(remainder8, step8, (byte)8);
|
|
remainder5 = AddMod(remainder5, step5, (byte)5);
|
|
remainder3 = AddMod(remainder3, step3, (byte)3);
|
|
remainder7 = AddMod(remainder7, step7, (byte)7);
|
|
remainder11 = AddMod(remainder11, step11, (byte)11);
|
|
}
|
|
|
|
factorCacheLease.Dispose();
|
|
factorCache = null; // The pooled dictionary was returned to the thread-static slot.
|
|
|
|
processedAll = divisor.CompareTo(limit) > 0;
|
|
divisor.High = 0UL;
|
|
divisor.Low = 0UL;
|
|
step.High = 0UL;
|
|
step.Low = 0UL;
|
|
limit.High = 0UL;
|
|
limit.Low = 0UL;
|
|
return false;
|
|
}
|
|
|
|
private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData)
|
|
{
|
|
ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData);
|
|
return residue == 1UL ? (byte)1 : (byte)0;
|
|
}
|
|
|
|
private static byte AddMod(byte value, byte delta, byte modulus)
|
|
{
|
|
int sum = value + delta;
|
|
if (sum >= modulus)
|
|
{
|
|
sum -= modulus;
|
|
}
|
|
|
|
return (byte)sum;
|
|
}
|
|
|
|
private void UpdateStatusUnsafe(ulong processedCount)
|
|
{
|
|
if (processedCount == 0UL)
|
|
{
|
|
return;
|
|
}
|
|
|
|
ulong interval = PerfectNumberConstants.ConsoleInterval;
|
|
if (interval == 0UL)
|
|
{
|
|
_lastStatusDivisor = 0UL;
|
|
return;
|
|
}
|
|
|
|
ulong total = _lastStatusDivisor + processedCount;
|
|
// TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
|
|
// status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly.
|
|
_lastStatusDivisor = total % interval;
|
|
}
|
|
|
|
private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
|
|
{
|
|
// The by-divisor CPU configuration only supplies primes greater than 1, so the guard below never trips.
|
|
// if (maxPrime <= 1UL)
|
|
// {
|
|
// return 0UL;
|
|
// }
|
|
if (maxPrime - 1UL >= 64UL)
|
|
{
|
|
return ulong.MaxValue;
|
|
}
|
|
|
|
return (1UL << (int)(maxPrime - 1UL)) - 1UL;
|
|
}
|
|
|
|
private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit)
|
|
{
|
|
// Production by-divisor scans only handle primes, so inputs never fall below 2.
|
|
// if (prime <= 1UL)
|
|
// {
|
|
// return 0UL;
|
|
// }
|
|
if (prime - 1UL >= 64UL)
|
|
{
|
|
return divisorLimit;
|
|
}
|
|
|
|
return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
|
|
}
|
|
|
|
|
|
}
|