AscendingdivisorRemainderStepper used in PrimeTesterByLastDigit

This commit is contained in:
mlemanczyk committed 2025-12-05 14:17:24 +01:00
1 parent b3a325521a
commit 928e847f1f
9 files changed
+210 -268

No files matched your search

@@ -17,7 +17,7 @@ public class ULongExtensionsTests
[InlineData(ulong.MaxValue)]
public void Mod10_returns_same_as_operator(ulong value)
{
ulong expected = value % 10UL;
byte expected = (byte)(value % 10UL);
ULongExtensions.Mod10(value).Should().Be(expected);
}
@@ -4,33 +4,34 @@ namespace PerfectNumbers.Core;
public struct AscendingDivisorRemainderStepper
{
private readonly ulong dividend;
private ulong maxDivisor;
private readonly ulong _dividend;
private ulong _maxDivisor;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public AscendingDivisorRemainderStepper(ulong value)
{
dividend = value;
_dividend = value;
ulong root = (ulong)Math.Sqrt(value);
if ((UInt128)root * root < value)
{
root++;
}
maxDivisor = root;
_maxDivisor = root;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public readonly bool ShouldContinue(ulong divisor) => divisor <= maxDivisor;
public readonly bool ShouldContinue(ulong divisor) => divisor <= _maxDivisor;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public bool Divides(ulong divisor)
{
ulong quotient = dividend / divisor;
ulong remainder = dividend - (quotient * divisor);
if (quotient < maxDivisor)
ulong remainder = _dividend;
ulong quotient = remainder / divisor;
remainder -= quotient * divisor;
if (quotient < _maxDivisor)
{
maxDivisor = quotient;
_maxDivisor = quotient;
}
return remainder == 0UL;
@@ -27,7 +27,7 @@ internal static class IncrementalKernels
ulong idxMod3 = idx % 3UL;
ulong idxMod5 = idx % 5UL;
// residue automaton
ulong step10 = (exponent.Mod10() << 1).Mod10();
ulong step10 = ((ulong)(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
@@ -30,7 +30,7 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u
KernelContainer kernels = GpuKernelPool.GetOrAddKernels(acceleratorIndex, stream, KernelType.Pow2ModKernelScan | KernelType.IncrementalKernelScan);
var pow2Kernel = kernels.Pow2Mod!;
var incKernel = kernels.Incremental!;
ulong step10 = (exponent.Mod10() << 1).Mod10();
ulong step10 = ((ulong)(exponent.Mod10() << 1)).Mod10();
ulong step8 = ((exponent & 7UL) << 1) & 7UL;
ulong step3 = ((exponent % 3UL) << 1) % 3UL;
ulong step5 = ((exponent % 5UL) << 1) % 5UL;
@@ -124,13 +124,13 @@ public static class MersenneNumberDivisorByDivisorTester
{
ulong candidate = primesSpan[index];
// The by-divisor CPU scan only operates on primes greater than 138,000,000, so the guard below never triggers.
// if (candidate <= 1UL)
// {
// markComposite();
// printResult(candidate, false, false, false);
// continue;
// }
// The by-divisor CPU scan only operates on primes greater than 138,000,000, so the guard below never triggers.
// if (candidate <= 1UL)
// {
// markComposite();
// printResult(candidate, false, false, false);
// continue;
// }
// This implementation is terribly slow, while this method expect prime p given as --filter-p input already.
// We don't need to additionally check it.
@@ -206,13 +206,13 @@ public static class MersenneNumberDivisorByDivisorTester
ulong prime = primeSpan[i];
ulong allowedMax = allowedMaxSpan[i];
// Primes in the production by-divisor flow yield massive divisor limits, so the short-circuit below never applies.
// if (allowedMax < 3UL)
// {
// clearComposite();
// printResult(prime, true, true, true);
// continue;
// }
// Primes in the production by-divisor flow yield massive divisor limits, so the short-circuit below never applies.
// if (allowedMax < 3UL)
// {
// clearComposite();
// printResult(prime, true, true, true);
// continue;
// }
filteredPrimes.Add(prime);
}
@@ -269,48 +269,56 @@ public static class MersenneNumberDivisorByDivisorTester
void ProcessPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime)
{
Console.WriteLine($"Processing {prime}");
IMersenneNumberDivisorByDivisorTester testerForPrime = CreateTester(tester);
string stateFile = Path.Combine(PerfectNumberConstants.ByDivisorStateDirectory, prime.ToString(CultureInfo.InvariantCulture) + ".bin");
ulong resumeK = tester.MinK;
if (File.Exists(stateFile))
try
{
ulong lastK = File.ReadLines(stateFile)
.Where(static line => !string.IsNullOrWhiteSpace(line))
.Select(static line => ulong.Parse(line, NumberStyles.None, CultureInfo.InvariantCulture))
.Max();
Console.WriteLine($"Processing {prime}");
IMersenneNumberDivisorByDivisorTester testerForPrime = CreateTester(tester);
string stateFile = Path.Combine(PerfectNumberConstants.ByDivisorStateDirectory, prime.ToString(CultureInfo.InvariantCulture) + ".bin");
resumeK = lastK + 1UL;
testerForPrime.ResumeFromState(lastK);
}
else
{
testerForPrime.ResumeFromState(0UL);
}
testerForPrime.MinK = resumeK;
testerForPrime.StateFilePath = stateFile;
testerForPrime.ResetStateTracking();
bool isPrime = testerForPrime.IsPrime(gpu, prime, out bool divisorsExhausted, out ulong divisor);
if (!isPrime)
{
if (!string.IsNullOrEmpty(stateFile) && File.Exists(stateFile))
ulong resumeK = tester.MinK;
if (File.Exists(stateFile))
{
File.Delete(stateFile);
ulong lastK = File.ReadLines(stateFile)
.Where(static line => !string.IsNullOrWhiteSpace(line))
.Select(static line => ulong.Parse(line, NumberStyles.None, CultureInfo.InvariantCulture))
.Max();
resumeK = lastK + 1UL;
testerForPrime.ResumeFromState(lastK);
}
else
{
testerForPrime.ResumeFromState(0UL);
}
markComposite();
printResult(prime, true, true, false, divisor);
Console.WriteLine($"Finished processing {prime}");
return;
}
testerForPrime.MinK = resumeK;
testerForPrime.StateFilePath = stateFile;
testerForPrime.ResetStateTracking();
clearComposite();
printResult(prime, true, divisorsExhausted, true, divisor);
Console.WriteLine($"Finished processing {prime}");
bool isPrime = testerForPrime.IsPrime(gpu, prime, out bool divisorsExhausted, out ulong divisor);
if (!isPrime)
{
if (!string.IsNullOrEmpty(stateFile) && File.Exists(stateFile))
{
File.Delete(stateFile);
}
markComposite();
printResult(prime, true, true, false, divisor);
Console.WriteLine($"Finished processing {prime}");
return;
}
clearComposite();
printResult(prime, true, divisorsExhausted, true, divisor);
Console.WriteLine($"Finished processing {prime}");
}
catch (Exception ex)
{
Console.WriteLine($"Error processing {ex.StackTrace}");
Environment.Exit(1);
}
}
if (workerCount == 1)
+56 -40
View File
@@ -9,13 +9,6 @@ namespace PerfectNumbers.Core;
public sealed class PrimeTester
{
public PrimeTester()
{
}
[ThreadStatic]
private static PrimeTester? _tester;
[ThreadStatic]
private static uint[]? _smallPrimeRemainders;
@@ -25,7 +18,6 @@ public sealed class PrimeTester
[ThreadStatic]
private static ulong _lastRemainderLimit;
public static PrimeTester Exclusive => _tester ??= new();
public static bool IsPrimeCpu(ulong n)
{
@@ -50,58 +42,82 @@ public sealed class PrimeTester
}
uint[] smallPrimeDivisors = PrimesGenerator.SmallPrimes;
ulong[] smallPrimeDivisorsPow2 = PrimesGenerator.SmallPrimesPow2;
int smallPrimeDivisorsLength = smallPrimeDivisors.Length;
uint[] remainders = _smallPrimeRemainders ??= new uint[smallPrimeDivisorsLength];
if (remainders.Length < smallPrimeDivisorsLength)
uint[]? remainders = _smallPrimeRemainders;
if (remainders == null)
{
remainders = _smallPrimeRemainders = new uint[smallPrimeDivisorsLength];
}
bool canIncrement = _lastTrackedValue != 0UL && n > _lastTrackedValue;
ulong delta = canIncrement ? n - _lastTrackedValue : 0UL;
ulong maxDivisor = (ulong)Math.Sqrt(n);
if ((UInt128)maxDivisor * maxDivisor < n)
{
maxDivisor++;
remainders = _smallPrimeRemainders = new uint[PrimesGenerator.SmallPrimes.Length];
}
for (int i = 1; i < smallPrimeDivisorsLength; i++)
ref ulong lastTrackedValue = ref _lastTrackedValue;
ref ulong lastRemainderLimit = ref _lastRemainderLimit;
bool canReusePrevious = lastTrackedValue != 0UL && n > lastTrackedValue;
ulong delta = canReusePrevious ? n - lastTrackedValue : 0UL;
canReusePrevious = canReusePrevious && lastRemainderLimit != 0UL && n < lastRemainderLimit;
ulong nextCutoff = 0UL;
uint divisor, remainder;
ulong divisorSquared;
int i;
if (canReusePrevious)
{
uint divisor = smallPrimeDivisors[i];
for (i = 1; i < smallPrimeDivisorsLength; i++)
{
divisorSquared = smallPrimeDivisorsPow2[i];
if (divisorSquared > n)
{
nextCutoff = divisorSquared;
break;
}
if (divisor > maxDivisor)
{
break;
}
uint remainder;
bool canReuse = canIncrement && divisor <= _lastRemainderLimit;
if (!canReuse)
{
remainder = (uint)(n % divisor);
remainders[i] = remainder;
}
else
{
uint previousRemainder = remainders[i];
divisor = smallPrimeDivisors[i];
uint deltaMod = delta >= divisor ? (uint)(delta % divisor) : (uint)delta;
remainder = previousRemainder + deltaMod;
remainder = remainders[i] + deltaMod;
if (remainder >= divisor)
{
remainder -= divisor;
}
remainders[i] = remainder;
if (remainder == 0UL)
{
lastTrackedValue = n;
lastRemainderLimit = nextCutoff == 0UL ? ulong.MaxValue : nextCutoff;
return false;
}
}
lastTrackedValue = n;
lastRemainderLimit = nextCutoff == 0UL ? ulong.MaxValue : nextCutoff;
return true;
}
for (i = 1; i < smallPrimeDivisorsLength; i++)
{
divisorSquared = smallPrimeDivisorsPow2[i];
if (divisorSquared > n)
{
nextCutoff = divisorSquared;
break;
}
divisor = smallPrimeDivisors[i];
remainder = (uint)(n % divisor);
remainders[i] = remainder;
if (remainder == 0UL)
{
_lastTrackedValue = n;
_lastRemainderLimit = maxDivisor;
lastTrackedValue = n;
lastRemainderLimit = nextCutoff == 0UL ? ulong.MaxValue : nextCutoff;
return false;
}
}
_lastTrackedValue = n;
_lastRemainderLimit = maxDivisor;
lastTrackedValue = n;
lastRemainderLimit = nextCutoff == 0UL ? ulong.MaxValue : nextCutoff;
return true;
}
+75 -160
View File
@@ -1,4 +1,3 @@
using System.Numerics;
using System.Runtime.CompilerServices;
using ILGPU.Runtime;
using PerfectNumbers.Core.Gpu.Accelerators;
@@ -7,15 +6,6 @@ namespace PerfectNumbers.Core;
public sealed class PrimeTesterByLastDigit
{
public PrimeTesterByLastDigit()
{
}
[ThreadStatic]
private static PrimeTesterByLastDigit? _tester;
public static PrimeTesterByLastDigit Exclusive => _tester ??= new();
public static bool IsPrimeCpu(ulong n)
{
// The below IFs never trigger in production code of EvenPerfectBitScanner on --mersenne=bydivisor path.
@@ -37,19 +27,19 @@ public sealed class PrimeTesterByLastDigit
uint[] smallPrimeDivisors;
ulong[] smallPrimeDivisorsMul;
byte nMod10 = n.Mod10();
ulong nMod10 = n.Mod10();
switch (nMod10)
{
case 1UL:
case 1:
smallPrimeDivisors = PrimesGenerator.SmallPrimesLastOne;
smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2LastOne;
break;
case 3UL:
case 3:
smallPrimeDivisors = DivisorGenerator.SmallPrimesLastThree;
smallPrimeDivisorsMul = DivisorGenerator.SmallPrimesPow2LastThree;
break;
case 7UL:
case 7:
smallPrimeDivisors = PrimesGenerator.SmallPrimesLastSeven;
smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2LastSeven;
break;
@@ -60,16 +50,18 @@ public sealed class PrimeTesterByLastDigit
}
int smallPrimeDivisorsLength = smallPrimeDivisors.Length;
AscendingDivisorRemainderStepper stepper = new(n);
for (int i = 0; i < smallPrimeDivisorsLength; i++)
{
ulong value = smallPrimeDivisorsMul[i];
if (value > n)
if (!stepper.ShouldContinue(value))
{
break;
}
value = n % smallPrimeDivisors[i];
if (value == 0)
uint divisor = smallPrimeDivisors[i];
if (stepper.Divides(divisor))
{
return false;
}
@@ -78,162 +70,85 @@ public sealed class PrimeTesterByLastDigit
return true;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsPrimeGpu(PrimeOrderCalculatorAccelerator gpu, ulong n)
{
byte flag = 0;
var inputView = gpu.InputView;
var outputView = gpu.OutputByteView;
// [MethodImpl(MethodImplOptions.AggressiveInlining)]
// public static bool IsPrimeGpu(PrimeOrderCalculatorAccelerator gpu, ulong n)
// {
// byte flag = 0;
// var inputView = gpu.InputView;
// var outputView = gpu.OutputByteView;
int acceleratorIndex = gpu.AcceleratorIndex;
AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex);
inputView.CopyFromCPU(stream, ref n, 1);
// int acceleratorIndex = gpu.AcceleratorIndex;
// AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex);
// inputView.CopyFromCPU(stream, ref n, 1);
var kernelLauncher = gpu.SmallPrimeSieveKernelLauncher;
// var kernelLauncher = gpu.SmallPrimeSieveKernelLauncher;
kernelLauncher(
stream,
1,
inputView,
gpu.DevicePrimesLastOne,
gpu.DevicePrimesLastSeven,
gpu.DevicePrimesLastThree,
gpu.DevicePrimesLastNine,
gpu.DevicePrimesPow2LastOne,
gpu.DevicePrimesPow2LastSeven,
gpu.DevicePrimesPow2LastThree,
gpu.DevicePrimesPow2LastNine,
outputView);
// kernelLauncher(
// stream,
// 1,
// inputView,
// gpu.DevicePrimesLastOne,
// gpu.DevicePrimesLastSeven,
// gpu.DevicePrimesLastThree,
// gpu.DevicePrimesLastNine,
// gpu.DevicePrimesPow2LastOne,
// gpu.DevicePrimesPow2LastSeven,
// gpu.DevicePrimesPow2LastThree,
// gpu.DevicePrimesPow2LastNine,
// outputView);
outputView.CopyToCPU(stream, ref flag, 1);
stream.Synchronize();
AcceleratorStreamPool.Return(acceleratorIndex, stream);
// outputView.CopyToCPU(stream, ref flag, 1);
// stream.Synchronize();
// AcceleratorStreamPool.Return(acceleratorIndex, stream);
return flag != 0;
}
// return flag != 0;
// }
public static int GpuBatchSize { get; set; } = 262_144;
// public static int GpuBatchSize { get; set; } = 262_144;
private static readonly object GpuWarmUpLock = new();
private static int WarmedGpuLeaseCount;
// [MethodImpl(MethodImplOptions.AggressiveInlining)]
// public static void IsPrimeBatchGpu(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> values, Span<byte> results)
// {
// int acceleratorIndex = gpu.AcceleratorIndex;
// int totalLength = values.Length;
// int batchSize = GpuBatchSize;
public static void WarmUpGpuKernels(int threadCount)
{
int target = threadCount >> 2;
if (target == 0)
{
target = threadCount;
}
// var inputView = gpu.InputView;
// var outputView = gpu.OutputByteView;
lock (GpuWarmUpLock)
{
if (target <= WarmedGpuLeaseCount)
{
return;
}
// int pos = 0;
// AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex);
// var kernelLauncher = gpu.SmallPrimeSieveKernelLauncher;
PrimeOrderCalculatorAccelerator.WarmUp();
WarmedGpuLeaseCount = target;
}
}
// while (pos < totalLength)
// {
// int remaining = totalLength - pos;
// int count = remaining > batchSize ? batchSize : remaining;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void IsPrimeBatchGpu(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> values, Span<byte> results)
{
// GpuPrimeWorkLimiter.Acquire();
int acceleratorIndex = gpu.AcceleratorIndex;
int totalLength = values.Length;
int batchSize = GpuBatchSize;
// var valueSlice = values.Slice(pos, count);
// inputView.CopyFromCPU(stream, valueSlice);
var inputView = gpu.InputView;
var outputView = gpu.OutputByteView;
// kernelLauncher(
// stream,
// count,
// inputView,
// gpu.DevicePrimesLastOne,
// gpu.DevicePrimesLastSeven,
// gpu.DevicePrimesLastThree,
// gpu.DevicePrimesLastNine,
// gpu.DevicePrimesPow2LastOne,
// gpu.DevicePrimesPow2LastSeven,
// gpu.DevicePrimesPow2LastThree,
// gpu.DevicePrimesPow2LastNine,
// outputView);
int pos = 0;
AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex);
var kernelLauncher = gpu.SmallPrimeSieveKernelLauncher;
// var resultSlice = results.Slice(pos, count);
// outputView.CopyToCPU(stream, resultSlice);
while (pos < totalLength)
{
int remaining = totalLength - pos;
int count = remaining > batchSize ? batchSize : remaining;
// pos += count;
// }
var valueSlice = values.Slice(pos, count);
inputView.CopyFromCPU(stream, valueSlice);
kernelLauncher(
stream,
count,
inputView,
gpu.DevicePrimesLastOne,
gpu.DevicePrimesLastSeven,
gpu.DevicePrimesLastThree,
gpu.DevicePrimesLastNine,
gpu.DevicePrimesPow2LastOne,
gpu.DevicePrimesPow2LastSeven,
gpu.DevicePrimesPow2LastThree,
gpu.DevicePrimesPow2LastNine,
outputView);
var resultSlice = results.Slice(pos, count);
outputView.CopyToCPU(stream, resultSlice);
pos += count;
}
stream.Synchronize();
AcceleratorStreamPool.Return(acceleratorIndex, stream);
// GpuPrimeWorkLimiter.Release();
}
// Expose cache clearing for accelerator disposal coordination
public static void ClearGpuCaches()
{
PrimeOrderCalculatorAccelerator.Clear();
}
internal static void DisposeGpuContexts()
{
PrimeOrderCalculatorAccelerator.DisposeAll();
lock (GpuWarmUpLock)
{
WarmedGpuLeaseCount = 0;
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static bool SharesFactorWithMaxExponent(ulong n)
{
// TODO: Replace this on-the-fly GCD probe with the cached factor table derived from
// ResidueComputationBenchmarks so divisor-cycle metadata can short-circuit the test
// instead of recomputing binary GCD for every candidate.
ulong m = (ulong)BitOperations.Log2(n);
return n.BinaryGcd(m) != 1UL;
}
internal static void SharesFactorWithMaxExponentBatch(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan<ulong> values, Span<byte> results)
{
// TODO: Route this batch helper through the shared GPU kernel pool from
// GpuUInt128BinaryGcdBenchmarks so we reuse cached kernels, pinned host buffers,
// and divisor-cycle staging instead of allocating new device buffers per call.
// Check in benchmarks, which implementation was the fastest, is compatible with GPU,
// and implement it.
int length = values.Length;
int acceleratorIndex = gpu.AcceleratorIndex;
gpu.EnsureCapacity(0, length);
var inputBufferView = gpu.InputView;
var resultBufferView = gpu.OutputByteView;
var kernelLauncher = gpu.SharesFactorKernelLauncher;
AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex);
inputBufferView.CopyFromCPU(stream, values);
kernelLauncher(stream, length, inputBufferView, resultBufferView);
resultBufferView.CopyToCPU(stream, in results);
stream.Synchronize();
AcceleratorStreamPool.Return(acceleratorIndex, stream);
}
// stream.Synchronize();
// AcceleratorStreamPool.Return(acceleratorIndex, stream);
// }
}
+6 -6
View File
@@ -434,13 +434,13 @@ public static partial class ULongExtensions
// Benchmarks (Mod5ULongBenchmarks) show the direct `% 5` is still cheaper (~0.26 ns vs 0.43 ns), so keep the modulo until a faster lookup is proven.
// (Mod8/Mod10 stay masked because they win; Mod5 currently does not.)
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod10(this ulong value) => (value % 5UL) switch
public static byte Mod10(this ulong value) => (byte)(value % 5UL) switch
{
0UL => 0UL,
1UL => 6UL,
2UL => 2UL,
3UL => 8UL,
_ => 4UL,
0 => 0,
1 => 6,
2 => 2,
3 => 8,
_ => 4,
};
[MethodImpl(MethodImplOptions.AggressiveInlining)]
+2
View File
@@ -57,6 +57,8 @@ This document explains **what the program does**, **how it works step-by-step**,
* `--mersenne=incremental`
* `--mersenne=lucas`
* `--mersenne=residue`
* `--mersenne=divisor`
* `--mersenne=bydivisor`
# EvenPerfectBitScanner – Configuration and Pipeline
Please run `EvenPerfectBitScanner.exe /?` to display all command line parameters with a short description. The description below describes `--mersenne=residue` mode.