From 86f20ed083dd5dd810ca39a965108c7268bfa525 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Fri, 5 Dec 2025 23:33:21 +0100 Subject: [PATCH] Support for q > 2^64 added in --mersenne=bydivisor mode --- .../CalculationResultHandler.cs | 52 +- EvenPerfectBitScanner/CliArguments.cs | 30 +- EvenPerfectBitScanner/Program.cs | 6 + ...rsenneNumberDivisorByDivisorTesterTests.cs | 19 +- .../MersenneNumberDivisorCpuTesterTests.cs | 19 +- .../MersenneNumberDivisorGpuTesterTests.cs | 31 +- ...MersenneNumberDivisorByDivisorCpuTester.cs | 499 +++++++++++++----- ...MersenneNumberDivisorByDivisorGpuTester.cs | 33 +- .../IMersenneNumberDivisorByDivisorTester.cs | 7 +- .../MersenneNumberDivisorByDivisorTester.cs | 67 +-- PerfectNumbers.Core/UInt128Extensions.cs | 1 + 11 files changed, 542 insertions(+), 222 deletions(-) diff --git a/EvenPerfectBitScanner/CalculationResultHandler.cs b/EvenPerfectBitScanner/CalculationResultHandler.cs index aacdf54..3fd8515 100644 --- a/EvenPerfectBitScanner/CalculationResultHandler.cs +++ b/EvenPerfectBitScanner/CalculationResultHandler.cs @@ -1,4 +1,7 @@ +using System; using System.Buffers.Text; +using System.Globalization; +using System.Numerics; using System.Runtime.CompilerServices; using System.Text; using PerfectNumbers.Core; @@ -93,19 +96,46 @@ internal static class CalculationResultHandler } [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static void HandleResult(ulong currentP, BigInteger divisor, bool searchedMersenne, bool detailedCheck, bool passedAllTests, bool lastWasComposite) + { + Span divisorBuffer = stackalloc char[256]; + if (!divisor.TryFormat(divisorBuffer, out int divisorChars, provider: CultureInfo.InvariantCulture)) + { + string divisorString = divisor.ToString(CultureInfo.InvariantCulture); + HandleResultInternal(currentP, divisorString.AsSpan(), searchedMersenne, detailedCheck, passedAllTests, lastWasComposite); + return; + } + + HandleResultInternal(currentP, divisorBuffer[..divisorChars], searchedMersenne, detailedCheck, passedAllTests, lastWasComposite); + } + internal static void HandleResult(ulong currentP, ulong divisor, bool searchedMersenne, bool detailedCheck, bool passedAllTests, bool lastWasComposite) + { + Span divisorBuffer = stackalloc char[32]; + if (!divisor.TryFormat(divisorBuffer, out int divisorChars, provider: CultureInfo.InvariantCulture)) + { + string divisorString = divisor.ToString(CultureInfo.InvariantCulture); + HandleResultInternal(currentP, divisorString.AsSpan(), searchedMersenne, detailedCheck, passedAllTests, lastWasComposite); + return; + } + + HandleResultInternal(currentP, divisorBuffer[..divisorChars], searchedMersenne, detailedCheck, passedAllTests, lastWasComposite); + } + + private static void HandleResultInternal(ulong currentP, ReadOnlySpan divisorText, bool searchedMersenne, bool detailedCheck, bool passedAllTests, bool lastWasComposite) { // Program flow ensures Initialize(...) and InitializeOutputBuffer() ran before this call. Do not add initialization guards here. - Span stackBuffer = stackalloc byte[128]; + Span stackBuffer = stackalloc byte[256]; byte[]? rentedBuffer = null; Span utf8Span = stackBuffer; - int byteCount = TryFormatResultLine(utf8Span, currentP, divisor, searchedMersenne, detailedCheck, passedAllTests); + int byteCount = TryFormatResultLine(utf8Span, currentP, divisorText, searchedMersenne, detailedCheck, passedAllTests); if (byteCount < 0) { - rentedBuffer = ThreadStaticPools.BytePool.Rent(256); + int minimumSize = Math.Max((divisorText.Length + 5) * 4, 256); + rentedBuffer = ThreadStaticPools.BytePool.Rent(minimumSize); utf8Span = rentedBuffer; - byteCount = TryFormatResultLine(utf8Span, currentP, divisor, searchedMersenne, detailedCheck, passedAllTests); + byteCount = TryFormatResultLine(utf8Span, currentP, divisorText, searchedMersenne, detailedCheck, passedAllTests); if (byteCount < 0) { ThreadStaticPools.BytePool.Return(rentedBuffer); @@ -114,7 +144,7 @@ internal static class CalculationResultHandler } ReadOnlySpan utf8Line = utf8Span[..byteCount]; - Span charBuffer = stackalloc char[128]; + Span charBuffer = stackalloc char[256]; int charCount = Utf8Encoding.GetChars(utf8Line, charBuffer); ReadOnlySpan recordSpan = charBuffer[..charCount]; @@ -229,7 +259,7 @@ internal static class CalculationResultHandler } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static int TryFormatResultLine(Span destination, ulong currentP, ulong divisor, bool searchedMersenne, bool detailedCheck, bool passedAllTests) + private static int TryFormatResultLine(Span destination, ulong currentP, ReadOnlySpan divisorText, bool searchedMersenne, bool detailedCheck, bool passedAllTests) { int offset = 0; @@ -246,16 +276,14 @@ internal static class CalculationResultHandler destination[offset++] = (byte)','; - if (!Utf8Formatter.TryFormat(divisor, destination[offset..], out written)) + int divisorByteCount = Utf8Encoding.GetByteCount(divisorText); + if (offset + divisorByteCount >= destination.Length) { return -1; } - offset += written; - if (offset >= destination.Length) - { - return -1; - } + divisorByteCount = Utf8Encoding.GetBytes(divisorText, destination[offset..]); + offset += divisorByteCount; destination[offset++] = (byte)','; diff --git a/EvenPerfectBitScanner/CliArguments.cs b/EvenPerfectBitScanner/CliArguments.cs index 9e8bb85..0bc37bb 100644 --- a/EvenPerfectBitScanner/CliArguments.cs +++ b/EvenPerfectBitScanner/CliArguments.cs @@ -1,4 +1,5 @@ using System.Globalization; +using System.Numerics; using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; @@ -50,7 +51,7 @@ internal readonly struct CliArguments internal readonly double ZeroFractionHard; internal readonly double ZeroFractionConjecture; internal readonly int MaxZeroConjecture; - internal readonly ulong MinK; + internal readonly BigInteger MinK; internal readonly int GpuRatio; private CliArguments( @@ -96,7 +97,7 @@ internal readonly struct CliArguments double zeroFractionHard, double zeroFractionConjecture, int maxZeroConjecture, - ulong minK, + BigInteger minK, int gpuRatio) { StartPrime = startPrime; @@ -191,7 +192,7 @@ internal readonly struct CliArguments double zeroFractionHard = -1.0; double zeroFractionConjecture = -1.0; int maxZeroConjecture = -1; - ulong minK = 1UL; + BigInteger minK = BigInteger.One; int gpuRatio = PerfectNumberConstants.GpuRatio; foreach (string argument in args) @@ -314,14 +315,31 @@ internal readonly struct CliArguments if (argument.StartsWith("--min-k=", StringComparison.OrdinalIgnoreCase)) { - if (Utf8CliParser.TryParseUInt64(argument.AsSpan("--min-k=".Length), out ulong parsedMinK) && parsedMinK > 0UL) + ReadOnlySpan minKSpan = argument.AsSpan("--min-k=".Length); + if (Utf8CliParser.TryParseUInt64(minKSpan, out ulong parsedMinK) && parsedMinK > 0UL) { minK = parsedMinK; } else { - Console.WriteLine("Invalid value for --min-k."); - return default; + try + { + BigInteger parsedMinKBig = BigInteger.Parse(minKSpan, CultureInfo.InvariantCulture); + if (parsedMinKBig > BigInteger.Zero) + { + minK = parsedMinKBig; + } + else + { + Console.WriteLine("Invalid value for --min-k."); + return default; + } + } + catch + { + Console.WriteLine("Invalid value for --min-k."); + return default; + } } continue; diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 6f8ec51..55adfb2 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -1,6 +1,7 @@ using System.Buffers; using System.Diagnostics; using System.Globalization; +using System.Numerics; using System.Runtime.CompilerServices; using EvenPerfectBitScanner.Candidates; using EvenPerfectBitScanner.IO; @@ -241,6 +242,11 @@ internal static class Program : new MersenneNumberDivisorByDivisorCpuTester(orderDevice); _byDivisorTester.MinK = _cliArguments.MinK; _byDivisorTester.BatchSize = scanBatchSize; + if (_byDivisorTester is MersenneNumberDivisorByDivisorGpuTester && _cliArguments.MinK > ulong.MaxValue) + { + Console.WriteLine("--min-k values above UInt64 are only supported on the CPU path for --mersenne=bydivisor."); + return; + } } // Load RLE blacklist (optional) diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs index 3f28334..9664adf 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs @@ -1,3 +1,4 @@ +using System.Numerics; using FluentAssertions; using PerfectNumbers.Core.Gpu.Accelerators; using Xunit; @@ -12,7 +13,7 @@ public class MersenneNumberDivisorByDivisorTesterTests { var candidates = new List { 5UL, 7UL, 11UL, 13UL }; var tester = new FakeTester(); - var results = new List<(ulong Prime, bool Searched, bool Detailed, bool Passed, ulong Divisor)>(); + var results = new List<(ulong Prime, bool Searched, bool Detailed, bool Passed, BigInteger Divisor)>(); int compositeMarks = 0; int clearedMarks = 0; @@ -31,12 +32,12 @@ public class MersenneNumberDivisorByDivisorTesterTests compositeMarks.Should().Be(2); clearedMarks.Should().Be(2); results.Should().Equal( - new[] + new (ulong Prime, bool Searched, bool Detailed, bool Passed, BigInteger Divisor)[] { - (5UL, true, true, true, 0UL), - (7UL, true, true, false, 7UL), - (11UL, true, true, false, 11UL), - (13UL, true, true, true, 0UL), + (5UL, true, true, true, BigInteger.Zero), + (7UL, true, true, false, new BigInteger(7)), + (11UL, true, true, false, new BigInteger(11)), + (13UL, true, true, true, BigInteger.Zero), }); } @@ -46,7 +47,7 @@ public class MersenneNumberDivisorByDivisorTesterTests public int BatchSize { get; set; } - public ulong MinK { get; set; } + public BigInteger MinK { get; set; } public string? StateFilePath { get; set; } @@ -54,7 +55,7 @@ public class MersenneNumberDivisorByDivisorTesterTests { } - public void ResumeFromState(ulong lastSavedK) + public void ResumeFromState(BigInteger lastSavedK) { } @@ -70,7 +71,7 @@ public class MersenneNumberDivisorByDivisorTesterTests return DivisorLimit; } - public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out ulong divisor) + public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out BigInteger divisor) { divisorsExhausted = true; if (prime % 7UL == 0UL) diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs index dddebee..b84444b 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs @@ -1,3 +1,4 @@ +using System.Numerics; using FluentAssertions; using PerfectNumbers.Core.Cpu; using PerfectNumbers.Core.Gpu.Accelerators; @@ -17,17 +18,17 @@ public class MersenneNumberDivisorCpuTesterTests var gpu = PrimeOrderCalculatorAccelerator.Rent(1); try { - tester.IsPrime(gpu, 5UL, out bool divisorsExhausted, out ulong divisor).Should().BeTrue(); - divisorsExhausted.Should().BeTrue(); - divisor.Should().Be(0UL); + tester.IsPrime(gpu, 5UL, out bool divisorsExhausted, out BigInteger divisor).Should().BeTrue(); + divisorsExhausted.Should().BeTrue(); + divisor.Should().Be(BigInteger.Zero); - tester.IsPrime(gpu, 7UL, out divisorsExhausted, out divisor).Should().BeTrue(); - divisorsExhausted.Should().BeTrue(); - divisor.Should().Be(0UL); + tester.IsPrime(gpu, 7UL, out divisorsExhausted, out divisor).Should().BeTrue(); + divisorsExhausted.Should().BeTrue(); + divisor.Should().Be(BigInteger.Zero); - tester.IsPrime(gpu, 11UL, out divisorsExhausted, out divisor).Should().BeFalse(); - divisorsExhausted.Should().BeTrue(); - divisor.Should().BeGreaterThan(0UL); + tester.IsPrime(gpu, 11UL, out divisorsExhausted, out divisor).Should().BeFalse(); + divisorsExhausted.Should().BeTrue(); + divisor.Should().BeGreaterThan(BigInteger.Zero); } finally { diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index 2248e5f..fa887ab 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -1,3 +1,4 @@ +using System.Numerics; using System.Reflection; using FluentAssertions; using PerfectNumbers.Core.Cpu; @@ -102,17 +103,17 @@ public class MersenneNumberDivisorGpuTesterTests var gpu = PrimeOrderCalculatorAccelerator.Rent(1); try { - tester.IsPrime(gpu, 31UL, out bool exhausted, out ulong divisor).Should().BeTrue(); + tester.IsPrime(gpu, 31UL, out bool exhausted, out BigInteger divisor).Should().BeTrue(); exhausted.Should().BeTrue(); - divisor.Should().Be(0UL); + divisor.Should().Be(BigInteger.Zero); tester.IsPrime(gpu, 37UL, out exhausted, out divisor).Should().BeFalse(); exhausted.Should().BeTrue(); - divisor.Should().BeGreaterThan(0UL); + divisor.Should().BeGreaterThan(BigInteger.Zero); tester.IsPrime(gpu, 43UL, out exhausted, out divisor).Should().BeFalse(); exhausted.Should().BeTrue(); - divisor.Should().BeGreaterThan(0UL); + divisor.Should().BeGreaterThan(BigInteger.Zero); } finally { @@ -145,15 +146,15 @@ public class MersenneNumberDivisorGpuTesterTests var gpu = PrimeOrderCalculatorAccelerator.Rent(1); try { - tester.IsPrime(gpu, 31UL, out _, out ulong divisor).Should().BeTrue(); - divisor.Should().Be(0UL); + tester.IsPrime(gpu, 31UL, out _, out BigInteger divisor).Should().BeTrue(); + divisor.Should().Be(BigInteger.Zero); tester.IsPrime(gpu, 37UL, out bool exhausted, out divisor).Should().BeFalse(); exhausted.Should().BeTrue(); - divisor.Should().BeGreaterThan(0UL); + divisor.Should().BeGreaterThan(BigInteger.Zero); tester.IsPrime(gpu, 41UL, out bool divisorsExhausted, out divisor).Should().BeFalse(); divisorsExhausted.Should().BeTrue(); - divisor.Should().BeGreaterThan(0UL); + divisor.Should().BeGreaterThan(BigInteger.Zero); } finally { @@ -175,16 +176,16 @@ public class MersenneNumberDivisorGpuTesterTests var gpu = PrimeOrderCalculatorAccelerator.Rent(1); try { - tester.IsPrime(gpu, exponent, out bool exhausted, out ulong divisor).Should().Be(expectedPrime); + tester.IsPrime(gpu, exponent, out bool exhausted, out BigInteger divisor).Should().Be(expectedPrime); exhausted.Should().BeTrue(); if (expectedPrime) { - divisor.Should().Be(0UL); + divisor.Should().Be(BigInteger.Zero); } else { - divisor.Should().BeGreaterThan(0UL); + divisor.Should().BeGreaterThan(BigInteger.Zero); } } finally @@ -207,9 +208,9 @@ public class MersenneNumberDivisorGpuTesterTests var gpu = PrimeOrderCalculatorAccelerator.Rent(1); try { - tester.IsPrime(gpu, 41UL, out bool divisorsExhausted, out ulong divisor).Should().BeFalse(); + tester.IsPrime(gpu, 41UL, out bool divisorsExhausted, out BigInteger divisor).Should().BeFalse(); divisorsExhausted.Should().BeTrue(); - divisor.Should().BeGreaterThan(0UL); + divisor.Should().BeGreaterThan(BigInteger.Zero); } finally { @@ -397,9 +398,9 @@ public class MersenneNumberDivisorGpuTesterTests .GetField("_divisorLimit", BindingFlags.NonPublic | BindingFlags.Instance)! .SetValue(tester, 200UL); - tester.IsPrime(gpu, 19UL, out bool divisorsExhausted, out ulong divisor).Should().BeTrue(); + tester.IsPrime(gpu, 19UL, out bool divisorsExhausted, out BigInteger divisor).Should().BeTrue(); divisorsExhausted.Should().BeTrue(); - divisor.Should().Be(0UL); + divisor.Should().Be(BigInteger.Zero); } finally { diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index b693b10..f34b688 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,3 +1,4 @@ +using System.Numerics; using System.Runtime.CompilerServices; using System.Globalization; using System.IO; @@ -8,15 +9,15 @@ namespace PerfectNumbers.Core.Cpu; public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester { - private ulong _divisorLimit; - private int _batchSize = 1_024; - private ulong _minK = 1UL; - private readonly ComputationDevice _orderDevice; - private readonly TryCycleLengthDelegate _tryCalculateCycleLengthForExponent; - private readonly CalculateCycleLengthDelegate _calculateCycleLength; - private string? _stateFilePath; - private int _stateCounter; - private ulong _lastSavedK; + private BigInteger _divisorLimit; + private int _batchSize = 1_024; + private BigInteger _minK = BigInteger.One; + private readonly ComputationDevice _orderDevice; + private readonly TryCycleLengthDelegate _tryCalculateCycleLengthForExponent; + private readonly CalculateCycleLengthDelegate _calculateCycleLength; + private string? _stateFilePath; + private int _stateCounter; + private BigInteger _lastSavedK; public ComputationDevice OrderDevice => _orderDevice; private delegate bool TryCycleLengthDelegate( @@ -67,11 +68,11 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv set => _batchSize = Math.Max(1, value); } - public ulong MinK - { - get => _minK; - set => _minK = value < 1UL ? 1UL : value; - } + public BigInteger MinK + { + get => _minK; + set => _minK = value < BigInteger.One ? BigInteger.One : value; + } public string? StateFilePath { @@ -84,37 +85,42 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv _stateCounter = 0; } - public void ResumeFromState(ulong lastSavedK) - { - _lastSavedK = lastSavedK; - _minK = lastSavedK + 1UL; - _stateCounter = 0; - } + public void ResumeFromState(BigInteger lastSavedK) + { + _lastSavedK = lastSavedK; + _minK = lastSavedK + BigInteger.One; + _stateCounter = 0; + } - public void ConfigureFromMaxPrime(ulong maxPrime) - { - _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); - } + public void ConfigureFromMaxPrime(ulong maxPrime) + { + _divisorLimit = ComputeDivisorLimitFromMaxPrimeBig(maxPrime); + } - public ulong DivisorLimit - { - [MethodImpl(MethodImplOptions.AggressiveInlining)] - get - { - return _divisorLimit; - } - } + public ulong DivisorLimit + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + if (_divisorLimit.IsZero) + { + return ulong.MaxValue; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public ulong GetAllowedMaxDivisor(ulong prime) => ComputeAllowedMaxDivisor(prime, _divisorLimit); + return _divisorLimit > ulong.MaxValue ? ulong.MaxValue : (ulong)_divisorLimit; + } + } - public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out ulong divisor) - { - ulong allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ulong GetAllowedMaxDivisor(ulong prime) => ComputeAllowedMaxDivisor(prime, _divisorLimit); - // The CPU by-divisor run always hands us primes with enormous divisor limits, so the fallback below never executes. - // if (allowedMax < 3UL) - // { + public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out BigInteger divisor) + { + BigInteger allowedMax = ComputeAllowedMaxDivisorBig(prime, _divisorLimit); + + // The CPU by-divisor run always hands us primes with enormous divisor limits, so the fallback below never executes. + // if (allowedMax < 3UL) + // { // // EvenPerfectBitScanner routes primes below the small-divisor cutoff to the GPU path, so the CPU path still sees // // trivial candidates during targeted tests. Short-circuit here to keep those runs aligned with the production flow. // divisorsExhausted = true; @@ -123,34 +129,34 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv bool processedAll; - bool composite = CheckDivisors( - gpu, - prime, - allowedMax, - _minK, - out processedAll, + bool composite = CheckDivisors( + gpu, + prime, + allowedMax, + _minK, + out processedAll, out divisor); if (composite) - { - divisorsExhausted = true; - return false; - } + { + divisorsExhausted = true; + return false; + } - divisorsExhausted = processedAll || composite; - divisor = 0UL; - return true; - } + divisorsExhausted = processedAll || composite; + divisor = BigInteger.Zero; + return true; + } - public void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues) - { - ulong divisorLimit = _divisorLimit; - int length = primes.Length; - for (int index = 0; index < length; index++) - { - allowedMaxValues[index] = ComputeAllowedMaxDivisor(primes[index], divisorLimit); - } - } + public void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues) + { + ulong divisorLimit = GetAllowedDivisorLimitForSpan(_divisorLimit); + int length = primes.Length; + for (int index = 0; index < length; index++) + { + allowedMaxValues[index] = ComputeAllowedMaxDivisor(primes[index], divisorLimit); + } + } public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu) { @@ -165,17 +171,84 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return new MersenneCpuDivisorScanSession(gpu, _orderDevice); } - private bool CheckDivisors( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong allowedMax, - ulong minK, - out bool processedAll, - out ulong foundDivisor) - { - foundDivisor = 0UL; - processedAll = true; - ulong currentK = minK < 1UL ? 1UL : minK; + private bool CheckDivisors( + PrimeOrderCalculatorAccelerator gpu, + ulong prime, + BigInteger allowedMax, + BigInteger minK, + out bool processedAll, + out BigInteger foundDivisor) + { + BigInteger normalizedMinK = minK < BigInteger.One ? BigInteger.One : minK; + BigInteger step = ((BigInteger)prime) << 1; + BigInteger firstDivisor = (step * normalizedMinK) + BigInteger.One; + bool unlimited = allowedMax.IsZero; + + bool fits64 = normalizedMinK <= ulong.MaxValue && firstDivisor <= ulong.MaxValue; + bool limitFits64 = unlimited || allowedMax <= ulong.MaxValue; + if (fits64 && limitFits64) + { + ulong allowedMax64 = unlimited ? ulong.MaxValue : (ulong)allowedMax; + bool composite64 = CheckDivisors64Bit( + gpu, + prime, + allowedMax64, + (ulong)normalizedMinK, + out processedAll, + out ulong foundDivisor64); + foundDivisor = foundDivisor64; + return composite64; + } + + if (fits64 && allowedMax > ulong.MaxValue) + { + ulong allowedMax64 = ulong.MaxValue; + bool composite64 = CheckDivisors64Bit( + gpu, + prime, + allowedMax64, + (ulong)normalizedMinK, + out processedAll, + out ulong foundDivisor64); + if (composite64) + { + foundDivisor = foundDivisor64; + return true; + } + + BigInteger iterations = ((BigInteger)(allowedMax64 - (ulong)firstDivisor) / (ulong)step) + BigInteger.One; + BigInteger nextK = normalizedMinK + iterations; + return CheckDivisorsLarge( + gpu, + prime, + allowedMax, + nextK, + step, + out processedAll, + out foundDivisor); + } + + return CheckDivisorsLarge( + gpu, + prime, + allowedMax, + normalizedMinK, + step, + out processedAll, + out foundDivisor); + } + + private bool CheckDivisors64Bit( + PrimeOrderCalculatorAccelerator gpu, + ulong prime, + ulong allowedMax, + ulong minK, + out bool processedAll, + out ulong foundDivisor) + { + foundDivisor = 0UL; + processedAll = true; + ulong currentK = minK < 1UL ? 1UL : minK; // The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs. // Keeping the guard commented out documents the reasoning for benchmarks and tests. @@ -405,16 +478,105 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv remainder19 = AddMod19(remainder19, step19); } - processedAll = true; - foundDivisor = 0UL; - return false; - } + processedAll = true; + foundDivisor = 0UL; + return false; + } - private bool CheckDivisors64Range( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong step, - ulong allowedMax, + private bool CheckDivisorsLarge( + PrimeOrderCalculatorAccelerator gpu, + ulong prime, + BigInteger allowedMax, + BigInteger minK, + BigInteger step, + out bool processedAll, + out BigInteger foundDivisor) + { + processedAll = true; + foundDivisor = BigInteger.Zero; + BigInteger currentK = minK < BigInteger.One ? BigInteger.One : minK; + BigInteger divisor = (step * currentK) + BigInteger.One; + _ = gpu; + if (allowedMax.IsZero) + { + processedAll = false; + return false; + } + if (divisor > allowedMax) + { + return false; + } + + LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; + ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); + + byte step10 = (byte)(step % 10); + byte step8 = (byte)(step % 8); + byte step3 = (byte)(step % 3); + byte step7 = (byte)(step % 7); + byte step11 = (byte)(step % 11); + byte step13 = (byte)(step % 13); + byte step17 = (byte)(step % 17); + byte step19 = (byte)(step % 19); + + byte remainder10 = (byte)(divisor % 10); + byte remainder8 = (byte)(divisor % 8); + byte remainder3 = (byte)(divisor % 3); + byte remainder7 = (byte)(divisor % 7); + byte remainder11 = (byte)(divisor % 11); + byte remainder13 = (byte)(divisor % 13); + byte remainder17 = (byte)(divisor % 17); + byte remainder19 = (byte)(divisor % 19); + + while (divisor <= allowedMax) + { + bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0; + if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) + { + if (IsProbablePrimeBigInteger(divisor)) + { + BigInteger powResult = BigInteger.ModPow(2, prime, divisor); + RecordState(currentK); + if (powResult.IsOne) + { + foundDivisor = divisor; + processedAll = true; + return true; + } + } + else + { + RecordState(currentK); + } + } + + currentK += BigInteger.One; + divisor += step; + remainder10 = AddMod10(remainder10, step10); + remainder8 = AddMod8(remainder8, step8); + remainder3 = AddMod3(remainder3, step3); + remainder7 = AddMod7(remainder7, step7); + remainder11 = AddMod11(remainder11, step11); + remainder13 = AddMod13(remainder13, step13); + remainder17 = AddMod17(remainder17, step17); + remainder19 = AddMod19(remainder19, step19); + + if (divisor > allowedMax) + { + break; + } + } + + processedAll = true; + foundDivisor = BigInteger.Zero; + return false; + } + + private bool CheckDivisors64Range( + PrimeOrderCalculatorAccelerator gpu, + ulong prime, + ulong step, + ulong allowedMax, ulong startK, ulong endK, ushort decimalMask, @@ -647,18 +809,18 @@ private bool CheckDivisors64( return residue == 1UL ? (byte)1 : (byte)0; } - private void RecordState(ulong k) - { - string? path = _stateFilePath; - if (string.IsNullOrEmpty(path)) - { - return; + private void RecordState(BigInteger k) + { + string? path = _stateFilePath; + if (string.IsNullOrEmpty(path)) + { + return; } - if (k <= _lastSavedK) - { - return; - } + if (k <= _lastSavedK) + { + return; + } int next = _stateCounter + 1; if (next >= PerfectNumberConstants.ByDivisorStateSaveInterval) @@ -674,15 +836,70 @@ private bool CheckDivisors64( _lastSavedK = k; } else - { - _stateCounter = next; - } - } + { + _stateCounter = next; + } + } - private static bool TryCalculateCycleLengthForExponentCpu( - ulong divisor, - PrimeOrderCalculatorAccelerator gpu, - ulong exponent, + private static bool IsProbablePrimeBigInteger(BigInteger value) + { + if (value <= 3) + { + return value >= 2; + } + + if (value.IsEven) + { + return false; + } + + BigInteger d = value - 1; + int s = 0; + while ((d & 1) == 0) + { + d >>= 1; + s++; + } + + ReadOnlySpan bases = stackalloc int[] { 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 ((BigInteger)baseValue >= value) + { + continue; + } + + BigInteger x = BigInteger.ModPow(baseValue, d, value); + if (x == 1 || x == value - 1) + { + continue; + } + + bool witnessFound = true; + for (int r = 1; r < s; r++) + { + x = BigInteger.ModPow(x, 2, value); + if (x == value - 1) + { + witnessFound = false; + break; + } + } + + if (witnessFound) + { + return false; + } + } + + return true; + } + + private static bool TryCalculateCycleLengthForExponentCpu( + ulong divisor, + PrimeOrderCalculatorAccelerator gpu, + ulong exponent, in MontgomeryDivisorData divisorData, out ulong cycleLength, out bool primeOrderFailed) => MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( @@ -854,35 +1071,63 @@ private bool CheckDivisors64( return (byte)((value + delta) & 7); } - private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime) - { - // The by-divisor CPU configuration only supplies primes greater than 1, so the guard below never trips. - // if (maxPrime <= 1UL) - // { - // return 0UL; - // } - if (maxPrime - 1UL >= 64UL) - { - return ulong.MaxValue; - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong GetAllowedDivisorLimitForSpan(BigInteger divisorLimit) + { + if (divisorLimit.IsZero) + { + return ulong.MaxValue; + } - return (1UL << (int)(maxPrime - 1UL)) - 1UL; - } + return divisorLimit > ulong.MaxValue ? ulong.MaxValue : (ulong)divisorLimit; + } - private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit) - { - // Production by-divisor scans only handle primes, so inputs never fall below 2. - // if (prime <= 1UL) - // { - // return 0UL; - // } - if (prime - 1UL >= 64UL) - { - return divisorLimit; - } + private static BigInteger ComputeDivisorLimitFromMaxPrimeBig(ulong maxPrime) + { + if (maxPrime <= 1UL) + { + return BigInteger.Zero; + } - return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); - } + if (maxPrime - 1UL >= 64UL) + { + return (BigInteger.One << 256) - BigInteger.One; + } + + return (BigInteger.One << (int)(maxPrime - 1UL)) - BigInteger.One; + } + + private static BigInteger ComputeAllowedMaxDivisorBig(ulong prime, BigInteger divisorLimit) + { + if (prime <= 1UL) + { + return BigInteger.Zero; + } + + if (divisorLimit.IsZero || prime - 1UL >= 64UL) + { + return divisorLimit; + } + + BigInteger candidateLimit = (BigInteger.One << (int)(prime - 1UL)) - BigInteger.One; + if (divisorLimit.IsZero) + { + return candidateLimit; + } + + return BigInteger.Min(candidateLimit, divisorLimit); + } + + private static ulong ComputeAllowedMaxDivisor(ulong prime, BigInteger divisorLimit) + { + BigInteger allowed = ComputeAllowedMaxDivisorBig(prime, divisorLimit); + if (allowed.IsZero) + { + return ulong.MaxValue; + } + + return allowed > ulong.MaxValue ? ulong.MaxValue : (ulong)allowed; + } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index acf086b..458993c 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -1,4 +1,6 @@ +using System; using System.Collections.Concurrent; +using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using ILGPU; @@ -18,7 +20,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN // so the synchronization fields from the previous implementation remain commented out here. // private readonly object _sync = new(); private ulong _divisorLimit; - private ulong _minK = 1UL; + private BigInteger _minK = BigInteger.One; // private bool _isConfigured; private readonly ConcurrentBag[] _resourcePool = [.. AcceleratorPool.Shared.Accelerators.Select(x => new ConcurrentBag())]; @@ -37,10 +39,10 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN set => GpuBatchSize = value; } - public ulong MinK + public BigInteger MinK { get => _minK; - set => _minK = value < 1UL ? 1UL : value; + set => _minK = value < BigInteger.One ? BigInteger.One : value; } public string? StateFilePath { get; set; } @@ -49,9 +51,14 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN { } - public void ResumeFromState(ulong lastSavedK) + public void ResumeFromState(BigInteger lastSavedK) { - _minK = lastSavedK + 1UL; + if (lastSavedK > ulong.MaxValue) + { + throw new NotSupportedException("GPU by-divisor tester does not support MinK above UInt64."); + } + + _minK = lastSavedK + BigInteger.One; } public void ConfigureFromMaxPrime(ulong maxPrime) @@ -68,7 +75,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN // _isConfigured = true; } - public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out ulong divisor) + public bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out BigInteger divisor) { ulong allowedMax; int batchCapacity; @@ -113,7 +120,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN gpu, prime, allowedMax, - _minK, + GetMinKOrThrow(), resources.CheckDivisorKernel, resources.DivisorDataBuffer, resources.OffsetBuffer, @@ -144,10 +151,20 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN } divisorsExhausted = coveredRange; - divisor = 0UL; + divisor = BigInteger.Zero; return true; } + private ulong GetMinKOrThrow() + { + if (_minK > ulong.MaxValue) + { + throw new NotSupportedException("GPU by-divisor tester does not support MinK above UInt64."); + } + + return (ulong)_minK; + } + public void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues) { // The scanner always supplies matching spans, so the previous validation remains commented out. diff --git a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs index 942864b..03071d5 100644 --- a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs @@ -1,3 +1,4 @@ +using System.Numerics; using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core; @@ -6,7 +7,7 @@ public interface IMersenneNumberDivisorByDivisorTester { int BatchSize { get; set; } - ulong MinK { get; set; } + BigInteger MinK { get; set; } void ConfigureFromMaxPrime(ulong maxPrime); @@ -14,7 +15,7 @@ public interface IMersenneNumberDivisorByDivisorTester ulong GetAllowedMaxDivisor(ulong prime); - bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out ulong divisor); + bool IsPrime(PrimeOrderCalculatorAccelerator gpu, ulong prime, out bool divisorsExhausted, out BigInteger divisor); IDivisorScanSession CreateDivisorSession(PrimeOrderCalculatorAccelerator gpu); @@ -22,7 +23,7 @@ public interface IMersenneNumberDivisorByDivisorTester void ResetStateTracking(); - void ResumeFromState(ulong lastSavedK); + void ResumeFromState(BigInteger lastSavedK); void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues); diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index 0961d9a..710bb39 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -1,5 +1,6 @@ using System.Buffers; using System.Globalization; +using System.Numerics; using System.Runtime.InteropServices; using PerfectNumbers.Core.Gpu; using PerfectNumbers.Core.Gpu.Accelerators; @@ -13,12 +14,12 @@ public static class MersenneNumberDivisorByDivisorTester List candidates, IMersenneNumberDivisorByDivisorTester tester, Dictionary? previousResults, - ulong startPrime, - Action markComposite, - Action clearComposite, - Action printResult, - int threadCount, - int primesPerTask) + ulong startPrime, + Action markComposite, + Action clearComposite, + Action printResult, + int threadCount, + int primesPerTask) { if (candidates.Count == 0) { @@ -243,15 +244,15 @@ public static class MersenneNumberDivisorByDivisorTester IMersenneNumberDivisorByDivisorTester CreateTester(IMersenneNumberDivisorByDivisorTester prototype) { - IMersenneNumberDivisorByDivisorTester tester = prototype switch - { - MersenneNumberDivisorByDivisorCpuTester cpu => new MersenneNumberDivisorByDivisorCpuTester(cpu.OrderDevice) - { - BatchSize = prototype.BatchSize, - MinK = prototype.MinK, - }, - _ => CreateGpuTester(prototype), - }; + IMersenneNumberDivisorByDivisorTester tester = prototype switch + { + MersenneNumberDivisorByDivisorCpuTester cpu => new MersenneNumberDivisorByDivisorCpuTester(cpu.OrderDevice) + { + BatchSize = prototype.BatchSize, + MinK = prototype.MinK, + }, + _ => CreateGpuTester(prototype), + }; tester.ConfigureFromMaxPrime(maxPrime); return tester; @@ -275,27 +276,27 @@ public static class MersenneNumberDivisorByDivisorTester IMersenneNumberDivisorByDivisorTester testerForPrime = CreateTester(tester); string stateFile = Path.Combine(PerfectNumberConstants.ByDivisorStateDirectory, prime.ToString(CultureInfo.InvariantCulture) + ".bin"); - ulong resumeK = tester.MinK; - if (File.Exists(stateFile)) - { - ulong lastK = File.ReadLines(stateFile) - .Where(static line => !string.IsNullOrWhiteSpace(line)) - .Select(static line => ulong.Parse(line, NumberStyles.None, CultureInfo.InvariantCulture)) - .Max(); + BigInteger resumeK = tester.MinK; + if (File.Exists(stateFile)) + { + BigInteger lastK = File.ReadLines(stateFile) + .Where(static line => !string.IsNullOrWhiteSpace(line)) + .Select(static line => BigInteger.Parse(line, NumberStyles.None, CultureInfo.InvariantCulture)) + .Max(); - resumeK = lastK + 1UL; - testerForPrime.ResumeFromState(lastK); - } - else - { - testerForPrime.ResumeFromState(0UL); - } + resumeK = lastK + BigInteger.One; + testerForPrime.ResumeFromState(lastK); + } + else + { + testerForPrime.ResumeFromState(BigInteger.Zero); + } - testerForPrime.MinK = resumeK; - testerForPrime.StateFilePath = stateFile; - testerForPrime.ResetStateTracking(); + testerForPrime.MinK = resumeK; + testerForPrime.StateFilePath = stateFile; + testerForPrime.ResetStateTracking(); - bool isPrime = testerForPrime.IsPrime(gpu, prime, out bool divisorsExhausted, out ulong divisor); + bool isPrime = testerForPrime.IsPrime(gpu, prime, out bool divisorsExhausted, out BigInteger divisor); if (!isPrime) { diff --git a/PerfectNumbers.Core/UInt128Extensions.cs b/PerfectNumbers.Core/UInt128Extensions.cs index 5152383..521c462 100644 --- a/PerfectNumbers.Core/UInt128Extensions.cs +++ b/PerfectNumbers.Core/UInt128Extensions.cs @@ -537,6 +537,7 @@ public static class UInt128Extensions int i = 0, primesLength = smallPrimes.Length; UInt128 cycle = MersenneDivisorCycles.GetCycleCpu(q); + // TODO: If the cache lacks this cycle, immediately schedule the configured device // (GPU by default) to compute it on the fly and skip inserting it into the cache so // wide-order factoring can still leverage cycle stepping without breaching the