PrimeOrderCalculator converted to a template.

This commit is contained in:
mlemanczyk committed 2025-12-21 18:02:59 +01:00
1 parent 48df3df290
commit 2bf9ff233c
19 files changed
+3395 -4320

No files matched your search

@@ -619,11 +619,11 @@ public class Pow2ModBenchmarks
return MersenneDivisorCyclesCpu.CalculateCycleLength(modulus, divisorData); return MersenneDivisorCyclesCpu.CalculateCycleLength(modulus, divisorData);
} }
ulong order = PrimeOrderCalculator.CalculateCpu( ulong order = PrimeOrderCalculatorCpu.Calculate(
modulus, modulus,
previousPrimeOrder, previousPrimeOrder,
divisorData, divisorData,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
if (order != 0UL) if (order != 0UL)
{ {
@@ -642,10 +642,10 @@ public class Pow2ModBenchmarks
return MersenneDivisorCyclesCpu.GetCycle(modulus); return MersenneDivisorCyclesCpu.GetCycle(modulus);
} }
UInt128 order = PrimeOrderCalculator.CalculateCpu( UInt128 order = PrimeOrderCalculatorCpu.Calculate(
modulus, modulus,
previousWidePrimeOrder, previousWidePrimeOrder,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
if (order != UInt128.Zero) if (order != UInt128.Zero)
{ {
@@ -689,11 +689,11 @@ public class Pow2ModBenchmarks
#if DEBUG #if DEBUG
Console.WriteLine("Trying heuristic. Prime order calculation"); Console.WriteLine("Trying heuristic. Prime order calculation");
#endif #endif
ulong order = PrimeOrderCalculator.CalculateCpu( ulong order = PrimeOrderCalculatorCpu.Calculate(
modulus, modulus,
_previousPrimeOrder, _previousPrimeOrder,
divisorData, divisorData,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
if (order != 0UL) if (order != 0UL)
{ {
@@ -718,12 +718,12 @@ public class Pow2ModBenchmarks
#if DEBUG #if DEBUG
Console.WriteLine("Trying heuristic. Prime order calculation"); Console.WriteLine("Trying heuristic. Prime order calculation");
#endif #endif
ulong order = PrimeOrderCalculator.CalculateGpu( ulong order = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
modulus, modulus,
_previousPrimeOrder, _previousPrimeOrder,
divisorData, divisorData,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
if (order != 0UL) if (order != 0UL)
{ {
@@ -747,10 +747,10 @@ public class Pow2ModBenchmarks
#if DEBUG #if DEBUG
Console.WriteLine("Trying heuristic. Wide prime order calculation"); Console.WriteLine("Trying heuristic. Wide prime order calculation");
#endif #endif
UInt128 order = PrimeOrderCalculator.CalculateCpu( UInt128 order = PrimeOrderCalculatorCpu.Calculate(
modulus, modulus,
_previousWidePrimeOrder, _previousWidePrimeOrder,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
if (order != UInt128.Zero) if (order != UInt128.Zero)
{ {
@@ -774,11 +774,11 @@ public class Pow2ModBenchmarks
#if DEBUG #if DEBUG
Console.WriteLine("Trying heuristic. Wide prime order calculation"); Console.WriteLine("Trying heuristic. Wide prime order calculation");
#endif #endif
UInt128 order = PrimeOrderCalculator.CalculateGpu( UInt128 order = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
modulus, modulus,
_previousWidePrimeOrder, _previousWidePrimeOrder,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
if (order != UInt128.Zero) if (order != UInt128.Zero)
{ {
@@ -398,7 +398,7 @@ public class PrimeOrderGpuHeuristicsTests
PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0; PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0;
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
ulong result = PrimeOrderCalculator.CalculateGpu(gpu, prime, null, divisorData, PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); ulong result = PrimeOrderCalculatorGpu.Calculate(gpu, prime, null, divisorData, PrimeOrderCalculatorConfig.HeuristicDefault);
result.Should().Be(3UL); result.Should().Be(3UL);
} }
@@ -45,12 +45,12 @@ public class PrimeOrderCalculatorTests
ulong expected = ComputeOrderByDoubling(prime); ulong expected = ComputeOrderByDoubling(prime);
var gpu = PrimeOrderCalculatorAccelerator.Rent(1); var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
ulong result = PrimeOrderCalculator.CalculateGpu( ulong result = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
prime, prime,
previousOrder: null, previousOrder: null,
MontgomeryDivisorData.FromModulus(prime), MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
PrimeOrderCalculatorAccelerator.Return(gpu); PrimeOrderCalculatorAccelerator.Return(gpu);
@@ -65,19 +65,19 @@ public class PrimeOrderCalculatorTests
ulong previousOrder = ComputeOrderByDoubling(previousPrime); ulong previousOrder = ComputeOrderByDoubling(previousPrime);
var gpu = PrimeOrderCalculatorAccelerator.Rent(1); var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
ulong heuristic = PrimeOrderCalculator.CalculateGpu( ulong heuristic = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
prime, prime,
previousOrder, previousOrder,
MontgomeryDivisorData.FromModulus(prime), MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
ulong strict = PrimeOrderCalculator.CalculateGpu( ulong strict = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
prime, prime,
previousOrder: null, previousOrder: null,
MontgomeryDivisorData.FromModulus(prime), MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderCalculatorConfig.StrictDefault); PrimeOrderCalculatorConfig.StrictDefault);
PrimeOrderCalculatorAccelerator.Return(gpu); PrimeOrderCalculatorAccelerator.Return(gpu);
@@ -91,12 +91,12 @@ public class PrimeOrderCalculatorTests
public void Calculate_handles_trivial_primes(ulong prime, ulong expectedOrder) public void Calculate_handles_trivial_primes(ulong prime, ulong expectedOrder)
{ {
var gpu = PrimeOrderCalculatorAccelerator.Rent(1); var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
ulong result = PrimeOrderCalculator.CalculateGpu( ulong result = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
prime, prime,
previousOrder: null, previousOrder: null,
MontgomeryDivisorData.FromModulus(prime), MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
PrimeOrderCalculatorAccelerator.Return(gpu); PrimeOrderCalculatorAccelerator.Return(gpu);
result.Should().Be(expectedOrder); result.Should().Be(expectedOrder);
@@ -106,14 +106,14 @@ public class PrimeOrderCalculatorTests
[Trait("Category", "Fast")] [Trait("Category", "Fast")]
public void Calculate_heuristic_mode_returns_unresolved_when_phi_cannot_be_factored() public void Calculate_heuristic_mode_returns_unresolved_when_phi_cannot_be_factored()
{ {
var config = new PrimeOrderCalculator.PrimeOrderCalculatorConfig( var config = new PrimeOrderCalculatorConfig(
smallFactorLimit: 1, smallFactorLimit: 1,
pollardRhoMilliseconds: 0, pollardRhoMilliseconds: 0,
maxPowChecks: 8, maxPowChecks: 8,
strictMode: false); strictMode: false);
var gpu = PrimeOrderCalculatorAccelerator.Rent(1); var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
ulong result = PrimeOrderCalculator.CalculateGpu( ulong result = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
prime: 13UL, prime: 13UL,
previousOrder: null, previousOrder: null,
@@ -128,7 +128,7 @@ public class PrimeOrderCalculatorTests
[Trait("Category", "Fast")] [Trait("Category", "Fast")]
public void Calculate_heuristic_mode_attempts_candidates_before_falling_back_to_strict_mode() public void Calculate_heuristic_mode_attempts_candidates_before_falling_back_to_strict_mode()
{ {
var heuristicConfig = new PrimeOrderCalculator.PrimeOrderCalculatorConfig( var heuristicConfig = new PrimeOrderCalculatorConfig(
smallFactorLimit: 2, smallFactorLimit: 2,
pollardRhoMilliseconds: 0, pollardRhoMilliseconds: 0,
maxPowChecks: 1, maxPowChecks: 1,
@@ -138,7 +138,7 @@ public class PrimeOrderCalculatorTests
try try
{ {
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL); MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL);
ulong heuristic = PrimeOrderCalculator.CalculateGpu( ulong heuristic = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
prime: 239UL, prime: 239UL,
previousOrder: null, previousOrder: null,
@@ -147,13 +147,13 @@ public class PrimeOrderCalculatorTests
heuristic.Should().Be(ComputeOrderByDoubling(239UL)); heuristic.Should().Be(ComputeOrderByDoubling(239UL));
var strictConfig = new PrimeOrderCalculator.PrimeOrderCalculatorConfig( var strictConfig = new PrimeOrderCalculatorConfig(
smallFactorLimit: 2, smallFactorLimit: 2,
pollardRhoMilliseconds: 0, pollardRhoMilliseconds: 0,
maxPowChecks: 1, maxPowChecks: 1,
strictMode: true); strictMode: true);
ulong strict = PrimeOrderCalculator.CalculateGpu( ulong strict = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
prime: 239UL, prime: 239UL,
previousOrder: null, previousOrder: null,
@@ -175,11 +175,11 @@ public class PrimeOrderCalculatorTests
UInt128 prime = UInt128.Parse("18446744073709641691"); UInt128 prime = UInt128.Parse("18446744073709641691");
var gpu = PrimeOrderCalculatorAccelerator.Rent(1); var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
UInt128 result = PrimeOrderCalculator.CalculateGpu( UInt128 result = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
prime, prime,
previousOrder: null, previousOrder: null,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
PrimeOrderCalculatorAccelerator.Return(gpu); PrimeOrderCalculatorAccelerator.Return(gpu);
result.Should().Be((UInt128)1229782938247309446UL); result.Should().Be((UInt128)1229782938247309446UL);
@@ -10,5 +10,9 @@ public static class PrimeOrderConstants
public const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1); public const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1);
public const int Pow2WindowSizeBits = 8; public const int Pow2WindowSizeBits = 8;
public const int WideStackThreshold = 12; public const int WideStackThreshold = 12;
public const int MaxCompositeSlots = 256;
public const int ExponentSlotLimit = 256;
public const int PrimeSlotLimit = 256;
public const int ScratchSlots = 64;
} }
@@ -253,7 +253,7 @@ internal static partial class PrimeOrderGpuHeuristics
PrimeOrderCalculatorAccelerator gpu, PrimeOrderCalculatorAccelerator gpu,
ulong prime, ulong prime,
ulong? previousOrder, ulong? previousOrder,
PrimeOrderCalculator.PrimeOrderCalculatorConfig config, PrimeOrderCalculatorConfig config,
in MontgomeryDivisorData divisorData, in MontgomeryDivisorData divisorData,
out ulong order) out ulong order)
{ {
@@ -259,22 +259,23 @@ public struct MersenneDivisorCyclesTemplate
} }
#if DEVICE_GPU #if DEVICE_GPU
UInt128 wideOrder = PrimeOrderCalculator.CalculateGpu( UInt128 wideOrder = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
divisor, divisor,
previousOrder: null, previousOrder: null,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
#elif DEVICE_HYBRID #elif DEVICE_HYBRID
UInt128 wideOrder = PrimeOrderCalculator.CalculateHybrid( UInt128 wideOrder = PrimeOrderCalculatorHybrid.Calculate(
gpu, gpu,
divisor, divisor,
previousOrder: null, previousOrder: null,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); MontgomeryDivisorData.FromModulus((ulong)divisor),
PrimeOrderCalculatorConfig.HeuristicDefault);
#else #else
UInt128 wideOrder = PrimeOrderCalculator.CalculateCpu( UInt128 wideOrder = PrimeOrderCalculatorCpu.Calculate(
divisor, divisor,
previousOrder: null, previousOrder: null,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
#endif #endif
if (wideOrder != UInt128.Zero) if (wideOrder != UInt128.Zero)
{ {
@@ -329,25 +330,25 @@ public struct MersenneDivisorCyclesTemplate
// } // }
#if DEVICE_GPU #if DEVICE_GPU
ulong computedOrder = PrimeOrderCalculator.CalculateGpu( ulong computedOrder = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
divisor, divisor,
previousOrder: null, previousOrder: null,
divisorData, divisorData,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
#elif DEVICE_HYBRID #elif DEVICE_HYBRID
ulong computedOrder = PrimeOrderCalculator.CalculateHybrid( ulong computedOrder = PrimeOrderCalculatorHybrid.Calculate(
gpu, gpu,
divisor, divisor,
previousOrder: null, previousOrder: null,
divisorData, divisorData,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
#else #else
ulong computedOrder = PrimeOrderCalculator.CalculateCpu( ulong computedOrder = PrimeOrderCalculatorCpu.Calculate(
divisor, divisor,
previousOrder: null, previousOrder: null,
divisorData, divisorData,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
#endif #endif
if (computedOrder != 0UL) if (computedOrder != 0UL)
{ {
@@ -1013,12 +1014,12 @@ public struct MersenneDivisorCyclesTemplate
#if DEVICE_GPU #if DEVICE_GPU
if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.IsPrimeGpu(gpu, divisor)) if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.IsPrimeGpu(gpu, divisor))
{ {
ulong computedOrder = PrimeOrderCalculator.CalculateGpu( ulong computedOrder = PrimeOrderCalculatorGpu.Calculate(
gpu, gpu,
divisor, divisor,
previousOrder: null, previousOrder: null,
divisorData, divisorData,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
if (computedOrder != 0UL) if (computedOrder != 0UL)
{ {
cycleLength = computedOrder; cycleLength = computedOrder;
@@ -1028,12 +1029,12 @@ public struct MersenneDivisorCyclesTemplate
#elif DEVICE_HYBRID #elif DEVICE_HYBRID
if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.IsPrimeCpu(divisor)) if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.IsPrimeCpu(divisor))
{ {
ulong computedOrder = PrimeOrderCalculator.CalculateHybrid( ulong computedOrder = PrimeOrderCalculatorHybrid.Calculate(
gpu, gpu,
divisor, divisor,
previousOrder: null, previousOrder: null,
divisorData, divisorData,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
if (computedOrder != 0UL) if (computedOrder != 0UL)
{ {
cycleLength = computedOrder; cycleLength = computedOrder;
@@ -1043,11 +1044,11 @@ public struct MersenneDivisorCyclesTemplate
#else #else
if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.IsPrimeCpu(divisor)) if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.IsPrimeCpu(divisor))
{ {
ulong computedOrder = PrimeOrderCalculator.CalculateCpu( ulong computedOrder = PrimeOrderCalculatorCpu.Calculate(
divisor, divisor,
previousOrder: null, previousOrder: null,
divisorData, divisorData,
PrimeOrderCalculator.PrimeOrderCalculatorConfig.HeuristicDefault); PrimeOrderCalculatorConfig.HeuristicDefault);
if (computedOrder != 0UL) if (computedOrder != 0UL)
{ {
cycleLength = computedOrder; cycleLength = computedOrder;
@@ -23,6 +23,7 @@ public static class PerfectNumberConstants
public const int ThreadsByAccelerator = 1; public const int ThreadsByAccelerator = 1;
public const int ByDivisorStateSaveInterval = 200_000; public const int ByDivisorStateSaveInterval = 200_000;
public const string ByDivisorStateDirectory = "Checks"; public const string ByDivisorStateDirectory = "Checks";
public const int PollardRhoMilliseconds = 100;
// TODO: Load these limits from the benchmark-driven configuration so CPU and GPU scans stay aligned with the optimal // 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. // divisor-cycle datasets we generate offline.
// TODO: Promote these magic numbers into a runtime profile derived from EvenPerfectBitScanner.Benchmarks so we can retune // TODO: Promote these magic numbers into a runtime profile derived from EvenPerfectBitScanner.Benchmarks so we can retune
@@ -5,6 +5,7 @@
<ImplicitUsings>enable</ImplicitUsings> <ImplicitUsings>enable</ImplicitUsings>
<Nullable>enable</Nullable> <Nullable>enable</Nullable>
<LangVersion>12.0</LangVersion> <LangVersion>12.0</LangVersion>
<DefineConstants>$(DefineConstants);DEVICE_CPU</DefineConstants>
<EmitCompilerGeneratedFiles>true</EmitCompilerGeneratedFiles> <EmitCompilerGeneratedFiles>true</EmitCompilerGeneratedFiles>
<CompilerGeneratedFilesOutputPath>$(IntermediateOutputPath)generated</CompilerGeneratedFilesOutputPath> <CompilerGeneratedFilesOutputPath>$(IntermediateOutputPath)generated</CompilerGeneratedFilesOutputPath>
</PropertyGroup> </PropertyGroup>
@@ -0,0 +1,6 @@
namespace PerfectNumbers.Core;
public static class PrimeNumberConstants
{
public const ulong GpuSmallPrimeBound = 1_000_000UL;
}
@@ -1,39 +0,0 @@
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core;
internal static partial class PrimeOrderCalculator
{
internal readonly struct PrimeOrderCalculatorConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, bool strictMode)
{
private static PrimeOrderCalculatorConfig _heuristicDefault = CreateHeuristicDefault(UnboundedTaskScheduler.ConfiguredThreadCount);
public static PrimeOrderCalculatorConfig HeuristicDefault
{
[MethodImpl(MethodImplOptions.AggressiveInlining)]
get => _heuristicDefault;
}
public static readonly PrimeOrderCalculatorConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, true);
public static void ConfigureHeuristicDefault(int threadCount)
{
_heuristicDefault = CreateHeuristicDefault(threadCount);
}
public static PrimeOrderCalculatorConfig CreateHeuristicDefault(int threadCount)
{
int normalizedThreadCount = threadCount < 1 ? 1 : threadCount;
return new(smallFactorLimit: 100_000, pollardRhoMilliseconds: normalizedThreadCount * 24, maxPowChecks: 24, false);
}
public readonly uint SmallFactorLimit = smallFactorLimit;
public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds;
public readonly int MaxPowChecks = maxPowChecks;
public readonly int MaxPowChecksCapacity = maxPowChecks << 2;
public readonly bool StrictMode = strictMode;
}
// TODO: Remove branching on the CPU / GPU paths
// TODO: Split big non-static functions into smaller / extract static code to limit JIT time
}
File diff suppressed because it is too large. Load diff
File diff suppressed because it is too large. Load diff
@@ -1,651 +0,0 @@
using System.Buffers;
using System.Diagnostics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
using ILGPU;
using PerfectNumbers.Core.Gpu;
using PerfectNumbers.Core.Gpu.Accelerators;
namespace PerfectNumbers.Core;
internal static partial class PrimeOrderCalculator
{
public static ulong CalculateHybrid(
PrimeOrderCalculatorAccelerator gpu,
ulong prime,
ulong? previousOrder,
in MontgomeryDivisorData divisorData,
in PrimeOrderCalculatorConfig config)
{
// TODO: Is this condition ever met on EvenPerfectBitScanner's execution path? If not, we can add a clarification comment and comment out the entire block. We want to support p candidates at least greater or equal to 31.
if (prime <= 3UL)
{
return prime == 3UL ? 2UL : 1UL;
}
ulong phi = prime - 1UL;
PartialFactorResult phiFactors = PartialFactorHybrid(gpu, phi, divisorData, config);
ulong result;
if (phiFactors.Factors is null)
{
result = CalculateByFactorizationCpu(prime, divisorData);
phiFactors.Dispose();
return result;
}
result = RunHeuristicPipelineHybrid(gpu, prime, previousOrder, config, divisorData, phi, phiFactors);
phiFactors.Dispose();
return result;
}
public static UInt128 CalculateHybrid(
PrimeOrderCalculatorAccelerator gpu,
in UInt128 prime,
in UInt128? previousOrder,
in PrimeOrderCalculatorConfig config)
{
MontgomeryDivisorData divisorData;
UInt128 result;
if (prime <= ulong.MaxValue)
{
ulong? previous = null;
if (previousOrder.HasValue)
{
UInt128 previousValue = previousOrder.Value;
if (previousValue <= ulong.MaxValue)
{
previous = (ulong)previousValue;
}
else
{
previous = ulong.MaxValue;
}
}
ulong prime64 = (ulong)prime;
divisorData = MontgomeryDivisorData.FromModulus(prime64);
ulong order64 = CalculateHybrid(gpu, prime64, previous, divisorData, config);
result = order64 == 0UL ? UInt128.Zero : (UInt128)order64;
}
else
{
divisorData = MontgomeryDivisorData.Empty;
result = CalculateWideInternalCpu(prime, previousOrder, divisorData, config);
}
return result;
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static void CollectFactorsHybrid(PrimeOrderCalculatorAccelerator gpu, Span<ulong> primeSlots, Span<int> exponentSlots, ref int factorCount, List<PartialFactorPendingEntry> pending, FixedCapacityStack<ulong> compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached)
{
while (compositeStack.Count > 0)
{
ulong composite = compositeStack.Pop();
bool isPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite);
if (isPrime)
{
AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite);
continue;
}
pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp;
if (pollardRhoDeadlineReached)
{
pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true));
continue;
}
if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor))
{
pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true));
continue;
}
ulong quotient = composite / factor;
compositeStack.Push(factor);
compositeStack.Push(quotient);
}
pollardRhoDeadlineReached = false;
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static ulong InitializeStartingOrderHybrid(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong phi, in MontgomeryDivisorData divisorData)
{
ulong order = phi;
if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL)
{
ulong half = phi >> 1;
if (Pow2EqualsOneHybrid(gpu, half, prime, divisorData))
{
order = half;
}
}
return order;
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static PartialFactorResult PartialFactorHybrid(PrimeOrderCalculatorAccelerator gpu, ulong value, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config)
{
if (value <= 1UL)
{
return PartialFactorResult.Empty;
}
// stackalloc is faster than pooling
var result = PartialFactorResult.Rent(divisorData);
ulong[] factors = result.Factors;
Span<ulong> primeSlots = new(factors);
int[] exponents = result.Exponents;
Span<int> exponentSlots = new(exponents);
// We don't need to worry about leftovers, because we always use indexes within the calculated counts
// primeSlots.Clear();
// exponentSlots.Clear();
uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
List<PartialFactorPendingEntry> pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2);
PopulateSmallPrimeFactorsCpu(
value,
limit,
primeSlots,
exponentSlots,
out int factorCount,
out ulong remaining);
bool limitReached;
if (remaining > 1UL)
{
if (config.PollardRhoMilliseconds > 0)
{
long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds);
FixedCapacityStack<ulong> compositeStack = ThreadStaticPools.RentUlongStack(4);
compositeStack.Push(remaining);
CollectFactorsHybrid(gpu, primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached);
if (limitReached)
{
while (compositeStack.Count > 0)
{
pending.Add(new PartialFactorPendingEntry(compositeStack.Pop(), knownComposite: false));
}
}
compositeStack.Clear();
ThreadStaticPools.ReturnUlongStack(compositeStack);
}
else
{
pending.Add(new PartialFactorPendingEntry(remaining, knownComposite: false));
}
}
ulong cofactor = 1UL;
bool cofactorContainsComposite = false;
int pendingCount = pending.Count;
int index = 0;
for (; index < pendingCount; index++)
{
PartialFactorPendingEntry entry = pending[index];
ulong composite = entry.Value;
if (entry.KnownComposite)
{
cofactor = checked(cofactor * composite);
cofactorContainsComposite = true;
continue;
}
if (!entry.HasKnownPrimality)
{
bool isPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite);
entry.WithPrimality(isPrime);
pending[index] = entry;
}
if (entry.IsPrime)
{
AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite);
}
// composite is never smaller on the execution path
// else if (composite > 1UL)
else
{
cofactor = checked(cofactor * composite);
cofactorContainsComposite = true;
}
}
bool cofactorIsPrime;
if (cofactor <= 1UL)
{
cofactorIsPrime = false;
}
else if (cofactorContainsComposite)
{
cofactorIsPrime = false;
}
else
{
cofactorIsPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(cofactor) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, cofactor);
}
// This will never happen in production code. We'll always get at least 1 factor
// if (factorCount == 0)
// {
// if (cofactor == value)
// {
// result = PartialFactorResult.Rent(null, value, false, 0);
// goto ReturnResult;
// }
// result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
// goto ReturnResult;
// }
Array.Sort(factors, exponents, 0, factorCount);
result.Cofactor = cofactor;
result.FullyFactored = cofactor == 1UL;
result.Count = factorCount;
result.CofactorIsPrime = cofactorIsPrime;
pending.Clear();
ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending);
return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
private static bool Pow2EqualsOneHybrid(PrimeOrderCalculatorAccelerator gpu, ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
{
// if (exponent >= PerfectNumberConstants.MaxQForDivisorCycles || prime >= PerfectNumberConstants.MaxQForDivisorCycles)
// {
// // Atomic.Add(ref _pow2GpuHits, 1);
// // Console.WriteLine($"pow2 GPU Hits {_pow2GpuHits}");
// ulong remainder = exponent.Pow2MontgomeryModWindowedConvertGpu(gpu, divisorData);
// // GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData);
// // if (status == GpuPow2ModStatus.Success)
// // {
// return remainder == 1UL;
// // }
// }
// Atomic.Add(ref _pow2CpuHits, 1);
// Console.WriteLine($"pow2 CPU Hits {_pow2CpuHits}");
return exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL;
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static ulong RunHeuristicPipelineHybrid(
PrimeOrderCalculatorAccelerator gpu,
ulong prime,
ulong? previousOrder,
in PrimeOrderCalculatorConfig config,
in MontgomeryDivisorData divisorData,
ulong phi,
PartialFactorResult phiFactors)
{
if (phiFactors.FullyFactored && TrySpecialMaxHybrid(gpu, phi, prime, phiFactors, divisorData))
{
return phi;
}
ulong candidateOrder = InitializeStartingOrderHybrid(gpu, prime, phi, divisorData);
candidateOrder = ExponentLoweringCpu(candidateOrder, phiFactors, divisorData);
if (TryConfirmOrderHybrid(gpu, prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors))
{
return candidateOrder;
}
if (config.StrictMode)
{
orderFactors?.Dispose();
return CalculateByFactorizationCpu(prime, divisorData);
}
if (TryHeuristicFinishHybrid(gpu, prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, out ulong order))
{
return order;
}
return candidateOrder;
}
private static int _cpuCount;
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static bool RunOnCpu()
{
int gpuRatio = EnvironmentConfiguration.GpuRatio;
if (gpuRatio == 0)
{
return false;
}
if (gpuRatio == 1)
{
while (true)
{
if (Interlocked.CompareExchange(ref _cpuCount, 0, 1) == 0)
{
return true;
}
if (Interlocked.CompareExchange(ref _cpuCount, 1, 0) == 1)
{
return false;
}
Thread.Yield();
}
}
int cpuCount = Interlocked.Increment(ref _cpuCount);
if (cpuCount == gpuRatio)
{
Interlocked.Add(ref _cpuCount, -gpuRatio);
}
else if (cpuCount > gpuRatio)
{
cpuCount -= gpuRatio;
}
return cpuCount != gpuRatio;
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static bool TryConfirmCandidateHybrid(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget)
{
PartialFactorResult factorization = PartialFactorHybrid(gpu, candidate, divisorData, config);
if (factorization.Factors is null)
{
factorization.Dispose();
return false;
}
if (!factorization.FullyFactored)
{
if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime)
{
factorization.Dispose();
return false;
}
factorization.WithAdditionalPrime(factorization.Cofactor);
}
ReadOnlySpan<ulong> factorsSpan = factorization.Factors;
ReadOnlySpan<int> exponentsSpan = factorization.Exponents;
int length = factorization.Count;
const int StackExponentCapacity = 32;
const int ExponentHardLimit = 256;
Span<ulong> stackBuffer = stackalloc ulong[StackExponentCapacity];
ulong[]? heapCandidateArray = null;
ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
bool violates = false;
for (int i = 0; i < length; i++)
{
ulong primeFactor = factorsSpan[i];
int exponent = exponentsSpan[i];
// Factorization entries for candidate confirmation are always positive here.
// if (exponent <= 0)
// {
// continue;
// }
if (exponent > ExponentHardLimit)
{
throw new InvalidOperationException($"Candidate factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}.");
}
if (exponent <= StackExponentCapacity)
{
if (CheckCandidateViolation(stackBuffer[..exponent], primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper))
{
violates = true;
break;
}
continue;
}
heapCandidateArray ??= ThreadStaticPools.UlongPool.Rent(ExponentHardLimit);
if (CheckCandidateViolation(heapCandidateArray.AsSpan(0, exponent), primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper))
{
violates = true;
break;
}
}
ThreadStaticPools.ReturnExponentStepperCpu(stepper);
if (heapCandidateArray is not null)
{
ThreadStaticPools.UlongPool.Return(heapCandidateArray);
}
factorization.Dispose();
return !violates;
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static bool TryConfirmOrderHybrid(
PrimeOrderCalculatorAccelerator gpu,
ulong prime,
ulong order,
in MontgomeryDivisorData divisorData,
in PrimeOrderCalculatorConfig config,
out PartialFactorResult? reusableFactorization)
{
reusableFactorization = null;
if (order == 0UL)
{
return false;
}
if (!Pow2EqualsOneHybrid(gpu, order, prime, divisorData))
{
return false;
}
PartialFactorResult factorization = PartialFactorHybrid(gpu, order, divisorData, config);
if (factorization.Factors is null)
{
factorization.Dispose();
return false;
}
if (!factorization.FullyFactored)
{
if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime)
{
factorization.Dispose();
return false;
}
factorization.WithAdditionalPrime(factorization.Cofactor);
}
if (!ValidateOrderAgainstFactors(prime, order, divisorData, factorization))
{
reusableFactorization = factorization;
return false;
}
factorization.Dispose();
return true;
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static bool TryHeuristicFinishHybrid(
PrimeOrderCalculatorAccelerator gpu,
ulong prime,
ulong order,
ulong? previousOrder,
in MontgomeryDivisorData divisorData,
in PrimeOrderCalculatorConfig config,
PartialFactorResult phiFactors,
PartialFactorResult? cachedOrderFactors,
out ulong result)
{
result = 0UL;
if (order <= 1UL)
{
cachedOrderFactors?.Dispose();
return false;
}
// Reuse the partial factorization from TryConfirmOrderCpu when available.
PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorHybrid(gpu, order, divisorData, config);
try
{
if (orderFactors.Factors is null)
{
return false;
}
if (!orderFactors.FullyFactored)
{
if (orderFactors.Cofactor <= 1UL || !orderFactors.CofactorIsPrime)
{
return false;
}
orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
}
int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity;
List<ulong> candidates = ThreadStaticPools.RentUlongList(capacity);
candidates.Clear();
ulong[] factorArray = orderFactors.Factors!;
int[] exponentsArray = orderFactors.Exponents!;
BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity);
int candidateCount = candidates.Count;
if (candidateCount == 0)
{
ThreadStaticPools.ReturnUlongList(candidates);
return false;
}
SortCandidates(prime, previousOrder, candidates);
int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks;
int powUsed = 0;
Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
bool powStepperInitialized = false;
int index = 0;
Span<ulong> stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount];
ArrayPool<ulong> pool = ThreadStaticPools.UlongPool;
while (index < candidateCount && powUsed < powBudget)
{
int remaining = candidateCount - index;
int budgetRemaining = powBudget - powUsed;
int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, PrimeOrderConstants.MaxGpuBatchSize));
if (batchSize <= 0)
{
break;
}
ReadOnlySpan<ulong> batch = candidateSpan.Slice(index, batchSize);
int startIndex = 0;
if (!powStepperInitialized)
{
ulong candidate = batch[0];
powUsed++;
bool equalsOne = powStepper.InitializeCpuIsUnity(candidate);
powStepperInitialized = true;
if (equalsOne && TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget))
{
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
result = candidate;
return true;
}
startIndex = 1;
}
for (int i = startIndex; i < batchSize && powUsed < powBudget; i++)
{
ulong candidate = batch[i];
powUsed++;
bool equalsOne = powStepper.ComputeNextIsUnity(candidate);
if (!equalsOne)
{
continue;
}
if (!TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget))
{
continue;
}
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
result = candidate;
return true;
}
index += batchSize;
}
candidates.Clear();
ThreadStaticPools.ReturnUlongList(candidates);
ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
return false;
}
finally
{
orderFactors.Dispose();
}
}
[MethodImpl(MethodImplOptions.AggressiveOptimization)]
private static bool TrySpecialMaxHybrid(PrimeOrderCalculatorAccelerator gpu, ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
{
int length = factors.Count;
// The phi factorization on the scanner path always yields at least one entry for phi >= 2,
// preventing the zero-length case from occurring outside synthetic tests.
// if (length == 0)
// {
// return true;
// }
ReadOnlySpan<ulong> factorSpan = new(factors.Factors, 0, length);
Span<ulong> stackBuffer = stackalloc ulong[length];
if (length <= 8)
{
return EvaluateSpecialMaxCandidatesCpu(stackBuffer, factorSpan, phi, prime, divisorData);
}
return EvaluateSpecialMaxCandidatesGpu(gpu, factorSpan, phi, prime, divisorData);
}
}
@@ -1,631 +0,0 @@
using System.Buffers;
using System.Diagnostics;
using System.Numerics;
namespace PerfectNumbers.Core;
internal static partial class PrimeOrderCalculator
{
private static void SortCandidatesWide(UInt128 prime, in UInt128? previousOrder, List<UInt128> candidates)
{
UInt128 previous = previousOrder ?? UInt128.Zero;
int previousGroup = previousOrder.HasValue ? GetGroupWide(previousOrder.Value, prime) : 1;
bool hasPrevious = previousOrder.HasValue;
candidates.Sort((x, y) =>
{
CandidateKey128 keyX = BuildKeyWide(x, prime, previous, previousGroup, hasPrevious);
CandidateKey128 keyY = BuildKeyWide(y, prime, previous, previousGroup, hasPrevious);
int primary = keyX.Primary.CompareTo(keyY.Primary);
if (primary != 0)
{
return primary;
}
int secondary = CompareComponents(keyX.SecondaryDescending, keyX.Secondary, keyY.SecondaryDescending, keyY.Secondary);
if (secondary != 0)
{
return secondary;
}
return CompareComponents(keyX.TertiaryDescending, keyX.Tertiary, keyY.TertiaryDescending, keyY.Tertiary);
});
}
private static CandidateKey128 BuildKeyWide(
in UInt128 value,
in UInt128 prime,
in UInt128 previous,
int previousGroup,
bool hasPrevious)
{
int group = GetGroupWide(value, prime);
if (group == 0)
{
return new CandidateKey128(int.MaxValue, false, UInt128.Zero, false, UInt128.Zero);
}
bool isGe = !hasPrevious || value >= previous;
int primary = ComputePrimary(group, isGe, previousGroup);
if (group == 3)
{
return new CandidateKey128(primary, true, value, true, value);
}
UInt128 distance = hasPrevious ? (value > previous ? value - previous : previous - value) : value;
return new CandidateKey128(primary, false, distance, false, value);
}
private static int CompareComponents(bool descendingX, in UInt128 valueX, bool descendingY, in UInt128 valueY)
{
if (descendingX == descendingY)
{
if (valueX == valueY)
{
return 0;
}
if (descendingX)
{
return valueX > valueY ? -1 : 1;
}
return valueX > valueY ? 1 : -1;
}
if (valueX == valueY)
{
return 0;
}
return valueX > valueY ? 1 : -1;
}
private static int GetGroupWide(in UInt128 value, in UInt128 prime)
{
UInt128 threshold1 = prime >> 3;
if (value <= threshold1)
{
return 1;
}
UInt128 threshold2 = prime >> 2;
if (value <= threshold2)
{
return 2;
}
UInt128 threshold3 = (UInt128)(((BigInteger)prime * 3) >> 3);
if (value <= threshold3)
{
return 3;
}
return 0;
}
private static void BuildCandidatesWide(in UInt128 order, FactorEntry128[] factors, int count, List<UInt128> candidates, int limit)
{
if (count == 0)
{
return;
}
Span<FactorEntry128> buffer = factors.AsSpan(0, count);
buffer.Sort(static (a, b) => a.Value.CompareTo(b.Value));
BuildCandidatesRecursiveWide(order, buffer, 0, UInt128.One, candidates, limit);
}
private static void BuildCandidatesRecursiveWide(
in UInt128 order,
in ReadOnlySpan<FactorEntry128> factors,
int index,
in UInt128 divisorProduct,
List<UInt128> candidates,
int limit)
{
if (candidates.Count >= limit)
{
return;
}
if (index >= factors.Length)
{
if (divisorProduct == UInt128.One || divisorProduct == order)
{
return;
}
UInt128 candidate = order / divisorProduct;
if (candidate > UInt128.One && candidate < order)
{
candidates.Add(candidate);
}
return;
}
FactorEntry128 factor = factors[index];
UInt128 primeFactor = factor.Value;
UInt128 contribution = UInt128.One;
for (int exponent = 0; exponent <= factor.Exponent; exponent++)
{
UInt128 nextDivisor = divisorProduct * contribution;
if (nextDivisor > order)
{
break;
}
BuildCandidatesRecursiveWide(order, factors, index + 1, nextDivisor, candidates, limit);
if (candidates.Count >= limit)
{
return;
}
if (exponent == factor.Exponent)
{
break;
}
if (contribution > order / primeFactor)
{
break;
}
contribution *= primeFactor;
}
}
private static PartialFactorResult128 PartialFactorWide(in UInt128 value, in PrimeOrderCalculatorConfig config)
{
if (value <= UInt128.One)
{
return PartialFactorResult128.Empty;
}
Dictionary<UInt128, int> counts = new(capacity: 8);
uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
UInt128 remaining;
// if (IsGpuPow2Allowed && value <= ulong.MaxValue)
// {
// Dictionary<ulong, int> narrowCounts = new(capacity: 8);
// if (TryPopulateSmallPrimeFactorsGpu((ulong)value, limit, narrowCounts, out ulong narrowRemaining))
// {
// foreach (KeyValuePair<ulong, int> entry in narrowCounts)
// {
// counts[(UInt128)entry.Key] = entry.Value;
// }
// remaining = (UInt128)narrowRemaining;
// }
// else
// {
// remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts);
// }
// }
// else
// {
remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts);
// }
List<UInt128> pending = new();
if (remaining > UInt128.One)
{
pending.Add(remaining);
}
if (config.PollardRhoMilliseconds > 0 && pending.Count > 0)
{
long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds);
Stack<UInt128> stack = new();
stack.Push(remaining);
pending.Clear();
long timestamp = 0L; // reused for deadline checks.
while (stack.Count > 0)
{
UInt128 composite = stack.Pop();
if (composite <= UInt128.One)
{
continue;
}
if (IsPrimeWide(composite))
{
AddFactor(counts, composite, 1);
continue;
}
timestamp = Stopwatch.GetTimestamp();
if (timestamp > deadlineTimestamp)
{
pending.Add(composite);
continue;
}
if (!TryPollardRhoWide(composite, deadlineTimestamp, out UInt128 factor))
{
pending.Add(composite);
continue;
}
UInt128 quotient = composite / factor;
stack.Push(factor);
stack.Push(quotient);
}
}
UInt128 cofactor = UInt128.One;
int pendingCount = pending.Count;
for (int i = 0; i < pendingCount; i++)
{
UInt128 composite = pending[i];
if (IsPrimeWide(composite))
{
AddFactor(counts, composite, 1);
}
else
{
BigInteger product = (BigInteger)cofactor * (BigInteger)composite;
cofactor = (UInt128)product;
}
}
if (counts.Count == 0 && cofactor == value)
{
return new PartialFactorResult128(null, value, false, 0);
}
ArrayPool<FactorEntry128>? pool = null;
FactorEntry128[]? rented = counts.Count > 0 ? (pool = ThreadStaticPools.FactorEntry128Pool).Rent(counts.Count) : null;
int index = 0;
if (rented is not null)
{
foreach (KeyValuePair<UInt128, int> entry in counts)
{
rented[index] = new FactorEntry128(entry.Key, entry.Value);
index++;
}
}
FactorEntry128[]? resultArray = null;
if (pool is not null)
{
Span<FactorEntry128> span = rented.AsSpan(0, index);
span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
resultArray = new FactorEntry128[index];
span.CopyTo(resultArray);
pool.Return(rented!, clearArray: false);
}
bool fullyFactored = cofactor == UInt128.One;
return new PartialFactorResult128(resultArray, cofactor, fullyFactored, index);
}
private static UInt128 PopulateSmallPrimeFactorsCpuWide(in UInt128 value, uint limit, Dictionary<UInt128, int> counts)
{
UInt128 remaining = value;
uint[] primes = PrimesGenerator.SmallPrimes;
ulong[] squares = PrimesGenerator.SmallPrimesPow2;
int primeCount = primes.Length;
for (int i = 0; i < primeCount; i++)
{
uint primeCandidate = primes[i];
if (primeCandidate > limit)
{
break;
}
UInt128 primeSquare = (UInt128)squares[i];
if (primeSquare > remaining)
{
break;
}
UInt128 primeValue = primeCandidate;
if ((remaining % primeValue) != UInt128.Zero)
{
continue;
}
int exponent = 0;
do
{
remaining /= primeValue;
exponent++;
}
while ((remaining % primeValue) == UInt128.Zero);
counts[primeValue] = exponent;
}
return remaining;
}
private static bool TryPollardRhoWide(in UInt128 n, long deadlineTimestamp, out UInt128 factor)
{
factor = UInt128.Zero;
if ((n & UInt128.One) == UInt128.Zero)
{
factor = (UInt128)2UL;
return true;
}
UInt128 c = UInt128.One;
UInt128 x = (UInt128)2UL;
UInt128 y = x;
long timestamp = 0L; // reused for deadline checks.
while (true)
{
timestamp = Stopwatch.GetTimestamp();
if (timestamp > deadlineTimestamp)
{
return false;
}
x = AdvancePolynomialWide(x, c, n);
y = AdvancePolynomialWide(y, c, n);
y = AdvancePolynomialWide(y, c, n);
UInt128 diff = x > y ? x - y : y - x;
UInt128 d = BinaryGcdWide(diff, n);
if (d == n)
{
c += UInt128.One;
x = (UInt128)2UL;
y = x;
continue;
}
if (d > UInt128.One)
{
factor = d;
return true;
}
}
}
private static UInt128 AdvancePolynomialWide(in UInt128 x, in UInt128 c, in UInt128 modulus)
{
BigInteger value = (BigInteger)x;
value = (value * value + (BigInteger)c) % (BigInteger)modulus;
return (UInt128)value;
}
private static UInt128 BinaryGcdWide(UInt128 a, UInt128 b)
{
if (a == UInt128.Zero)
{
return b;
}
if (b == UInt128.Zero)
{
return a;
}
int shift = TrailingZeroCountWide(a | b);
a >>= TrailingZeroCountWide(a);
while (true)
{
b >>= TrailingZeroCountWide(b);
if (a > b)
{
(a, b) = (b, a);
}
b -= a;
if (b == UInt128.Zero)
{
return a << shift;
}
}
}
private static int TrailingZeroCountWide(in UInt128 value)
{
if (value == UInt128.Zero)
{
return 128;
}
ulong low = (ulong)value;
if (low != 0UL)
{
return BitOperations.TrailingZeroCount(low);
}
ulong high = (ulong)(value >> 64);
return 64 + BitOperations.TrailingZeroCount(high);
}
private static void FactorCompletelyWide(in UInt128 value, Dictionary<UInt128, int> counts)
{
if (value <= UInt128.One)
{
return;
}
if (IsPrimeWide(value))
{
AddFactor(counts, value, 1);
return;
}
UInt128 factor = PollardRhoWide(value);
FactorCompletelyWide(factor, counts);
FactorCompletelyWide(value / factor, counts);
}
private static UInt128 PollardRhoWide(in UInt128 n)
{
if ((n & UInt128.One) == UInt128.Zero)
{
return (UInt128)2UL;
}
UInt128 c = UInt128.One;
while (true)
{
UInt128 x = (UInt128)2UL;
UInt128 y = x;
UInt128 d = UInt128.One;
while (d == UInt128.One)
{
x = AdvancePolynomialWide(x, c, n);
y = AdvancePolynomialWide(y, c, n);
y = AdvancePolynomialWide(y, c, n);
UInt128 diff = x > y ? x - y : y - x;
d = BinaryGcdWide(diff, n);
}
if (d != n)
{
return d;
}
c += UInt128.One;
}
}
private static bool IsPrimeWide(in UInt128 value)
{
if (value <= ulong.MaxValue)
{
// return HeuristicPrimeTester.Exclusive.IsPrimeCpu((ulong)value, CancellationToken.None);
return PrimeTester.IsPrimeCpu((ulong)value);
}
if ((value & UInt128.One) == UInt128.Zero)
{
return false;
}
uint[] smallPrimes = PrimesGenerator.SmallPrimes;
int length = smallPrimes.Length;
for (int i = 0; i < length; i++)
{
uint prime = smallPrimes[i];
UInt128 primeValue = prime;
UInt128 primeSquare = primeValue * primeValue;
if (primeSquare > value)
{
break;
}
if ((value % primeValue) == UInt128.Zero)
{
return false;
}
}
BigInteger n = (BigInteger)value;
BigInteger d = n - BigInteger.One;
int s = 0;
while ((d & 1) == 0)
{
d >>= 1;
s++;
}
int[] bases = new[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
for (int i = 0; i < bases.Length; i++)
{
int baseValue = bases[i];
if (baseValue <= 1)
{
continue;
}
BigInteger a = baseValue;
if (a >= n)
{
continue;
}
BigInteger x = BigInteger.ModPow(a, d, n);
if (x == BigInteger.One || x == n - BigInteger.One)
{
continue;
}
bool witnessFound = true;
for (int r = 1; r < s; r++)
{
x = BigInteger.ModPow(x, 2, n);
if (x == n - BigInteger.One)
{
witnessFound = false;
break;
}
}
if (witnessFound)
{
return false;
}
}
return true;
}
private static void AddFactor(Dictionary<UInt128, int> counts, in UInt128 prime, int exponent)
{
if (counts.TryGetValue(prime, out int existing))
{
counts[prime] = existing + exponent;
}
else
{
counts[prime] = exponent;
}
}
private readonly struct CandidateKey128(int primary, bool secondaryDescending, in UInt128 secondary, bool tertiaryDescending, in UInt128 tertiary)
{
public int Primary { get; } = primary;
public bool SecondaryDescending { get; } = secondaryDescending;
public UInt128 Secondary { get; } = secondary;
public bool TertiaryDescending { get; } = tertiaryDescending;
public UInt128 Tertiary { get; } = tertiary;
}
private readonly struct PartialFactorResult128(FactorEntry128[]? factors, UInt128 cofactor, bool fullyFactored, int count)
{
public FactorEntry128[]? Factors { get; } = factors;
public UInt128 Cofactor { get; } = cofactor;
public bool FullyFactored { get; } = fullyFactored;
public int Count { get; } = count;
public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0);
public PartialFactorResult128 WithAdditionalPrime(UInt128 prime)
{
if (Factors is null)
{
FactorEntry128[] local = new FactorEntry128[1];
local[0] = new FactorEntry128(prime, 1);
return new PartialFactorResult128(local, UInt128.One, true, 1);
}
FactorEntry128[] extended = new FactorEntry128[Count + 1];
Array.Copy(Factors, extended, Count);
extended[Count] = new FactorEntry128(prime, 1);
return new PartialFactorResult128(extended, UInt128.One, true, Count + 1);
}
}
}
@@ -0,0 +1,34 @@
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core;
public readonly struct PrimeOrderCalculatorConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, bool strictMode, int maxSeconds = 0)
{
private static PrimeOrderCalculatorConfig _heuristicDefault = CreateHeuristicDefault(UnboundedTaskScheduler.ConfiguredThreadCount);
public static PrimeOrderCalculatorConfig HeuristicDefault
{
[MethodImpl(MethodImplOptions.AggressiveInlining)]
get => _heuristicDefault;
}
public static readonly PrimeOrderCalculatorConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, true);
public static void ConfigureHeuristicDefault(int threadCount)
{
_heuristicDefault = CreateHeuristicDefault(threadCount);
}
public static PrimeOrderCalculatorConfig CreateHeuristicDefault(int threadCount)
{
int normalizedThreadCount = threadCount < 1 ? 1 : threadCount;
return new(smallFactorLimit: 100_000, pollardRhoMilliseconds: normalizedThreadCount * 24, maxPowChecks: 24, false);
}
public readonly uint SmallFactorLimit = smallFactorLimit;
public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds;
public readonly int MaxPowChecks = maxPowChecks;
public readonly int MaxPowChecksCapacity = maxPowChecks << 2;
public readonly bool StrictMode = strictMode;
public readonly int MaxSeconds = maxSeconds;
}
File diff suppressed because it is too large. Load diff
+25
View File
@@ -100,6 +100,31 @@ namespace PerfectNumbers.Core
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static void Return(this FixedCapacityStack<Dictionary<UInt128, int>> pool, Dictionary<UInt128, int> dictionary) => pool.Push(dictionary); public static void Return(this FixedCapacityStack<Dictionary<UInt128, int>> pool, Dictionary<UInt128, int> dictionary) => pool.Push(dictionary);
[ThreadStatic]
private static FixedCapacityStack<Dictionary<ulong, int>>? _ulongIntDictionaryPool;
public static FixedCapacityStack<Dictionary<ulong, int>> UlongIntDictionaryPool
{
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
get => _ulongIntDictionaryPool ??= new(PerfectNumberConstants.DefaultPoolCapacity);
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static Dictionary<ulong, int> Rent(this FixedCapacityStack<Dictionary<ulong, int>> pool, int capacity)
{
if (pool.Pop() is { } pooled)
{
pooled.EnsureCapacity(capacity);
pooled.Clear();
return pooled;
}
return new(capacity);
}
[MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
public static void Return(this FixedCapacityStack<Dictionary<ulong, int>> pool, Dictionary<ulong, int> dictionary) => pool.Push(dictionary);
[ThreadStatic] [ThreadStatic]
private static ArrayPool<FactorEntry>? _factorEntryPool; private static ArrayPool<FactorEntry>? _factorEntryPool;
+1
View File
@@ -5,6 +5,7 @@ namespace PerfectNumbers.Core
public static class UInt128Numbers public static class UInt128Numbers
{ {
public static readonly UInt128 Two = 2UL; public static readonly UInt128 Two = 2UL;
public static readonly UInt128 Three = 3UL;
public static readonly UInt128 Four = 4UL; public static readonly UInt128 Four = 4UL;
public static readonly UInt128 SixtyFour = 64UL; public static readonly UInt128 SixtyFour = 64UL;
public static readonly UInt128 OneHundredTwentySeven = 127UL; public static readonly UInt128 OneHundredTwentySeven = 127UL;
@@ -28,7 +28,7 @@ public sealed class UnboundedTaskScheduler : TaskScheduler
public static void ConfigureThreadCount(int threadCount) public static void ConfigureThreadCount(int threadCount)
{ {
_configuredThreadCount = threadCount; _configuredThreadCount = threadCount;
PrimeOrderCalculator.PrimeOrderCalculatorConfig.ConfigureHeuristicDefault(threadCount); PrimeOrderCalculatorConfig.ConfigureHeuristicDefault(threadCount);
// On EvenPerfectBitScanner's execution path the scheduler instance is created after configuration, // On EvenPerfectBitScanner's execution path the scheduler instance is created after configuration,
// so updating an existing pool would never run. Leave the update path disabled for now. // so updating an existing pool would never run. Leave the update path disabled for now.