Restore cycle-aware divisor scans and background generation

This commit is contained in:
Marek Lemańczyk committed 2025-09-25 02:12:40 +02:00
1 parent bc4e526946
commit b798970ce9
6 files changed
+1501 -1747

No files matched your search

@@ -1,5 +1,4 @@
using System;
using System.Buffers;
using System.Collections.Concurrent;
using System.Runtime.CompilerServices;
using PerfectNumbers.Core;
@@ -8,556 +7,368 @@ using MontgomeryDivisorData = PerfectNumbers.Core.Gpu.MersenneNumberDivisorByDiv
namespace PerfectNumbers.Core.Cpu;
// TODO: The current implementation is over-complicated for CPU. We shouldn't need all the hassle with intermediate buffers, checking the status there etc.
public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester
{
private readonly object _sync = new();
private readonly ConcurrentBag<DivisorScanSession> _sessionPool = new();
private ulong _divisorLimit;
private ulong _lastStatusDivisor;
private bool _isConfigured;
private bool _useDivisorCycles;
private int _batchSize = 1_024;
public bool UseDivisorCycles
{
get => _useDivisorCycles;
set => _useDivisorCycles = value;
}
public int BatchSize
{
get => _batchSize;
set => _batchSize = Math.Max(1, value);
}
public void ConfigureFromMaxPrime(ulong maxPrime)
{
lock (_sync)
{
_divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
_lastStatusDivisor = 0UL;
_isConfigured = true;
}
}
public ulong DivisorLimit
{
get
{
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
return _divisorLimit;
}
}
}
public ulong GetAllowedMaxDivisor(ulong prime)
{
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
return ComputeAllowedMaxDivisor(prime, _divisorLimit);
}
}
public bool IsPrime(ulong prime, out bool divisorsExhausted)
{
ulong allowedMax;
bool useCycles;
int batchCapacity;
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit);
useCycles = _useDivisorCycles;
batchCapacity = _batchSize;
}
if (allowedMax < 3UL)
{
divisorsExhausted = true;
return true;
}
ulong[] divisors = ArrayPool<ulong>.Shared.Rent(batchCapacity);
byte[] hits = ArrayPool<byte>.Shared.Rent(batchCapacity);
MontgomeryDivisorData[] divisorData = ArrayPool<MontgomeryDivisorData>.Shared.Rent(batchCapacity);
ulong processedCount = 0UL;
ulong lastProcessed = 0UL;
bool composite = false;
bool processedAll = false;
try
{
composite = CheckDivisors(
prime,
allowedMax,
useCycles,
divisors,
hits,
divisorData,
out lastProcessed,
out processedAll,
out processedCount);
}
finally
{
ArrayPool<ulong>.Shared.Return(divisors, clearArray: false);
ArrayPool<byte>.Shared.Return(hits, clearArray: false);
ArrayPool<MontgomeryDivisorData>.Shared.Return(divisorData, clearArray: false);
}
if (processedCount > 0UL)
{
lock (_sync)
{
UpdateStatusUnsafe(lastProcessed, processedCount);
}
}
if (composite)
{
divisorsExhausted = true;
return false;
}
divisorsExhausted = processedAll || composite;
return true;
}
public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession()
{
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
if (_sessionPool.TryTake(out DivisorScanSession? session))
{
session.Reset();
return session;
}
return new DivisorScanSession(this);
}
}
private void ReturnSession(DivisorScanSession session)
{
_sessionPool.Add(session);
}
private bool CheckDivisors(
ulong prime,
ulong allowedMax,
bool useCycles,
ulong[] divisors,
byte[] hits,
MontgomeryDivisorData[] divisorData,
out ulong lastProcessed,
out bool processedAll,
out ulong processedCount)
{
lastProcessed = 0UL;
processedCount = 0UL;
processedAll = false;
if (allowedMax < 3UL)
{
return false;
}
int batchCapacity = _batchSize;
ulong divisor = 3UL;
bool composite = false;
while (divisor <= allowedMax)
{
int batchSize = 0;
bool reachedEndInBatch = false;
while (batchSize < batchCapacity && divisor <= allowedMax)
{
ulong currentDivisor = divisor;
ulong nextDivisor = currentDivisor + 2UL;
divisors[batchSize++] = currentDivisor;
processedCount++;
lastProcessed = currentDivisor;
if (nextDivisor <= currentDivisor)
{
reachedEndInBatch = true;
break;
}
if (nextDivisor > allowedMax)
{
divisor = nextDivisor;
reachedEndInBatch = true;
break;
}
divisor = nextDivisor;
}
if (batchSize == 0)
{
if (reachedEndInBatch)
{
processedAll = true;
}
break;
}
Span<MontgomeryDivisorData> divisorDataSpan = divisorData.AsSpan(0, batchSize);
Span<ulong> divisorSpan = divisors.AsSpan(0, batchSize);
Span<byte> hitsSpan = hits.AsSpan(0, batchSize);
for (int i = 0; i < batchSize; i++)
{
divisorDataSpan[i] = CreateMontgomeryDivisorData(divisorSpan[i]);
}
for (int i = 0; i < batchSize; i++)
{
hitsSpan[i] = CheckDivisor(prime, useCycles, divisorSpan[i], divisorDataSpan[i]);
if (hitsSpan[i] != 0)
{
composite = true;
lastProcessed = divisorSpan[i];
break;
}
}
if (!composite)
{
lastProcessed = divisorSpan[batchSize - 1];
}
if (composite)
{
break;
}
if (reachedEndInBatch)
{
processedAll = true;
break;
}
}
processedAll = processedAll || divisor > allowedMax;
return composite;
}
private static byte CheckDivisor(ulong prime, bool useCycles, ulong divisor, in MontgomeryDivisorData divisorData)
{
ulong modulus = divisorData.Modulus;
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
return 0;
}
ulong exponent = prime;
if (useCycles)
{
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(divisor);
if (cycle == 0UL)
{
return 0;
}
ulong remainder = prime % cycle;
if (remainder != 0UL)
{
return 0;
}
exponent = remainder;
}
return exponent.Pow2MontgomeryMod(divisorData) == 1UL ? (byte)1 : (byte)0;
}
private void UpdateStatusUnsafe(ulong lastProcessed, ulong processedCount)
{
if (processedCount == 0UL)
{
return;
}
ulong interval = PerfectNumberConstants.ConsoleInterval;
if (interval == 0UL)
{
_lastStatusDivisor = 0UL;
return;
}
ulong total = _lastStatusDivisor + processedCount;
_lastStatusDivisor = total % interval;
}
private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
{
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)
{
if (prime <= 1UL)
{
return 0UL;
}
if (prime - 1UL >= 64UL)
{
return divisorLimit;
}
return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
}
private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
{
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL);
}
return new MontgomeryDivisorData(
modulus,
ComputeMontgomeryNPrime(modulus),
ComputeMontgomeryResidue(1UL, modulus),
ComputeMontgomeryResidue(2UL, modulus));
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
private static ulong ComputeMontgomeryNPrime(ulong modulus)
{
ulong inv = modulus;
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
return unchecked(0UL - inv);
}
private sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
{
private readonly MersenneNumberDivisorByDivisorCpuTester _owner;
private ulong[] _primeBuffer;
private int[] _positionBuffer;
private ulong[] _exponentBuffer;
private int _capacity;
private bool _disposed;
internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner)
{
_owner = owner;
_capacity = Math.Max(1, owner._batchSize);
_primeBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
_positionBuffer = ArrayPool<int>.Shared.Rent(_capacity);
_exponentBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
}
internal void Reset()
{
_disposed = false;
}
public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits)
{
if (_disposed)
{
throw new ObjectDisposedException(nameof(DivisorScanSession));
}
int primesLength = primes.Length;
if (primesLength == 0)
{
return;
}
MersenneNumberDivisorByDivisorCpuTester owner = _owner;
int batchSize = owner._batchSize;
EnsureCapacity(batchSize);
ulong modulus = divisor;
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
return;
}
bool cycleEnabled = owner._useDivisorCycles && divisorCycle != 0UL;
ulong cycle = divisorCycle;
MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor);
int offset = 0;
Span<ulong> primeSpan = _primeBuffer.AsSpan();
Span<int> positionSpan = _positionBuffer.AsSpan();
Span<ulong> exponentSpan = _exponentBuffer.AsSpan();
ulong previousResidue = 0UL;
ulong previousPrime = 0UL;
bool hasPreviousResidue = false;
while (offset < primesLength)
{
int sliceLength = Math.Min(batchSize, primesLength - offset);
ReadOnlySpan<ulong> primesSlice = primes.Slice(offset, sliceLength);
Span<byte> hitsSlice = hits.Slice(offset, sliceLength);
hitsSlice.Clear();
int computeCount = 0;
if (cycleEnabled)
{
for (int i = 0; i < sliceLength; i++)
{
ulong prime = primesSlice[i];
if (prime % cycle != 0UL)
{
continue;
}
primeSpan[computeCount] = prime;
positionSpan[computeCount] = i;
computeCount++;
}
}
else
{
for (int i = 0; i < sliceLength; i++)
{
primeSpan[computeCount] = primesSlice[i];
positionSpan[computeCount] = i;
computeCount++;
}
}
if (computeCount > 0)
{
Span<ulong> exponentSlice = exponentSpan[..computeCount];
ulong currentPrime = hasPreviousResidue ? previousPrime : 0UL;
bool hasDeltaSource = hasPreviousResidue;
for (int i = 0; i < computeCount; i++)
{
ulong prime = primeSpan[i];
ulong exponentValue;
if (hasDeltaSource)
{
exponentValue = prime - currentPrime;
}
else
{
exponentValue = prime;
hasDeltaSource = true;
}
if (cycleEnabled)
{
exponentValue %= cycle;
}
exponentSlice[i] = exponentValue;
currentPrime = prime;
}
for (int i = 0; i < computeCount; i++)
{
exponentSlice[i] = exponentSlice[i].Pow2MontgomeryMod(divisorData);
}
for (int i = 0; i < computeCount; i++)
{
int position = positionSpan[i];
ulong stepResidue = exponentSlice[i];
ulong prime = primeSpan[i];
ulong residue = hasPreviousResidue ? MultiplyMod(previousResidue, stepResidue, modulus) : stepResidue;
hitsSlice[position] = residue == 1UL ? (byte)1 : (byte)0;
previousResidue = residue;
previousPrime = prime;
hasPreviousResidue = true;
}
}
offset += sliceLength;
}
}
public void Dispose()
{
if (_disposed)
{
return;
}
_disposed = true;
_owner.ReturnSession(this);
}
private void EnsureCapacity(int requiredCapacity)
{
if (requiredCapacity <= _capacity)
{
return;
}
ArrayPool<ulong>.Shared.Return(_primeBuffer, clearArray: false);
ArrayPool<int>.Shared.Return(_positionBuffer, clearArray: false);
ArrayPool<ulong>.Shared.Return(_exponentBuffer, clearArray: false);
_capacity = requiredCapacity;
_primeBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
_positionBuffer = ArrayPool<int>.Shared.Rent(_capacity);
_exponentBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
}
private static ulong MultiplyMod(ulong left, ulong right, ulong modulus)
{
if (modulus == 0UL)
{
return 0UL;
}
UInt128 product = (UInt128)left * right;
return (ulong)(product % modulus);
}
}
private readonly object _sync = new();
private readonly ConcurrentBag<DivisorScanSession> _sessionPool = new();
private ulong _divisorLimit;
private ulong _lastStatusDivisor;
private bool _isConfigured;
private bool _useDivisorCycles;
private int _batchSize = 1_024;
public bool UseDivisorCycles
{
get => _useDivisorCycles;
set => _useDivisorCycles = value;
}
public int BatchSize
{
get => _batchSize;
set => _batchSize = Math.Max(1, value);
}
public void ConfigureFromMaxPrime(ulong maxPrime)
{
lock (_sync)
{
_divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
_lastStatusDivisor = 0UL;
_isConfigured = true;
}
}
public ulong DivisorLimit
{
get
{
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
return _divisorLimit;
}
}
}
public ulong GetAllowedMaxDivisor(ulong prime)
{
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
return ComputeAllowedMaxDivisor(prime, _divisorLimit);
}
}
public bool IsPrime(ulong prime, out bool divisorsExhausted)
{
ulong allowedMax;
bool useCycles;
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit);
useCycles = _useDivisorCycles;
}
if (allowedMax < 3UL)
{
divisorsExhausted = true;
return true;
}
ulong processedCount = 0UL;
ulong lastProcessed = 0UL;
bool processedAll = false;
bool composite = CheckDivisors(
prime,
allowedMax,
useCycles,
out lastProcessed,
out processedAll,
out processedCount);
if (processedCount > 0UL)
{
lock (_sync)
{
UpdateStatusUnsafe(lastProcessed, processedCount);
}
}
if (composite)
{
divisorsExhausted = true;
return false;
}
divisorsExhausted = processedAll || composite;
return true;
}
public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession()
{
lock (_sync)
{
if (!_isConfigured)
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
}
if (_sessionPool.TryTake(out DivisorScanSession? session))
{
session.Reset();
return session;
}
return new DivisorScanSession(this);
}
}
private void ReturnSession(DivisorScanSession session)
{
_sessionPool.Add(session);
}
private bool CheckDivisors(
ulong prime,
ulong allowedMax,
bool useCycles,
out ulong lastProcessed,
out bool processedAll,
out ulong processedCount)
{
lastProcessed = 0UL;
processedCount = 0UL;
processedAll = false;
if (allowedMax < 3UL)
{
return false;
}
ulong divisor = 3UL;
while (divisor <= allowedMax)
{
processedCount++;
lastProcessed = divisor;
MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor);
byte hit;
if (useCycles)
{
using DivisorCycleCache.Lease cycleLease = DivisorCycleCache.Shared.Acquire(divisor);
ulong divisorCycle = cycleLease.GetCycle(divisor);
hit = CheckDivisor(prime, divisorCycle != 0UL, divisorCycle, divisorData);
}
else
{
hit = CheckDivisor(prime, false, 0UL, divisorData);
}
if (hit != 0)
{
processedAll = true;
return true;
}
if (divisor >= allowedMax - 1UL)
{
processedAll = true;
break;
}
if (divisor > ulong.MaxValue - 2UL)
{
break;
}
divisor += 2UL;
}
if (!processedAll)
{
processedAll = divisor > allowedMax;
}
return false;
}
private static byte CheckDivisor(ulong prime, bool useCycles, ulong divisorCycle, in MontgomeryDivisorData divisorData)
{
ulong modulus = divisorData.Modulus;
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
return 0;
}
ulong residue = useCycles && divisorCycle != 0UL
? prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData)
: prime.Pow2MontgomeryMod(divisorData);
return residue == 1UL ? (byte)1 : (byte)0;
}
private void UpdateStatusUnsafe(ulong lastProcessed, ulong processedCount)
{
if (processedCount == 0UL)
{
return;
}
ulong interval = PerfectNumberConstants.ConsoleInterval;
if (interval == 0UL)
{
_lastStatusDivisor = 0UL;
return;
}
ulong total = _lastStatusDivisor + processedCount;
_lastStatusDivisor = total % interval;
}
private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
{
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)
{
if (prime <= 1UL)
{
return 0UL;
}
if (prime - 1UL >= 64UL)
{
return divisorLimit;
}
return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
}
private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
{
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL);
}
return new MontgomeryDivisorData(
modulus,
ComputeMontgomeryNPrime(modulus),
ComputeMontgomeryResidue(1UL, modulus),
ComputeMontgomeryResidue(2UL, modulus));
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
private static ulong ComputeMontgomeryNPrime(ulong modulus)
{
ulong inv = modulus;
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
return unchecked(0UL - inv);
}
private sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
{
private readonly MersenneNumberDivisorByDivisorCpuTester _owner;
private bool _disposed;
internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner)
{
_owner = owner;
}
internal void Reset()
{
_disposed = false;
}
public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits)
{
if (_disposed)
{
throw new ObjectDisposedException(nameof(DivisorScanSession));
}
int length = primes.Length;
if (length == 0)
{
return;
}
ulong modulus = divisor;
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
hits.Clear();
return;
}
bool cycleEnabled = _owner._useDivisorCycles && divisorCycle != 0UL;
MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor);
for (int i = 0; i < length; i++)
{
ulong prime = primes[i];
if (cycleEnabled)
{
ulong remainder = prime % divisorCycle;
if (remainder != 0UL)
{
hits[i] = 0;
continue;
}
}
ulong residue = cycleEnabled
? prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData)
: prime.Pow2MontgomeryMod(divisorData);
hits[i] = residue == 1UL ? (byte)1 : (byte)0;
}
}
public void Dispose()
{
if (_disposed)
{
return;
}
_disposed = true;
_owner.ReturnSession(this);
}
}
}
+370 -420
View File
@@ -1,428 +1,378 @@
using System;
using System.Buffers;
using System.Collections.Concurrent;
using System.Runtime.CompilerServices;
using System.Threading.Tasks;
using ILGPU;
using ILGPU.Runtime;
using PerfectNumbers.Core.Gpu;
namespace PerfectNumbers.Core;
public sealed class DivisorCycleCache
{
internal sealed class CycleBlock
{
// TODO: This can be removed. See other comments.
private int _referenceCount;
internal CycleBlock(int index, ulong start, ulong[] cycles)
{
Index = index;
Start = start;
Cycles = cycles;
End = start + (ulong)cycles.Length - 1UL;
}
// TODO: Let's use fields instead of the properties below for better performance
internal int Index { get; }
internal ulong Start { get; }
internal ulong End { get; }
internal ulong[] Cycles { get; }
// TODO: We don't need full reference tracking. It's enough to assign divisor block to a task, instead. We can remove excessive blocks as required.
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal void AddRef()
{
Interlocked.Increment(ref _referenceCount);
}
// TODO: We don't need full reference tracking. It's enough to assign divisor block to a task, instead. We can remove excessive blocks as required.
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal void Release()
{
Interlocked.Decrement(ref _referenceCount);
}
// TODO: We don't need full reference tracking. It's enough to assign divisor block to a task, instead. This can be removed
internal bool IsInUse => Volatile.Read(ref _referenceCount) > 0;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal ulong GetCycle(ulong divisor)
{
// TODO: The production code should guarantee that this will never occur and we'll use divisors with their corresponding blocks, only. This can be removed and we can simplify this method to =>
if (divisor < Start || divisor > End)
{
return 0UL;
}
return Cycles[(int)(divisor - Start)];
}
}
public struct Lease : IDisposable
{
// TODO: We can probably remove _owner, if we don't need full cycle management which will help GC to clear these out
private readonly DivisorCycleCache? _owner;
private CycleBlock? _block;
internal Lease(DivisorCycleCache owner, CycleBlock block)
{
_owner = owner;
_block = block;
_block.AddRef();
}
// TODO: We're highly over-complicating this. _block is always assigned by the constructor. We don't need to make it nullable. This can be removed.
public readonly bool IsValid => _block is not null;
// TODO: Once we make _block non-nullable we can simplify these to just fields assigned during the creation for the best performance.
public readonly ulong Start => _block?.Start ?? 0UL;
public readonly ulong End => _block?.End ?? 0UL;
public readonly ulong[]? Values => _block?.Cycles;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public readonly ulong GetCycle(ulong divisor) => _block?.GetCycle(divisor) ?? 0UL;
public void Dispose()
{
// TODO: We don't need full cycle block management. We can rely on GC to do its job. Thanks to that we could possibly make the struct read-only.
CycleBlock? block = _block;
if (block is null)
{
return;
}
block.Release();
_owner?.Release(block);
_block = null;
}
}
private readonly object _sync = new();
private readonly ConcurrentDictionary<int, Task<CycleBlock>> _pending = new();
// TODO: _baseBlock should be like no other block. Why do we need a separate property for it? It should become current upon initialization
private CycleBlock _baseBlock = null!;
private CycleBlock? _previous;
private CycleBlock? _current;
private CycleBlock? _next;
private bool _initialized;
// TODO: Let's not use lazy but directly create the instance, instead for better performance. Then we can remove _lazy and just assign the instance to Shared.
private static readonly Lazy<DivisorCycleCache> _lazy = new(() => new DivisorCycleCache());
public static DivisorCycleCache Shared => _lazy.Value;
private DivisorCycleCache()
{
ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
Initialize(snapshot);
}
public void ReloadFromCurrentSnapshot()
{
ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
Initialize(snapshot);
}
private void Initialize(ulong[] snapshot)
{
lock (_sync)
{
_pending.Clear();
_baseBlock = new CycleBlock(index: 0, start: 0UL, snapshot);
_previous = _baseBlock;
_current = null;
_next = null;
_initialized = true;
StartPrefetchLocked(1);
}
}
public Lease Acquire(ulong divisor)
{
// TODO: We don't need _initialized. The constructor calls Initialize and makes sure it's always initialized, before use. This check can be removed, too.
if (!_initialized)
{
throw new InvalidOperationException("DivisorCycleCache must be initialized before use.");
}
// TODO: Do we need Lease instances, if we don't care who is using our blocks once they were given and just remove them internally from future use?
if (divisor <= _baseBlock.End)
{
return new Lease(this, _baseBlock);
}
CycleBlock block = AcquireDynamicBlock(divisor);
return new Lease(this, block);
}
private CycleBlock AcquireDynamicBlock(ulong divisor)
{
int blockIndex = GetBlockIndex(divisor);
while (true)
{
Task<CycleBlock>? taskToAwait = null;
lock (_sync)
{
if (TryGetCachedBlockLocked(blockIndex, out CycleBlock? cached))
{
return cached!;
}
if (!_pending.TryGetValue(blockIndex, out taskToAwait))
{
taskToAwait = Task.Run(() => GenerateBlock(blockIndex));
_pending[blockIndex] = taskToAwait;
}
}
CycleBlock block = taskToAwait.GetAwaiter().GetResult();
lock (_sync)
{
_pending.TryRemove(blockIndex, out _);
if (TryGetCachedBlockLocked(blockIndex, out CycleBlock? cachedAfterAwait))
{
return cachedAfterAwait!;
}
PromoteBlockLocked(block);
return block;
}
}
}
private bool TryGetCachedBlockLocked(int blockIndex, out CycleBlock? block)
{
if (_current is { Index: var currentIndex } current && currentIndex == blockIndex)
{
block = current;
return true;
}
if (_previous is { Index: var previousIndex } previous && previousIndex == blockIndex)
{
block = previous;
return true;
}
if (_next is { Index: var nextIndex } next && nextIndex == blockIndex)
{
PromoteNextLocked();
block = _current!;
return true;
}
block = null;
return false;
}
private void PromoteNextLocked()
{
if (_next is null)
{
return;
}
CycleBlock? oldPrevious = _previous;
_previous = _current ?? _baseBlock;
_current = _next;
_next = null;
StartPrefetchLocked(_current!.Index + 1);
TryReleaseBlockLocked(oldPrevious);
}
private void PromoteBlockLocked(CycleBlock block)
{
if (_current is null)
{
_previous = _baseBlock;
_current = block;
StartPrefetchLocked(block.Index + 1);
return;
}
if (block.Index > _current.Index)
{
CycleBlock? oldPrevious = _previous;
_previous = _current;
_current = block;
StartPrefetchLocked(block.Index + 1);
TryReleaseBlockLocked(oldPrevious);
return;
}
if (_previous is null || block.Index >= _previous.Index)
{
CycleBlock? oldPrevious = _previous;
_previous = block;
TryReleaseBlockLocked(oldPrevious);
}
}
private void StartPrefetchLocked(int blockIndex)
{
// TODO: Why do we need this check? This should never happen in production code.
if (blockIndex <= 0)
{
return;
}
// TODO: Why do we need this check? This should never happen in production code.
if (_next is { Index: var nextIndex } next && nextIndex == blockIndex)
{
return;
}
if (_pending.TryGetValue(blockIndex, out Task<CycleBlock>? existing))
{
if (existing.IsCompletedSuccessfully)
{
if (_next is null || _next.Index < blockIndex)
{
_next = existing.Result;
}
_pending.TryRemove(blockIndex, out _);
}
return;
}
Task<CycleBlock> task = Task.Run(() => GenerateBlock(blockIndex));
_pending[blockIndex] = task;
task.ContinueWith(t => OnPrefetchCompleted(blockIndex, t), TaskScheduler.Default);
}
private void OnPrefetchCompleted(int blockIndex, Task<CycleBlock> task)
{
if (!task.IsCompletedSuccessfully)
{
_pending.TryRemove(blockIndex, out _);
return;
}
lock (_sync)
{
if (_next is null || _next.Index < blockIndex)
{
_next = task.Result;
}
_pending.TryRemove(blockIndex, out _);
TrimCacheLocked();
}
}
private void TrimCacheLocked()
{
if (_previous is null)
{
return;
}
if (ReferenceEquals(_previous, _baseBlock))
{
return;
}
if (_previous.IsInUse)
{
return;
}
if (_next is not null && _current is not null)
{
_previous = null;
}
}
private void TryReleaseBlockLocked(CycleBlock? block)
{
if (block is null || ReferenceEquals(block, _baseBlock))
{
return;
}
if (ReferenceEquals(block, _current) || ReferenceEquals(block, _next))
{
return;
}
if (block.IsInUse)
{
return;
}
if (ReferenceEquals(block, _previous))
{
_previous = null;
}
}
private void Release(CycleBlock block)
{
if (!block.IsInUse)
{
lock (_sync)
{
TryReleaseBlockLocked(block);
}
}
}
private CycleBlock GenerateBlock(int blockIndex)
{
if (blockIndex == 0)
{
return _baseBlock;
}
ulong start = GetBlockStart(blockIndex);
int length = GetBlockLength(blockIndex);
ulong[] cycles = new ulong[length];
ulong divisor = start;
for (int i = 0; i < length; i++)
{
cycles[i] = MersenneDivisorCycles.CalculateCycleLength(divisor);
divisor++;
}
return new CycleBlock(blockIndex, start, cycles);
}
private int GetBlockLength(int blockIndex)
{
// TODO: We're over-complicating this again. We should just remember and use the _baseBlock.Cycles.Length. Remember it in a field for better performance during initialization / creation
if (blockIndex == 0)
{
return _baseBlock.Cycles.Length;
}
return PerfectNumberConstants.MaxQForDivisorCycles;
}
private ulong GetBlockStart(int blockIndex)
{
// TODO: We're over-complicating this. There should be only 3 blocks kept. We should simply check all 3 to return the correct one, instead of expensive multiplications.
if (blockIndex == 0)
{
return 0UL;
}
ulong baseEnd = _baseBlock.End;
return baseEnd + 1UL + (ulong)(blockIndex - 1) * (ulong)PerfectNumberConstants.MaxQForDivisorCycles;
}
private int GetBlockIndex(ulong divisor)
{
// TODO: We're over-complicating this. There should be only 3 blocks kept. We should simply check all 3 to return the correct one, instead of expensive multiplications.
if (divisor <= _baseBlock.End)
{
return 0;
}
ulong offset = divisor - (_baseBlock.End + 1UL);
return (int)(offset / (ulong)PerfectNumberConstants.MaxQForDivisorCycles) + 1;
}
internal sealed class CycleBlock
{
internal CycleBlock(int index, ulong start, ulong[] cycles)
{
Index = index;
Start = start;
Cycles = cycles;
End = start + (ulong)cycles.Length - 1UL;
}
internal int Index { get; }
internal ulong Start { get; }
internal ulong End { get; }
internal ulong[] Cycles { get; }
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal ulong GetCycle(ulong divisor)
{
if (divisor < Start || divisor > End)
{
return 0UL;
}
return Cycles[(int)(divisor - Start)];
}
}
public struct Lease : IDisposable
{
private CycleBlock? _block;
internal Lease(CycleBlock block)
{
_block = block;
}
public readonly bool IsValid => _block is not null;
public readonly ulong Start => _block?.Start ?? 0UL;
public readonly ulong End => _block?.End ?? 0UL;
public readonly ulong[]? Values => _block?.Cycles;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public readonly ulong GetCycle(ulong divisor) => _block?.GetCycle(divisor) ?? 0UL;
public void Dispose()
{
_block = null;
}
}
private readonly object _sync = new();
private readonly ConcurrentDictionary<int, Task<CycleBlock>> _pending = new();
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>> _gpuKernelCache = new(AcceleratorReferenceComparer.Instance);
private CycleBlock _baseBlock = null!;
private CycleBlock? _activeBlock;
private CycleBlock? _prefetchedBlock;
private bool _useGpuGeneration = true;
public static DivisorCycleCache Shared { get; } = new DivisorCycleCache();
private DivisorCycleCache()
{
ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
Initialize(snapshot);
}
public void ReloadFromCurrentSnapshot()
{
ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
Initialize(snapshot);
}
public void ConfigureGeneratorDevice(bool useGpu)
{
lock (_sync)
{
_useGpuGeneration = useGpu;
}
}
private void Initialize(ulong[] snapshot)
{
lock (_sync)
{
_pending.Clear();
_baseBlock = new CycleBlock(index: 0, start: 0UL, snapshot);
_activeBlock = null;
_prefetchedBlock = null;
StartPrefetchLocked(1);
}
}
public Lease Acquire(ulong divisor)
{
if (divisor <= _baseBlock.End)
{
return new Lease(_baseBlock);
}
CycleBlock block = AcquireDynamicBlock(divisor);
return new Lease(block);
}
private CycleBlock AcquireDynamicBlock(ulong divisor)
{
int blockIndex = GetBlockIndex(divisor);
while (true)
{
Task<CycleBlock>? pendingTask = null;
lock (_sync)
{
if (TryGetCachedBlockLocked(blockIndex, out CycleBlock? cached))
{
return cached!;
}
if (!_pending.TryGetValue(blockIndex, out pendingTask))
{
pendingTask = Task.Run(() => GenerateBlock(blockIndex));
_pending[blockIndex] = pendingTask;
}
}
CycleBlock block = pendingTask.GetAwaiter().GetResult();
lock (_sync)
{
_pending.TryRemove(blockIndex, out _);
if (TryGetCachedBlockLocked(blockIndex, out CycleBlock? cachedAfterAwait))
{
return cachedAfterAwait!;
}
SetActiveBlockLocked(block);
return block;
}
}
}
private bool TryGetCachedBlockLocked(int blockIndex, out CycleBlock? block)
{
CycleBlock? active = _activeBlock;
if (active is not null && active.Index == blockIndex)
{
block = active;
return true;
}
CycleBlock? prefetched = _prefetchedBlock;
if (prefetched is not null && prefetched.Index == blockIndex)
{
PromotePrefetchedLocked();
block = _activeBlock!;
return true;
}
block = null;
return false;
}
private void PromotePrefetchedLocked()
{
CycleBlock? prefetched = _prefetchedBlock;
if (prefetched is null)
{
return;
}
_activeBlock = prefetched;
_prefetchedBlock = null;
StartPrefetchLocked(prefetched.Index + 1);
}
private void SetActiveBlockLocked(CycleBlock block)
{
_activeBlock = block;
CycleBlock? prefetched = _prefetchedBlock;
if (prefetched is not null && prefetched.Index != block.Index + 1)
{
_prefetchedBlock = null;
}
StartPrefetchLocked(block.Index + 1);
}
private void StartPrefetchLocked(int blockIndex)
{
if (blockIndex <= 0)
{
return;
}
CycleBlock? prefetched = _prefetchedBlock;
if (prefetched is not null && prefetched.Index == blockIndex)
{
return;
}
if (_pending.TryGetValue(blockIndex, out Task<CycleBlock>? existing))
{
if (existing.IsCompletedSuccessfully)
{
_pending.TryRemove(blockIndex, out _);
AssignPrefetchedBlockLocked(blockIndex, existing.Result);
}
return;
}
Task<CycleBlock> task = Task.Run(() => GenerateBlock(blockIndex));
_pending[blockIndex] = task;
task.ContinueWith(t => OnPrefetchCompleted(blockIndex, t), TaskScheduler.Default);
}
private void OnPrefetchCompleted(int blockIndex, Task<CycleBlock> task)
{
if (!task.IsCompletedSuccessfully)
{
_pending.TryRemove(blockIndex, out _);
return;
}
lock (_sync)
{
_pending.TryRemove(blockIndex, out _);
AssignPrefetchedBlockLocked(blockIndex, task.Result);
}
}
private void AssignPrefetchedBlockLocked(int blockIndex, CycleBlock block)
{
CycleBlock? active = _activeBlock;
if (active is not null && blockIndex != active.Index + 1)
{
return;
}
if (_prefetchedBlock is null || _prefetchedBlock.Index < blockIndex)
{
_prefetchedBlock = block;
}
}
private CycleBlock GenerateBlock(int blockIndex)
{
if (blockIndex == 0)
{
return _baseBlock;
}
ulong start = GetBlockStart(blockIndex);
int length = _baseBlock.Cycles.Length;
ulong[] cycles = new ulong[length];
if (_useGpuGeneration)
{
ComputeCyclesGpu(start, cycles);
}
else
{
ComputeCyclesCpu(start, cycles);
}
return new CycleBlock(blockIndex, start, cycles);
}
private int GetBlockIndex(ulong divisor)
{
if (divisor <= _baseBlock.End)
{
return 0;
}
ulong offset = divisor - (_baseBlock.End + 1UL);
return (int)(offset / (ulong)PerfectNumberConstants.MaxQForDivisorCycles) + 1;
}
private ulong GetBlockStart(int blockIndex)
{
if (blockIndex == 0)
{
return 0UL;
}
ulong baseEnd = _baseBlock.End;
return baseEnd + 1UL + (ulong)(blockIndex - 1) * (ulong)PerfectNumberConstants.MaxQForDivisorCycles;
}
private void ComputeCyclesCpu(ulong start, ulong[] destination)
{
ulong divisor = start;
for (int i = 0; i < destination.Length; i++)
{
destination[i] = MersenneDivisorCycles.CalculateCycleLength(divisor);
divisor++;
}
}
private void ComputeCyclesGpu(ulong start, ulong[] destination)
{
var gpuLease = GpuKernelPool.GetKernel(_useGpuGeneration);
var execution = gpuLease.EnterExecutionScope();
try
{
Accelerator accelerator = gpuLease.Accelerator;
Action<Index1D, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>> kernel = _gpuKernelCache.GetOrAdd(accelerator, LoadKernel);
int length = destination.Length;
ulong[] divisors = ArrayPool<ulong>.Shared.Rent(length);
try
{
for (int i = 0; i < length; i++)
{
divisors[i] = start + (ulong)i;
}
using MemoryBuffer1D<ulong, Stride1D.Dense> divisorBuffer = accelerator.Allocate1D<ulong>(length);
using MemoryBuffer1D<ulong, Stride1D.Dense> resultBuffer = accelerator.Allocate1D<ulong>(length);
divisorBuffer.View.CopyFromCPU(ref divisors[0], length);
kernel(length, divisorBuffer.View, resultBuffer.View);
accelerator.Synchronize();
resultBuffer.View.CopyToCPU(ref destination[0], length);
}
finally
{
ArrayPool<ulong>.Shared.Return(divisors, clearArray: false);
}
}
finally
{
execution.Dispose();
gpuLease.Dispose();
}
}
private static Action<Index1D, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>> LoadKernel(Accelerator accelerator)
{
return accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(GpuDivisorCycleKernel);
}
private static void GpuDivisorCycleKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> divisors, ArrayView1D<ulong, Stride1D.Dense> cycles)
{
cycles[index] = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[index]);
}
private sealed class AcceleratorReferenceComparer : IEqualityComparer<Accelerator>
{
internal static AcceleratorReferenceComparer Instance { get; } = new();
public bool Equals(Accelerator? x, Accelerator? y) => ReferenceEquals(x, y);
public int GetHashCode(Accelerator obj) => RuntimeHelpers.GetHashCode(obj);
}
}
File diff suppressed because it is too large. Load diff
+17
View File
@@ -246,6 +246,23 @@ public static class ULongExtensions
return result.MontgomeryMultiply(1UL, modulus, nPrime);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor)
{
if (cycleLength == 0UL)
{
return exponent.Pow2MontgomeryMod(divisor);
}
ulong rotationCount = exponent % cycleLength;
if (rotationCount == 0UL)
{
return 1UL;
}
return rotationCount.Pow2MontgomeryMod(divisor);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 PowMod(this ulong exponent, UInt128 modulus)
{