535 lines
20 KiB
C#
535 lines
20 KiB
C#
using System.Runtime.CompilerServices;
|
|
using PerfectNumbers.Core.Gpu.Accelerators;
|
|
using UInt128 = System.UInt128;
|
|
|
|
namespace PerfectNumbers.Core;
|
|
|
|
public enum ResidueModel
|
|
{
|
|
Identity, // residue of N modulo d
|
|
Mersenne, // residue of (2^N - 1) modulo d
|
|
}
|
|
|
|
/// Tracks residues for a monotonically increasing sequence of numbers against
|
|
/// a dynamic set of divisors. Designed to support repeated divisor checks
|
|
/// while the tested number grows step-by-step without ever decreasing.
|
|
public sealed class ModResidueTracker
|
|
{
|
|
public static readonly ModResidueTracker Shared = new(ResidueModel.Mersenne, initialNumber: PerfectNumberConstants.BiggestKnownEvenPerfectP, initialized: true);
|
|
|
|
private readonly ResidueModel _model;
|
|
private UInt128 _currentNumber;
|
|
private bool _initialized;
|
|
// Per-divisor sorted state at _currentNumber. Aligned arrays for speed.
|
|
private readonly List<UInt128> _divisors = new(1024); // sorted ascending
|
|
private readonly List<UInt128> _residues = new(1024); // same index as _divisors
|
|
private int _cursor; // merge pointer for ascending candidate stream
|
|
private bool _hasLastAppended;
|
|
private UInt128 _lastAppended;
|
|
|
|
public ModResidueTracker(ResidueModel model, UInt128 initialNumber = default, bool initialized = false)
|
|
{
|
|
_model = model;
|
|
_currentNumber = initialNumber;
|
|
_initialized = initialized;
|
|
}
|
|
|
|
/// Returns true if the provided number (either p or M_p depending on model)
|
|
/// is divisible by the divisor. Advances internal state forward if needed
|
|
/// (number must be non-decreasing across calls). Adds and initializes new
|
|
/// divisors on first use.
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
public bool IsDivisibleCpu(UInt128 number, UInt128 divisor)
|
|
{
|
|
ArgumentOutOfRangeException.ThrowIfEqual(divisor, UInt128.Zero);
|
|
|
|
EnsureAt(number);
|
|
|
|
List<UInt128> divisors = _divisors;
|
|
int idx = LowerBound(divisors, divisor);
|
|
if (idx == divisors.Count || divisors[idx] != divisor)
|
|
{
|
|
UInt128 residue = ComputeInitialResidueCpu(number, divisor);
|
|
divisors.Insert(idx, divisor);
|
|
_residues.Insert(idx, residue);
|
|
if (_hasLastAppended && idx == divisors.Count - 1 && divisor >= _lastAppended)
|
|
{
|
|
_lastAppended = divisor;
|
|
}
|
|
else if (!_hasLastAppended)
|
|
{
|
|
_hasLastAppended = true;
|
|
_lastAppended = divisor;
|
|
}
|
|
|
|
return residue == UInt128.Zero;
|
|
}
|
|
|
|
return _residues[idx] == UInt128.Zero;
|
|
}
|
|
|
|
/// Returns true if the provided number (either p or M_p depending on model)
|
|
/// is divisible by the divisor. Advances internal state forward if needed
|
|
/// (number must be non-decreasing across calls). Adds and initializes new
|
|
/// divisors on first use.
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
public bool IsDivisibleGpu(PrimeOrderCalculatorAccelerator gpu, UInt128 number, UInt128 divisor)
|
|
{
|
|
ArgumentOutOfRangeException.ThrowIfEqual(divisor, UInt128.Zero);
|
|
|
|
EnsureAt(number);
|
|
|
|
List<UInt128> divisors = _divisors;
|
|
int idx = LowerBound(divisors, divisor);
|
|
if (idx == divisors.Count || divisors[idx] != divisor)
|
|
{
|
|
UInt128 residue = ComputeInitialResidueGpu(gpu, number, divisor);
|
|
divisors.Insert(idx, divisor);
|
|
_residues.Insert(idx, residue);
|
|
if (_hasLastAppended && idx == divisors.Count - 1 && divisor >= _lastAppended)
|
|
{
|
|
_lastAppended = divisor;
|
|
}
|
|
else if (!_hasLastAppended)
|
|
{
|
|
_hasLastAppended = true;
|
|
_lastAppended = divisor;
|
|
}
|
|
|
|
return residue == UInt128.Zero;
|
|
}
|
|
|
|
return _residues[idx] == UInt128.Zero;
|
|
}
|
|
|
|
/// Returns true if the provided number (either p or M_p depending on model)
|
|
/// is divisible by the divisor. Advances internal state forward if needed
|
|
/// (number must be non-decreasing across calls). Adds and initializes new
|
|
/// divisors on first use.
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
public bool IsDivisibleHybrid(PrimeOrderCalculatorAccelerator gpu, UInt128 number, UInt128 divisor)
|
|
{
|
|
ArgumentOutOfRangeException.ThrowIfEqual(divisor, UInt128.Zero);
|
|
|
|
EnsureAt(number);
|
|
|
|
List<UInt128> divisors = _divisors;
|
|
int idx = LowerBound(divisors, divisor);
|
|
if (idx == divisors.Count || divisors[idx] != divisor)
|
|
{
|
|
UInt128 residue = ComputeInitialResidueHybrid(gpu, number, divisor);
|
|
divisors.Insert(idx, divisor);
|
|
_residues.Insert(idx, residue);
|
|
if (_hasLastAppended && idx == divisors.Count - 1 && divisor >= _lastAppended)
|
|
{
|
|
_lastAppended = divisor;
|
|
}
|
|
else if (!_hasLastAppended)
|
|
{
|
|
_hasLastAppended = true;
|
|
_lastAppended = divisor;
|
|
}
|
|
|
|
return residue == UInt128.Zero;
|
|
}
|
|
|
|
return _residues[idx] == UInt128.Zero;
|
|
}
|
|
|
|
// Merge-walk support: prepare for increasing sequence of candidate divisors.
|
|
public void BeginMerge(UInt128 number)
|
|
{
|
|
EnsureAt(number);
|
|
_cursor = 0;
|
|
}
|
|
|
|
// Merge or append current candidate divisor (ascending order required across calls).
|
|
// Returns true and sets 'divisible' based on residue at 'number'.
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
public bool MergeOrAppendCpu(UInt128 number, UInt128 divisor, out bool divisible)
|
|
{
|
|
// Assume number is non-decreasing; advance residues lazily when number changes.
|
|
List<UInt128> divisors = _divisors;
|
|
if (number != _currentNumber)
|
|
{
|
|
EnsureAt(number);
|
|
// After EnsureAt we can keep _cursor (still valid monotonic pointer)
|
|
if (_cursor > divisors.Count) _cursor = divisors.Count;
|
|
}
|
|
|
|
UInt128 residue;
|
|
// Fast path: append to tail when candidates are strictly increasing.
|
|
List<UInt128> residues = _residues;
|
|
if (divisors.Count == 0 || (_hasLastAppended && divisor >= _lastAppended && divisors[^1] <= divisor))
|
|
{
|
|
residue = ComputeInitialResidueCpu(number, divisor);
|
|
divisors.Add(divisor);
|
|
residues.Add(residue);
|
|
_hasLastAppended = true;
|
|
_lastAppended = divisor;
|
|
_cursor = divisors.Count; // position after last
|
|
divisible = residue == UInt128.Zero;
|
|
return true;
|
|
}
|
|
|
|
// Merge pointer: advance until current >= divisor
|
|
var cursor = _cursor;
|
|
while (cursor < divisors.Count && divisors[cursor] < divisor)
|
|
{
|
|
cursor++;
|
|
}
|
|
|
|
if (cursor < divisors.Count && divisors[cursor] == divisor)
|
|
{
|
|
divisible = residues[cursor] == UInt128.Zero;
|
|
_cursor = cursor + 1; // position after matched
|
|
return true;
|
|
}
|
|
|
|
|
|
// Insert at cursor to keep order
|
|
residue = ComputeInitialResidueCpu(number, divisor);
|
|
divisors.Insert(cursor, divisor);
|
|
residues.Insert(cursor, residue);
|
|
divisible = residue == UInt128.Zero;
|
|
_cursor = cursor + 1;
|
|
|
|
// update last appended hint
|
|
if (!_hasLastAppended || divisor >= _lastAppended)
|
|
{
|
|
_hasLastAppended = true;
|
|
_lastAppended = divisor;
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
// Merge or append current candidate divisor (ascending order required across calls).
|
|
// Returns true and sets 'divisible' based on residue at 'number'.
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
public bool MergeOrAppendGpu(PrimeOrderCalculatorAccelerator gpu, UInt128 number, UInt128 divisor, out bool divisible)
|
|
{
|
|
// Assume number is non-decreasing; advance residues lazily when number changes.
|
|
List<UInt128> divisors = _divisors;
|
|
if (number != _currentNumber)
|
|
{
|
|
EnsureAt(number);
|
|
// After EnsureAt we can keep _cursor (still valid monotonic pointer)
|
|
if (_cursor > divisors.Count) _cursor = divisors.Count;
|
|
}
|
|
|
|
UInt128 residue;
|
|
// Fast path: append to tail when candidates are strictly increasing.
|
|
List<UInt128> residues = _residues;
|
|
if (divisors.Count == 0 || (_hasLastAppended && divisor >= _lastAppended && divisors[^1] <= divisor))
|
|
{
|
|
residue = ComputeInitialResidueGpu(gpu, number, divisor);
|
|
divisors.Add(divisor);
|
|
residues.Add(residue);
|
|
_hasLastAppended = true;
|
|
_lastAppended = divisor;
|
|
_cursor = divisors.Count; // position after last
|
|
divisible = residue == UInt128.Zero;
|
|
return true;
|
|
}
|
|
|
|
// Merge pointer: advance until current >= divisor
|
|
var cursor = _cursor;
|
|
while (cursor < divisors.Count && divisors[cursor] < divisor)
|
|
{
|
|
cursor++;
|
|
}
|
|
|
|
if (cursor < divisors.Count && divisors[cursor] == divisor)
|
|
{
|
|
divisible = residues[cursor] == UInt128.Zero;
|
|
_cursor = cursor + 1; // position after matched
|
|
return true;
|
|
}
|
|
|
|
|
|
// Insert at cursor to keep order
|
|
residue = ComputeInitialResidueGpu(gpu, number, divisor);
|
|
divisors.Insert(cursor, divisor);
|
|
residues.Insert(cursor, residue);
|
|
divisible = residue == UInt128.Zero;
|
|
_cursor = cursor + 1;
|
|
|
|
// update last appended hint
|
|
if (!_hasLastAppended || divisor >= _lastAppended)
|
|
{
|
|
_hasLastAppended = true;
|
|
_lastAppended = divisor;
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
// Merge or append current candidate divisor (ascending order required across calls).
|
|
// Returns true and sets 'divisible' based on residue at 'number'.
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
public bool MergeOrAppendHybrid(PrimeOrderCalculatorAccelerator gpu, UInt128 number, UInt128 divisor, out bool divisible)
|
|
{
|
|
// Assume number is non-decreasing; advance residues lazily when number changes.
|
|
List<UInt128> divisors = _divisors;
|
|
if (number != _currentNumber)
|
|
{
|
|
EnsureAt(number);
|
|
// After EnsureAt we can keep _cursor (still valid monotonic pointer)
|
|
if (_cursor > divisors.Count) _cursor = divisors.Count;
|
|
}
|
|
|
|
UInt128 residue;
|
|
// Fast path: append to tail when candidates are strictly increasing.
|
|
List<UInt128> residues = _residues;
|
|
if (divisors.Count == 0 || (_hasLastAppended && divisor >= _lastAppended && divisors[^1] <= divisor))
|
|
{
|
|
residue = ComputeInitialResidueHybrid(gpu, number, divisor);
|
|
divisors.Add(divisor);
|
|
residues.Add(residue);
|
|
_hasLastAppended = true;
|
|
_lastAppended = divisor;
|
|
_cursor = divisors.Count; // position after last
|
|
divisible = residue == UInt128.Zero;
|
|
return true;
|
|
}
|
|
|
|
// Merge pointer: advance until current >= divisor
|
|
var cursor = _cursor;
|
|
while (cursor < divisors.Count && divisors[cursor] < divisor)
|
|
{
|
|
cursor++;
|
|
}
|
|
|
|
if (cursor < divisors.Count && divisors[cursor] == divisor)
|
|
{
|
|
divisible = residues[cursor] == UInt128.Zero;
|
|
_cursor = cursor + 1; // position after matched
|
|
return true;
|
|
}
|
|
|
|
|
|
// Insert at cursor to keep order
|
|
residue = ComputeInitialResidueHybrid(gpu, number, divisor);
|
|
divisors.Insert(cursor, divisor);
|
|
residues.Insert(cursor, residue);
|
|
divisible = residue == UInt128.Zero;
|
|
_cursor = cursor + 1;
|
|
|
|
// update last appended hint
|
|
if (!_hasLastAppended || divisor >= _lastAppended)
|
|
{
|
|
_hasLastAppended = true;
|
|
_lastAppended = divisor;
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
/// Ensures the internal state corresponds to the provided number by moving
|
|
/// forward from the current number (never backwards). O(Count(divisors)).
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
public void EnsureAt(UInt128 number)
|
|
{
|
|
if (!_initialized)
|
|
{
|
|
_currentNumber = number;
|
|
// Initialize residues lazily upon first check for each divisor.
|
|
_initialized = true;
|
|
return;
|
|
}
|
|
|
|
if (number == _currentNumber)
|
|
{
|
|
return;
|
|
}
|
|
|
|
if (number < _currentNumber)
|
|
{
|
|
throw new InvalidOperationException("Number must be non-decreasing.");
|
|
}
|
|
|
|
UInt128 delta = number - _currentNumber;
|
|
List<UInt128> divisors = _divisors;
|
|
int divisorsCount = divisors.Count;
|
|
if (divisorsCount == 0)
|
|
{
|
|
_currentNumber = number;
|
|
return;
|
|
}
|
|
|
|
// Advance all tracked residues from _currentNumber to 'number'.
|
|
// For Identity: r' = (r + delta) mod d
|
|
// For Mersenne: r' = (2^delta * r + (2^delta - 1)) mod d
|
|
// Optimize small deltas with simple forms when possible.
|
|
int i = 0;
|
|
UInt128 d, value;
|
|
List<UInt128> residues = _residues;
|
|
|
|
if (_model == ResidueModel.Identity)
|
|
{
|
|
// Avoid 128-bit division: delta is small in scanner (<= 6) so use add/sub.
|
|
for (; i < divisorsCount; i++)
|
|
{
|
|
d = divisors[i];
|
|
value = residues[i] + delta;
|
|
// TODO: Replace this subtraction loop with the divisor-cycle stepping helper so identity
|
|
// residues reuse cached cycle deltas, computing the missing cycle on the configured
|
|
// device without caching the result or scheduling extra blocks when the snapshot lacks
|
|
// the divisor.
|
|
while (value >= d)
|
|
{
|
|
value -= d;
|
|
}
|
|
|
|
residues[i] = value;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
for (; i < divisors.Count; i++)
|
|
{
|
|
d = divisors[i];
|
|
// The pow-delta term now calls the windowed Pow2 Montgomery helper so residue updates reuse the fast ladder.
|
|
UInt128 powDelta = delta.Pow2MontgomeryModWindowed(d);
|
|
// r' = pow2Delta*r + (pow2Delta - 1) mod d
|
|
value = powDelta.MulMod(residues[i], d);
|
|
value = value.AddMod(powDelta, d);
|
|
value = value.SubtractOneMod(d);
|
|
|
|
residues[i] = value;
|
|
}
|
|
}
|
|
|
|
_currentNumber = number;
|
|
}
|
|
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
private static int LowerBound(List<UInt128> list, UInt128 value)
|
|
{
|
|
int lo = 0, hi = list.Count, mid;
|
|
while (lo < hi)
|
|
{
|
|
mid = (int)(((uint)lo + (uint)hi) >> 1);
|
|
if (list[mid] < value)
|
|
{
|
|
lo = mid + 1;
|
|
}
|
|
else
|
|
{
|
|
hi = mid;
|
|
}
|
|
}
|
|
|
|
return lo;
|
|
}
|
|
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
private UInt128 ComputeInitialResidueCpu(UInt128 number, UInt128 divisor)
|
|
{
|
|
if (_model == ResidueModel.Identity)
|
|
{
|
|
// Fast path for 64-bit inputs
|
|
if (number <= ulong.MaxValue && divisor <= ulong.MaxValue)
|
|
{
|
|
// TODO: Swap this modulo fallback for the shared divisor-cycle remainder helper so 64-bit
|
|
// residues reuse cached cycles instead of recomputing `%` on every lookup.
|
|
return (ulong)number % (ulong)divisor;
|
|
}
|
|
|
|
// TODO: Use the Montgomery folding helper measured faster in the Pow2Montgomery benchmarks so this
|
|
// wide residue path avoids general `%` on UInt128 inputs.
|
|
return number % divisor;
|
|
}
|
|
|
|
// residue of Mersenne number M_number = 2^number - 1 modulo divisor
|
|
// TODO: Consult MersenneDivisorCycles.Shared (or similar caches) here so we reuse precomputed cycle lengths
|
|
// instead of recomputing powmods for every divisor; the divisor-cycle benchmarks showed large wins once the
|
|
// cached orders were used across scans, and a miss should trigger the configured device to compute the cycle
|
|
// immediately without storing it back or requesting additional cache blocks.
|
|
UInt128 pow = number.Pow2MontgomeryModWindowed(divisor);
|
|
|
|
// q divides M_p ⇔ 2^p ≡ 1 (mod q) AND ord_q(2) | p
|
|
if (pow != UInt128.One)
|
|
return (pow - UInt128.One) % divisor;
|
|
|
|
ulong ord = divisor.CalculateOrderCpu();
|
|
if (ord == 0UL || number % (UInt128)ord != UInt128.Zero)
|
|
return UInt128.One; // does NOT divide
|
|
|
|
return UInt128.Zero; // q really divides M_p
|
|
}
|
|
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
private UInt128 ComputeInitialResidueGpu(PrimeOrderCalculatorAccelerator gpu, UInt128 number, UInt128 divisor)
|
|
{
|
|
if (_model == ResidueModel.Identity)
|
|
{
|
|
// Fast path for 64-bit inputs
|
|
if (number <= ulong.MaxValue && divisor <= ulong.MaxValue)
|
|
{
|
|
// TODO: Swap this modulo fallback for the shared divisor-cycle remainder helper so 64-bit
|
|
// residues reuse cached cycles instead of recomputing `%` on every lookup.
|
|
return (ulong)number % (ulong)divisor;
|
|
}
|
|
|
|
// TODO: Use the Montgomery folding helper measured faster in the Pow2Montgomery benchmarks so this
|
|
// wide residue path avoids general `%` on UInt128 inputs.
|
|
return number % divisor;
|
|
}
|
|
|
|
// residue of Mersenne number M_number = 2^number - 1 modulo divisor
|
|
// TODO: Consult MersenneDivisorCycles.Shared (or similar caches) here so we reuse precomputed cycle lengths
|
|
// instead of recomputing powmods for every divisor; the divisor-cycle benchmarks showed large wins once the
|
|
// cached orders were used across scans, and a miss should trigger the configured device to compute the cycle
|
|
// immediately without storing it back or requesting additional cache blocks.
|
|
UInt128 pow = number.Pow2MontgomeryModWindowed(divisor);
|
|
|
|
// q divides M_p ⇔ 2^p ≡ 1 (mod q) AND ord_q(2) | p
|
|
if (pow != UInt128.One)
|
|
return (pow - UInt128.One) % divisor;
|
|
|
|
ulong ord = divisor.CalculateOrderGpu(gpu);
|
|
if (ord == 0UL || number % (UInt128)ord != UInt128.Zero)
|
|
return UInt128.One; // does NOT divide
|
|
|
|
return UInt128.Zero; // q really divides M_p
|
|
}
|
|
|
|
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
|
private UInt128 ComputeInitialResidueHybrid(PrimeOrderCalculatorAccelerator gpu, UInt128 number, UInt128 divisor)
|
|
{
|
|
if (_model == ResidueModel.Identity)
|
|
{
|
|
// Fast path for 64-bit inputs
|
|
if (number <= ulong.MaxValue && divisor <= ulong.MaxValue)
|
|
{
|
|
// TODO: Swap this modulo fallback for the shared divisor-cycle remainder helper so 64-bit
|
|
// residues reuse cached cycles instead of recomputing `%` on every lookup.
|
|
return (ulong)number % (ulong)divisor;
|
|
}
|
|
|
|
// TODO: Use the Montgomery folding helper measured faster in the Pow2Montgomery benchmarks so this
|
|
// wide residue path avoids general `%` on UInt128 inputs.
|
|
return number % divisor;
|
|
}
|
|
|
|
// residue of Mersenne number M_number = 2^number - 1 modulo divisor
|
|
// TODO: Consult MersenneDivisorCycles.Shared (or similar caches) here so we reuse precomputed cycle lengths
|
|
// instead of recomputing powmods for every divisor; the divisor-cycle benchmarks showed large wins once the
|
|
// cached orders were used across scans, and a miss should trigger the configured device to compute the cycle
|
|
// immediately without storing it back or requesting additional cache blocks.
|
|
UInt128 pow = number.Pow2MontgomeryModWindowed(divisor);
|
|
|
|
// q divides M_p ⇔ 2^p ≡ 1 (mod q) AND ord_q(2) | p
|
|
if (pow != UInt128.One)
|
|
return (pow - UInt128.One) % divisor;
|
|
|
|
ulong ord = divisor.CalculateOrderHybrid(gpu);
|
|
if (ord == 0UL || number % (UInt128)ord != UInt128.Zero)
|
|
return UInt128.One; // does NOT divide
|
|
|
|
return UInt128.Zero; // q really divides M_p
|
|
}
|
|
}
|