From cd1debe9d69b4b96c512dde8df38af1d6c2d7f13 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Tue, 7 Oct 2025 22:14:34 +0200 Subject: [PATCH] Various micro- and medium- optimizations --- .../Pow2MontgomeryModBenchmarks.cs | 2 +- ...MontgomeryModCycleComputationBenchmarks.cs | 6 +- .../Pow2OddPowerTableBenchmarks.cs | 2482 +++++++++++++++++ EvenPerfectBitScanner.Benchmarks/Program.cs | 16 +- .../Gpu/PrimeOrderGpuHeuristicsTests.cs | 56 +- .../MersenneDivisorCyclesGenerateGpuTests.cs | 6 +- .../MersenneDivisorCyclesTests.cs | 3 +- .../MersenneNumberDivisorCpuTesterTests.cs | 8 +- .../MersenneNumberDivisorGpuTesterTests.cs | 16 +- .../PrimeOrderCalculatorTests.cs | 8 + ...MersenneNumberDivisorByDivisorCpuTester.cs | 7 +- PerfectNumbers.Core/DivisorCycleCache.cs | 6 +- PerfectNumbers.Core/Gpu/GpuKernelPool.cs | 511 ++-- PerfectNumbers.Core/Gpu/GpuUInt128.cs | 802 +----- ...MersenneNumberDivisorByDivisorGpuTester.cs | 5 +- PerfectNumbers.Core/Gpu/Pow2OddPowerTable.cs | 681 +++++ .../Gpu/PrimeOrderGpuHeuristics.cs | 1208 ++------ .../IMersenneNumberDivisorByDivisorTester.cs | 2 +- PerfectNumbers.Core/MersenneDivisorCycles.cs | 961 +++---- PerfectNumbers.Core/MontgomeryDivisorData.cs | 47 +- PerfectNumbers.Core/PerfectNumberConstants.cs | 1 + .../PrimeOrderCalculator.Wide.cs | 71 +- PerfectNumbers.Core/PrimeOrderCalculator.cs | 20 +- PerfectNumbers.Core/UInt128Numbers.cs | 4 + PerfectNumbers.Core/ULongExtensions.cs | 684 +---- 25 files changed, 4403 insertions(+), 3210 deletions(-) create mode 100644 EvenPerfectBitScanner.Benchmarks/Pow2OddPowerTableBenchmarks.cs create mode 100644 PerfectNumbers.Core/Gpu/Pow2OddPowerTable.cs diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs index 2d6731e..48c2050 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs @@ -39,7 +39,7 @@ public class Pow2MontgomeryModBenchmarks ulong smallModulus = NextSmallOddModulus(); _smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus); _smallExponents[i] = NextSmallExponent(); - _smallCycles[i] = MersenneDivisorCycles.CalculateCycleLength(smallModulus); + _smallCycles[i] = MersenneDivisorCycles.CalculateCycleLength(smallModulus, MontgomeryDivisorData.FromModulus(smallModulus)); (ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle(); _largeDivisors[i] = CreateMontgomeryDivisorData(largeModulus); diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs index 6497724..8c773fb 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs @@ -261,9 +261,10 @@ public class Pow2MontgomeryModCycleComputationBenchmarks private ulong CalculateCycleLengthWithHeuristics(ulong modulus, ulong fallbackOrder = 0UL) { + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus); if (modulus <= 1UL || (modulus & 1UL) == 0UL) { - return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus); + return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData); } // When the below is uncommented, the benchmark never completes due to the primality test cost on very large random set. @@ -281,6 +282,7 @@ public class Pow2MontgomeryModCycleComputationBenchmarks PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate( modulus, _previousPrimeOrder, + divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault); if (orderResult.Order != 0UL) @@ -292,7 +294,7 @@ public class Pow2MontgomeryModCycleComputationBenchmarks _previousPrimeOrder = null; Console.WriteLine($"Heuristic failed for {modulus}, falling back to full cycle calculation"); - return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus); + return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData); } private UInt128 CalculateCycleLengthWithHeuristics(UInt128 modulus) diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2OddPowerTableBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2OddPowerTableBenchmarks.cs new file mode 100644 index 0000000..9c9e1bf --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/Pow2OddPowerTableBenchmarks.cs @@ -0,0 +1,2482 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Runtime.CompilerServices; +using System.Threading.Tasks; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks +{ + public readonly struct Pow2OddPowerTableDelegate + { + private static readonly Action _emptyAction = () => { }; + + public Pow2OddPowerTableDelegate(GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount) + { + Element0 = baseValue; + GpuUInt128 current; + Action transform; + if (oddPowerCount > 1) + { + // Reusing baseValue to hold base^2 for the shared odd-power ladder that follows. + baseValue.MulMod(baseValue, modulus); + // TODO: We can calculate baseValue % modulus before loop and use it to increase ladderEntry calculation speed - we'll reuse the base for incremental calculations. + current = baseValue; + transform = () => current.MulMod(baseValue, modulus); + } + else + { + current = GpuUInt128.Zero; + transform = _emptyAction; + } + + // We're manually assigning each field to prevent the compiler to initialize each field twice due to auto-initialization. We're using the action to lower the code base size. + transform(); + Element1 = current; + transform(); + Element2 = current; + transform(); + Element3 = current; + transform(); + Element4 = current; + transform(); + Element5 = current; + transform(); + Element6 = current; + transform(); + Element7 = current; + transform(); + Element8 = current; + transform(); + Element9 = current; + transform(); + Element10 = current; + transform(); + Element11 = current; + transform(); + Element12 = current; + transform(); + Element13 = current; + transform(); + Element14 = current; + transform(); + Element15 = current; + transform(); + Element16 = current; + transform(); + Element17 = current; + transform(); + Element18 = current; + transform(); + Element19 = current; + transform(); + Element20 = current; + transform(); + Element21 = current; + transform(); + Element22 = current; + transform(); + Element23 = current; + transform(); + Element24 = current; + transform(); + Element25 = current; + transform(); + Element26 = current; + transform(); + Element27 = current; + transform(); + Element28 = current; + transform(); + Element29 = current; + transform(); + Element30 = current; + transform(); + Element31 = current; + transform(); + Element32 = current; + transform(); + Element33 = current; + transform(); + Element34 = current; + transform(); + Element35 = current; + transform(); + Element36 = current; + transform(); + Element37 = current; + transform(); + Element38 = current; + transform(); + Element39 = current; + transform(); + Element40 = current; + transform(); + Element41 = current; + transform(); + Element42 = current; + transform(); + Element43 = current; + transform(); + Element44 = current; + transform(); + Element45 = current; + transform(); + Element46 = current; + transform(); + Element47 = current; + transform(); + Element48 = current; + transform(); + Element49 = current; + transform(); + Element50 = current; + transform(); + Element51 = current; + transform(); + Element52 = current; + transform(); + Element53 = current; + transform(); + Element54 = current; + transform(); + Element55 = current; + transform(); + Element56 = current; + transform(); + Element57 = current; + transform(); + Element58 = current; + transform(); + Element59 = current; + transform(); + Element60 = current; + transform(); + Element61 = current; + transform(); + Element62 = current; + transform(); + Element63 = current; + transform(); + Element64 = current; + transform(); + Element65 = current; + transform(); + Element66 = current; + transform(); + Element67 = current; + transform(); + Element68 = current; + transform(); + Element69 = current; + transform(); + Element70 = current; + transform(); + Element71 = current; + transform(); + Element72 = current; + transform(); + Element73 = current; + transform(); + Element74 = current; + transform(); + Element75 = current; + transform(); + Element76 = current; + transform(); + Element77 = current; + transform(); + Element78 = current; + transform(); + Element79 = current; + transform(); + Element80 = current; + transform(); + Element81 = current; + transform(); + Element82 = current; + transform(); + Element83 = current; + transform(); + Element84 = current; + transform(); + Element85 = current; + transform(); + Element86 = current; + transform(); + Element87 = current; + transform(); + Element88 = current; + transform(); + Element89 = current; + transform(); + Element90 = current; + transform(); + Element91 = current; + transform(); + Element92 = current; + transform(); + Element93 = current; + transform(); + Element94 = current; + transform(); + Element95 = current; + transform(); + Element96 = current; + transform(); + Element97 = current; + transform(); + Element98 = current; + transform(); + Element99 = current; + transform(); + Element100 = current; + transform(); + Element101 = current; + transform(); + Element102 = current; + transform(); + Element103 = current; + transform(); + Element104 = current; + transform(); + Element105 = current; + transform(); + Element106 = current; + transform(); + Element107 = current; + transform(); + Element108 = current; + transform(); + Element109 = current; + transform(); + Element110 = current; + transform(); + Element111 = current; + transform(); + Element112 = current; + transform(); + Element113 = current; + transform(); + Element114 = current; + transform(); + Element115 = current; + transform(); + Element116 = current; + transform(); + Element117 = current; + transform(); + Element118 = current; + transform(); + Element119 = current; + transform(); + Element120 = current; + transform(); + Element121 = current; + transform(); + Element122 = current; + transform(); + Element123 = current; + transform(); + Element124 = current; + transform(); + Element125 = current; + transform(); + Element126 = current; + transform(); + Element127 = current; + } + + public readonly GpuUInt128 Element0; + public readonly GpuUInt128 Element1; + public readonly GpuUInt128 Element2; + public readonly GpuUInt128 Element3; + public readonly GpuUInt128 Element4; + public readonly GpuUInt128 Element5; + public readonly GpuUInt128 Element6; + public readonly GpuUInt128 Element7; + public readonly GpuUInt128 Element8; + public readonly GpuUInt128 Element9; + public readonly GpuUInt128 Element10; + public readonly GpuUInt128 Element11; + public readonly GpuUInt128 Element12; + public readonly GpuUInt128 Element13; + public readonly GpuUInt128 Element14; + public readonly GpuUInt128 Element15; + public readonly GpuUInt128 Element16; + public readonly GpuUInt128 Element17; + public readonly GpuUInt128 Element18; + public readonly GpuUInt128 Element19; + public readonly GpuUInt128 Element20; + public readonly GpuUInt128 Element21; + public readonly GpuUInt128 Element22; + public readonly GpuUInt128 Element23; + public readonly GpuUInt128 Element24; + public readonly GpuUInt128 Element25; + public readonly GpuUInt128 Element26; + public readonly GpuUInt128 Element27; + public readonly GpuUInt128 Element28; + public readonly GpuUInt128 Element29; + public readonly GpuUInt128 Element30; + public readonly GpuUInt128 Element31; + public readonly GpuUInt128 Element32; + public readonly GpuUInt128 Element33; + public readonly GpuUInt128 Element34; + public readonly GpuUInt128 Element35; + public readonly GpuUInt128 Element36; + public readonly GpuUInt128 Element37; + public readonly GpuUInt128 Element38; + public readonly GpuUInt128 Element39; + public readonly GpuUInt128 Element40; + public readonly GpuUInt128 Element41; + public readonly GpuUInt128 Element42; + public readonly GpuUInt128 Element43; + public readonly GpuUInt128 Element44; + public readonly GpuUInt128 Element45; + public readonly GpuUInt128 Element46; + public readonly GpuUInt128 Element47; + public readonly GpuUInt128 Element48; + public readonly GpuUInt128 Element49; + public readonly GpuUInt128 Element50; + public readonly GpuUInt128 Element51; + public readonly GpuUInt128 Element52; + public readonly GpuUInt128 Element53; + public readonly GpuUInt128 Element54; + public readonly GpuUInt128 Element55; + public readonly GpuUInt128 Element56; + public readonly GpuUInt128 Element57; + public readonly GpuUInt128 Element58; + public readonly GpuUInt128 Element59; + public readonly GpuUInt128 Element60; + public readonly GpuUInt128 Element61; + public readonly GpuUInt128 Element62; + public readonly GpuUInt128 Element63; + public readonly GpuUInt128 Element64; + public readonly GpuUInt128 Element65; + public readonly GpuUInt128 Element66; + public readonly GpuUInt128 Element67; + public readonly GpuUInt128 Element68; + public readonly GpuUInt128 Element69; + public readonly GpuUInt128 Element70; + public readonly GpuUInt128 Element71; + public readonly GpuUInt128 Element72; + public readonly GpuUInt128 Element73; + public readonly GpuUInt128 Element74; + public readonly GpuUInt128 Element75; + public readonly GpuUInt128 Element76; + public readonly GpuUInt128 Element77; + public readonly GpuUInt128 Element78; + public readonly GpuUInt128 Element79; + public readonly GpuUInt128 Element80; + public readonly GpuUInt128 Element81; + public readonly GpuUInt128 Element82; + public readonly GpuUInt128 Element83; + public readonly GpuUInt128 Element84; + public readonly GpuUInt128 Element85; + public readonly GpuUInt128 Element86; + public readonly GpuUInt128 Element87; + public readonly GpuUInt128 Element88; + public readonly GpuUInt128 Element89; + public readonly GpuUInt128 Element90; + public readonly GpuUInt128 Element91; + public readonly GpuUInt128 Element92; + public readonly GpuUInt128 Element93; + public readonly GpuUInt128 Element94; + public readonly GpuUInt128 Element95; + public readonly GpuUInt128 Element96; + public readonly GpuUInt128 Element97; + public readonly GpuUInt128 Element98; + public readonly GpuUInt128 Element99; + public readonly GpuUInt128 Element100; + public readonly GpuUInt128 Element101; + public readonly GpuUInt128 Element102; + public readonly GpuUInt128 Element103; + public readonly GpuUInt128 Element104; + public readonly GpuUInt128 Element105; + public readonly GpuUInt128 Element106; + public readonly GpuUInt128 Element107; + public readonly GpuUInt128 Element108; + public readonly GpuUInt128 Element109; + public readonly GpuUInt128 Element110; + public readonly GpuUInt128 Element111; + public readonly GpuUInt128 Element112; + public readonly GpuUInt128 Element113; + public readonly GpuUInt128 Element114; + public readonly GpuUInt128 Element115; + public readonly GpuUInt128 Element116; + public readonly GpuUInt128 Element117; + public readonly GpuUInt128 Element118; + public readonly GpuUInt128 Element119; + public readonly GpuUInt128 Element120; + public readonly GpuUInt128 Element121; + public readonly GpuUInt128 Element122; + public readonly GpuUInt128 Element123; + public readonly GpuUInt128 Element124; + public readonly GpuUInt128 Element125; + public readonly GpuUInt128 Element126; + public readonly GpuUInt128 Element127; + + public readonly GpuUInt128 this[int index] + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return index switch + { + 0 => Element0, + 1 => Element1, + 2 => Element2, + 3 => Element3, + 4 => Element4, + 5 => Element5, + 6 => Element6, + 7 => Element7, + 8 => Element8, + 9 => Element9, + 10 => Element10, + 11 => Element11, + 12 => Element12, + 13 => Element13, + 14 => Element14, + 15 => Element15, + 16 => Element16, + 17 => Element17, + 18 => Element18, + 19 => Element19, + 20 => Element20, + 21 => Element21, + 22 => Element22, + 23 => Element23, + 24 => Element24, + 25 => Element25, + 26 => Element26, + 27 => Element27, + 28 => Element28, + 29 => Element29, + 30 => Element30, + 31 => Element31, + 32 => Element32, + 33 => Element33, + 34 => Element34, + 35 => Element35, + 36 => Element36, + 37 => Element37, + 38 => Element38, + 39 => Element39, + 40 => Element40, + 41 => Element41, + 42 => Element42, + 43 => Element43, + 44 => Element44, + 45 => Element45, + 46 => Element46, + 47 => Element47, + 48 => Element48, + 49 => Element49, + 50 => Element50, + 51 => Element51, + 52 => Element52, + 53 => Element53, + 54 => Element54, + 55 => Element55, + 56 => Element56, + 57 => Element57, + 58 => Element58, + 59 => Element59, + 60 => Element60, + 61 => Element61, + 62 => Element62, + 63 => Element63, + 64 => Element64, + 65 => Element65, + 66 => Element66, + 67 => Element67, + 68 => Element68, + 69 => Element69, + 70 => Element70, + 71 => Element71, + 72 => Element72, + 73 => Element73, + 74 => Element74, + 75 => Element75, + 76 => Element76, + 77 => Element77, + 78 => Element78, + 79 => Element79, + 80 => Element80, + 81 => Element81, + 82 => Element82, + 83 => Element83, + 84 => Element84, + 85 => Element85, + 86 => Element86, + 87 => Element87, + 88 => Element88, + 89 => Element89, + 90 => Element90, + 91 => Element91, + 92 => Element92, + 93 => Element93, + 94 => Element94, + 95 => Element95, + 96 => Element96, + 97 => Element97, + 98 => Element98, + 99 => Element99, + 100 => Element100, + 101 => Element101, + 102 => Element102, + 103 => Element103, + 104 => Element104, + 105 => Element105, + 106 => Element106, + 107 => Element107, + 108 => Element108, + 109 => Element109, + 110 => Element110, + 111 => Element111, + 112 => Element112, + 113 => Element113, + 114 => Element114, + 115 => Element115, + 116 => Element116, + 117 => Element117, + 118 => Element118, + 119 => Element119, + 120 => Element120, + 121 => Element121, + 122 => Element122, + 123 => Element123, + 124 => Element124, + 125 => Element125, + 126 => Element126, + 127 => Element127, + _ => Element0, + }; + } + } + } + + public readonly struct Pow2OddPowerTableReturn + { + public Pow2OddPowerTableReturn(GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount) + { + Element0 = baseValue; + if (oddPowerCount == 1) + { + return; + } + + // Reusing baseValue to hold base^2 for the shared odd-power ladder that follows. + baseValue.MulMod(baseValue, modulus); + // TODO: We can calculate baseValue % modulus before loop and use it to increase ladderEntry calculation speed - we'll reuse the base for incremental calculations. + GpuUInt128 current = baseValue; + + // We're manually assigning each field to prevent the compiler to initialize each field twice due to auto-initialization. We're using the action to lower the code base size. + current.MulMod(baseValue, modulus); + Element1 = current; + current.MulMod(baseValue, modulus); + Element2 = current; + current.MulMod(baseValue, modulus); + Element3 = current; + current.MulMod(baseValue, modulus); + Element4 = current; + current.MulMod(baseValue, modulus); + Element5 = current; + current.MulMod(baseValue, modulus); + Element6 = current; + current.MulMod(baseValue, modulus); + Element7 = current; + current.MulMod(baseValue, modulus); + Element8 = current; + current.MulMod(baseValue, modulus); + Element9 = current; + current.MulMod(baseValue, modulus); + Element10 = current; + current.MulMod(baseValue, modulus); + Element11 = current; + current.MulMod(baseValue, modulus); + Element12 = current; + current.MulMod(baseValue, modulus); + Element13 = current; + current.MulMod(baseValue, modulus); + Element14 = current; + current.MulMod(baseValue, modulus); + Element15 = current; + current.MulMod(baseValue, modulus); + Element16 = current; + current.MulMod(baseValue, modulus); + Element17 = current; + current.MulMod(baseValue, modulus); + Element18 = current; + current.MulMod(baseValue, modulus); + Element19 = current; + current.MulMod(baseValue, modulus); + Element20 = current; + current.MulMod(baseValue, modulus); + Element21 = current; + current.MulMod(baseValue, modulus); + Element22 = current; + current.MulMod(baseValue, modulus); + Element23 = current; + current.MulMod(baseValue, modulus); + Element24 = current; + current.MulMod(baseValue, modulus); + Element25 = current; + current.MulMod(baseValue, modulus); + Element26 = current; + current.MulMod(baseValue, modulus); + Element27 = current; + current.MulMod(baseValue, modulus); + Element28 = current; + current.MulMod(baseValue, modulus); + Element29 = current; + current.MulMod(baseValue, modulus); + Element30 = current; + current.MulMod(baseValue, modulus); + Element31 = current; + current.MulMod(baseValue, modulus); + Element32 = current; + current.MulMod(baseValue, modulus); + Element33 = current; + current.MulMod(baseValue, modulus); + Element34 = current; + current.MulMod(baseValue, modulus); + Element35 = current; + current.MulMod(baseValue, modulus); + Element36 = current; + current.MulMod(baseValue, modulus); + Element37 = current; + current.MulMod(baseValue, modulus); + Element38 = current; + current.MulMod(baseValue, modulus); + Element39 = current; + current.MulMod(baseValue, modulus); + Element40 = current; + current.MulMod(baseValue, modulus); + Element41 = current; + current.MulMod(baseValue, modulus); + Element42 = current; + current.MulMod(baseValue, modulus); + Element43 = current; + current.MulMod(baseValue, modulus); + Element44 = current; + current.MulMod(baseValue, modulus); + Element45 = current; + current.MulMod(baseValue, modulus); + Element46 = current; + current.MulMod(baseValue, modulus); + Element47 = current; + current.MulMod(baseValue, modulus); + Element48 = current; + current.MulMod(baseValue, modulus); + Element49 = current; + current.MulMod(baseValue, modulus); + Element50 = current; + current.MulMod(baseValue, modulus); + Element51 = current; + current.MulMod(baseValue, modulus); + Element52 = current; + current.MulMod(baseValue, modulus); + Element53 = current; + current.MulMod(baseValue, modulus); + Element54 = current; + current.MulMod(baseValue, modulus); + Element55 = current; + current.MulMod(baseValue, modulus); + Element56 = current; + current.MulMod(baseValue, modulus); + Element57 = current; + current.MulMod(baseValue, modulus); + Element58 = current; + current.MulMod(baseValue, modulus); + Element59 = current; + current.MulMod(baseValue, modulus); + Element60 = current; + current.MulMod(baseValue, modulus); + Element61 = current; + current.MulMod(baseValue, modulus); + Element62 = current; + current.MulMod(baseValue, modulus); + Element63 = current; + current.MulMod(baseValue, modulus); + Element64 = current; + current.MulMod(baseValue, modulus); + Element65 = current; + current.MulMod(baseValue, modulus); + Element66 = current; + current.MulMod(baseValue, modulus); + Element67 = current; + current.MulMod(baseValue, modulus); + Element68 = current; + current.MulMod(baseValue, modulus); + Element69 = current; + current.MulMod(baseValue, modulus); + Element70 = current; + current.MulMod(baseValue, modulus); + Element71 = current; + current.MulMod(baseValue, modulus); + Element72 = current; + current.MulMod(baseValue, modulus); + Element73 = current; + current.MulMod(baseValue, modulus); + Element74 = current; + current.MulMod(baseValue, modulus); + Element75 = current; + current.MulMod(baseValue, modulus); + Element76 = current; + current.MulMod(baseValue, modulus); + Element77 = current; + current.MulMod(baseValue, modulus); + Element78 = current; + current.MulMod(baseValue, modulus); + Element79 = current; + current.MulMod(baseValue, modulus); + Element80 = current; + current.MulMod(baseValue, modulus); + Element81 = current; + current.MulMod(baseValue, modulus); + Element82 = current; + current.MulMod(baseValue, modulus); + Element83 = current; + current.MulMod(baseValue, modulus); + Element84 = current; + current.MulMod(baseValue, modulus); + Element85 = current; + current.MulMod(baseValue, modulus); + Element86 = current; + current.MulMod(baseValue, modulus); + Element87 = current; + current.MulMod(baseValue, modulus); + Element88 = current; + current.MulMod(baseValue, modulus); + Element89 = current; + current.MulMod(baseValue, modulus); + Element90 = current; + current.MulMod(baseValue, modulus); + Element91 = current; + current.MulMod(baseValue, modulus); + Element92 = current; + current.MulMod(baseValue, modulus); + Element93 = current; + current.MulMod(baseValue, modulus); + Element94 = current; + current.MulMod(baseValue, modulus); + Element95 = current; + current.MulMod(baseValue, modulus); + Element96 = current; + current.MulMod(baseValue, modulus); + Element97 = current; + current.MulMod(baseValue, modulus); + Element98 = current; + current.MulMod(baseValue, modulus); + Element99 = current; + current.MulMod(baseValue, modulus); + Element100 = current; + current.MulMod(baseValue, modulus); + Element101 = current; + current.MulMod(baseValue, modulus); + Element102 = current; + current.MulMod(baseValue, modulus); + Element103 = current; + current.MulMod(baseValue, modulus); + Element104 = current; + current.MulMod(baseValue, modulus); + Element105 = current; + current.MulMod(baseValue, modulus); + Element106 = current; + current.MulMod(baseValue, modulus); + Element107 = current; + current.MulMod(baseValue, modulus); + Element108 = current; + current.MulMod(baseValue, modulus); + Element109 = current; + current.MulMod(baseValue, modulus); + Element110 = current; + current.MulMod(baseValue, modulus); + Element111 = current; + current.MulMod(baseValue, modulus); + Element112 = current; + current.MulMod(baseValue, modulus); + Element113 = current; + current.MulMod(baseValue, modulus); + Element114 = current; + current.MulMod(baseValue, modulus); + Element115 = current; + current.MulMod(baseValue, modulus); + Element116 = current; + current.MulMod(baseValue, modulus); + Element117 = current; + current.MulMod(baseValue, modulus); + Element118 = current; + current.MulMod(baseValue, modulus); + Element119 = current; + current.MulMod(baseValue, modulus); + Element120 = current; + current.MulMod(baseValue, modulus); + Element121 = current; + current.MulMod(baseValue, modulus); + Element122 = current; + current.MulMod(baseValue, modulus); + Element123 = current; + current.MulMod(baseValue, modulus); + Element124 = current; + current.MulMod(baseValue, modulus); + Element125 = current; + current.MulMod(baseValue, modulus); + Element126 = current; + current.MulMod(baseValue, modulus); + Element127 = current; + } + + public readonly GpuUInt128 Element0; + public readonly GpuUInt128 Element1; + public readonly GpuUInt128 Element2; + public readonly GpuUInt128 Element3; + public readonly GpuUInt128 Element4; + public readonly GpuUInt128 Element5; + public readonly GpuUInt128 Element6; + public readonly GpuUInt128 Element7; + public readonly GpuUInt128 Element8; + public readonly GpuUInt128 Element9; + public readonly GpuUInt128 Element10; + public readonly GpuUInt128 Element11; + public readonly GpuUInt128 Element12; + public readonly GpuUInt128 Element13; + public readonly GpuUInt128 Element14; + public readonly GpuUInt128 Element15; + public readonly GpuUInt128 Element16; + public readonly GpuUInt128 Element17; + public readonly GpuUInt128 Element18; + public readonly GpuUInt128 Element19; + public readonly GpuUInt128 Element20; + public readonly GpuUInt128 Element21; + public readonly GpuUInt128 Element22; + public readonly GpuUInt128 Element23; + public readonly GpuUInt128 Element24; + public readonly GpuUInt128 Element25; + public readonly GpuUInt128 Element26; + public readonly GpuUInt128 Element27; + public readonly GpuUInt128 Element28; + public readonly GpuUInt128 Element29; + public readonly GpuUInt128 Element30; + public readonly GpuUInt128 Element31; + public readonly GpuUInt128 Element32; + public readonly GpuUInt128 Element33; + public readonly GpuUInt128 Element34; + public readonly GpuUInt128 Element35; + public readonly GpuUInt128 Element36; + public readonly GpuUInt128 Element37; + public readonly GpuUInt128 Element38; + public readonly GpuUInt128 Element39; + public readonly GpuUInt128 Element40; + public readonly GpuUInt128 Element41; + public readonly GpuUInt128 Element42; + public readonly GpuUInt128 Element43; + public readonly GpuUInt128 Element44; + public readonly GpuUInt128 Element45; + public readonly GpuUInt128 Element46; + public readonly GpuUInt128 Element47; + public readonly GpuUInt128 Element48; + public readonly GpuUInt128 Element49; + public readonly GpuUInt128 Element50; + public readonly GpuUInt128 Element51; + public readonly GpuUInt128 Element52; + public readonly GpuUInt128 Element53; + public readonly GpuUInt128 Element54; + public readonly GpuUInt128 Element55; + public readonly GpuUInt128 Element56; + public readonly GpuUInt128 Element57; + public readonly GpuUInt128 Element58; + public readonly GpuUInt128 Element59; + public readonly GpuUInt128 Element60; + public readonly GpuUInt128 Element61; + public readonly GpuUInt128 Element62; + public readonly GpuUInt128 Element63; + public readonly GpuUInt128 Element64; + public readonly GpuUInt128 Element65; + public readonly GpuUInt128 Element66; + public readonly GpuUInt128 Element67; + public readonly GpuUInt128 Element68; + public readonly GpuUInt128 Element69; + public readonly GpuUInt128 Element70; + public readonly GpuUInt128 Element71; + public readonly GpuUInt128 Element72; + public readonly GpuUInt128 Element73; + public readonly GpuUInt128 Element74; + public readonly GpuUInt128 Element75; + public readonly GpuUInt128 Element76; + public readonly GpuUInt128 Element77; + public readonly GpuUInt128 Element78; + public readonly GpuUInt128 Element79; + public readonly GpuUInt128 Element80; + public readonly GpuUInt128 Element81; + public readonly GpuUInt128 Element82; + public readonly GpuUInt128 Element83; + public readonly GpuUInt128 Element84; + public readonly GpuUInt128 Element85; + public readonly GpuUInt128 Element86; + public readonly GpuUInt128 Element87; + public readonly GpuUInt128 Element88; + public readonly GpuUInt128 Element89; + public readonly GpuUInt128 Element90; + public readonly GpuUInt128 Element91; + public readonly GpuUInt128 Element92; + public readonly GpuUInt128 Element93; + public readonly GpuUInt128 Element94; + public readonly GpuUInt128 Element95; + public readonly GpuUInt128 Element96; + public readonly GpuUInt128 Element97; + public readonly GpuUInt128 Element98; + public readonly GpuUInt128 Element99; + public readonly GpuUInt128 Element100; + public readonly GpuUInt128 Element101; + public readonly GpuUInt128 Element102; + public readonly GpuUInt128 Element103; + public readonly GpuUInt128 Element104; + public readonly GpuUInt128 Element105; + public readonly GpuUInt128 Element106; + public readonly GpuUInt128 Element107; + public readonly GpuUInt128 Element108; + public readonly GpuUInt128 Element109; + public readonly GpuUInt128 Element110; + public readonly GpuUInt128 Element111; + public readonly GpuUInt128 Element112; + public readonly GpuUInt128 Element113; + public readonly GpuUInt128 Element114; + public readonly GpuUInt128 Element115; + public readonly GpuUInt128 Element116; + public readonly GpuUInt128 Element117; + public readonly GpuUInt128 Element118; + public readonly GpuUInt128 Element119; + public readonly GpuUInt128 Element120; + public readonly GpuUInt128 Element121; + public readonly GpuUInt128 Element122; + public readonly GpuUInt128 Element123; + public readonly GpuUInt128 Element124; + public readonly GpuUInt128 Element125; + public readonly GpuUInt128 Element126; + public readonly GpuUInt128 Element127; + + public readonly GpuUInt128 this[int index] + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return index switch + { + 0 => Element0, + 1 => Element1, + 2 => Element2, + 3 => Element3, + 4 => Element4, + 5 => Element5, + 6 => Element6, + 7 => Element7, + 8 => Element8, + 9 => Element9, + 10 => Element10, + 11 => Element11, + 12 => Element12, + 13 => Element13, + 14 => Element14, + 15 => Element15, + 16 => Element16, + 17 => Element17, + 18 => Element18, + 19 => Element19, + 20 => Element20, + 21 => Element21, + 22 => Element22, + 23 => Element23, + 24 => Element24, + 25 => Element25, + 26 => Element26, + 27 => Element27, + 28 => Element28, + 29 => Element29, + 30 => Element30, + 31 => Element31, + 32 => Element32, + 33 => Element33, + 34 => Element34, + 35 => Element35, + 36 => Element36, + 37 => Element37, + 38 => Element38, + 39 => Element39, + 40 => Element40, + 41 => Element41, + 42 => Element42, + 43 => Element43, + 44 => Element44, + 45 => Element45, + 46 => Element46, + 47 => Element47, + 48 => Element48, + 49 => Element49, + 50 => Element50, + 51 => Element51, + 52 => Element52, + 53 => Element53, + 54 => Element54, + 55 => Element55, + 56 => Element56, + 57 => Element57, + 58 => Element58, + 59 => Element59, + 60 => Element60, + 61 => Element61, + 62 => Element62, + 63 => Element63, + 64 => Element64, + 65 => Element65, + 66 => Element66, + 67 => Element67, + 68 => Element68, + 69 => Element69, + 70 => Element70, + 71 => Element71, + 72 => Element72, + 73 => Element73, + 74 => Element74, + 75 => Element75, + 76 => Element76, + 77 => Element77, + 78 => Element78, + 79 => Element79, + 80 => Element80, + 81 => Element81, + 82 => Element82, + 83 => Element83, + 84 => Element84, + 85 => Element85, + 86 => Element86, + 87 => Element87, + 88 => Element88, + 89 => Element89, + 90 => Element90, + 91 => Element91, + 92 => Element92, + 93 => Element93, + 94 => Element94, + 95 => Element95, + 96 => Element96, + 97 => Element97, + 98 => Element98, + 99 => Element99, + 100 => Element100, + 101 => Element101, + 102 => Element102, + 103 => Element103, + 104 => Element104, + 105 => Element105, + 106 => Element106, + 107 => Element107, + 108 => Element108, + 109 => Element109, + 110 => Element110, + 111 => Element111, + 112 => Element112, + 113 => Element113, + 114 => Element114, + 115 => Element115, + 116 => Element116, + 117 => Element117, + 118 => Element118, + 119 => Element119, + 120 => Element120, + 121 => Element121, + 122 => Element122, + 123 => Element123, + 124 => Element124, + 125 => Element125, + 126 => Element126, + 127 => Element127, + _ => Element0, + }; + } + } + } + + public readonly struct Pow2OddPowerTableLoop + { + public Pow2OddPowerTableLoop(GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount) + { + Element0 = baseValue; + if (oddPowerCount == 1) + { + return; + } + + // Reusing baseValue to hold base^2 for the shared odd-power ladder that follows. + baseValue.MulMod(baseValue, modulus); + // TODO: We can calculate baseValue % modulus before loop and use it to increase ladderEntry calculation speed - we'll reuse the base for incremental calculations. + GpuUInt128 current = baseValue; + + // We're manually assigning each field to prevent the compiler to initialize each field twice due to auto-initialization. We're using the action to lower the code base size. + for (int i = 1; i < oddPowerCount; i++) + { + current.MulMod(baseValue, modulus); + switch (i) + { + case 1: + Element1 = current; + break; + case 2: + Element2 = current; + break; + case 3: + Element3 = current; + break; + case 4: + Element4 = current; + break; + case 5: + Element5 = current; + break; + case 6: + Element6 = current; + break; + case 7: + Element7 = current; + break; + case 8: + Element8 = current; + break; + case 9: + Element9 = current; + break; + case 10: + Element10 = current; + break; + case 11: + Element11 = current; + break; + case 12: + Element12 = current; + break; + case 13: + Element13 = current; + break; + case 14: + Element14 = current; + break; + case 15: + Element15 = current; + break; + case 16: + Element16 = current; + break; + case 17: + Element17 = current; + break; + case 18: + Element18 = current; + break; + case 19: + Element19 = current; + break; + case 20: + Element20 = current; + break; + case 21: + Element21 = current; + break; + case 22: + Element22 = current; + break; + case 23: + Element23 = current; + break; + case 24: + Element24 = current; + break; + case 25: + Element25 = current; + break; + case 26: + Element26 = current; + break; + case 27: + Element27 = current; + break; + case 28: + Element28 = current; + break; + case 29: + Element29 = current; + break; + case 30: + Element30 = current; + break; + case 31: + Element31 = current; + break; + case 32: + Element32 = current; + break; + case 33: + Element33 = current; + break; + case 34: + Element34 = current; + break; + case 35: + Element35 = current; + break; + case 36: + Element36 = current; + break; + case 37: + Element37 = current; + break; + case 38: + Element38 = current; + break; + case 39: + Element39 = current; + break; + case 40: + Element40 = current; + break; + case 41: + Element41 = current; + break; + case 42: + Element42 = current; + break; + case 43: + Element43 = current; + break; + case 44: + Element44 = current; + break; + case 45: + Element45 = current; + break; + case 46: + Element46 = current; + break; + case 47: + Element47 = current; + break; + case 48: + Element48 = current; + break; + case 49: + Element49 = current; + break; + case 50: + Element50 = current; + break; + case 51: + Element51 = current; + break; + case 52: + Element52 = current; + break; + case 53: + Element53 = current; + break; + case 54: + Element54 = current; + break; + case 55: + Element55 = current; + break; + case 56: + Element56 = current; + break; + case 57: + Element57 = current; + break; + case 58: + Element58 = current; + break; + case 59: + Element59 = current; + break; + case 60: + Element60 = current; + break; + case 61: + Element61 = current; + break; + case 62: + Element62 = current; + break; + case 63: + Element63 = current; + break; + case 64: + Element64 = current; + break; + case 65: + Element65 = current; + break; + case 66: + Element66 = current; + break; + case 67: + Element67 = current; + break; + case 68: + Element68 = current; + break; + case 69: + Element69 = current; + break; + case 70: + Element70 = current; + break; + case 71: + Element71 = current; + break; + case 72: + Element72 = current; + break; + case 73: + Element73 = current; + break; + case 74: + Element74 = current; + break; + case 75: + Element75 = current; + break; + case 76: + Element76 = current; + break; + case 77: + Element77 = current; + break; + case 78: + Element78 = current; + break; + case 79: + Element79 = current; + break; + case 80: + Element80 = current; + break; + case 81: + Element81 = current; + break; + case 82: + Element82 = current; + break; + case 83: + Element83 = current; + break; + case 84: + Element84 = current; + break; + case 85: + Element85 = current; + break; + case 86: + Element86 = current; + break; + case 87: + Element87 = current; + break; + case 88: + Element88 = current; + break; + case 89: + Element89 = current; + break; + case 90: + Element90 = current; + break; + case 91: + Element91 = current; + break; + case 92: + Element92 = current; + break; + case 93: + Element93 = current; + break; + case 94: + Element94 = current; + break; + case 95: + Element95 = current; + break; + case 96: + Element96 = current; + break; + case 97: + Element97 = current; + break; + case 98: + Element98 = current; + break; + case 99: + Element99 = current; + break; + case 100: + Element100 = current; + break; + case 101: + Element101 = current; + break; + case 102: + Element102 = current; + break; + case 103: + Element103 = current; + break; + case 104: + Element104 = current; + break; + case 105: + Element105 = current; + break; + case 106: + Element106 = current; + break; + case 107: + Element107 = current; + break; + case 108: + Element108 = current; + break; + case 109: + Element109 = current; + break; + case 110: + Element110 = current; + break; + case 111: + Element111 = current; + break; + case 112: + Element112 = current; + break; + case 113: + Element113 = current; + break; + case 114: + Element114 = current; + break; + case 115: + Element115 = current; + break; + case 116: + Element116 = current; + break; + case 117: + Element117 = current; + break; + case 118: + Element118 = current; + break; + case 119: + Element119 = current; + break; + case 120: + Element120 = current; + break; + case 121: + Element121 = current; + break; + case 122: + Element122 = current; + break; + case 123: + Element123 = current; + break; + case 124: + Element124 = current; + break; + case 125: + Element125 = current; + break; + case 126: + Element126 = current; + break; + case 127: + Element127 = current; + break; + } + } + } + + public readonly GpuUInt128 Element0; + public readonly GpuUInt128 Element1; + public readonly GpuUInt128 Element2; + public readonly GpuUInt128 Element3; + public readonly GpuUInt128 Element4; + public readonly GpuUInt128 Element5; + public readonly GpuUInt128 Element6; + public readonly GpuUInt128 Element7; + public readonly GpuUInt128 Element8; + public readonly GpuUInt128 Element9; + public readonly GpuUInt128 Element10; + public readonly GpuUInt128 Element11; + public readonly GpuUInt128 Element12; + public readonly GpuUInt128 Element13; + public readonly GpuUInt128 Element14; + public readonly GpuUInt128 Element15; + public readonly GpuUInt128 Element16; + public readonly GpuUInt128 Element17; + public readonly GpuUInt128 Element18; + public readonly GpuUInt128 Element19; + public readonly GpuUInt128 Element20; + public readonly GpuUInt128 Element21; + public readonly GpuUInt128 Element22; + public readonly GpuUInt128 Element23; + public readonly GpuUInt128 Element24; + public readonly GpuUInt128 Element25; + public readonly GpuUInt128 Element26; + public readonly GpuUInt128 Element27; + public readonly GpuUInt128 Element28; + public readonly GpuUInt128 Element29; + public readonly GpuUInt128 Element30; + public readonly GpuUInt128 Element31; + public readonly GpuUInt128 Element32; + public readonly GpuUInt128 Element33; + public readonly GpuUInt128 Element34; + public readonly GpuUInt128 Element35; + public readonly GpuUInt128 Element36; + public readonly GpuUInt128 Element37; + public readonly GpuUInt128 Element38; + public readonly GpuUInt128 Element39; + public readonly GpuUInt128 Element40; + public readonly GpuUInt128 Element41; + public readonly GpuUInt128 Element42; + public readonly GpuUInt128 Element43; + public readonly GpuUInt128 Element44; + public readonly GpuUInt128 Element45; + public readonly GpuUInt128 Element46; + public readonly GpuUInt128 Element47; + public readonly GpuUInt128 Element48; + public readonly GpuUInt128 Element49; + public readonly GpuUInt128 Element50; + public readonly GpuUInt128 Element51; + public readonly GpuUInt128 Element52; + public readonly GpuUInt128 Element53; + public readonly GpuUInt128 Element54; + public readonly GpuUInt128 Element55; + public readonly GpuUInt128 Element56; + public readonly GpuUInt128 Element57; + public readonly GpuUInt128 Element58; + public readonly GpuUInt128 Element59; + public readonly GpuUInt128 Element60; + public readonly GpuUInt128 Element61; + public readonly GpuUInt128 Element62; + public readonly GpuUInt128 Element63; + public readonly GpuUInt128 Element64; + public readonly GpuUInt128 Element65; + public readonly GpuUInt128 Element66; + public readonly GpuUInt128 Element67; + public readonly GpuUInt128 Element68; + public readonly GpuUInt128 Element69; + public readonly GpuUInt128 Element70; + public readonly GpuUInt128 Element71; + public readonly GpuUInt128 Element72; + public readonly GpuUInt128 Element73; + public readonly GpuUInt128 Element74; + public readonly GpuUInt128 Element75; + public readonly GpuUInt128 Element76; + public readonly GpuUInt128 Element77; + public readonly GpuUInt128 Element78; + public readonly GpuUInt128 Element79; + public readonly GpuUInt128 Element80; + public readonly GpuUInt128 Element81; + public readonly GpuUInt128 Element82; + public readonly GpuUInt128 Element83; + public readonly GpuUInt128 Element84; + public readonly GpuUInt128 Element85; + public readonly GpuUInt128 Element86; + public readonly GpuUInt128 Element87; + public readonly GpuUInt128 Element88; + public readonly GpuUInt128 Element89; + public readonly GpuUInt128 Element90; + public readonly GpuUInt128 Element91; + public readonly GpuUInt128 Element92; + public readonly GpuUInt128 Element93; + public readonly GpuUInt128 Element94; + public readonly GpuUInt128 Element95; + public readonly GpuUInt128 Element96; + public readonly GpuUInt128 Element97; + public readonly GpuUInt128 Element98; + public readonly GpuUInt128 Element99; + public readonly GpuUInt128 Element100; + public readonly GpuUInt128 Element101; + public readonly GpuUInt128 Element102; + public readonly GpuUInt128 Element103; + public readonly GpuUInt128 Element104; + public readonly GpuUInt128 Element105; + public readonly GpuUInt128 Element106; + public readonly GpuUInt128 Element107; + public readonly GpuUInt128 Element108; + public readonly GpuUInt128 Element109; + public readonly GpuUInt128 Element110; + public readonly GpuUInt128 Element111; + public readonly GpuUInt128 Element112; + public readonly GpuUInt128 Element113; + public readonly GpuUInt128 Element114; + public readonly GpuUInt128 Element115; + public readonly GpuUInt128 Element116; + public readonly GpuUInt128 Element117; + public readonly GpuUInt128 Element118; + public readonly GpuUInt128 Element119; + public readonly GpuUInt128 Element120; + public readonly GpuUInt128 Element121; + public readonly GpuUInt128 Element122; + public readonly GpuUInt128 Element123; + public readonly GpuUInt128 Element124; + public readonly GpuUInt128 Element125; + public readonly GpuUInt128 Element126; + public readonly GpuUInt128 Element127; + + public readonly GpuUInt128 this[int index] + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get => index switch + { + 0 => Element0, + 1 => Element1, + 2 => Element2, + 3 => Element3, + 4 => Element4, + 5 => Element5, + 6 => Element6, + 7 => Element7, + 8 => Element8, + 9 => Element9, + 10 => Element10, + 11 => Element11, + 12 => Element12, + 13 => Element13, + 14 => Element14, + 15 => Element15, + 16 => Element16, + 17 => Element17, + 18 => Element18, + 19 => Element19, + 20 => Element20, + 21 => Element21, + 22 => Element22, + 23 => Element23, + 24 => Element24, + 25 => Element25, + 26 => Element26, + 27 => Element27, + 28 => Element28, + 29 => Element29, + 30 => Element30, + 31 => Element31, + 32 => Element32, + 33 => Element33, + 34 => Element34, + 35 => Element35, + 36 => Element36, + 37 => Element37, + 38 => Element38, + 39 => Element39, + 40 => Element40, + 41 => Element41, + 42 => Element42, + 43 => Element43, + 44 => Element44, + 45 => Element45, + 46 => Element46, + 47 => Element47, + 48 => Element48, + 49 => Element49, + 50 => Element50, + 51 => Element51, + 52 => Element52, + 53 => Element53, + 54 => Element54, + 55 => Element55, + 56 => Element56, + 57 => Element57, + 58 => Element58, + 59 => Element59, + 60 => Element60, + 61 => Element61, + 62 => Element62, + 63 => Element63, + 64 => Element64, + 65 => Element65, + 66 => Element66, + 67 => Element67, + 68 => Element68, + 69 => Element69, + 70 => Element70, + 71 => Element71, + 72 => Element72, + 73 => Element73, + 74 => Element74, + 75 => Element75, + 76 => Element76, + 77 => Element77, + 78 => Element78, + 79 => Element79, + 80 => Element80, + 81 => Element81, + 82 => Element82, + 83 => Element83, + 84 => Element84, + 85 => Element85, + 86 => Element86, + 87 => Element87, + 88 => Element88, + 89 => Element89, + 90 => Element90, + 91 => Element91, + 92 => Element92, + 93 => Element93, + 94 => Element94, + 95 => Element95, + 96 => Element96, + 97 => Element97, + 98 => Element98, + 99 => Element99, + 100 => Element100, + 101 => Element101, + 102 => Element102, + 103 => Element103, + 104 => Element104, + 105 => Element105, + 106 => Element106, + 107 => Element107, + 108 => Element108, + 109 => Element109, + 110 => Element110, + 111 => Element111, + 112 => Element112, + 113 => Element113, + 114 => Element114, + 115 => Element115, + 116 => Element116, + 117 => Element117, + 118 => Element118, + 119 => Element119, + 120 => Element120, + 121 => Element121, + 122 => Element122, + 123 => Element123, + 124 => Element124, + 125 => Element125, + 126 => Element126, + 127 => Element127, + _ => Element0, + }; + } + } + + public readonly struct Pow2OddPowerTableSetter + { + public Pow2OddPowerTableSetter(GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount) + { + Element0 = baseValue; + if (oddPowerCount == 1) + { + return; + } + + // Reusing baseValue to hold base^2 for the shared odd-power ladder that follows. + baseValue.MulMod(baseValue, modulus); + // TODO: We can calculate baseValue % modulus before loop and use it to increase ladderEntry calculation speed - we'll reuse the base for incremental calculations. + GpuUInt128 current = baseValue; + + // We're manually assigning each field to prevent the compiler to initialize each field twice due to auto-initialization. We're using the action to lower the code base size. + for (int i = 1; i < oddPowerCount; i++) + { + current.MulMod(baseValue, modulus); + this[i] = current; + } + } + + public readonly GpuUInt128 Element0; + public readonly GpuUInt128 Element1; + public readonly GpuUInt128 Element2; + public readonly GpuUInt128 Element3; + public readonly GpuUInt128 Element4; + public readonly GpuUInt128 Element5; + public readonly GpuUInt128 Element6; + public readonly GpuUInt128 Element7; + public readonly GpuUInt128 Element8; + public readonly GpuUInt128 Element9; + public readonly GpuUInt128 Element10; + public readonly GpuUInt128 Element11; + public readonly GpuUInt128 Element12; + public readonly GpuUInt128 Element13; + public readonly GpuUInt128 Element14; + public readonly GpuUInt128 Element15; + public readonly GpuUInt128 Element16; + public readonly GpuUInt128 Element17; + public readonly GpuUInt128 Element18; + public readonly GpuUInt128 Element19; + public readonly GpuUInt128 Element20; + public readonly GpuUInt128 Element21; + public readonly GpuUInt128 Element22; + public readonly GpuUInt128 Element23; + public readonly GpuUInt128 Element24; + public readonly GpuUInt128 Element25; + public readonly GpuUInt128 Element26; + public readonly GpuUInt128 Element27; + public readonly GpuUInt128 Element28; + public readonly GpuUInt128 Element29; + public readonly GpuUInt128 Element30; + public readonly GpuUInt128 Element31; + public readonly GpuUInt128 Element32; + public readonly GpuUInt128 Element33; + public readonly GpuUInt128 Element34; + public readonly GpuUInt128 Element35; + public readonly GpuUInt128 Element36; + public readonly GpuUInt128 Element37; + public readonly GpuUInt128 Element38; + public readonly GpuUInt128 Element39; + public readonly GpuUInt128 Element40; + public readonly GpuUInt128 Element41; + public readonly GpuUInt128 Element42; + public readonly GpuUInt128 Element43; + public readonly GpuUInt128 Element44; + public readonly GpuUInt128 Element45; + public readonly GpuUInt128 Element46; + public readonly GpuUInt128 Element47; + public readonly GpuUInt128 Element48; + public readonly GpuUInt128 Element49; + public readonly GpuUInt128 Element50; + public readonly GpuUInt128 Element51; + public readonly GpuUInt128 Element52; + public readonly GpuUInt128 Element53; + public readonly GpuUInt128 Element54; + public readonly GpuUInt128 Element55; + public readonly GpuUInt128 Element56; + public readonly GpuUInt128 Element57; + public readonly GpuUInt128 Element58; + public readonly GpuUInt128 Element59; + public readonly GpuUInt128 Element60; + public readonly GpuUInt128 Element61; + public readonly GpuUInt128 Element62; + public readonly GpuUInt128 Element63; + public readonly GpuUInt128 Element64; + public readonly GpuUInt128 Element65; + public readonly GpuUInt128 Element66; + public readonly GpuUInt128 Element67; + public readonly GpuUInt128 Element68; + public readonly GpuUInt128 Element69; + public readonly GpuUInt128 Element70; + public readonly GpuUInt128 Element71; + public readonly GpuUInt128 Element72; + public readonly GpuUInt128 Element73; + public readonly GpuUInt128 Element74; + public readonly GpuUInt128 Element75; + public readonly GpuUInt128 Element76; + public readonly GpuUInt128 Element77; + public readonly GpuUInt128 Element78; + public readonly GpuUInt128 Element79; + public readonly GpuUInt128 Element80; + public readonly GpuUInt128 Element81; + public readonly GpuUInt128 Element82; + public readonly GpuUInt128 Element83; + public readonly GpuUInt128 Element84; + public readonly GpuUInt128 Element85; + public readonly GpuUInt128 Element86; + public readonly GpuUInt128 Element87; + public readonly GpuUInt128 Element88; + public readonly GpuUInt128 Element89; + public readonly GpuUInt128 Element90; + public readonly GpuUInt128 Element91; + public readonly GpuUInt128 Element92; + public readonly GpuUInt128 Element93; + public readonly GpuUInt128 Element94; + public readonly GpuUInt128 Element95; + public readonly GpuUInt128 Element96; + public readonly GpuUInt128 Element97; + public readonly GpuUInt128 Element98; + public readonly GpuUInt128 Element99; + public readonly GpuUInt128 Element100; + public readonly GpuUInt128 Element101; + public readonly GpuUInt128 Element102; + public readonly GpuUInt128 Element103; + public readonly GpuUInt128 Element104; + public readonly GpuUInt128 Element105; + public readonly GpuUInt128 Element106; + public readonly GpuUInt128 Element107; + public readonly GpuUInt128 Element108; + public readonly GpuUInt128 Element109; + public readonly GpuUInt128 Element110; + public readonly GpuUInt128 Element111; + public readonly GpuUInt128 Element112; + public readonly GpuUInt128 Element113; + public readonly GpuUInt128 Element114; + public readonly GpuUInt128 Element115; + public readonly GpuUInt128 Element116; + public readonly GpuUInt128 Element117; + public readonly GpuUInt128 Element118; + public readonly GpuUInt128 Element119; + public readonly GpuUInt128 Element120; + public readonly GpuUInt128 Element121; + public readonly GpuUInt128 Element122; + public readonly GpuUInt128 Element123; + public readonly GpuUInt128 Element124; + public readonly GpuUInt128 Element125; + public readonly GpuUInt128 Element126; + public readonly GpuUInt128 Element127; + + public readonly GpuUInt128 this[int index] + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get => index switch + { + 0 => Element0, + 1 => Element1, + 2 => Element2, + 3 => Element3, + 4 => Element4, + 5 => Element5, + 6 => Element6, + 7 => Element7, + 8 => Element8, + 9 => Element9, + 10 => Element10, + 11 => Element11, + 12 => Element12, + 13 => Element13, + 14 => Element14, + 15 => Element15, + 16 => Element16, + 17 => Element17, + 18 => Element18, + 19 => Element19, + 20 => Element20, + 21 => Element21, + 22 => Element22, + 23 => Element23, + 24 => Element24, + 25 => Element25, + 26 => Element26, + 27 => Element27, + 28 => Element28, + 29 => Element29, + 30 => Element30, + 31 => Element31, + 32 => Element32, + 33 => Element33, + 34 => Element34, + 35 => Element35, + 36 => Element36, + 37 => Element37, + 38 => Element38, + 39 => Element39, + 40 => Element40, + 41 => Element41, + 42 => Element42, + 43 => Element43, + 44 => Element44, + 45 => Element45, + 46 => Element46, + 47 => Element47, + 48 => Element48, + 49 => Element49, + 50 => Element50, + 51 => Element51, + 52 => Element52, + 53 => Element53, + 54 => Element54, + 55 => Element55, + 56 => Element56, + 57 => Element57, + 58 => Element58, + 59 => Element59, + 60 => Element60, + 61 => Element61, + 62 => Element62, + 63 => Element63, + 64 => Element64, + 65 => Element65, + 66 => Element66, + 67 => Element67, + 68 => Element68, + 69 => Element69, + 70 => Element70, + 71 => Element71, + 72 => Element72, + 73 => Element73, + 74 => Element74, + 75 => Element75, + 76 => Element76, + 77 => Element77, + 78 => Element78, + 79 => Element79, + 80 => Element80, + 81 => Element81, + 82 => Element82, + 83 => Element83, + 84 => Element84, + 85 => Element85, + 86 => Element86, + 87 => Element87, + 88 => Element88, + 89 => Element89, + 90 => Element90, + 91 => Element91, + 92 => Element92, + 93 => Element93, + 94 => Element94, + 95 => Element95, + 96 => Element96, + 97 => Element97, + 98 => Element98, + 99 => Element99, + 100 => Element100, + 101 => Element101, + 102 => Element102, + 103 => Element103, + 104 => Element104, + 105 => Element105, + 106 => Element106, + 107 => Element107, + 108 => Element108, + 109 => Element109, + 110 => Element110, + 111 => Element111, + 112 => Element112, + 113 => Element113, + 114 => Element114, + 115 => Element115, + 116 => Element116, + 117 => Element117, + 118 => Element118, + 119 => Element119, + 120 => Element120, + 121 => Element121, + 122 => Element122, + 123 => Element123, + 124 => Element124, + 125 => Element125, + 126 => Element126, + 127 => Element127, + _ => Element0, + }; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + init + { + switch (index) + { + case 0: + Element0 = value; + return; + case 1: + Element1 = value; + return; + case 2: + Element2 = value; + return; + case 3: + Element3 = value; + return; + case 4: + Element4 = value; + return; + case 5: + Element5 = value; + return; + case 6: + Element6 = value; + return; + case 7: + Element7 = value; + return; + case 8: + Element8 = value; + return; + case 9: + Element9 = value; + return; + case 10: + Element10 = value; + return; + case 11: + Element11 = value; + return; + case 12: + Element12 = value; + return; + case 13: + Element13 = value; + return; + case 14: + Element14 = value; + return; + case 15: + Element15 = value; + return; + case 16: + Element16 = value; + return; + case 17: + Element17 = value; + return; + case 18: + Element18 = value; + return; + case 19: + Element19 = value; + return; + case 20: + Element20 = value; + return; + case 21: + Element21 = value; + return; + case 22: + Element22 = value; + return; + case 23: + Element23 = value; + return; + case 24: + Element24 = value; + return; + case 25: + Element25 = value; + return; + case 26: + Element26 = value; + return; + case 27: + Element27 = value; + return; + case 28: + Element28 = value; + return; + case 29: + Element29 = value; + return; + case 30: + Element30 = value; + return; + case 31: + Element31 = value; + return; + case 32: + Element32 = value; + return; + case 33: + Element33 = value; + return; + case 34: + Element34 = value; + return; + case 35: + Element35 = value; + return; + case 36: + Element36 = value; + return; + case 37: + Element37 = value; + return; + case 38: + Element38 = value; + return; + case 39: + Element39 = value; + return; + case 40: + Element40 = value; + return; + case 41: + Element41 = value; + return; + case 42: + Element42 = value; + return; + case 43: + Element43 = value; + return; + case 44: + Element44 = value; + return; + case 45: + Element45 = value; + return; + case 46: + Element46 = value; + return; + case 47: + Element47 = value; + return; + case 48: + Element48 = value; + return; + case 49: + Element49 = value; + return; + case 50: + Element50 = value; + return; + case 51: + Element51 = value; + return; + case 52: + Element52 = value; + return; + case 53: + Element53 = value; + return; + case 54: + Element54 = value; + return; + case 55: + Element55 = value; + return; + case 56: + Element56 = value; + return; + case 57: + Element57 = value; + return; + case 58: + Element58 = value; + return; + case 59: + Element59 = value; + return; + case 60: + Element60 = value; + return; + case 61: + Element61 = value; + return; + case 62: + Element62 = value; + return; + case 63: + Element63 = value; + return; + case 64: + Element64 = value; + return; + case 65: + Element65 = value; + return; + case 66: + Element66 = value; + return; + case 67: + Element67 = value; + return; + case 68: + Element68 = value; + return; + case 69: + Element69 = value; + return; + case 70: + Element70 = value; + return; + case 71: + Element71 = value; + return; + case 72: + Element72 = value; + return; + case 73: + Element73 = value; + return; + case 74: + Element74 = value; + return; + case 75: + Element75 = value; + return; + case 76: + Element76 = value; + return; + case 77: + Element77 = value; + return; + case 78: + Element78 = value; + return; + case 79: + Element79 = value; + return; + case 80: + Element80 = value; + return; + case 81: + Element81 = value; + return; + case 82: + Element82 = value; + return; + case 83: + Element83 = value; + return; + case 84: + Element84 = value; + return; + case 85: + Element85 = value; + return; + case 86: + Element86 = value; + return; + case 87: + Element87 = value; + return; + case 88: + Element88 = value; + return; + case 89: + Element89 = value; + return; + case 90: + Element90 = value; + return; + case 91: + Element91 = value; + return; + case 92: + Element92 = value; + return; + case 93: + Element93 = value; + return; + case 94: + Element94 = value; + return; + case 95: + Element95 = value; + return; + case 96: + Element96 = value; + return; + case 97: + Element97 = value; + return; + case 98: + Element98 = value; + return; + case 99: + Element99 = value; + return; + case 100: + Element100 = value; + return; + case 101: + Element101 = value; + return; + case 102: + Element102 = value; + return; + case 103: + Element103 = value; + return; + case 104: + Element104 = value; + return; + case 105: + Element105 = value; + return; + case 106: + Element106 = value; + return; + case 107: + Element107 = value; + return; + case 108: + Element108 = value; + return; + case 109: + Element109 = value; + return; + case 110: + Element110 = value; + return; + case 111: + Element111 = value; + return; + case 112: + Element112 = value; + return; + case 113: + Element113 = value; + return; + case 114: + Element114 = value; + return; + case 115: + Element115 = value; + return; + case 116: + Element116 = value; + return; + case 117: + Element117 = value; + return; + case 118: + Element118 = value; + return; + case 119: + Element119 = value; + return; + case 120: + Element120 = value; + return; + case 121: + Element121 = value; + return; + case 122: + Element122 = value; + return; + case 123: + Element123 = value; + return; + case 124: + Element124 = value; + return; + case 125: + Element125 = value; + return; + case 126: + Element126 = value; + return; + case 127: + Element127 = value; + return; + default: + // ILGPU kernels cannot throw exceptions, and callers guarantee the index range. + return; + } + } + } + } + + [MemoryDiagnoser] + [SimpleJob(RuntimeMoniker.Net80)] + public class Pow2OddPowerTableBenchmarks + { + private static readonly GpuUInt128 BaseValue = new GpuUInt128(3, 123456789); + private static readonly GpuUInt128 Modulus = new GpuUInt128(17, 987654321); + private static readonly int OddPowerCount = 128; + + [Benchmark(Baseline = true)] + public Pow2OddPowerTableLoop CreateLoop() => new Pow2OddPowerTableLoop(BaseValue, Modulus, OddPowerCount); + + [Benchmark] + public Pow2OddPowerTableSetter CreateSetter() => new Pow2OddPowerTableSetter(BaseValue, Modulus, OddPowerCount); + + [Benchmark] + public Pow2OddPowerTableDelegate CreateDelegate() => new Pow2OddPowerTableDelegate(BaseValue, Modulus, OddPowerCount); + + [Benchmark] + public Pow2OddPowerTableReturn CreateReturn() => new Pow2OddPowerTableReturn(BaseValue, Modulus, OddPowerCount); + } +} \ No newline at end of file diff --git a/EvenPerfectBitScanner.Benchmarks/Program.cs b/EvenPerfectBitScanner.Benchmarks/Program.cs index 892c277..b4e8c81 100644 --- a/EvenPerfectBitScanner.Benchmarks/Program.cs +++ b/EvenPerfectBitScanner.Benchmarks/Program.cs @@ -1,5 +1,6 @@ using BenchmarkDotNet.Configs; using BenchmarkDotNet.Running; +using PerfectNumbers.Core.Gpu; namespace EvenPerfectBitScanner.Benchmarks; @@ -7,15 +8,16 @@ public static class Program { public static void Main(string[] args) { - var config = DefaultConfig.Instance.WithOptions(ConfigOptions.JoinSummary); - BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly).Run(args, config); + // var config = DefaultConfig.Instance.WithOptions(ConfigOptions.JoinSummary | ConfigOptions.DisableOptimizationsValidator); + // BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly).Run(args, config); - // Pow2MontgomeryModCycleComputationBenchmarks benchmarks = new() - // { - // Scale = Pow2MontgomeryModCycleComputationBenchmarks.InputScale.Large - // }; + Pow2MontgomeryModCycleComputationBenchmarks benchmarks = new() + { + Scale = Pow2MontgomeryModCycleComputationBenchmarks.InputScale.Large + }; - // benchmarks.Setup(); + GpuPrimeWorkLimiter.SetLimit(1); + benchmarks.Setup(); // Console.WriteLine("Press ENTER to run"); // Console.ReadLine(); diff --git a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs index aca4483..3338cd0 100644 --- a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs @@ -25,8 +25,8 @@ public class PrimeOrderGpuHeuristicsTests try { PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0; - - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); remainder.Should().Be(0UL); @@ -40,7 +40,9 @@ public class PrimeOrderGpuHeuristicsTests [Fact] public void TryPow2Mod_returns_success_for_supported_prime() { - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, 193UL, out ulong remainder); + const ulong prime = 193UL; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder, divisorData); status.Should().Be(GpuPow2ModStatus.Success); remainder.Should().Be(2UL); @@ -54,7 +56,8 @@ public class PrimeOrderGpuHeuristicsTests PrimeOrderGpuHeuristics.OverrideCapabilitiesForTesting(capability); try { - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(1UL, prime, out ulong remainder, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); remainder.Should().Be(0UL); @@ -73,12 +76,14 @@ public class PrimeOrderGpuHeuristicsTests PrimeOrderGpuHeuristics.OverrideCapabilitiesForTesting(capability); try { - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(16UL, 193UL, out ulong remainder); + const ulong prime = 193UL; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(16UL, prime, out ulong remainder, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); remainder.Should().Be(0UL); - GpuPow2ModStatus followUp = PrimeOrderGpuHeuristics.TryPow2Mod(4UL, 193UL, out ulong followUpRemainder); + GpuPow2ModStatus followUp = PrimeOrderGpuHeuristics.TryPow2Mod(4UL, prime, out ulong followUpRemainder, divisorData); followUp.Should().Be(GpuPow2ModStatus.Success); followUpRemainder.Should().Be(16UL); } @@ -98,7 +103,8 @@ public class PrimeOrderGpuHeuristicsTests ulong[] remainders = new ulong[2]; ulong[] exponents = { 1UL, 2UL }; - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); remainders[0].Should().Be(0UL); @@ -115,7 +121,9 @@ public class PrimeOrderGpuHeuristicsTests { ulong[] remainders = { 123UL, 456UL }; ulong[] exponents = { 1UL, 2UL }; - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, 257UL, remainders); + const ulong prime = 257UL; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData); status.Should().Be(GpuPow2ModStatus.Success); remainders[0].Should().Be(2UL); @@ -133,7 +141,8 @@ public class PrimeOrderGpuHeuristicsTests ulong[] exponents = { 1UL, 2UL }; ulong[] remainders = { 5UL, 6UL }; - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); remainders[0].Should().Be(0UL); @@ -156,7 +165,9 @@ public class PrimeOrderGpuHeuristicsTests ulong[] exponents = { 2UL, 16UL }; ulong[] remainders = { 9UL, 11UL }; - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, 193UL, remainders); + const ulong prime = 193UL; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData); status.Should().Be(GpuPow2ModStatus.Overflow); remainders[0].Should().Be(0UL); @@ -322,16 +333,22 @@ public class PrimeOrderGpuHeuristicsTests } } - [Fact] - public void TryPow2ModBatch_throws_when_remainder_span_is_too_small() - { - ulong[] remainders = new ulong[1]; - ulong[] exponents = { 1UL, 2UL }; + // This will never happen in production code because the caller ensures the remainder span is large enough. + // [Fact] + // public void TryPow2ModBatch_throws_when_remainder_span_is_too_small() + // { + // ulong[] remainders = new ulong[1]; + // ulong[] exponents = { 1UL, 2UL }; - Action act = () => PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, 5UL, remainders); + // Action act = () => + // { + // const ulong prime = 5UL; + // MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + // PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders, divisorData); + // }; - act.Should().Throw(); - } + // act.Should().Throw(); + // } [Fact] public void Calculate_falls_back_to_cpu_when_gpu_overflow_is_marked() @@ -341,7 +358,8 @@ public class PrimeOrderGpuHeuristicsTests { PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0; - PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(prime, null, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime); + PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(prime, null, divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault); result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found); result.Order.Should().Be(3UL); diff --git a/PerfectNumbers.Core.Tests/MersenneDivisorCyclesGenerateGpuTests.cs b/PerfectNumbers.Core.Tests/MersenneDivisorCyclesGenerateGpuTests.cs index 008f7b0..af63152 100644 --- a/PerfectNumbers.Core.Tests/MersenneDivisorCyclesGenerateGpuTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneDivisorCyclesGenerateGpuTests.cs @@ -38,7 +38,8 @@ public class MersenneDivisorCyclesGenerateGpuTests foreach (var pair in pairs) { expected.Should().Contain(pair.divisor); - pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLength(pair.divisor)); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(pair.divisor); + pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLength(pair.divisor, divisorData)); } } finally @@ -86,7 +87,8 @@ public class MersenneDivisorCyclesGenerateGpuTests foreach (var pair in pairs) { uniqueDivisors.Add(pair.divisor).Should().BeTrue(); - pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLength(pair.divisor)); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(pair.divisor); + pair.cycle.Should().Be(MersenneDivisorCycles.CalculateCycleLength(pair.divisor, divisorData)); } uniqueDivisors.Should().BeEquivalentTo(expected); diff --git a/PerfectNumbers.Core.Tests/MersenneDivisorCyclesTests.cs b/PerfectNumbers.Core.Tests/MersenneDivisorCyclesTests.cs index 05461b4..80c844a 100644 --- a/PerfectNumbers.Core.Tests/MersenneDivisorCyclesTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneDivisorCyclesTests.cs @@ -18,7 +18,8 @@ public class MersenneDivisorCyclesTests [MemberData(nameof(Divisors))] public void CalculateCycleLength_returns_expected_value(ulong divisor, ulong expected) { - MersenneDivisorCycles.CalculateCycleLength(divisor).Should().Be(expected); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData).Should().Be(expected); MersenneDivisorCycles.CalculateCycleLengthGpu(divisor).Should().Be(expected); } } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs index eaf7f71..fdd85ea 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs @@ -35,13 +35,13 @@ public class MersenneNumberDivisorCpuTesterTests ulong[] primes = { 5UL, 7UL, 11UL, 13UL }; byte[] hits = new byte[primes.Length]; - ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL); - session.CheckDivisor(23UL, cycle23, primes, hits); + ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL)); + session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits); hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 }); Array.Fill(hits, (byte)0); - ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL); - session.CheckDivisor(31UL, cycle31, primes, hits); + ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL)); + session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits); hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 }); } } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index 0aa5657..765349d 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -146,13 +146,13 @@ public class MersenneNumberDivisorGpuTesterTests ulong[] primes = { 31UL, 37UL, 41UL, 43UL }; byte[] hits = new byte[primes.Length]; - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL); - session.CheckDivisor(223UL, cycle223, primes, hits); + ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL)); + session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits); hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0 }); Array.Fill(hits, (byte)0); - ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLength(13367UL); - session.CheckDivisor(13367UL, cycle13367, primes, hits); + ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLength(13367UL, MontgomeryDivisorData.FromModulus(13367UL)); + session.CheckDivisor(13367UL, MontgomeryDivisorData.FromModulus(13367UL), cycle13367, primes, hits); hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 }); } @@ -171,8 +171,8 @@ public class MersenneNumberDivisorGpuTesterTests ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL }; byte[] hits = new byte[primes.Length]; - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL); - session.CheckDivisor(223UL, cycle223, primes, hits); + ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL)); + session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits); hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 }); } @@ -192,8 +192,8 @@ public class MersenneNumberDivisorGpuTesterTests ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL }; byte[] hits = new byte[primes.Length]; - ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL); - session.CheckDivisor(223UL, cycle223, primes, hits); + ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL)); + session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits); hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 }); } diff --git a/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs b/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs index da6b4e9..d7dd522 100644 --- a/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs +++ b/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs @@ -48,6 +48,7 @@ public class PrimeOrderCalculatorTests PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate( prime, previousOrder: null, + MontgomeryDivisorData.FromModulus(prime), PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault); result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found); @@ -64,11 +65,13 @@ public class PrimeOrderCalculatorTests PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate( prime, previousOrder, + MontgomeryDivisorData.FromModulus(prime), PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault); PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate( prime, previousOrder: null, + MontgomeryDivisorData.FromModulus(prime), PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault); heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found); @@ -85,6 +88,7 @@ public class PrimeOrderCalculatorTests PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate( prime, previousOrder: null, + MontgomeryDivisorData.FromModulus(prime), PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault); result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found); @@ -104,6 +108,7 @@ public class PrimeOrderCalculatorTests PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate( prime: 13UL, previousOrder: null, + MontgomeryDivisorData.FromModulus(13UL), config); result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved); @@ -120,9 +125,11 @@ public class PrimeOrderCalculatorTests maxPowChecks: 1, mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic); + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL); PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate( prime: 239UL, previousOrder: null, + divisorData, heuristicConfig); heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved); @@ -137,6 +144,7 @@ public class PrimeOrderCalculatorTests PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate( prime: 239UL, previousOrder: null, + divisorData, strictConfig); strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found); diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 8cec408..5921ecf 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -194,7 +194,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv processedCount++; lastProcessed = candidate; - if (MersenneDivisorCycles.CycleEqualsExponentForMersenneCandidate(candidate, prime)) + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); + if (MersenneDivisorCycles.CycleEqualsExponentForMersenneCandidate(candidate, divisorData, prime)) { processedAll = true; return true; @@ -272,7 +273,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv _disposed = false; } - public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits) + public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan primes, Span hits) { if (_disposed) { @@ -288,7 +289,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv for (int i = 0; i < length; i++) { ulong prime = primes[i]; - bool divides = MersenneDivisorCycles.CycleEqualsExponentForMersenneCandidate(divisor, prime); + bool divides = MersenneDivisorCycles.CycleEqualsExponentForMersenneCandidate(divisor, divisorData, prime); hits[i] = divides ? (byte)1 : (byte)0; } } diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index 0da555a..b01e5e2 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -109,14 +109,16 @@ public sealed class DivisorCycleCache ArrayPool.Shared.Return(rentedMissing, clearArray: false); } } - + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void ComputeCyclesCpu(ReadOnlySpan divisors, Span cycles, ReadOnlySpan indices) { for (int i = 0; i < indices.Length; i++) { int targetIndex = indices[i]; - cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLength(divisors[targetIndex]); + ulong divisor = divisors[targetIndex]; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + cycles[targetIndex] = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData); } } diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index e96e290..d39cf06 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -25,16 +25,16 @@ public readonly struct ResidueAutomatonArgs public sealed class KernelContainer { - // Serializes first-time initialization of kernels/buffers per accelerator. - public Action, ArrayView>? Order; - public Action, ArrayView1D>? Incremental; + // Serializes first-time initialization of kernels/buffers per accelerator. + public Action, ArrayView>? Order; + public Action, ArrayView1D>? Incremental; public Action, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>? Pow2Mod; - public Action, ArrayView1D>? IncrementalOrder; + public Action, ArrayView1D>? IncrementalOrder; public Action, ArrayView1D>? Pow2ModOrder; @@ -45,22 +45,22 @@ public sealed class KernelContainer public MemoryBuffer1D? SmallPrimesLastSeven; public MemoryBuffer1D? SmallPrimesPow2LastSeven; - public static T InitOnce(ref T? slot, Func factory) where T : class + public static T InitOnce(ref T? slot, Func factory) where T : class + { + var current = Volatile.Read(ref slot); + if (current is not null) { - var current = Volatile.Read(ref slot); - if (current is not null) - { - return current; - } - - current = Volatile.Read(ref slot); - if (current is null) - { - current = factory(); - Volatile.Write(ref slot, current); - } - return current; + return current; } + + current = Volatile.Read(ref slot); + if (current is null) + { + current = factory(); + Volatile.Write(ref slot, current); + } + return current; + } } public readonly struct ResiduePrimeViews( @@ -77,122 +77,120 @@ public readonly struct ResiduePrimeViews( public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels) { - private bool disposedValue; - private IDisposable? _limiter = limiter; - private GpuContextLease? _gpu = gpu; - private KernelContainer _kernels = kernels; - private AcceleratorStream? _stream; - private object? _executionLock = gpu.ExecutionLock; + private bool disposedValue; + private IDisposable? _limiter = limiter; + private GpuContextLease? _gpu = gpu; + private KernelContainer _kernels = kernels; + private AcceleratorStream? _stream; + private object? _executionLock = gpu.ExecutionLock; - public readonly Accelerator Accelerator => _gpu?.Accelerator ?? throw new NullReferenceException("GPU context is already released"); + public readonly Accelerator Accelerator => _gpu?.Accelerator ?? throw new NullReferenceException("GPU context is already released"); - public AcceleratorStream Stream + public AcceleratorStream Stream + { + get { - get + if (_stream is { } stream) + { + return stream; + } + + // It feels that we can combine these into 1 return statement, but that increases the max stack frame size to 3 from 2 + stream = Accelerator.CreateStream(); + _stream = stream; + return stream; + } + } + + public readonly ExecutionScope EnterExecutionScope() => new ExecutionScope(_executionLock); + + public readonly Action, ArrayView> OrderKernel + { + get + { + var accel = Accelerator; // avoid capturing 'this' in lambda + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.Order, () => { - if (_stream is { } stream) - { - return stream; - } - - stream = Accelerator.CreateStream(); - _stream = stream; - return stream; - } - } - - public readonly ExecutionScope EnterExecutionScope() - { - return new ExecutionScope(_executionLock); - } - - public readonly Action, ArrayView> OrderKernel - { - get - { - var accel = Accelerator; // avoid capturing 'this' in lambda - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.Order, () => - { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView>(OrderKernelScan); - var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, ArrayView>>(); - }); - } - } + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView>(OrderKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView>>(); + }); + } } + } public readonly Action, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> Pow2ModKernel + { + get { - get + var accel = Accelerator; + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () => { - var accel = Accelerator; - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () => - { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(Pow2ModKernelScan); - var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); - }); - } - } + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(Pow2ModKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); + }); + } } + } public readonly Action, ArrayView1D> IncrementalKernel + { + get { - get + var accel = Accelerator; + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.Incremental, () => { - var accel = Accelerator; - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.Incremental, () => - { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalKernelScan); - var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, ArrayView1D>>(); - }); - } - } + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView1D>>(); + }); + } } + } public readonly Action, ArrayView1D> IncrementalOrderKernel + { + get { - get + var accel = Accelerator; + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.IncrementalOrder, () => { - var accel = Accelerator; - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.IncrementalOrder, () => - { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalOrderKernelScan); - var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, ArrayView1D>>(); - }); - } - } + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(IncrementalOrderKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView1D>>(); + }); + } } + } public readonly Action, ArrayView1D> Pow2ModOrderKernel + { + get { - get + var accel = Accelerator; + lock (_gpu!.Value.KernelInitLock) + { + return KernelContainer.InitOnce(ref _kernels.Pow2ModOrder, () => { - var accel = Accelerator; - lock (_gpu!.Value.KernelInitLock) - { - return KernelContainer.InitOnce(ref _kernels.Pow2ModOrder, () => - { - var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(Pow2ModOrderKernelScan); - var kernel = KernelUtil.GetKernel(loaded); - return kernel.CreateLauncherDelegate, ArrayView1D>>(); - }); - } - } + var loaded = accel.LoadAutoGroupedStreamKernel, ArrayView1D>(Pow2ModOrderKernelScan); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, ArrayView1D>>(); + }); + } } + } - // TODO: Plumb the small-cycles device buffer into all kernels that can benefit - // (some already accept it). Consider a compact type (byte/ushort) for memory footprint. + // TODO: Plumb the small-cycles device buffer into all kernels that can benefit + // (some already accept it). Consider a compact type (byte/ushort) for memory footprint. private static void IncrementalKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong divMul, ulong q0m10, ulong q0m8, ulong q0m3, ulong q0m5, ArrayView orders, @@ -206,21 +204,21 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne // residue automaton ulong step10 = (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 - if (shouldCheck) - { + r10 -= (r10 >= 10UL) ? 10UL : 0UL; + bool shouldCheck = r10 != 5UL; // skip q ≡ 5 (mod 10); other residues handled by r8/r3/r5 + if (shouldCheck) + { ulong step8 = ((exponent & 7UL) << 1) & 7UL; ulong r8 = (q0m8 + ((step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } else { - // TODO: Swap these `%` filters to the shared Mod helpers so the residue automaton matches the benchmarked - // bitmask-based implementation for GPU workloads. - ulong step3 = ((exponent % 3UL) << 1) % 3UL; + // TODO: Swap these `%` filters to the shared Mod helpers so the residue automaton matches the benchmarked + // bitmask-based implementation for GPU workloads. + ulong step3 = ((exponent % 3UL) << 1) % 3UL; ulong r3 = q0m3 + (step3 * idxMod3); r3 -= (r3 >= 3UL) ? 3UL : 0UL; if (r3 == 0UL) { @@ -236,19 +234,20 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne } } } - } + } - if (!shouldCheck) - { - orders[index] = 0UL; - return; - } + if (!shouldCheck) + { + orders[index] = 0UL; + return; + } + // TODO: Is this expected that the value of kStart is modified? GpuUInt128 k = kStart + (GpuUInt128)idx; GpuUInt128 q = twoP; q.Mul64(k); q.Add(GpuUInt128.One); - + // Small-cycles in-kernel early rejection from device table if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) { @@ -260,7 +259,9 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne } } + // TODO: q will be lowered by GpuUInt128.One here. Is this expected? Maybe that's why we have issues with the incremental kernel? Review other places where we use GpuUInt128 operators. GpuUInt128 phi = q - GpuUInt128.One; + // TODO: Do we have the check here because we don't support high values or this is expected result? if (phi.High != 0UL) { orders[index] = 0UL; @@ -332,29 +333,29 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne orders[index] = exponent; } - private void Dispose(bool disposing) + private void Dispose(bool disposing) + { + if (!disposedValue) { - if (!disposedValue) - { - if (disposing) - { - _stream?.Dispose(); - _stream = null; + if (disposing) + { + _stream?.Dispose(); + _stream = null; - _gpu?.Dispose(); - _gpu = null; + _gpu?.Dispose(); + _gpu = null; - _executionLock = null; + _executionLock = null; - _limiter?.Dispose(); - _limiter = null; - } + _limiter?.Dispose(); + _limiter = null; + } - // free unmanaged resources (unmanaged objects) and override finalizer - // set large fields to null - disposedValue = true; - } - } + // free unmanaged resources (unmanaged objects) and override finalizer + // set large fields to null + disposedValue = true; + } + } private static void Pow2ModKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _, ResidueAutomatonArgs ra, ArrayView orders, @@ -404,13 +405,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne return; } - GpuUInt128 k = kStart + (GpuUInt128)idx; + // TODO: kStart is modified after this. Is this expected? + kStart.Add(idx); GpuUInt128 q = twoP; - q.Mul64(k); + q.Mul64(kStart); q.Add(GpuUInt128.One); if (q.High == 0UL && q.Low < (ulong)smallCycles.Length) { ulong cycle = smallCycles[(int)q.Low]; + // cycle should always be initialized if we're within array limit in production code if (cycle != 0UL && cycle <= exponent && (exponent % cycle) != 0UL) { orders[index] = 0UL; @@ -517,40 +520,40 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne ulong idxMod3 = idx % 3UL; ulong idxMod5 = idx % 5UL; ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; - bool shouldCheck = r10 != 5UL; - if (shouldCheck) - { - ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } - else - { - ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); - if (r3 >= 6UL) r3 -= 6UL; - if (r3 >= 3UL) r3 -= 3UL; - if (r3 == 0UL) - { - shouldCheck = false; - } - else - { - ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); - if (r5 >= 15UL) r5 -= 15UL; - if (r5 >= 10UL) r5 -= 10UL; - if (r5 >= 5UL) r5 -= 5UL; - if (r5 == 0UL) - { - shouldCheck = false; - } - } - } - } - if (!shouldCheck) - { - return; - } + bool shouldCheck = r10 != 5UL; + if (shouldCheck) + { + ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } + else + { + ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); + if (r3 >= 6UL) r3 -= 6UL; + if (r3 >= 3UL) r3 -= 3UL; + if (r3 == 0UL) + { + shouldCheck = false; + } + else + { + ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); + if (r5 >= 15UL) r5 -= 15UL; + if (r5 >= 10UL) r5 -= 10UL; + if (r5 >= 5UL) r5 -= 5UL; + if (r5 == 0UL) + { + shouldCheck = false; + } + } + } + } + if (!shouldCheck) + { + return; + } GpuUInt128 k = kStart + (GpuUInt128)idx; GpuUInt128 q = twoP; @@ -603,40 +606,40 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne ulong idxMod3 = idx % 3UL; ulong idxMod5 = idx % 5UL; ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; - bool shouldCheck = r10 != 5UL; - if (shouldCheck) - { - ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; - if (r8 != 1UL && r8 != 7UL) - { - shouldCheck = false; - } - else - { - ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); - if (r3 >= 6UL) r3 -= 6UL; - if (r3 >= 3UL) r3 -= 3UL; - if (r3 == 0UL) - { - shouldCheck = false; - } - else - { - ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); - if (r5 >= 15UL) r5 -= 15UL; - if (r5 >= 10UL) r5 -= 10UL; - if (r5 >= 5UL) r5 -= 5UL; - if (r5 == 0UL) - { - shouldCheck = false; - } - } - } - } - if (!shouldCheck) - { - return; - } + bool shouldCheck = r10 != 5UL; + if (shouldCheck) + { + ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL; + if (r8 != 1UL && r8 != 7UL) + { + shouldCheck = false; + } + else + { + ulong r3 = ra.Q0M3 + (ra.Step3 * idxMod3); + if (r3 >= 6UL) r3 -= 6UL; + if (r3 >= 3UL) r3 -= 3UL; + if (r3 == 0UL) + { + shouldCheck = false; + } + else + { + ulong r5 = ra.Q0M5 + (ra.Step5 * idxMod5); + if (r5 >= 15UL) r5 -= 15UL; + if (r5 >= 10UL) r5 -= 10UL; + if (r5 >= 5UL) r5 -= 5UL; + if (r5 == 0UL) + { + shouldCheck = false; + } + } + } + } + if (!shouldCheck) + { + return; + } GpuUInt128 k = kStart + (GpuUInt128)idx; GpuUInt128 q = twoP; @@ -661,35 +664,35 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne Atomic.Or(ref found[0], 1); } + public void Dispose() + { + // Do not change this code. Put cleanup code in 'Dispose(bool disposing)' method + Dispose(disposing: true); + // Enabled this when the struct is converted to a class + // GC.SuppressFinalize(this); + } + + public readonly struct ExecutionScope : IDisposable + { + private readonly object? _lock; + + public ExecutionScope(object? sync) + { + _lock = sync; + if (_lock is not null) + { + Monitor.Enter(_lock); + } + } + public void Dispose() { - // Do not change this code. Put cleanup code in 'Dispose(bool disposing)' method - Dispose(disposing: true); - // Enabled this when the struct is converted to a class - // GC.SuppressFinalize(this); - } - - public readonly struct ExecutionScope : IDisposable - { - private readonly object? _lock; - - public ExecutionScope(object? sync) - { - _lock = sync; - if (_lock is not null) - { - Monitor.Enter(_lock); - } - } - - public void Dispose() - { - if (_lock is not null) - { - Monitor.Exit(_lock); - } - } + if (_lock is not null) + { + Monitor.Exit(_lock); + } } + } } public class GpuKernelPool @@ -771,19 +774,19 @@ public class GpuKernelPool } } - public static GpuKernelLease GetKernel(bool useGpuOrder) - { - var limiter = GpuPrimeWorkLimiter.Acquire(); - var gpu = RentPreferred(preferCpu: !useGpuOrder); - var accelerator = gpu.Accelerator; - var kernels = GetKernels(accelerator); - return new(limiter, gpu, kernels); - } + public static GpuKernelLease GetKernel(bool useGpuOrder) + { + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = RentPreferred(preferCpu: !useGpuOrder); + var accelerator = gpu.Accelerator; + var kernels = GetKernels(accelerator); + return new(limiter, gpu, kernels); + } - /// - /// Runs a GPU action with an acquired accelerator and stream. - /// - /// Action to run with (Accelerator, Stream). + /// + /// Runs a GPU action with an acquired accelerator and stream. + /// + /// Action to run with (Accelerator, Stream). public static void Run(Action action) { var lease = GetKernel(useGpuOrder: true); diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs index c47b0bd..9bb712b 100644 --- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs +++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs @@ -8,7 +8,6 @@ namespace PerfectNumbers.Core.Gpu; public struct GpuUInt128 : IComparable, IEquatable { public ulong High; - public ulong Low; private const int NativeModuloChunkBits = 8; @@ -21,690 +20,8 @@ public struct GpuUInt128 : IComparable, IEquatable private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1); private const int Pow2WindowMaxExponent = (1 << Pow2WindowSizeBits) - 1; - private struct Pow2OddPowerTable - { - public GpuUInt128 Element0; - public GpuUInt128 Element1; - public GpuUInt128 Element2; - public GpuUInt128 Element3; - public GpuUInt128 Element4; - public GpuUInt128 Element5; - public GpuUInt128 Element6; - public GpuUInt128 Element7; - public GpuUInt128 Element8; - public GpuUInt128 Element9; - public GpuUInt128 Element10; - public GpuUInt128 Element11; - public GpuUInt128 Element12; - public GpuUInt128 Element13; - public GpuUInt128 Element14; - public GpuUInt128 Element15; - public GpuUInt128 Element16; - public GpuUInt128 Element17; - public GpuUInt128 Element18; - public GpuUInt128 Element19; - public GpuUInt128 Element20; - public GpuUInt128 Element21; - public GpuUInt128 Element22; - public GpuUInt128 Element23; - public GpuUInt128 Element24; - public GpuUInt128 Element25; - public GpuUInt128 Element26; - public GpuUInt128 Element27; - public GpuUInt128 Element28; - public GpuUInt128 Element29; - public GpuUInt128 Element30; - public GpuUInt128 Element31; - public GpuUInt128 Element32; - public GpuUInt128 Element33; - public GpuUInt128 Element34; - public GpuUInt128 Element35; - public GpuUInt128 Element36; - public GpuUInt128 Element37; - public GpuUInt128 Element38; - public GpuUInt128 Element39; - public GpuUInt128 Element40; - public GpuUInt128 Element41; - public GpuUInt128 Element42; - public GpuUInt128 Element43; - public GpuUInt128 Element44; - public GpuUInt128 Element45; - public GpuUInt128 Element46; - public GpuUInt128 Element47; - public GpuUInt128 Element48; - public GpuUInt128 Element49; - public GpuUInt128 Element50; - public GpuUInt128 Element51; - public GpuUInt128 Element52; - public GpuUInt128 Element53; - public GpuUInt128 Element54; - public GpuUInt128 Element55; - public GpuUInt128 Element56; - public GpuUInt128 Element57; - public GpuUInt128 Element58; - public GpuUInt128 Element59; - public GpuUInt128 Element60; - public GpuUInt128 Element61; - public GpuUInt128 Element62; - public GpuUInt128 Element63; - public GpuUInt128 Element64; - public GpuUInt128 Element65; - public GpuUInt128 Element66; - public GpuUInt128 Element67; - public GpuUInt128 Element68; - public GpuUInt128 Element69; - public GpuUInt128 Element70; - public GpuUInt128 Element71; - public GpuUInt128 Element72; - public GpuUInt128 Element73; - public GpuUInt128 Element74; - public GpuUInt128 Element75; - public GpuUInt128 Element76; - public GpuUInt128 Element77; - public GpuUInt128 Element78; - public GpuUInt128 Element79; - public GpuUInt128 Element80; - public GpuUInt128 Element81; - public GpuUInt128 Element82; - public GpuUInt128 Element83; - public GpuUInt128 Element84; - public GpuUInt128 Element85; - public GpuUInt128 Element86; - public GpuUInt128 Element87; - public GpuUInt128 Element88; - public GpuUInt128 Element89; - public GpuUInt128 Element90; - public GpuUInt128 Element91; - public GpuUInt128 Element92; - public GpuUInt128 Element93; - public GpuUInt128 Element94; - public GpuUInt128 Element95; - public GpuUInt128 Element96; - public GpuUInt128 Element97; - public GpuUInt128 Element98; - public GpuUInt128 Element99; - public GpuUInt128 Element100; - public GpuUInt128 Element101; - public GpuUInt128 Element102; - public GpuUInt128 Element103; - public GpuUInt128 Element104; - public GpuUInt128 Element105; - public GpuUInt128 Element106; - public GpuUInt128 Element107; - public GpuUInt128 Element108; - public GpuUInt128 Element109; - public GpuUInt128 Element110; - public GpuUInt128 Element111; - public GpuUInt128 Element112; - public GpuUInt128 Element113; - public GpuUInt128 Element114; - public GpuUInt128 Element115; - public GpuUInt128 Element116; - public GpuUInt128 Element117; - public GpuUInt128 Element118; - public GpuUInt128 Element119; - public GpuUInt128 Element120; - public GpuUInt128 Element121; - public GpuUInt128 Element122; - public GpuUInt128 Element123; - public GpuUInt128 Element124; - public GpuUInt128 Element125; - public GpuUInt128 Element126; - public GpuUInt128 Element127; - - public GpuUInt128 this[int index] - { - readonly get - { - return GetElement(index); - } - - set - { - SetElement(index, value); - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private readonly GpuUInt128 GetElement(int index) => index switch - { - 0 => Element0, - 1 => Element1, - 2 => Element2, - 3 => Element3, - 4 => Element4, - 5 => Element5, - 6 => Element6, - 7 => Element7, - 8 => Element8, - 9 => Element9, - 10 => Element10, - 11 => Element11, - 12 => Element12, - 13 => Element13, - 14 => Element14, - 15 => Element15, - 16 => Element16, - 17 => Element17, - 18 => Element18, - 19 => Element19, - 20 => Element20, - 21 => Element21, - 22 => Element22, - 23 => Element23, - 24 => Element24, - 25 => Element25, - 26 => Element26, - 27 => Element27, - 28 => Element28, - 29 => Element29, - 30 => Element30, - 31 => Element31, - 32 => Element32, - 33 => Element33, - 34 => Element34, - 35 => Element35, - 36 => Element36, - 37 => Element37, - 38 => Element38, - 39 => Element39, - 40 => Element40, - 41 => Element41, - 42 => Element42, - 43 => Element43, - 44 => Element44, - 45 => Element45, - 46 => Element46, - 47 => Element47, - 48 => Element48, - 49 => Element49, - 50 => Element50, - 51 => Element51, - 52 => Element52, - 53 => Element53, - 54 => Element54, - 55 => Element55, - 56 => Element56, - 57 => Element57, - 58 => Element58, - 59 => Element59, - 60 => Element60, - 61 => Element61, - 62 => Element62, - 63 => Element63, - 64 => Element64, - 65 => Element65, - 66 => Element66, - 67 => Element67, - 68 => Element68, - 69 => Element69, - 70 => Element70, - 71 => Element71, - 72 => Element72, - 73 => Element73, - 74 => Element74, - 75 => Element75, - 76 => Element76, - 77 => Element77, - 78 => Element78, - 79 => Element79, - 80 => Element80, - 81 => Element81, - 82 => Element82, - 83 => Element83, - 84 => Element84, - 85 => Element85, - 86 => Element86, - 87 => Element87, - 88 => Element88, - 89 => Element89, - 90 => Element90, - 91 => Element91, - 92 => Element92, - 93 => Element93, - 94 => Element94, - 95 => Element95, - 96 => Element96, - 97 => Element97, - 98 => Element98, - 99 => Element99, - 100 => Element100, - 101 => Element101, - 102 => Element102, - 103 => Element103, - 104 => Element104, - 105 => Element105, - 106 => Element106, - 107 => Element107, - 108 => Element108, - 109 => Element109, - 110 => Element110, - 111 => Element111, - 112 => Element112, - 113 => Element113, - 114 => Element114, - 115 => Element115, - 116 => Element116, - 117 => Element117, - 118 => Element118, - 119 => Element119, - 120 => Element120, - 121 => Element121, - 122 => Element122, - 123 => Element123, - 124 => Element124, - 125 => Element125, - 126 => Element126, - 127 => Element127, - _ => Element0, - }; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private void SetElement(int index, GpuUInt128 value) - { - switch (index) - { - case 0: - Element0 = value; - return; - case 1: - Element1 = value; - return; - case 2: - Element2 = value; - return; - case 3: - Element3 = value; - return; - case 4: - Element4 = value; - return; - case 5: - Element5 = value; - return; - case 6: - Element6 = value; - return; - case 7: - Element7 = value; - return; - case 8: - Element8 = value; - return; - case 9: - Element9 = value; - return; - case 10: - Element10 = value; - return; - case 11: - Element11 = value; - return; - case 12: - Element12 = value; - return; - case 13: - Element13 = value; - return; - case 14: - Element14 = value; - return; - case 15: - Element15 = value; - return; - case 16: - Element16 = value; - return; - case 17: - Element17 = value; - return; - case 18: - Element18 = value; - return; - case 19: - Element19 = value; - return; - case 20: - Element20 = value; - return; - case 21: - Element21 = value; - return; - case 22: - Element22 = value; - return; - case 23: - Element23 = value; - return; - case 24: - Element24 = value; - return; - case 25: - Element25 = value; - return; - case 26: - Element26 = value; - return; - case 27: - Element27 = value; - return; - case 28: - Element28 = value; - return; - case 29: - Element29 = value; - return; - case 30: - Element30 = value; - return; - case 31: - Element31 = value; - return; - case 32: - Element32 = value; - return; - case 33: - Element33 = value; - return; - case 34: - Element34 = value; - return; - case 35: - Element35 = value; - return; - case 36: - Element36 = value; - return; - case 37: - Element37 = value; - return; - case 38: - Element38 = value; - return; - case 39: - Element39 = value; - return; - case 40: - Element40 = value; - return; - case 41: - Element41 = value; - return; - case 42: - Element42 = value; - return; - case 43: - Element43 = value; - return; - case 44: - Element44 = value; - return; - case 45: - Element45 = value; - return; - case 46: - Element46 = value; - return; - case 47: - Element47 = value; - return; - case 48: - Element48 = value; - return; - case 49: - Element49 = value; - return; - case 50: - Element50 = value; - return; - case 51: - Element51 = value; - return; - case 52: - Element52 = value; - return; - case 53: - Element53 = value; - return; - case 54: - Element54 = value; - return; - case 55: - Element55 = value; - return; - case 56: - Element56 = value; - return; - case 57: - Element57 = value; - return; - case 58: - Element58 = value; - return; - case 59: - Element59 = value; - return; - case 60: - Element60 = value; - return; - case 61: - Element61 = value; - return; - case 62: - Element62 = value; - return; - case 63: - Element63 = value; - return; - case 64: - Element64 = value; - return; - case 65: - Element65 = value; - return; - case 66: - Element66 = value; - return; - case 67: - Element67 = value; - return; - case 68: - Element68 = value; - return; - case 69: - Element69 = value; - return; - case 70: - Element70 = value; - return; - case 71: - Element71 = value; - return; - case 72: - Element72 = value; - return; - case 73: - Element73 = value; - return; - case 74: - Element74 = value; - return; - case 75: - Element75 = value; - return; - case 76: - Element76 = value; - return; - case 77: - Element77 = value; - return; - case 78: - Element78 = value; - return; - case 79: - Element79 = value; - return; - case 80: - Element80 = value; - return; - case 81: - Element81 = value; - return; - case 82: - Element82 = value; - return; - case 83: - Element83 = value; - return; - case 84: - Element84 = value; - return; - case 85: - Element85 = value; - return; - case 86: - Element86 = value; - return; - case 87: - Element87 = value; - return; - case 88: - Element88 = value; - return; - case 89: - Element89 = value; - return; - case 90: - Element90 = value; - return; - case 91: - Element91 = value; - return; - case 92: - Element92 = value; - return; - case 93: - Element93 = value; - return; - case 94: - Element94 = value; - return; - case 95: - Element95 = value; - return; - case 96: - Element96 = value; - return; - case 97: - Element97 = value; - return; - case 98: - Element98 = value; - return; - case 99: - Element99 = value; - return; - case 100: - Element100 = value; - return; - case 101: - Element101 = value; - return; - case 102: - Element102 = value; - return; - case 103: - Element103 = value; - return; - case 104: - Element104 = value; - return; - case 105: - Element105 = value; - return; - case 106: - Element106 = value; - return; - case 107: - Element107 = value; - return; - case 108: - Element108 = value; - return; - case 109: - Element109 = value; - return; - case 110: - Element110 = value; - return; - case 111: - Element111 = value; - return; - case 112: - Element112 = value; - return; - case 113: - Element113 = value; - return; - case 114: - Element114 = value; - return; - case 115: - Element115 = value; - return; - case 116: - Element116 = value; - return; - case 117: - Element117 = value; - return; - case 118: - Element118 = value; - return; - case 119: - Element119 = value; - return; - case 120: - Element120 = value; - return; - case 121: - Element121 = value; - return; - case 122: - Element122 = value; - return; - case 123: - Element123 = value; - return; - case 124: - Element124 = value; - return; - case 125: - Element125 = value; - return; - case 126: - Element126 = value; - return; - case 127: - Element127 = value; - return; - default: - // ILGPU kernels cannot throw exceptions, and callers guarantee the index range. - return; - } - } - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static int GetScalarBitLength(ulong value) - { - if (value == 0UL) - { - return 0; - } - - return 64 - XMath.LeadingZeroCount(value); - } + private static int GetScalarBitLength(ulong value) => value == 0UL ? 0 : 64 - XMath.LeadingZeroCount(value); public static readonly GpuUInt128 Zero = new(); public static readonly GpuUInt128 One = new(1UL); @@ -743,23 +60,21 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public readonly int GetBitLength() { - if (High != 0UL) + ulong temp = High; + if (temp != 0UL) { - return 64 + (64 - BitOperations.LeadingZeroCount(High)); + return 64 + (64 - BitOperations.LeadingZeroCount(temp)); } - ulong low = Low; - if (low == 0UL) + temp = Low; + if (temp == 0UL) { return 0; } - return 64 - BitOperations.LeadingZeroCount(low); + return 64 - BitOperations.LeadingZeroCount(temp); } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static int GetBitLength(GpuUInt128 value) => value.GetBitLength(); - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static implicit operator GpuUInt128(UInt128 value) => new(value); @@ -770,6 +85,7 @@ public struct GpuUInt128 : IComparable, IEquatable public static implicit operator UInt128(GpuUInt128 value) => ((UInt128)value.High << 64) | value.Low; + // TODO: The operators modifying the state of left is intuitive and can easily lead to bugs. Change this back to return a new instance, but replace where possible with method calls, if we want to modify the state of the instance. [MethodImpl(MethodImplOptions.AggressiveInlining)] public static GpuUInt128 operator +(GpuUInt128 left, GpuUInt128 right) { @@ -809,16 +125,17 @@ public struct GpuUInt128 : IComparable, IEquatable ulong high = value.High; ulong low = value.Low; + int shiftedValue; if (shift >= 64) { - int longShift = shift - 64; - value.High = low << longShift; + shiftedValue = shift - 64; + value.High = low << shiftedValue; value.Low = 0UL; return value; } - int inverseShift = 64 - shift; - value.High = (high << shift) | (low >> inverseShift); + shiftedValue = 64 - shift; + value.High = (high << shift) | (low >> shiftedValue); value.Low = low << shift; return value; } @@ -835,16 +152,17 @@ public struct GpuUInt128 : IComparable, IEquatable ulong high = value.High; ulong low = value.Low; + int shiftedValue; if (shift >= 64) { - int longShift = shift - 64; - value.Low = longShift == 0 ? high : high >> longShift; + shiftedValue = shift - 64; + value.Low = shiftedValue == 0 ? high : high >> shiftedValue; value.High = 0UL; return value; } - int inverseShift = 64 - shift; - value.Low = (low >> shift) | (high << inverseShift); + shiftedValue = 64 - shift; + value.Low = (low >> shift) | (high << shiftedValue); value.High = high >> shift; return value; } @@ -913,17 +231,20 @@ public struct GpuUInt128 : IComparable, IEquatable public void Add(GpuUInt128 other) { ulong originalLow = Low; + // We're reusing originalLow for carry to avoid an extra local. ulong low = originalLow + other.Low; - ulong carry = low < originalLow ? 1UL : 0UL; - High = High + other.High + carry; + originalLow = low < originalLow ? 1UL : 0UL; Low = low; + // We're reusing low for the final sum to avoid an extra local and to keep max stack frame at 3 instead of 4 with one statement. + low = High + other.High; + High = low + originalLow; } [MethodImpl(MethodImplOptions.AggressiveInlining)] public void Add(ulong value) { ulong low = Low + value; - High += (low < Low ? 1UL : 0UL); + High += low < Low ? 1UL : 0UL; Low = low; } @@ -962,10 +283,11 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public static GpuUInt128 Pow2Minus1Mod(ulong exponent, GpuUInt128 modulus) { - if (modulus.IsZero) - { - return Zero; - } + // This should never happen in production code + // if (modulus.IsZero) + // { + // return Zero; + // } GpuUInt128 pow = Pow2Mod(exponent, modulus); if (pow.IsZero) @@ -1034,30 +356,34 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MulHigh(ulong x, ulong y) { + // TODO: Review if we can reuse any other existing variables to reduce registry pressure. ulong xLow = (uint)x; ulong yLow = (uint)y; ulong lowProduct = xLow * yLow; ulong yHigh = y >> 32; - ulong cross = xLow * yHigh; + // We're reusing xLow for the cross product to avoid an extra local. + xLow *= yHigh; ulong xHigh = x >> 32; - ulong result = xHigh * yHigh; - ulong temp = xHigh * yLow; + // We're reusing yHigh for the result to void an extra local. + yHigh *= xHigh; + // We're reusing yLow for the other cross product to avoid an extra local. + yLow *= xHigh; // Keeping the wide partial sum in a dedicated local prevents the JIT from // materialising it on the stack before the final carry propagation. The // additional store looks redundant in C#, but it shortens the generated // instruction sequence by avoiding an extra temporary and results in a // measurable throughput win in the MulHigh benchmarks. - result += cross >> 32; - result += temp >> 32; + yHigh += (xLow >> 32); + // We've more instructions with this but the max method stack frame is 3 instead of 4. + yHigh += (yLow >> 32); lowProduct >>= 32; - lowProduct += (uint)cross; - lowProduct += (uint)temp; - result += lowProduct >> 32; - return result; + lowProduct += (uint)xLow + (uint)yLow; + yHigh += lowProduct >> 32; + return yHigh; } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -1099,23 +425,25 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public static GpuUInt128 Pow2Mod(ulong exponent, GpuUInt128 modulus) { - if (modulus.IsZero || modulus == One) - { - return new(); - } + // This should never happen in production code. + // if (modulus.IsZero || modulus == One) + // { + // return new(); + // } - if (exponent == 0UL) - { - return new GpuUInt128(1UL); - } + // This should never happen in production code. + // if (exponent == 0UL) + // { + // return new GpuUInt128(1UL); + // } - Pow2OddPowerTable oddPowers = default; - GpuUInt128 power = new GpuUInt128(1UL); + GpuUInt128[] oddPowers = new GpuUInt128[PerfectNumberConstants.MaxOddPowersCount]; + GpuUInt128 power = new(1UL); int oddIndex = 0; for (int bit = 1; bit <= Pow2WindowMaxExponent; bit++) { - power += power; + power.Add(power); if (power.CompareTo(modulus) >= 0) { power.Sub(modulus); @@ -1128,7 +456,7 @@ public struct GpuUInt128 : IComparable, IEquatable } } - GpuUInt128 result = new GpuUInt128(1UL); + GpuUInt128 result = GpuUInt128.One; int bitLength = GetScalarBitLength(exponent); int index = bitLength - 1; @@ -2144,14 +1472,14 @@ public struct GpuUInt128 : IComparable, IEquatable } } - // TODO: Check if the TODO below is still relevant. +// TODO: Check if the TODO below is still relevant. - // TODO(MOD-OPT): Montgomery/Barrett integration plan - // - Introduce caches of modulus-dependent constants: - // * Montgomery: n' (-(n^{-1}) mod 2^64 or 2^128), R2 = (R^2 mod n) - // * Barrett: mu = floor(2^k / n) for k ∈ {128, 192, 256} - // - Add fast-path for 64-bit NTT primes (modulus.High == 0UL) using pure 64-bit ops. - // - Expose helpers to retrieve/calc constants once per modulus and reuse in kernels. - // - Wire these into MulMod and SquareMod hot paths under feature toggles. - // - Ensure ILGPU compatibility (no BigInteger, no % inside kernels). +// TODO(MOD-OPT): Montgomery/Barrett integration plan +// - Introduce caches of modulus-dependent constants: +// * Montgomery: n' (-(n^{-1}) mod 2^64 or 2^128), R2 = (R^2 mod n) +// * Barrett: mu = floor(2^k / n) for k ∈ {128, 192, 256} +// - Add fast-path for 64-bit NTT primes (modulus.High == 0UL) using pure 64-bit ops. +// - Expose helpers to retrieve/calc constants once per modulus and reuse in kernels. +// - Wire these into MulMod and SquareMod hot paths under feature toggles. +// - Ensure ILGPU compatibility (no BigInteger, no % inside kernels). diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index cb3eec2..ff5bd6d 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -375,7 +375,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv for (int i = 0; i < batchSize; i++) { - divisorDataSpan[i] = MontgomeryDivisorDataCache.Get(divisorSpan[i]); + divisorDataSpan[i] = MontgomeryDivisorData.FromModulus(divisorSpan[i]); } divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize); @@ -538,7 +538,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv _disposed = false; } - public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits) + public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan primes, Span hits) { if (_disposed) { @@ -551,7 +551,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return; } - MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL) { hits.Clear(); diff --git a/PerfectNumbers.Core/Gpu/Pow2OddPowerTable.cs b/PerfectNumbers.Core/Gpu/Pow2OddPowerTable.cs new file mode 100644 index 0000000..ffe3ca2 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/Pow2OddPowerTable.cs @@ -0,0 +1,681 @@ +using System.Runtime.CompilerServices; +namespace PerfectNumbers.Core.Gpu; + +internal static partial class PrimeOrderGpuHeuristics +{ + internal readonly struct Pow2OddPowerTableGpuUInt128 + { + public Pow2OddPowerTableGpuUInt128(GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount) + { + Element0 = baseValue; + if (oddPowerCount == 1) + { + return; + } + + // Reusing baseValue to hold base^2 for the shared odd-power ladder that follows. + baseValue.MulMod(baseValue, modulus); + // TODO: We can calculate baseValue % modulus before loop and use it to increase ladderEntry calculation speed - we'll reuse the base for incremental calculations. + GpuUInt128 current = baseValue; + + // We're manually assigning each field to prevent the compiler to initialize each field twice due to auto-initialization. We're using the action to lower the code base size. + for (int i = 1; i < oddPowerCount; i++) + { + current.MulMod(baseValue, modulus); + switch (i) + { + case 1: + Element1 = current; + break; + case 2: + Element2 = current; + break; + case 3: + Element3 = current; + break; + case 4: + Element4 = current; + break; + case 5: + Element5 = current; + break; + case 6: + Element6 = current; + break; + case 7: + Element7 = current; + break; + case 8: + Element8 = current; + break; + case 9: + Element9 = current; + break; + case 10: + Element10 = current; + break; + case 11: + Element11 = current; + break; + case 12: + Element12 = current; + break; + case 13: + Element13 = current; + break; + case 14: + Element14 = current; + break; + case 15: + Element15 = current; + break; + case 16: + Element16 = current; + break; + case 17: + Element17 = current; + break; + case 18: + Element18 = current; + break; + case 19: + Element19 = current; + break; + case 20: + Element20 = current; + break; + case 21: + Element21 = current; + break; + case 22: + Element22 = current; + break; + case 23: + Element23 = current; + break; + case 24: + Element24 = current; + break; + case 25: + Element25 = current; + break; + case 26: + Element26 = current; + break; + case 27: + Element27 = current; + break; + case 28: + Element28 = current; + break; + case 29: + Element29 = current; + break; + case 30: + Element30 = current; + break; + case 31: + Element31 = current; + break; + case 32: + Element32 = current; + break; + case 33: + Element33 = current; + break; + case 34: + Element34 = current; + break; + case 35: + Element35 = current; + break; + case 36: + Element36 = current; + break; + case 37: + Element37 = current; + break; + case 38: + Element38 = current; + break; + case 39: + Element39 = current; + break; + case 40: + Element40 = current; + break; + case 41: + Element41 = current; + break; + case 42: + Element42 = current; + break; + case 43: + Element43 = current; + break; + case 44: + Element44 = current; + break; + case 45: + Element45 = current; + break; + case 46: + Element46 = current; + break; + case 47: + Element47 = current; + break; + case 48: + Element48 = current; + break; + case 49: + Element49 = current; + break; + case 50: + Element50 = current; + break; + case 51: + Element51 = current; + break; + case 52: + Element52 = current; + break; + case 53: + Element53 = current; + break; + case 54: + Element54 = current; + break; + case 55: + Element55 = current; + break; + case 56: + Element56 = current; + break; + case 57: + Element57 = current; + break; + case 58: + Element58 = current; + break; + case 59: + Element59 = current; + break; + case 60: + Element60 = current; + break; + case 61: + Element61 = current; + break; + case 62: + Element62 = current; + break; + case 63: + Element63 = current; + break; + case 64: + Element64 = current; + break; + case 65: + Element65 = current; + break; + case 66: + Element66 = current; + break; + case 67: + Element67 = current; + break; + case 68: + Element68 = current; + break; + case 69: + Element69 = current; + break; + case 70: + Element70 = current; + break; + case 71: + Element71 = current; + break; + case 72: + Element72 = current; + break; + case 73: + Element73 = current; + break; + case 74: + Element74 = current; + break; + case 75: + Element75 = current; + break; + case 76: + Element76 = current; + break; + case 77: + Element77 = current; + break; + case 78: + Element78 = current; + break; + case 79: + Element79 = current; + break; + case 80: + Element80 = current; + break; + case 81: + Element81 = current; + break; + case 82: + Element82 = current; + break; + case 83: + Element83 = current; + break; + case 84: + Element84 = current; + break; + case 85: + Element85 = current; + break; + case 86: + Element86 = current; + break; + case 87: + Element87 = current; + break; + case 88: + Element88 = current; + break; + case 89: + Element89 = current; + break; + case 90: + Element90 = current; + break; + case 91: + Element91 = current; + break; + case 92: + Element92 = current; + break; + case 93: + Element93 = current; + break; + case 94: + Element94 = current; + break; + case 95: + Element95 = current; + break; + case 96: + Element96 = current; + break; + case 97: + Element97 = current; + break; + case 98: + Element98 = current; + break; + case 99: + Element99 = current; + break; + case 100: + Element100 = current; + break; + case 101: + Element101 = current; + break; + case 102: + Element102 = current; + break; + case 103: + Element103 = current; + break; + case 104: + Element104 = current; + break; + case 105: + Element105 = current; + break; + case 106: + Element106 = current; + break; + case 107: + Element107 = current; + break; + case 108: + Element108 = current; + break; + case 109: + Element109 = current; + break; + case 110: + Element110 = current; + break; + case 111: + Element111 = current; + break; + case 112: + Element112 = current; + break; + case 113: + Element113 = current; + break; + case 114: + Element114 = current; + break; + case 115: + Element115 = current; + break; + case 116: + Element116 = current; + break; + case 117: + Element117 = current; + break; + case 118: + Element118 = current; + break; + case 119: + Element119 = current; + break; + case 120: + Element120 = current; + break; + case 121: + Element121 = current; + break; + case 122: + Element122 = current; + break; + case 123: + Element123 = current; + break; + case 124: + Element124 = current; + break; + case 125: + Element125 = current; + break; + case 126: + Element126 = current; + break; + case 127: + Element127 = current; + break; + } + } + } + + public readonly GpuUInt128 Element0; + public readonly GpuUInt128 Element1; + public readonly GpuUInt128 Element2; + public readonly GpuUInt128 Element3; + public readonly GpuUInt128 Element4; + public readonly GpuUInt128 Element5; + public readonly GpuUInt128 Element6; + public readonly GpuUInt128 Element7; + public readonly GpuUInt128 Element8; + public readonly GpuUInt128 Element9; + public readonly GpuUInt128 Element10; + public readonly GpuUInt128 Element11; + public readonly GpuUInt128 Element12; + public readonly GpuUInt128 Element13; + public readonly GpuUInt128 Element14; + public readonly GpuUInt128 Element15; + public readonly GpuUInt128 Element16; + public readonly GpuUInt128 Element17; + public readonly GpuUInt128 Element18; + public readonly GpuUInt128 Element19; + public readonly GpuUInt128 Element20; + public readonly GpuUInt128 Element21; + public readonly GpuUInt128 Element22; + public readonly GpuUInt128 Element23; + public readonly GpuUInt128 Element24; + public readonly GpuUInt128 Element25; + public readonly GpuUInt128 Element26; + public readonly GpuUInt128 Element27; + public readonly GpuUInt128 Element28; + public readonly GpuUInt128 Element29; + public readonly GpuUInt128 Element30; + public readonly GpuUInt128 Element31; + public readonly GpuUInt128 Element32; + public readonly GpuUInt128 Element33; + public readonly GpuUInt128 Element34; + public readonly GpuUInt128 Element35; + public readonly GpuUInt128 Element36; + public readonly GpuUInt128 Element37; + public readonly GpuUInt128 Element38; + public readonly GpuUInt128 Element39; + public readonly GpuUInt128 Element40; + public readonly GpuUInt128 Element41; + public readonly GpuUInt128 Element42; + public readonly GpuUInt128 Element43; + public readonly GpuUInt128 Element44; + public readonly GpuUInt128 Element45; + public readonly GpuUInt128 Element46; + public readonly GpuUInt128 Element47; + public readonly GpuUInt128 Element48; + public readonly GpuUInt128 Element49; + public readonly GpuUInt128 Element50; + public readonly GpuUInt128 Element51; + public readonly GpuUInt128 Element52; + public readonly GpuUInt128 Element53; + public readonly GpuUInt128 Element54; + public readonly GpuUInt128 Element55; + public readonly GpuUInt128 Element56; + public readonly GpuUInt128 Element57; + public readonly GpuUInt128 Element58; + public readonly GpuUInt128 Element59; + public readonly GpuUInt128 Element60; + public readonly GpuUInt128 Element61; + public readonly GpuUInt128 Element62; + public readonly GpuUInt128 Element63; + public readonly GpuUInt128 Element64; + public readonly GpuUInt128 Element65; + public readonly GpuUInt128 Element66; + public readonly GpuUInt128 Element67; + public readonly GpuUInt128 Element68; + public readonly GpuUInt128 Element69; + public readonly GpuUInt128 Element70; + public readonly GpuUInt128 Element71; + public readonly GpuUInt128 Element72; + public readonly GpuUInt128 Element73; + public readonly GpuUInt128 Element74; + public readonly GpuUInt128 Element75; + public readonly GpuUInt128 Element76; + public readonly GpuUInt128 Element77; + public readonly GpuUInt128 Element78; + public readonly GpuUInt128 Element79; + public readonly GpuUInt128 Element80; + public readonly GpuUInt128 Element81; + public readonly GpuUInt128 Element82; + public readonly GpuUInt128 Element83; + public readonly GpuUInt128 Element84; + public readonly GpuUInt128 Element85; + public readonly GpuUInt128 Element86; + public readonly GpuUInt128 Element87; + public readonly GpuUInt128 Element88; + public readonly GpuUInt128 Element89; + public readonly GpuUInt128 Element90; + public readonly GpuUInt128 Element91; + public readonly GpuUInt128 Element92; + public readonly GpuUInt128 Element93; + public readonly GpuUInt128 Element94; + public readonly GpuUInt128 Element95; + public readonly GpuUInt128 Element96; + public readonly GpuUInt128 Element97; + public readonly GpuUInt128 Element98; + public readonly GpuUInt128 Element99; + public readonly GpuUInt128 Element100; + public readonly GpuUInt128 Element101; + public readonly GpuUInt128 Element102; + public readonly GpuUInt128 Element103; + public readonly GpuUInt128 Element104; + public readonly GpuUInt128 Element105; + public readonly GpuUInt128 Element106; + public readonly GpuUInt128 Element107; + public readonly GpuUInt128 Element108; + public readonly GpuUInt128 Element109; + public readonly GpuUInt128 Element110; + public readonly GpuUInt128 Element111; + public readonly GpuUInt128 Element112; + public readonly GpuUInt128 Element113; + public readonly GpuUInt128 Element114; + public readonly GpuUInt128 Element115; + public readonly GpuUInt128 Element116; + public readonly GpuUInt128 Element117; + public readonly GpuUInt128 Element118; + public readonly GpuUInt128 Element119; + public readonly GpuUInt128 Element120; + public readonly GpuUInt128 Element121; + public readonly GpuUInt128 Element122; + public readonly GpuUInt128 Element123; + public readonly GpuUInt128 Element124; + public readonly GpuUInt128 Element125; + public readonly GpuUInt128 Element126; + public readonly GpuUInt128 Element127; + + public readonly GpuUInt128 this[int index] + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return index switch + { + 0 => Element0, + 1 => Element1, + 2 => Element2, + 3 => Element3, + 4 => Element4, + 5 => Element5, + 6 => Element6, + 7 => Element7, + 8 => Element8, + 9 => Element9, + 10 => Element10, + 11 => Element11, + 12 => Element12, + 13 => Element13, + 14 => Element14, + 15 => Element15, + 16 => Element16, + 17 => Element17, + 18 => Element18, + 19 => Element19, + 20 => Element20, + 21 => Element21, + 22 => Element22, + 23 => Element23, + 24 => Element24, + 25 => Element25, + 26 => Element26, + 27 => Element27, + 28 => Element28, + 29 => Element29, + 30 => Element30, + 31 => Element31, + 32 => Element32, + 33 => Element33, + 34 => Element34, + 35 => Element35, + 36 => Element36, + 37 => Element37, + 38 => Element38, + 39 => Element39, + 40 => Element40, + 41 => Element41, + 42 => Element42, + 43 => Element43, + 44 => Element44, + 45 => Element45, + 46 => Element46, + 47 => Element47, + 48 => Element48, + 49 => Element49, + 50 => Element50, + 51 => Element51, + 52 => Element52, + 53 => Element53, + 54 => Element54, + 55 => Element55, + 56 => Element56, + 57 => Element57, + 58 => Element58, + 59 => Element59, + 60 => Element60, + 61 => Element61, + 62 => Element62, + 63 => Element63, + 64 => Element64, + 65 => Element65, + 66 => Element66, + 67 => Element67, + 68 => Element68, + 69 => Element69, + 70 => Element70, + 71 => Element71, + 72 => Element72, + 73 => Element73, + 74 => Element74, + 75 => Element75, + 76 => Element76, + 77 => Element77, + 78 => Element78, + 79 => Element79, + 80 => Element80, + 81 => Element81, + 82 => Element82, + 83 => Element83, + 84 => Element84, + 85 => Element85, + 86 => Element86, + 87 => Element87, + 88 => Element88, + 89 => Element89, + 90 => Element90, + 91 => Element91, + 92 => Element92, + 93 => Element93, + 94 => Element94, + 95 => Element95, + 96 => Element96, + 97 => Element97, + 98 => Element98, + 99 => Element99, + 100 => Element100, + 101 => Element101, + 102 => Element102, + 103 => Element103, + 104 => Element104, + 105 => Element105, + 106 => Element106, + 107 => Element107, + 108 => Element108, + 109 => Element109, + 110 => Element110, + 111 => Element111, + 112 => Element112, + 113 => Element113, + 114 => Element114, + 115 => Element115, + 116 => Element116, + 117 => Element117, + 118 => Element118, + 119 => Element119, + 120 => Element120, + 121 => Element121, + 122 => Element122, + 123 => Element123, + 124 => Element124, + 125 => Element125, + 126 => Element126, + 127 => Element127, + _ => Element0, + }; + } + } + } +} diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 886296a..fa4251e 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -27,73 +27,37 @@ internal static partial class PrimeOrderGpuHeuristics private static readonly ConcurrentDictionary OrderKernelCache = new(); private static readonly ConcurrentDictionary SmallPrimeDeviceCaches = new(); - public readonly struct OrderKernelConfig + public readonly struct OrderKernelConfig(ulong previousOrder, byte hasPreviousOrder, uint smallFactorLimit, int maxPowChecks, int mode) { - public OrderKernelConfig(ulong previousOrder, byte hasPreviousOrder, uint smallFactorLimit, int maxPowChecks, int mode) - { - PreviousOrder = previousOrder; - HasPreviousOrder = hasPreviousOrder; - SmallFactorLimit = smallFactorLimit; - MaxPowChecks = maxPowChecks; - Mode = mode; - } - - public ulong PreviousOrder { get; } - - public byte HasPreviousOrder { get; } - - public uint SmallFactorLimit { get; } - - public int MaxPowChecks { get; } - - public int Mode { get; } + public readonly ulong PreviousOrder = previousOrder; + public readonly byte HasPreviousOrder = hasPreviousOrder; + public readonly uint SmallFactorLimit = smallFactorLimit; + public readonly int MaxPowChecks = maxPowChecks; + public readonly int Mode = mode; } - public readonly struct OrderKernelBuffers + public readonly struct OrderKernelBuffers( + ArrayView1D phiFactors, + ArrayView1D phiExponents, + ArrayView1D workFactors, + ArrayView1D workExponents, + ArrayView1D candidates, + ArrayView1D stackIndex, + ArrayView1D stackExponent, + ArrayView1D stackProduct, + ArrayView1D result, + ArrayView1D status) { - public OrderKernelBuffers( - ArrayView1D phiFactors, - ArrayView1D phiExponents, - ArrayView1D workFactors, - ArrayView1D workExponents, - ArrayView1D candidates, - ArrayView1D stackIndex, - ArrayView1D stackExponent, - ArrayView1D stackProduct, - ArrayView1D result, - ArrayView1D status) - { - PhiFactors = phiFactors; - PhiExponents = phiExponents; - WorkFactors = workFactors; - WorkExponents = workExponents; - Candidates = candidates; - StackIndex = stackIndex; - StackExponent = stackExponent; - StackProduct = stackProduct; - Result = result; - Status = status; - } - - public ArrayView1D PhiFactors { get; } - - public ArrayView1D PhiExponents { get; } - - public ArrayView1D WorkFactors { get; } - - public ArrayView1D WorkExponents { get; } - - public ArrayView1D Candidates { get; } - - public ArrayView1D StackIndex { get; } - - public ArrayView1D StackExponent { get; } - - public ArrayView1D StackProduct { get; } - - public ArrayView1D Result { get; } - - public ArrayView1D Status { get; } + public readonly ArrayView1D PhiFactors = phiFactors; + public readonly ArrayView1D PhiExponents = phiExponents; + public readonly ArrayView1D WorkFactors = workFactors; + public readonly ArrayView1D WorkExponents = workExponents; + public readonly ArrayView1D Candidates = candidates; + public readonly ArrayView1D StackIndex = stackIndex; + public readonly ArrayView1D StackExponent = stackExponent; + public readonly ArrayView1D StackProduct = stackProduct; + public readonly ArrayView1D Result = result; + public readonly ArrayView1D Status = status; } private delegate void OrderKernelLauncher( @@ -125,681 +89,6 @@ internal static partial class PrimeOrderGpuHeuristics private const int GpuSmallPrimeFactorSlots = 64; - - private struct Pow2OddPowerTable - { - public GpuUInt128 Element0; - public GpuUInt128 Element1; - public GpuUInt128 Element2; - public GpuUInt128 Element3; - public GpuUInt128 Element4; - public GpuUInt128 Element5; - public GpuUInt128 Element6; - public GpuUInt128 Element7; - public GpuUInt128 Element8; - public GpuUInt128 Element9; - public GpuUInt128 Element10; - public GpuUInt128 Element11; - public GpuUInt128 Element12; - public GpuUInt128 Element13; - public GpuUInt128 Element14; - public GpuUInt128 Element15; - public GpuUInt128 Element16; - public GpuUInt128 Element17; - public GpuUInt128 Element18; - public GpuUInt128 Element19; - public GpuUInt128 Element20; - public GpuUInt128 Element21; - public GpuUInt128 Element22; - public GpuUInt128 Element23; - public GpuUInt128 Element24; - public GpuUInt128 Element25; - public GpuUInt128 Element26; - public GpuUInt128 Element27; - public GpuUInt128 Element28; - public GpuUInt128 Element29; - public GpuUInt128 Element30; - public GpuUInt128 Element31; - public GpuUInt128 Element32; - public GpuUInt128 Element33; - public GpuUInt128 Element34; - public GpuUInt128 Element35; - public GpuUInt128 Element36; - public GpuUInt128 Element37; - public GpuUInt128 Element38; - public GpuUInt128 Element39; - public GpuUInt128 Element40; - public GpuUInt128 Element41; - public GpuUInt128 Element42; - public GpuUInt128 Element43; - public GpuUInt128 Element44; - public GpuUInt128 Element45; - public GpuUInt128 Element46; - public GpuUInt128 Element47; - public GpuUInt128 Element48; - public GpuUInt128 Element49; - public GpuUInt128 Element50; - public GpuUInt128 Element51; - public GpuUInt128 Element52; - public GpuUInt128 Element53; - public GpuUInt128 Element54; - public GpuUInt128 Element55; - public GpuUInt128 Element56; - public GpuUInt128 Element57; - public GpuUInt128 Element58; - public GpuUInt128 Element59; - public GpuUInt128 Element60; - public GpuUInt128 Element61; - public GpuUInt128 Element62; - public GpuUInt128 Element63; - public GpuUInt128 Element64; - public GpuUInt128 Element65; - public GpuUInt128 Element66; - public GpuUInt128 Element67; - public GpuUInt128 Element68; - public GpuUInt128 Element69; - public GpuUInt128 Element70; - public GpuUInt128 Element71; - public GpuUInt128 Element72; - public GpuUInt128 Element73; - public GpuUInt128 Element74; - public GpuUInt128 Element75; - public GpuUInt128 Element76; - public GpuUInt128 Element77; - public GpuUInt128 Element78; - public GpuUInt128 Element79; - public GpuUInt128 Element80; - public GpuUInt128 Element81; - public GpuUInt128 Element82; - public GpuUInt128 Element83; - public GpuUInt128 Element84; - public GpuUInt128 Element85; - public GpuUInt128 Element86; - public GpuUInt128 Element87; - public GpuUInt128 Element88; - public GpuUInt128 Element89; - public GpuUInt128 Element90; - public GpuUInt128 Element91; - public GpuUInt128 Element92; - public GpuUInt128 Element93; - public GpuUInt128 Element94; - public GpuUInt128 Element95; - public GpuUInt128 Element96; - public GpuUInt128 Element97; - public GpuUInt128 Element98; - public GpuUInt128 Element99; - public GpuUInt128 Element100; - public GpuUInt128 Element101; - public GpuUInt128 Element102; - public GpuUInt128 Element103; - public GpuUInt128 Element104; - public GpuUInt128 Element105; - public GpuUInt128 Element106; - public GpuUInt128 Element107; - public GpuUInt128 Element108; - public GpuUInt128 Element109; - public GpuUInt128 Element110; - public GpuUInt128 Element111; - public GpuUInt128 Element112; - public GpuUInt128 Element113; - public GpuUInt128 Element114; - public GpuUInt128 Element115; - public GpuUInt128 Element116; - public GpuUInt128 Element117; - public GpuUInt128 Element118; - public GpuUInt128 Element119; - public GpuUInt128 Element120; - public GpuUInt128 Element121; - public GpuUInt128 Element122; - public GpuUInt128 Element123; - public GpuUInt128 Element124; - public GpuUInt128 Element125; - public GpuUInt128 Element126; - public GpuUInt128 Element127; - - public GpuUInt128 this[int index] - { - readonly get - { - return GetElement(index); - } - - set - { - SetElement(index, value); - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private readonly GpuUInt128 GetElement(int index) => index switch - { - 0 => Element0, - 1 => Element1, - 2 => Element2, - 3 => Element3, - 4 => Element4, - 5 => Element5, - 6 => Element6, - 7 => Element7, - 8 => Element8, - 9 => Element9, - 10 => Element10, - 11 => Element11, - 12 => Element12, - 13 => Element13, - 14 => Element14, - 15 => Element15, - 16 => Element16, - 17 => Element17, - 18 => Element18, - 19 => Element19, - 20 => Element20, - 21 => Element21, - 22 => Element22, - 23 => Element23, - 24 => Element24, - 25 => Element25, - 26 => Element26, - 27 => Element27, - 28 => Element28, - 29 => Element29, - 30 => Element30, - 31 => Element31, - 32 => Element32, - 33 => Element33, - 34 => Element34, - 35 => Element35, - 36 => Element36, - 37 => Element37, - 38 => Element38, - 39 => Element39, - 40 => Element40, - 41 => Element41, - 42 => Element42, - 43 => Element43, - 44 => Element44, - 45 => Element45, - 46 => Element46, - 47 => Element47, - 48 => Element48, - 49 => Element49, - 50 => Element50, - 51 => Element51, - 52 => Element52, - 53 => Element53, - 54 => Element54, - 55 => Element55, - 56 => Element56, - 57 => Element57, - 58 => Element58, - 59 => Element59, - 60 => Element60, - 61 => Element61, - 62 => Element62, - 63 => Element63, - 64 => Element64, - 65 => Element65, - 66 => Element66, - 67 => Element67, - 68 => Element68, - 69 => Element69, - 70 => Element70, - 71 => Element71, - 72 => Element72, - 73 => Element73, - 74 => Element74, - 75 => Element75, - 76 => Element76, - 77 => Element77, - 78 => Element78, - 79 => Element79, - 80 => Element80, - 81 => Element81, - 82 => Element82, - 83 => Element83, - 84 => Element84, - 85 => Element85, - 86 => Element86, - 87 => Element87, - 88 => Element88, - 89 => Element89, - 90 => Element90, - 91 => Element91, - 92 => Element92, - 93 => Element93, - 94 => Element94, - 95 => Element95, - 96 => Element96, - 97 => Element97, - 98 => Element98, - 99 => Element99, - 100 => Element100, - 101 => Element101, - 102 => Element102, - 103 => Element103, - 104 => Element104, - 105 => Element105, - 106 => Element106, - 107 => Element107, - 108 => Element108, - 109 => Element109, - 110 => Element110, - 111 => Element111, - 112 => Element112, - 113 => Element113, - 114 => Element114, - 115 => Element115, - 116 => Element116, - 117 => Element117, - 118 => Element118, - 119 => Element119, - 120 => Element120, - 121 => Element121, - 122 => Element122, - 123 => Element123, - 124 => Element124, - 125 => Element125, - 126 => Element126, - 127 => Element127, - _ => Element0, - }; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private void SetElement(int index, GpuUInt128 value) - { - switch (index) - { - case 0: - Element0 = value; - return; - case 1: - Element1 = value; - return; - case 2: - Element2 = value; - return; - case 3: - Element3 = value; - return; - case 4: - Element4 = value; - return; - case 5: - Element5 = value; - return; - case 6: - Element6 = value; - return; - case 7: - Element7 = value; - return; - case 8: - Element8 = value; - return; - case 9: - Element9 = value; - return; - case 10: - Element10 = value; - return; - case 11: - Element11 = value; - return; - case 12: - Element12 = value; - return; - case 13: - Element13 = value; - return; - case 14: - Element14 = value; - return; - case 15: - Element15 = value; - return; - case 16: - Element16 = value; - return; - case 17: - Element17 = value; - return; - case 18: - Element18 = value; - return; - case 19: - Element19 = value; - return; - case 20: - Element20 = value; - return; - case 21: - Element21 = value; - return; - case 22: - Element22 = value; - return; - case 23: - Element23 = value; - return; - case 24: - Element24 = value; - return; - case 25: - Element25 = value; - return; - case 26: - Element26 = value; - return; - case 27: - Element27 = value; - return; - case 28: - Element28 = value; - return; - case 29: - Element29 = value; - return; - case 30: - Element30 = value; - return; - case 31: - Element31 = value; - return; - case 32: - Element32 = value; - return; - case 33: - Element33 = value; - return; - case 34: - Element34 = value; - return; - case 35: - Element35 = value; - return; - case 36: - Element36 = value; - return; - case 37: - Element37 = value; - return; - case 38: - Element38 = value; - return; - case 39: - Element39 = value; - return; - case 40: - Element40 = value; - return; - case 41: - Element41 = value; - return; - case 42: - Element42 = value; - return; - case 43: - Element43 = value; - return; - case 44: - Element44 = value; - return; - case 45: - Element45 = value; - return; - case 46: - Element46 = value; - return; - case 47: - Element47 = value; - return; - case 48: - Element48 = value; - return; - case 49: - Element49 = value; - return; - case 50: - Element50 = value; - return; - case 51: - Element51 = value; - return; - case 52: - Element52 = value; - return; - case 53: - Element53 = value; - return; - case 54: - Element54 = value; - return; - case 55: - Element55 = value; - return; - case 56: - Element56 = value; - return; - case 57: - Element57 = value; - return; - case 58: - Element58 = value; - return; - case 59: - Element59 = value; - return; - case 60: - Element60 = value; - return; - case 61: - Element61 = value; - return; - case 62: - Element62 = value; - return; - case 63: - Element63 = value; - return; - case 64: - Element64 = value; - return; - case 65: - Element65 = value; - return; - case 66: - Element66 = value; - return; - case 67: - Element67 = value; - return; - case 68: - Element68 = value; - return; - case 69: - Element69 = value; - return; - case 70: - Element70 = value; - return; - case 71: - Element71 = value; - return; - case 72: - Element72 = value; - return; - case 73: - Element73 = value; - return; - case 74: - Element74 = value; - return; - case 75: - Element75 = value; - return; - case 76: - Element76 = value; - return; - case 77: - Element77 = value; - return; - case 78: - Element78 = value; - return; - case 79: - Element79 = value; - return; - case 80: - Element80 = value; - return; - case 81: - Element81 = value; - return; - case 82: - Element82 = value; - return; - case 83: - Element83 = value; - return; - case 84: - Element84 = value; - return; - case 85: - Element85 = value; - return; - case 86: - Element86 = value; - return; - case 87: - Element87 = value; - return; - case 88: - Element88 = value; - return; - case 89: - Element89 = value; - return; - case 90: - Element90 = value; - return; - case 91: - Element91 = value; - return; - case 92: - Element92 = value; - return; - case 93: - Element93 = value; - return; - case 94: - Element94 = value; - return; - case 95: - Element95 = value; - return; - case 96: - Element96 = value; - return; - case 97: - Element97 = value; - return; - case 98: - Element98 = value; - return; - case 99: - Element99 = value; - return; - case 100: - Element100 = value; - return; - case 101: - Element101 = value; - return; - case 102: - Element102 = value; - return; - case 103: - Element103 = value; - return; - case 104: - Element104 = value; - return; - case 105: - Element105 = value; - return; - case 106: - Element106 = value; - return; - case 107: - Element107 = value; - return; - case 108: - Element108 = value; - return; - case 109: - Element109 = value; - return; - case 110: - Element110 = value; - return; - case 111: - Element111 = value; - return; - case 112: - Element112 = value; - return; - case 113: - Element113 = value; - return; - case 114: - Element114 = value; - return; - case 115: - Element115 = value; - return; - case 116: - Element116 = value; - return; - case 117: - Element117 = value; - return; - case 118: - Element118 = value; - return; - case 119: - Element119 = value; - return; - case 120: - Element120 = value; - return; - case 121: - Element121 = value; - return; - case 122: - Element122 = value; - return; - case 123: - Element123 = value; - return; - case 124: - Element124 = value; - return; - case 125: - Element125 = value; - return; - case 126: - Element126 = value; - return; - case 127: - Element127 = value; - return; - default: - // ILGPU kernels cannot throw exceptions, and callers guarantee the index range. - return; - } - } - } - internal static ConcurrentDictionary OverflowRegistry => OverflowedPrimes; internal static ConcurrentDictionary OverflowRegistryWide => OverflowedPrimesWide; @@ -1039,36 +328,44 @@ internal static partial class PrimeOrderGpuHeuristics remainingOut[0] = remainingLocal; fullyFactoredOut[0] = remainingLocal == 1UL ? (byte)1 : (byte)0; } - public static GpuPow2ModStatus TryPow2Mod(ulong exponent, ulong prime, out ulong remainder) + public static GpuPow2ModStatus TryPow2Mod(ulong exponent, ulong prime, out ulong remainder, in MontgomeryDivisorData divisorData) { Span exponents = stackalloc ulong[1]; Span remainders = stackalloc ulong[1]; exponents[0] = exponent; - GpuPow2ModStatus status = TryPow2ModBatch(exponents, prime, remainders); + GpuPow2ModStatus status = TryPow2ModBatch(exponents, prime, remainders, divisorData); remainder = remainders[0]; return status; } - public static GpuPow2ModStatus TryPow2ModBatch(ReadOnlySpan exponents, ulong prime, Span remainders) + public static GpuPow2ModStatus TryPow2ModBatch(ReadOnlySpan exponents, ulong prime, Span remainders, in MontgomeryDivisorData divisorData) { - if (exponents.Length == 0) - { - return GpuPow2ModStatus.Success; - } + // This will never occur in production code + // if (exponents.Length == 0) + // { + // return GpuPow2ModStatus.Success; + // } - if (remainders.Length < exponents.Length) - { - throw new ArgumentException("Remainder span is shorter than the exponent span.", nameof(remainders)); - } + // This will never occur in production code + // if (remainders.Length < exponents.Length) + // { + // throw new ArgumentException("Remainder span is shorter than the exponent span.", nameof(remainders)); + // } - Span target = remainders.Slice(0, exponents.Length); - target.Clear(); + // This will never occur in production code + // if (remainders.Length > exponents.Length) + // throw new ArgumentException("Remainder span is longer than the exponent span.", nameof(remainders)); - if (prime <= 1UL) - { - return GpuPow2ModStatus.Unavailable; - } + // We need to clear more, but we save on allocations + // Span target = remainders[..exponents.Length]; + remainders.Clear(); + + // This will never occur in production code + // if (prime <= 1UL) + // { + // return GpuPow2ModStatus.Unavailable; + // } ConcurrentDictionary overflowRegistry = OverflowedPrimes; @@ -1093,11 +390,11 @@ internal static partial class PrimeOrderGpuHeuristics } } - bool computed = TryComputeOnGpu(exponents, prime, target); + bool computed = TryComputeOnGpu(exponents, prime, divisorData, remainders); return computed ? GpuPow2ModStatus.Success : GpuPow2ModStatus.Unavailable; } - public static GpuPow2ModStatus TryPow2Mod(UInt128 exponent, UInt128 prime, out UInt128 remainder) + public static GpuPow2ModStatus TryPow2Mod(in UInt128 exponent, in UInt128 prime, out UInt128 remainder) { Span exponents = stackalloc UInt128[1]; Span remainders = stackalloc UInt128[1]; @@ -1122,104 +419,111 @@ internal static partial class PrimeOrderGpuHeuristics { result = default; - try + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var execution = lease.EnterExecutionScope(); + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; + + var kernel = GetOrderKernel(accelerator); + SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator); + + var phiFactorBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); + var phiExponentBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); + var workFactorBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); + var workExponentBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); + var candidateBuffer = accelerator.Allocate1D(HeuristicCandidateLimit); + var stackIndexBuffer = accelerator.Allocate1D(HeuristicStackCapacity); + var stackExponentBuffer = accelerator.Allocate1D(HeuristicStackCapacity); + var stackProductBuffer = accelerator.Allocate1D(HeuristicStackCapacity); + var resultBuffer = accelerator.Allocate1D(1); + var statusBuffer = accelerator.Allocate1D(1); + + phiFactorBuffer.MemSetToZero(); + phiExponentBuffer.MemSetToZero(); + workFactorBuffer.MemSetToZero(); + workExponentBuffer.MemSetToZero(); + candidateBuffer.MemSetToZero(); + stackIndexBuffer.MemSetToZero(); + stackExponentBuffer.MemSetToZero(); + stackProductBuffer.MemSetToZero(); + resultBuffer.MemSetToZero(); + statusBuffer.MemSetToZero(); + + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + byte hasPrevious = previousOrder.HasValue ? (byte)1 : (byte)0; + ulong previousValue = previousOrder ?? 0UL; + + var kernelConfig = new OrderKernelConfig(previousValue, hasPrevious, limit, config.MaxPowChecks, (int)config.Mode); + var buffers = new OrderKernelBuffers( + phiFactorBuffer.View, + phiExponentBuffer.View, + workFactorBuffer.View, + workExponentBuffer.View, + candidateBuffer.View, + stackIndexBuffer.View, + stackExponentBuffer.View, + stackProductBuffer.View, + resultBuffer.View, + statusBuffer.View); + + kernel( + stream, + 1, + prime, + kernelConfig, + divisorData, + cache.Primes!.View, + cache.Squares!.View, + cache.Count, + buffers); + + stream.Synchronize(); + + byte status = 0; + statusBuffer.View.CopyToCPU(ref status, 1); + + PrimeOrderKernelStatus kernelStatus = (PrimeOrderKernelStatus)status; + if (kernelStatus == PrimeOrderKernelStatus.Fallback) { - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); - try - { - using var execution = lease.EnterExecutionScope(); - Accelerator accelerator = lease.Accelerator; - AcceleratorStream stream = lease.Stream; - - var kernel = GetOrderKernel(accelerator); - SmallPrimeDeviceCache cache = GetSmallPrimeDeviceCache(accelerator); - - using var phiFactorBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); - using var phiExponentBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); - using var workFactorBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); - using var workExponentBuffer = accelerator.Allocate1D(GpuSmallPrimeFactorSlots); - using var candidateBuffer = accelerator.Allocate1D(HeuristicCandidateLimit); - using var stackIndexBuffer = accelerator.Allocate1D(HeuristicStackCapacity); - using var stackExponentBuffer = accelerator.Allocate1D(HeuristicStackCapacity); - using var stackProductBuffer = accelerator.Allocate1D(HeuristicStackCapacity); - using var resultBuffer = accelerator.Allocate1D(1); - using var statusBuffer = accelerator.Allocate1D(1); - - phiFactorBuffer.MemSetToZero(); - phiExponentBuffer.MemSetToZero(); - workFactorBuffer.MemSetToZero(); - workExponentBuffer.MemSetToZero(); - candidateBuffer.MemSetToZero(); - stackIndexBuffer.MemSetToZero(); - stackExponentBuffer.MemSetToZero(); - stackProductBuffer.MemSetToZero(); - resultBuffer.MemSetToZero(); - statusBuffer.MemSetToZero(); - - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - byte hasPrevious = previousOrder.HasValue ? (byte)1 : (byte)0; - ulong previousValue = previousOrder ?? 0UL; - - var kernelConfig = new OrderKernelConfig(previousValue, hasPrevious, limit, config.MaxPowChecks, (int)config.Mode); - var buffers = new OrderKernelBuffers( - phiFactorBuffer.View, - phiExponentBuffer.View, - workFactorBuffer.View, - workExponentBuffer.View, - candidateBuffer.View, - stackIndexBuffer.View, - stackExponentBuffer.View, - stackProductBuffer.View, - resultBuffer.View, - statusBuffer.View); - - kernel( - stream, - 1, - prime, - kernelConfig, - divisorData, - cache.Primes!.View, - cache.Squares!.View, - cache.Count, - buffers); - - stream.Synchronize(); - - byte status = 0; - statusBuffer.View.CopyToCPU(ref status, 1); - - PrimeOrderKernelStatus kernelStatus = (PrimeOrderKernelStatus)status; - if (kernelStatus == PrimeOrderKernelStatus.Fallback) - { - return false; - } - - if (kernelStatus == PrimeOrderKernelStatus.PollardOverflow) - { - throw new InvalidOperationException("GPU Pollard Rho stack overflow; increase HeuristicStackCapacity."); - } - - ulong order = 0UL; - resultBuffer.View.CopyToCPU(ref order, 1); - - PrimeOrderCalculator.PrimeOrderStatus finalStatus = kernelStatus == PrimeOrderKernelStatus.Found - ? PrimeOrderCalculator.PrimeOrderStatus.Found - : PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved; - - result = new PrimeOrderCalculator.PrimeOrderResult(finalStatus, order); - return true; - } - finally - { - lease.Dispose(); - } - } - catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) - { - result = default; + DisposeResources(); + lease.Dispose(); return false; } + + if (kernelStatus == PrimeOrderKernelStatus.PollardOverflow) + { + DisposeResources(); + lease.Dispose(); + throw new InvalidOperationException("GPU Pollard Rho stack overflow; increase HeuristicStackCapacity."); + } + + ulong order = 0UL; + resultBuffer.View.CopyToCPU(ref order, 1); + + PrimeOrderCalculator.PrimeOrderStatus finalStatus = kernelStatus == PrimeOrderKernelStatus.Found + ? PrimeOrderCalculator.PrimeOrderStatus.Found + : PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved; + + result = new PrimeOrderCalculator.PrimeOrderResult(finalStatus, order); + DisposeResources(); + lease.Dispose(); + + return true; + + void DisposeResources() + { + phiFactorBuffer.Dispose(); + phiExponentBuffer.Dispose(); + workFactorBuffer.Dispose(); + workExponentBuffer.Dispose(); + candidateBuffer.Dispose(); + stackIndexBuffer.Dispose(); + stackExponentBuffer.Dispose(); + stackProductBuffer.Dispose(); + resultBuffer.Dispose(); + statusBuffer.Dispose(); + execution.Dispose(); + } } private static OrderKernelLauncher GetOrderKernel(Accelerator accelerator) @@ -2356,49 +1660,48 @@ internal static partial class PrimeOrderGpuHeuristics return computed ? GpuPow2ModStatus.Success : GpuPow2ModStatus.Unavailable; } - private static bool TryComputeOnGpu(ReadOnlySpan exponents, ulong prime, Span results) + private static bool TryComputeOnGpu(ReadOnlySpan exponents, ulong prime, in MontgomeryDivisorData divisorData, Span results) { + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var execution = lease.EnterExecutionScope(); + Accelerator accelerator = lease.Accelerator; + AcceleratorStream stream = lease.Stream; + var kernel = GetPow2ModKernel(accelerator); + var exponentBuffer = accelerator.Allocate1D(exponents.Length); + var remainderBuffer = accelerator.Allocate1D(exponents.Length); + + exponentBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(exponents), exponents.Length); + remainderBuffer.MemSetToZero(); + try { - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); - try - { - using var execution = lease.EnterExecutionScope(); - Accelerator accelerator = lease.Accelerator; - AcceleratorStream stream = lease.Stream; - var kernel = GetPow2ModKernel(accelerator); - using var exponentBuffer = accelerator.Allocate1D(exponents.Length); - using var remainderBuffer = accelerator.Allocate1D(exponents.Length); - - exponentBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(exponents), exponents.Length); - remainderBuffer.MemSetToZero(); - - MontgomeryDivisorData divisor = MontgomeryDivisorDataCache.Get(prime); - kernel(stream, exponents.Length, exponentBuffer.View, divisor, remainderBuffer.View); - - stream.Synchronize(); - - remainderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(results), results.Length); - return true; - } - finally - { - lease.Dispose(); - } + kernel(stream, exponents.Length, exponentBuffer.View, divisorData, remainderBuffer.View); + stream.Synchronize(); } - catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) + catch (Exception) { - ComputePow2ModCpu(exponents, prime, results); - return true; + Console.WriteLine($"Exception for {prime} and exponents: {string.Join(",", exponents.ToArray())}."); + exponentBuffer.Dispose(); + remainderBuffer.Dispose(); + execution.Dispose(); + lease.Dispose(); + throw; } + + remainderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(results), exponents.Length); + exponentBuffer.Dispose(); + remainderBuffer.Dispose(); + execution.Dispose(); + lease.Dispose(); + return true; } - private static void ComputePow2ModCpu(ReadOnlySpan exponents, ulong prime, Span results) + private static void ComputePow2ModCpu(ReadOnlySpan exponents, ulong prime, in MontgomeryDivisorData divisorData, Span results) { int length = exponents.Length; for (int i = 0; i < length; i++) { - results[i] = Pow2ModCpu(exponents[i], prime); + results[i] = Pow2ModCpu(exponents[i], prime, divisorData); } } @@ -2416,12 +1719,12 @@ internal static partial class PrimeOrderGpuHeuristics try { - using var execution = lease.EnterExecutionScope(); + var execution = lease.EnterExecutionScope(); Accelerator accelerator = lease.Accelerator; AcceleratorStream stream = lease.Stream; var kernel = GetPow2ModWideKernel(accelerator); - using var exponentBuffer = accelerator.Allocate1D(length); - using var remainderBuffer = accelerator.Allocate1D(length); + var exponentBuffer = accelerator.Allocate1D(length); + var remainderBuffer = accelerator.Allocate1D(length); Span exponentSpan = length <= WideStackThreshold ? stackalloc GpuUInt128[length] @@ -2451,12 +1754,15 @@ internal static partial class PrimeOrderGpuHeuristics results[i] = (UInt128)resultSpan[i]; } + exponentBuffer.Dispose(); + remainderBuffer.Dispose(); + execution.Dispose(); return true; } - catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) + catch (Exception) { - ComputePow2ModCpuWide(exponents, prime, results); - return true; + Console.WriteLine($"Exception when computing {prime} with exponents: {string.Join(", ", exponents.ToArray())}"); + throw; } finally { @@ -2483,15 +1789,14 @@ internal static partial class PrimeOrderGpuHeuristics } } - private static ulong Pow2ModCpu(ulong exponent, ulong modulus) + private static ulong Pow2ModCpu(ulong exponent, ulong modulus, in MontgomeryDivisorData divisorData) { if (modulus <= 1UL) { return 0UL; } - MontgomeryDivisorData divisor = MontgomeryDivisorDataCache.Get(modulus); - return exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false); + return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false); } private static Action, MontgomeryDivisorData, ArrayView1D> GetPow2ModKernel(Accelerator accelerator) @@ -2526,39 +1831,68 @@ internal static partial class PrimeOrderGpuHeuristics remainders[index] = Pow2ModKernelCore(exponent, modulus); } + private static GpuUInt128[] InitializeOddPowersTable(GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount) + { + GpuUInt128[] result = new GpuUInt128[PerfectNumberConstants.MaxOddPowersCount]; + result[0] = baseValue; + if (oddPowerCount == 1) + { + return result; + } + + // Reusing baseValue to hold base^2 for the shared odd-power ladder that follows. + baseValue.MulMod(baseValue, modulus); + + // TODO: We can calculate baseValue % modulus before loop and use it to increase ladderEntry calculation speed - we'll reuse the base for incremental calculations. + GpuUInt128 current = baseValue; + + // We're manually assigning each field to prevent the compiler to initialize each field twice due to auto-initialization. We're using the action to lower the code base size. + for (int i = 1; i < oddPowerCount; i++) + { + current.MulMod(baseValue, modulus); + result[i] = current; + } + + return result; + } + private static GpuUInt128 Pow2ModKernelCore(GpuUInt128 exponent, GpuUInt128 modulus) { - if (modulus == GpuUInt128.One) - { - return GpuUInt128.Zero; - } + // This should never happen in production code. + // if (modulus == GpuUInt128.One) + // { + // return GpuUInt128.Zero; + // } - if (exponent.IsZero) - { - return GpuUInt128.One; - } + // This should never happen in production code. + // if (exponent.IsZero) + // { + // return GpuUInt128.One; + // } - GpuUInt128 baseValue = new GpuUInt128(2UL); - if (baseValue.CompareTo(modulus) >= 0) - { - baseValue.Sub(modulus); - } + GpuUInt128 baseValue = UInt128Numbers.TwoGpu; + + // This should never happen in production code - 2 should never be greater or equal modulus. + // if (baseValue.CompareTo(modulus) >= 0) + // { + // baseValue.Sub(modulus); + // } if (ShouldUseSingleBit(exponent)) { return Pow2MontgomeryModSingleBit(exponent, modulus, baseValue); } - int bitLength = exponent.GetBitLength(); - int windowSize = GetWindowSize(bitLength); + int index = exponent.GetBitLength(); + int windowSize = GetWindowSize(index); + index--; + int oddPowerCount = 1 << (windowSize - 1); - - Pow2OddPowerTable oddPowers = default; - InitializeOddPowers(ref oddPowers, baseValue, modulus, oddPowerCount); - + GpuUInt128[] oddPowers = InitializeOddPowersTable(baseValue, modulus, oddPowerCount); GpuUInt128 result = GpuUInt128.One; - int index = bitLength - 1; + int windowStart; + ulong windowValue; while (index >= 0) { if (!IsBitSet(exponent, index)) @@ -2568,7 +1902,7 @@ internal static partial class PrimeOrderGpuHeuristics continue; } - int windowStart = index - windowSize + 1; + windowStart = index - windowSize + 1; if (windowStart < 0) { windowStart = 0; @@ -2579,24 +1913,27 @@ internal static partial class PrimeOrderGpuHeuristics windowStart++; } - int windowBitCount = index - windowStart + 1; - for (int square = 0; square < windowBitCount; square++) + // We're reusing oddPowerCount as windowBitCount here to lower registry pressure & avoid additional allocation + oddPowerCount = index - windowStart + 1; + index = windowStart - 1; + windowValue = ExtractWindowValue(exponent, windowStart, oddPowerCount); + + // We're reusing windowStart as square here to lower registry pressure & avoid additional allocation + for (windowStart = 0; windowStart < oddPowerCount; windowStart++) { result.MulMod(result, modulus); } - ulong windowValue = ExtractWindowValue(exponent, windowStart, windowBitCount); - int tableIndex = (int)((windowValue - 1UL) >> 1); - GpuUInt128 factor = oddPowers[tableIndex]; - result.MulMod(factor, modulus); - - index = windowStart - 1; + // We're reusing windowStart as tableIndex here to lower registry pressure & avoid additional allocation + windowStart = (int)((windowValue - 1UL) >> 1); + result.MulMod(oddPowers[windowStart], modulus); } return result; } - private static bool ShouldUseSingleBit(in GpuUInt128 exponent) => exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool ShouldUseSingleBit(GpuUInt128 exponent) => exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold; private static int GetWindowSize(int bitLength) { @@ -2628,24 +1965,6 @@ internal static partial class PrimeOrderGpuHeuristics return Pow2WindowSizeBits; } - private static void InitializeOddPowers(ref Pow2OddPowerTable oddPowers, GpuUInt128 baseValue, GpuUInt128 modulus, int oddPowerCount) - { - oddPowers[0] = baseValue; - if (oddPowerCount == 1) - { - return; - } - - // Reusing baseValue to hold base^2 for the shared odd-power ladder that follows. - baseValue.MulMod(baseValue, modulus); - for (int i = 1; i < oddPowerCount; i++) - { - GpuUInt128 ladderEntry = oddPowers[i - 1]; - ladderEntry.MulMod(baseValue, modulus); - oddPowers[i] = ladderEntry; - } - } - private static GpuUInt128 Pow2MontgomeryModSingleBit(GpuUInt128 exponent, GpuUInt128 modulus, GpuUInt128 baseValue) { GpuUInt128 result = GpuUInt128.One; @@ -2670,7 +1989,7 @@ internal static partial class PrimeOrderGpuHeuristics return result; } - private static bool IsBitSet(in GpuUInt128 value, int bitIndex) + private static bool IsBitSet(GpuUInt128 value, int bitIndex) { if (bitIndex >= 64) { @@ -2680,14 +1999,13 @@ internal static partial class PrimeOrderGpuHeuristics return ((value.Low >> bitIndex) & 1UL) != 0UL; } - private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount) + private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount) { if (windowStart != 0) { - GpuUInt128 shifted = exponent; - shifted.ShiftRight(windowStart); + exponent.ShiftRight(windowStart); ulong mask = (1UL << windowBitCount) - 1UL; - return shifted.Low & mask; + return exponent.Low & mask; } ulong directMask = (1UL << windowBitCount) - 1UL; diff --git a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs index cfdbbd6..e700843 100644 --- a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs @@ -20,6 +20,6 @@ public interface IMersenneNumberDivisorByDivisorTester public interface IDivisorScanSession : IDisposable { - void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits); + void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan primes, Span hits); } } diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index c585dcc..cc825d5 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -14,55 +14,55 @@ namespace PerfectNumbers.Core; public class MersenneDivisorCycles { - private List<(ulong divisor, ulong cycleLength)> _table = []; - // Lightweight read-mostly cache for small divisors (<= 4,000,000). 0 => unknown - private ulong[]? _smallCycles; + private List<(ulong divisor, ulong cycleLength)> _table = []; + // Lightweight read-mostly cache for small divisors (<= 4,000,000). 0 => unknown + private ulong[]? _smallCycles; - public static MersenneDivisorCycles Shared { get; } = new MersenneDivisorCycles(); + public static MersenneDivisorCycles Shared { get; } = new MersenneDivisorCycles(); - public MersenneDivisorCycles() - { - } - - public void LoadFrom(string path) - { - EnsureSmallBuffer(); - List<(ulong divisor, ulong cycle)> cycles = []; - using Stream outputStream = new FileStream(path, FileMode.Open, FileAccess.Read, FileShare.Read, BufferSize10M, useAsync: false); - foreach (var (d, c) in EnumerateStream(outputStream)) - { - cycles.Add((d, c)); - if (d <= PerfectNumberConstants.MaxQForDivisorCycles) - { - _smallCycles![(int)d] = c == 0UL ? 1UL : c; - } - } - - cycles.Sort((a, b) => a.divisor < b.divisor ? -1 : (a.divisor > b.divisor ? 1 : 0)); - _table = cycles; - } - - public const int SmallDivisorsMax = PerfectNumberConstants.MaxQForDivisorCycles; - - // Provides a snapshot array of small cycles [0..PerfectNumberConstants.MaxQForDivisorCycles]. Index is the divisor. - // The returned array is a copy safe for use across threads and device uploads. - public ulong[] ExportSmallCyclesSnapshot() - { - EnsureSmallBuffer(); - ulong[] snapshot = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; - Array.Copy(_smallCycles!, snapshot, snapshot.Length); - return snapshot; - } - - public static bool CycleEqualsExponent(ulong divisor, ulong exponent) + public MersenneDivisorCycles() { - Dictionary? cache = null; - return CycleEqualsExponent(divisor, exponent, ref cache); } - public static bool CycleEqualsExponent(ulong divisor, ulong exponent, ref Dictionary? factorCache) + public void LoadFrom(string path) { - if (CycleEqualsExponentForMersenneCandidate(divisor, exponent)) + EnsureSmallBuffer(); + List<(ulong divisor, ulong cycle)> cycles = []; + using Stream outputStream = new FileStream(path, FileMode.Open, FileAccess.Read, FileShare.Read, BufferSize10M, useAsync: false); + foreach (var (d, c) in EnumerateStream(outputStream)) + { + cycles.Add((d, c)); + if (d <= PerfectNumberConstants.MaxQForDivisorCycles) + { + _smallCycles![(int)d] = c == 0UL ? 1UL : c; + } + } + + cycles.Sort((a, b) => a.divisor < b.divisor ? -1 : (a.divisor > b.divisor ? 1 : 0)); + _table = cycles; + } + + public const int SmallDivisorsMax = PerfectNumberConstants.MaxQForDivisorCycles; + + // Provides a snapshot array of small cycles [0..PerfectNumberConstants.MaxQForDivisorCycles]. Index is the divisor. + // The returned array is a copy safe for use across threads and device uploads. + public ulong[] ExportSmallCyclesSnapshot() + { + EnsureSmallBuffer(); + ulong[] snapshot = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; + Array.Copy(_smallCycles!, snapshot, snapshot.Length); + return snapshot; + } + + public static bool CycleEqualsExponent(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) + { + Dictionary? cache = null; + return CycleEqualsExponent(divisor, exponent, divisorData, ref cache); + } + + public static bool CycleEqualsExponent(ulong divisor, ulong exponent, in MontgomeryDivisorData divisorData, ref Dictionary? factorCache) + { + if (CycleEqualsExponentForMersenneCandidate(divisor, divisorData, exponent)) { return true; } @@ -85,7 +85,7 @@ public class MersenneDivisorCycles } } - if (TryCalculateCycleLengthHeuristic(divisor, out ulong heuristicCycle) && heuristicCycle != 0UL) + if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong heuristicCycle) && heuristicCycle != 0UL) { return heuristicCycle == exponent; } @@ -97,16 +97,15 @@ public class MersenneDivisorCycles factorCache = cache; } - if (TryCalculateCycleLengthForExponent(divisor, exponent, cache, out ulong cycleLength) && cycleLength != 0UL) + if (TryCalculateCycleLengthForExponent(divisor, exponent, divisorData, cache, out ulong cycleLength) && cycleLength != 0UL) { return cycleLength == exponent; } - MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL; } - public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, ulong exponent) + public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent) { if (!IsValidMersenneDivisorCandidate(divisor, exponent)) { @@ -127,7 +126,6 @@ public class MersenneDivisorCycles } } - MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL; } @@ -148,77 +146,77 @@ public class MersenneDivisorCycles return adjusted % step == UInt128.Zero; } - public static IEnumerable<(ulong divisor, ulong cycleLength)> EnumerateStream(Stream compressor) - { - // Binary pairs: (ulong divisor, ulong cycle) - using var reader = new BinaryReader(compressor, Encoding.UTF8, leaveOpen: true); - while (true) - { - ulong d, c; - try - { - d = reader.ReadUInt64(); - c = reader.ReadUInt64(); - } - catch (EndOfStreamException) - { - yield break; - } + public static IEnumerable<(ulong divisor, ulong cycleLength)> EnumerateStream(Stream compressor) + { + // Binary pairs: (ulong divisor, ulong cycle) + using var reader = new BinaryReader(compressor, Encoding.UTF8, leaveOpen: true); + while (true) + { + ulong d, c; + try + { + d = reader.ReadUInt64(); + c = reader.ReadUInt64(); + } + catch (EndOfStreamException) + { + yield break; + } - yield return (d, c); - } - } + yield return (d, c); + } + } - public ulong GetCycle(ulong divisor) - { - // Fast-path: in-memory array for small divisors - if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) - { - var arr = _smallCycles; - if (arr is not null) - { - ulong cached = arr[(int)divisor]; - if (cached != 0UL) - { - return cached; - } - } - } + public ulong GetCycle(ulong divisor, in MontgomeryDivisorData divisorData) + { + // Fast-path: in-memory array for small divisors + if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) + { + var arr = _smallCycles; + if (arr is not null) + { + ulong cached = arr[(int)divisor]; + if (cached != 0UL) + { + return cached; + } + } + } - // Binary search for divisor in sorted _table - int mid, right, left = 0; - ulong d; - List<(ulong divisor, ulong cycleLength)> cycles = _table; + // Binary search for divisor in sorted _table + int mid, right, left = 0; + ulong d; + List<(ulong divisor, ulong cycleLength)> cycles = _table; - right = cycles.Count - 1; - while (left <= right) - { - mid = left + ((right - left) >> 1); - d = cycles[mid].divisor; + right = cycles.Count - 1; + while (left <= right) + { + mid = left + ((right - left) >> 1); + d = cycles[mid].divisor; - if (d == divisor) - { - return cycles[mid].cycleLength; - } + if (d == divisor) + { + return cycles[mid].cycleLength; + } - if (d < divisor) - { - left = mid + 1; - } - else - { - right = mid - 1; - } - } + if (d < divisor) + { + left = mid + 1; + } + else + { + right = mid - 1; + } + } - // TODO: Replace this naive doubling fallback with the unrolled-hex generator from - // MersenneDivisorCycleLengthGpuBenchmarks so large divisors avoid the millions of - // iterations measured in the scalar loop. - // TODO: Route this miss to an ephemeral single-cycle computation on the device selected by the current - // configuration (and skip persisting the result) when the divisor falls outside the in-memory snapshot - // so we respect the large-p memory limits and avoid extra cache locks. - return CalculateCycleLength(divisor); - } + // TODO: Replace this naive doubling fallback with the unrolled-hex generator from + // MersenneDivisorCycleLengthGpuBenchmarks so large divisors avoid the millions of + // iterations measured in the scalar loop. + // TODO: Route this miss to an ephemeral single-cycle computation on the device selected by the current + // configuration (and skip persisting the result) when the divisor falls outside the in-memory snapshot + // so we respect the large-p memory limits and avoid extra cache locks. + return CalculateCycleLength(divisor, divisorData); + } public static UInt128 GetCycle(UInt128 divisor) { @@ -230,7 +228,9 @@ public class MersenneDivisorCycles if (divisor <= ulong.MaxValue) { - ulong cycle = CalculateCycleLength((ulong)divisor); + ulong divisor64 = (ulong)divisor; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor64); + ulong cycle = CalculateCycleLength(divisor64, divisorData); return (UInt128)cycle; } @@ -265,8 +265,9 @@ public class MersenneDivisorCycles } public static bool TryCalculateCycleLengthForExponent( - ulong divisor, + ulong divisor, ulong exponent, + in MontgomeryDivisorData divisorData, Dictionary? factorCache, out ulong cycleLength) { @@ -280,13 +281,14 @@ public class MersenneDivisorCycles if (divisor <= 3UL || exponent <= 1UL) { - cycleLength = CalculateCycleLength(divisor); + cycleLength = CalculateCycleLength(divisor, divisorData); return true; } PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate( divisor, previousOrder: null, + divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); if (orderResult.Order != 0UL) @@ -325,7 +327,7 @@ public class MersenneDivisorCycles return false; } - cycleLength = ReduceOrder(divisor, phi, factorCounts); + cycleLength = ReduceOrder(divisorData, phi, factorCounts); return true; } @@ -467,14 +469,13 @@ public class MersenneDivisorCycles } } - private static ulong ReduceOrder(ulong divisor, ulong initialOrder, Dictionary factorCounts) + private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary factorCounts) { if (factorCounts.Count == 0) { return initialOrder; } - MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); int count = factorCounts.Count; ulong[] primes = ArrayPool.Shared.Rent(count); try @@ -611,376 +612,378 @@ public class MersenneDivisorCycles public int Count { get; } } - public static (long nextPosition, long completeCount) FindLast(string path) - { - using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false); - long count = 0; - foreach (var _ in EnumerateStream(outputStream)) - { - count++; - } - return (outputStream.Position, count); - } + public static (long nextPosition, long completeCount) FindLast(string path) + { + using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false); + long count = 0; + foreach (var _ in EnumerateStream(outputStream)) + { + count++; + } + return (outputStream.Position, count); + } - private const int BufferSize10M = 10 * 1024 * 1024; + private const int BufferSize10M = 10 * 1024 * 1024; - public static void Generate(string path, ulong maxDivisor, int threads = 16) - { - // Binary-only generation of (divisor, cycle) pairs. CSV scaffolding removed. - ulong start = 2; - ulong blockSize = maxDivisor / (ulong)threads + 1UL; - Task[] tasks = new Task[threads]; - using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false); - // Start fresh to avoid mixing formats - outputStream.SetLength(0L); - outputStream.Position = 0L; - using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false); - for (var taskIndex = 0; taskIndex < threads; taskIndex++) - { - ulong threadStart = start + (ulong)taskIndex * blockSize; - if (threadStart > maxDivisor) - { - tasks[taskIndex] = Task.CompletedTask; - continue; - } + public static void Generate(string path, ulong maxDivisor, int threads = 16) + { + // Binary-only generation of (divisor, cycle) pairs. CSV scaffolding removed. + ulong start = 2; + ulong blockSize = maxDivisor / (ulong)threads + 1UL; + Task[] tasks = new Task[threads]; + using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false); + // Start fresh to avoid mixing formats + outputStream.SetLength(0L); + outputStream.Position = 0L; + using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false); + for (var taskIndex = 0; taskIndex < threads; taskIndex++) + { + ulong threadStart = start + (ulong)taskIndex * blockSize; + if (threadStart > maxDivisor) + { + tasks[taskIndex] = Task.CompletedTask; + continue; + } - ulong threadEnd = threadStart + blockSize - 1UL; - if (threadEnd > maxDivisor) - { - threadEnd = maxDivisor; - } + ulong threadEnd = threadStart + blockSize - 1UL; + if (threadEnd > maxDivisor) + { + threadEnd = maxDivisor; + } - ulong localStart = threadStart; - ulong localEnd = threadEnd; - tasks[taskIndex] = Task.Run(() => - { - ulong rangeLength = localEnd - localStart + 1UL; - (ulong divisor, ulong cycle)[] localCycles = ArrayPool<(ulong, ulong)>.Shared.Rent(checked((int)rangeLength)); - int localCycleIndex = 0; + ulong localStart = threadStart; + ulong localEnd = threadEnd; + tasks[taskIndex] = Task.Run(() => + { + ulong rangeLength = localEnd - localStart + 1UL; + (ulong divisor, ulong cycle)[] localCycles = ArrayPool<(ulong, ulong)>.Shared.Rent(checked((int)rangeLength)); + int localCycleIndex = 0; - try - { - for (ulong divisor = localStart; divisor <= localEnd; divisor++) - { - if ((divisor & 1UL) == 0UL) - { - continue; - } - - // TODO: Replace this modulo-based filter with the cached Mod3/Mod5/Mod7/Mod11 - // helpers so the CPU generator stops paying for `%` in the small-prime sieve - // and lines up with the divisor-cycle pipeline. - if ((divisor % 3UL) == 0UL || (divisor % 5UL) == 0UL || (divisor % 7UL) == 0UL || (divisor % 11UL) == 0UL) - { - continue; - } - - // TODO: Migrate this generation path to the shared - // unrolled-hex calculator so the CPU generator matches - // the GPU benchmark leader for large divisors. - ulong cycle = CalculateCycleLength(divisor); - localCycles[localCycleIndex++] = (divisor, cycle); - } - - if (localCycleIndex == 0) - { - return; - } - - lock (writer) - { - writer.BaseStream.Position = writer.BaseStream.Length; - for (int i = 0; i < localCycleIndex; i++) - { - (ulong divisor, ulong cycle) = localCycles[i]; - writer.Write(divisor); - writer.Write(cycle); - } - - writer.Flush(); - } - } - finally - { - ArrayPool<(ulong, ulong)>.Shared.Return(localCycles, clearArray: false); - } - }); - } - - Task.WaitAll(tasks); - } - - public static void GenerateGpu(string path, ulong maxDivisor, int batchSize = 1_000_000, long skipCount = 0L, long nextPosition = 0L) - { - // Prepare output file with header - - ulong start = 3UL; - var locker = new object(); - - var pool = ArrayPool.Shared; - ulong batchSizeUL = (ulong)batchSize, d, end; - int count = (int)Math.Min(batchSizeUL, maxDivisor), i, idx; - ulong[] divisors, outCycles, validDivisors; - - divisors = pool.Rent(count); - outCycles = pool.Rent(count); - var lease = GpuKernelPool.GetKernel(useGpuOrder: true); - var accelerator = lease.Accelerator; - var stream = accelerator.CreateStream(); - var kernel = accelerator.LoadAutoGroupedStreamKernel< - Index1D, - ArrayView1D, - ArrayView1D>(GpuDivisorCycleKernel); - - - using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: true); - if (nextPosition > 0L) - { - outputStream.Position = nextPosition; - } - else - { - outputStream.SetLength(0L); - outputStream.Position = 0L; - } - - using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false); - try - { - while (skipCount > 0L && start <= maxDivisor) + try + { + for (ulong divisor = localStart; divisor <= localEnd; divisor++) + { + if ((divisor & 1UL) == 0UL) { - // TODO: Swap these modulo checks for the cached Mod3/Mod5/Mod7/Mod11 helpers so - // range initialization avoids repeated `%` calls before the unrolled-hex pipeline - // takes over. - if ((start & 1UL) == 0UL || (start % 3UL) == 0UL || (start % 5UL) == 0UL || (start % 7UL) == 0UL || (start % 11UL) == 0UL) - { - start++; - continue; - } - - start++; - skipCount--; - } - - while (start <= maxDivisor) - { - end = Math.Min(start + batchSizeUL - 1UL, maxDivisor); - count = checked((int)(end - start + 1UL)); - - idx = 0; - for (d = start; d <= end; d++) - { - if ((d & 1UL) == 0UL) - continue; - - divisors[idx++] = d; - } - - validDivisors = divisors[..idx]; - - // Use GpuKernelPool to get a kernel and context - // GpuKernelPool.Run((accelerator, stream) => - var bufferDiv = accelerator.Allocate1D(validDivisors); - var bufferCycle = accelerator.Allocate1D(idx); - - kernel( - idx, - bufferDiv.View, - bufferCycle.View); - - accelerator.Synchronize(); - - bufferCycle.View.CopyToCPU(ref outCycles[0], idx); - bufferCycle.Dispose(); - bufferDiv.Dispose(); - - // Collect results in order - for (i = 0; i < idx; i++) - { - d = divisors[i]; - - // TODO: Replace this modulo sieve with the cached Mod helpers once the divisor-cycle cache - // is mandatory so cycle enumeration can drop the slower `%` operations in this hot loop. - if ((d % 3UL) == 0UL || (d % 5UL) == 0UL || (d % 7UL) == 0UL || (d % 11UL) == 0UL) - continue; - - writer.Write(d); - writer.Write(outCycles[i]); - } - - writer.Flush(); - Console.WriteLine($"...processed divisor = {d}"); - start = end + 1UL; - } - } - finally - { - pool.Return(divisors, clearArray: false); - pool.Return(outCycles, clearArray: false); - } - - stream.Dispose(); - lease.Dispose(); - - } - - const byte ByteZero = 0; - const byte ByteOne = 1; - - // GPU kernel for divisor cycle calculation - static void GpuDivisorCycleKernel( - Index1D index, - ArrayView1D divisors, - ArrayView1D outCycles) - { - int i = index.X; - outCycles[i] = CalculateCycleLengthGpu(divisors[i]); - } - - // GPU-friendly version of cycle length calculation - /// - /// GPU-friendly cycle calculator that unrolls sixteen doubling steps; it wins the 8,388,607 benchmark and stays within ~2% of the octo loop at divisor 131,071. - /// - public static ulong CalculateCycleLengthGpu(ulong divisor) - { - if ((divisor & (divisor - 1UL)) == 0UL) - return 1UL; - - ulong order = 1UL; - ulong pow = 2UL; - - while (true) - { - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - - if (GpuStep(ref pow, divisor, ref order)) - return order; - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool GpuStep(ref ulong pow, ulong divisor, ref ulong order) - { - pow += pow; - if (pow >= divisor) - pow -= divisor; - - order++; - return pow == 1UL; - } - - public static ulong CalculateCycleLength(ulong divisor) - { - if (TryCalculateCycleLengthHeuristic(divisor, out ulong cycleLength)) - return cycleLength; - - return CalculateCycleLengthFallback(divisor); - } - - internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, out ulong cycleLength) - { - if ((divisor & (divisor - 1UL)) == 0UL) - { - cycleLength = 1UL; - return true; - } - - if (divisor <= 3UL) - { - cycleLength = CalculateCycleLengthFallback(divisor); - return true; - } - - if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None)) - { - PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate( - divisor, - previousOrder: null, - PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, - PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); - if (orderResult.Order != 0UL) - { - cycleLength = orderResult.Order; - return true; + continue; } + + // TODO: Replace this modulo-based filter with the cached Mod3/Mod5/Mod7/Mod11 + // helpers so the CPU generator stops paying for `%` in the small-prime sieve + // and lines up with the divisor-cycle pipeline. + if ((divisor % 3UL) == 0UL || (divisor % 5UL) == 0UL || (divisor % 7UL) == 0UL || (divisor % 11UL) == 0UL) + { + continue; + } + + // TODO: Migrate this generation path to the shared + // unrolled-hex calculator so the CPU generator matches + // the GPU benchmark leader for large divisors. + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + ulong cycle = CalculateCycleLength(divisor, divisorData); + localCycles[localCycleIndex++] = (divisor, cycle); + } + + if (localCycleIndex == 0) + { + return; + } + + lock (writer) + { + writer.BaseStream.Position = writer.BaseStream.Length; + for (int i = 0; i < localCycleIndex; i++) + { + (ulong divisor, ulong cycle) = localCycles[i]; + writer.Write(divisor); + writer.Write(cycle); + } + + writer.Flush(); + } } - - cycleLength = 0UL; - return false; - } - - private static ulong CalculateCycleLengthFallback(ulong divisor) - { - // Otherwise, find order of 2 mod divisor - // TODO: Switch this scalar fallback to the unrolled-hex stepping sequence once - // the generator is shared with CPU callers; the benchmark shows the unrolled - // variant winning decisively for divisors >= 131,071. - // TODO: Expose a GPU-first branch here so high divisors leverage the ProcessEightBitWindows - // kernel measured fastest in CycleLengthGpuVsCpuBenchmarks, returning the result without - // storing it in the shared cache. - ulong order = 1UL; - ulong pow = 2UL; - while (pow != 1UL) - { - pow <<= 1; - if (pow >= divisor) - pow -= divisor; - - order++; - } - - return order; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private void EnsureSmallBuffer() - { - if (_smallCycles is null) + finally { - _smallCycles = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; + ArrayPool<(ulong, ulong)>.Shared.Return(localCycles, clearArray: false); } + }); } + + Task.WaitAll(tasks); + } + + public static void GenerateGpu(string path, ulong maxDivisor, int batchSize = 1_000_000, long skipCount = 0L, long nextPosition = 0L) + { + // Prepare output file with header + + ulong start = 3UL; + var locker = new object(); + + var pool = ArrayPool.Shared; + ulong batchSizeUL = (ulong)batchSize, d, end; + int count = (int)Math.Min(batchSizeUL, maxDivisor), i, idx; + ulong[] divisors, outCycles, validDivisors; + + divisors = pool.Rent(count); + outCycles = pool.Rent(count); + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var accelerator = lease.Accelerator; + var stream = accelerator.CreateStream(); + var kernel = accelerator.LoadAutoGroupedStreamKernel< + Index1D, + ArrayView1D, + ArrayView1D>(GpuDivisorCycleKernel); + + + using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: true); + if (nextPosition > 0L) + { + outputStream.Position = nextPosition; + } + else + { + outputStream.SetLength(0L); + outputStream.Position = 0L; + } + + using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false); + try + { + while (skipCount > 0L && start <= maxDivisor) + { + // TODO: Swap these modulo checks for the cached Mod3/Mod5/Mod7/Mod11 helpers so + // range initialization avoids repeated `%` calls before the unrolled-hex pipeline + // takes over. + if ((start & 1UL) == 0UL || (start % 3UL) == 0UL || (start % 5UL) == 0UL || (start % 7UL) == 0UL || (start % 11UL) == 0UL) + { + start++; + continue; + } + + start++; + skipCount--; + } + + while (start <= maxDivisor) + { + end = Math.Min(start + batchSizeUL - 1UL, maxDivisor); + count = checked((int)(end - start + 1UL)); + + idx = 0; + for (d = start; d <= end; d++) + { + if ((d & 1UL) == 0UL) + continue; + + divisors[idx++] = d; + } + + validDivisors = divisors[..idx]; + + // Use GpuKernelPool to get a kernel and context + // GpuKernelPool.Run((accelerator, stream) => + var bufferDiv = accelerator.Allocate1D(validDivisors); + var bufferCycle = accelerator.Allocate1D(idx); + + kernel( + idx, + bufferDiv.View, + bufferCycle.View); + + accelerator.Synchronize(); + + bufferCycle.View.CopyToCPU(ref outCycles[0], idx); + bufferCycle.Dispose(); + bufferDiv.Dispose(); + + // Collect results in order + for (i = 0; i < idx; i++) + { + d = divisors[i]; + + // TODO: Replace this modulo sieve with the cached Mod helpers once the divisor-cycle cache + // is mandatory so cycle enumeration can drop the slower `%` operations in this hot loop. + if ((d % 3UL) == 0UL || (d % 5UL) == 0UL || (d % 7UL) == 0UL || (d % 11UL) == 0UL) + continue; + + writer.Write(d); + writer.Write(outCycles[i]); + } + + writer.Flush(); + Console.WriteLine($"...processed divisor = {d}"); + start = end + 1UL; + } + } + finally + { + pool.Return(divisors, clearArray: false); + pool.Return(outCycles, clearArray: false); + } + + stream.Dispose(); + lease.Dispose(); + + } + + const byte ByteZero = 0; + const byte ByteOne = 1; + + // GPU kernel for divisor cycle calculation + static void GpuDivisorCycleKernel( + Index1D index, + ArrayView1D divisors, + ArrayView1D outCycles) + { + int i = index.X; + outCycles[i] = CalculateCycleLengthGpu(divisors[i]); + } + + // GPU-friendly version of cycle length calculation + /// + /// GPU-friendly cycle calculator that unrolls sixteen doubling steps; it wins the 8,388,607 benchmark and stays within ~2% of the octo loop at divisor 131,071. + /// + public static ulong CalculateCycleLengthGpu(ulong divisor) + { + if ((divisor & (divisor - 1UL)) == 0UL) + return 1UL; + + ulong order = 1UL; + ulong pow = 2UL; + + while (true) + { + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + + if (GpuStep(ref pow, divisor, ref order)) + return order; + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool GpuStep(ref ulong pow, ulong divisor, ref ulong order) + { + pow += pow; + if (pow >= divisor) + pow -= divisor; + + order++; + return pow == 1UL; + } + + public static ulong CalculateCycleLength(ulong divisor, in MontgomeryDivisorData divisorData) + { + if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong cycleLength)) + return cycleLength; + + return CalculateCycleLengthFallback(divisor); + } + + internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength) + { + if ((divisor & (divisor - 1UL)) == 0UL) + { + cycleLength = 1UL; + return true; + } + + if (divisor <= 3UL) + { + cycleLength = CalculateCycleLengthFallback(divisor); + return true; + } + + if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None)) + { + PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate( + divisor, + previousOrder: null, + divisorData, + PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, + PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu); + if (orderResult.Order != 0UL) + { + cycleLength = orderResult.Order; + return true; + } + } + + cycleLength = 0UL; + return false; + } + + private static ulong CalculateCycleLengthFallback(ulong divisor) + { + // Otherwise, find order of 2 mod divisor + // TODO: Switch this scalar fallback to the unrolled-hex stepping sequence once + // the generator is shared with CPU callers; the benchmark shows the unrolled + // variant winning decisively for divisors >= 131,071. + // TODO: Expose a GPU-first branch here so high divisors leverage the ProcessEightBitWindows + // kernel measured fastest in CycleLengthGpuVsCpuBenchmarks, returning the result without + // storing it in the shared cache. + ulong order = 1UL; + ulong pow = 2UL; + while (pow != 1UL) + { + pow <<= 1; + if (pow >= divisor) + pow -= divisor; + + order++; + } + + return order; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private void EnsureSmallBuffer() + { + if (_smallCycles is null) + { + _smallCycles = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; + } + } } diff --git a/PerfectNumbers.Core/MontgomeryDivisorData.cs b/PerfectNumbers.Core/MontgomeryDivisorData.cs index 0027afb..95ea18b 100644 --- a/PerfectNumbers.Core/MontgomeryDivisorData.cs +++ b/PerfectNumbers.Core/MontgomeryDivisorData.cs @@ -14,29 +14,6 @@ public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong public readonly ulong MontgomeryTwo = montgomeryTwo; public readonly ulong MontgomeryTwoSquared = montgomeryTwoSquared; -} - -internal static class MontgomeryDivisorDataCache -{ - private static readonly ConcurrentDictionary Cache = new(); - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static MontgomeryDivisorData Get(ulong modulus) => Cache.GetOrAdd(modulus, static m => Create(m)); - - private static MontgomeryDivisorData Create(ulong modulus) - { - ulong nPrime = ComputeMontgomeryNPrime(modulus); - ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus); - ulong montgomeryTwo = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64x2, modulus); - ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime); - - return new MontgomeryDivisorData( - modulus, - nPrime, - montgomeryOne, - montgomeryTwo, - montgomeryTwoSquared); - } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong ComputeMontgomeryResidue(UInt128 value, ulong modulus) => (ulong)(value % modulus); @@ -53,4 +30,28 @@ internal static class MontgomeryDivisorDataCache inv *= unchecked(2UL - modulus * inv); return unchecked(0UL - inv); } + + public static MontgomeryDivisorData FromModulus(ulong modulus) + { + ulong nPrime = ComputeMontgomeryNPrime(modulus); + ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus); + ulong montgomeryTwo = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64x2, modulus); + ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime); + + return new MontgomeryDivisorData( + modulus, + nPrime, + montgomeryOne, + montgomeryTwo, + montgomeryTwoSquared); + } + +} + +internal static class MontgomeryDivisorDataCache +{ + private static readonly ConcurrentDictionary Cache = new(); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static MontgomeryDivisorData Get(ulong modulus) => Cache.GetOrAdd(modulus, static m => MontgomeryDivisorData.FromModulus(m)); } diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index b425b11..4eb6e93 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -7,6 +7,7 @@ public static class PerfectNumberConstants public const ulong ExtraDivisorCycleSearchLimit = 64UL; public const int MaxQForDivisorCycles = 4_000_000; public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; + public const int MaxOddPowersCount = 128; // 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. // TODO: Promote these magic numbers into a runtime profile derived from EvenPerfectBitScanner.Benchmarks so we can retune diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs index 3d618ba..9e1fd27 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs @@ -32,6 +32,7 @@ internal static partial class PrimeOrderCalculator PrimeOrderHeuristicDevice device) { using var scope = UsePow2Mode(device); + MontgomeryDivisorData divisorData; if (prime <= ulong.MaxValue) { ulong? previous = null; @@ -48,14 +49,20 @@ internal static partial class PrimeOrderCalculator } } - PrimeOrderResult result = Calculate((ulong)prime, previous, config, device); + ulong prime64 = (ulong)prime; + divisorData = MontgomeryDivisorData.FromModulus(prime64); + PrimeOrderResult result = Calculate(prime64, previous, divisorData, config, device); return new PrimeOrderResultWide(result.Status, result.Order); } + else + { + divisorData = default; + } - return CalculateWideInternal(prime, previousOrder, config); + return CalculateWideInternal(prime, previousOrder, divisorData, config); } - private static PrimeOrderResultWide CalculateWideInternal(UInt128 prime, UInt128? previousOrder, PrimeOrderSearchConfig config) + private static PrimeOrderResultWide CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config) { if (prime <= UInt128.One) { @@ -72,53 +79,54 @@ internal static partial class PrimeOrderCalculator PartialFactorResult128 phiFactors = PartialFactorWide(phi, config); if (phiFactors.Factors is null) { - return FinishStrictlyWide(prime, config.Mode); + return FinishStrictlyWide(prime, divisorData, config.Mode); } - return RunHeuristicPipelineWide(prime, previousOrder, config, phi, phiFactors); + return RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors); } private static PrimeOrderResultWide RunHeuristicPipelineWide( UInt128 prime, UInt128? previousOrder, + in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config, UInt128 phi, PartialFactorResult128 phiFactors) { - if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, phiFactors)) + if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, divisorData, phiFactors)) { return new PrimeOrderResultWide(PrimeOrderStatus.Found, phi); } - UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi); - candidateOrder = ExponentLoweringWide(candidateOrder, prime, phiFactors); + UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi, divisorData); + candidateOrder = ExponentLoweringWide(candidateOrder, prime, divisorData, phiFactors); - if (TryConfirmOrderWide(prime, candidateOrder, config)) + if (TryConfirmOrderWide(prime, candidateOrder, divisorData, config)) { return new PrimeOrderResultWide(PrimeOrderStatus.Found, candidateOrder); } if (config.Mode == PrimeOrderMode.Strict) { - return FinishStrictlyWide(prime, PrimeOrderMode.Strict); + return FinishStrictlyWide(prime, divisorData, PrimeOrderMode.Strict); } - if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, config, out UInt128 order)) + if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order)) { return new PrimeOrderResultWide(PrimeOrderStatus.Found, order); } - return FinishStrictlyWide(prime, config.Mode); + return FinishStrictlyWide(prime, divisorData, config.Mode); } - private static PrimeOrderResultWide FinishStrictlyWide(UInt128 prime, PrimeOrderMode mode) + private static PrimeOrderResultWide FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData, PrimeOrderMode mode) { - UInt128 strictOrder = CalculateByFactorizationWide(prime); + UInt128 strictOrder = CalculateByFactorizationWide(prime, divisorData); PrimeOrderStatus status = mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved; return new PrimeOrderResultWide(status, strictOrder); } - private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, PartialFactorResult128 factors) + private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors) { ReadOnlySpan factorSpan = factors.Factors; int length = factors.Count; @@ -126,7 +134,7 @@ internal static partial class PrimeOrderCalculator { UInt128 factor = factorSpan[i].Value; UInt128 reduced = phi / factor; - if (Pow2ModWide(reduced, prime) == UInt128.One) + if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One) { return false; } @@ -135,14 +143,14 @@ internal static partial class PrimeOrderCalculator return true; } - private static UInt128 InitializeStartingOrderWide(UInt128 prime, UInt128 phi) + private static UInt128 InitializeStartingOrderWide(UInt128 prime, UInt128 phi, in MontgomeryDivisorData divisorData) { UInt128 order = phi; UInt128 mod8 = prime & (UInt128)7UL; if (mod8 == UInt128.One || mod8 == (UInt128)7UL) { UInt128 half = phi >> 1; - if (Pow2ModWide(half, prime) == UInt128.One) + if (Pow2ModWide(half, prime, divisorData) == UInt128.One) { order = half; } @@ -151,7 +159,7 @@ internal static partial class PrimeOrderCalculator return order; } - private static UInt128 ExponentLoweringWide(UInt128 order, UInt128 prime, PartialFactorResult128 factors) + private static UInt128 ExponentLoweringWide(UInt128 order, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors) { FactorEntry128[]? tempArray = null; try @@ -184,7 +192,7 @@ internal static partial class PrimeOrderCalculator if ((order % primeFactor) == UInt128.Zero) { UInt128 reduced = order / primeFactor; - if (Pow2ModWide(reduced, prime) == UInt128.One) + if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One) { order = reduced; continue; @@ -206,14 +214,14 @@ internal static partial class PrimeOrderCalculator } } - private static bool TryConfirmOrderWide(UInt128 prime, UInt128 order, PrimeOrderSearchConfig config) + private static bool TryConfirmOrderWide(UInt128 prime, UInt128 order, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config) { if (order == UInt128.Zero) { return false; } - if (Pow2ModWide(order, prime) != UInt128.One) + if (Pow2ModWide(order, prime, divisorData) != UInt128.One) { return false; } @@ -253,7 +261,7 @@ internal static partial class PrimeOrderCalculator } reduced /= primeFactor; - if (Pow2ModWide(reduced, prime) == UInt128.One) + if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One) { return false; } @@ -267,6 +275,7 @@ internal static partial class PrimeOrderCalculator UInt128 prime, UInt128 order, UInt128? previousOrder, + in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config, out UInt128 result) { @@ -322,12 +331,12 @@ internal static partial class PrimeOrderCalculator UInt128 candidate = candidates[i]; powUsed++; - if (Pow2ModWide(candidate, prime) != UInt128.One) + if (Pow2ModWide(candidate, prime, divisorData) != UInt128.One) { continue; } - if (!TryConfirmCandidateWide(prime, candidate, config, ref powUsed, powBudget)) + if (!TryConfirmCandidateWide(prime, candidate, config, ref powUsed, powBudget, divisorData)) { continue; } @@ -511,7 +520,7 @@ internal static partial class PrimeOrderCalculator } } - private static bool TryConfirmCandidateWide(UInt128 prime, UInt128 candidate, PrimeOrderSearchConfig config, ref int powUsed, int powBudget) + private static bool TryConfirmCandidateWide(UInt128 prime, UInt128 candidate, PrimeOrderSearchConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData divisorData) { PartialFactorResult128 factorization = PartialFactorWide(candidate, config); if (factorization.Factors is null) @@ -554,7 +563,7 @@ internal static partial class PrimeOrderCalculator } powUsed++; - if (Pow2ModWide(reduced, prime) == UInt128.One) + if (Pow2ModWide(reduced, prime, divisorData) == UInt128.One) { return false; } @@ -826,7 +835,7 @@ internal static partial class PrimeOrderCalculator return 64 + BitOperations.TrailingZeroCount(high); } - private static UInt128 CalculateByFactorizationWide(UInt128 prime) + private static UInt128 CalculateByFactorizationWide(in UInt128 prime, in MontgomeryDivisorData divisorData) { UInt128 phi = prime - UInt128.One; Dictionary counts = new(capacity: 8); @@ -850,7 +859,7 @@ internal static partial class PrimeOrderCalculator if ((order % primeFactor) != UInt128.Zero) { UInt128 candidate = order / primeFactor; - if (Pow2ModWide(candidate, prime) == UInt128.One) + if (Pow2ModWide(candidate, prime, divisorData) == UInt128.One) { order = candidate; } @@ -996,7 +1005,7 @@ internal static partial class PrimeOrderCalculator return true; } - private static UInt128 Pow2ModWide(UInt128 exponent, UInt128 modulus) + private static UInt128 Pow2ModWide(UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData divisorData) { if (modulus == UInt128.One) { @@ -1009,7 +1018,7 @@ internal static partial class PrimeOrderCalculator { ulong prime64 = (ulong)modulus; ulong exponent64 = (ulong)exponent; - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent64, prime64, out ulong remainder); + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent64, prime64, out ulong remainder, divisorData); if (status == GpuPow2ModStatus.Success) { return remainder; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index 7541403..f079c42 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -146,20 +146,21 @@ internal static partial class PrimeOrderCalculator } } - public static PrimeOrderResult Calculate(ulong prime, ulong? previousOrder, PrimeOrderSearchConfig config) - => Calculate(prime, previousOrder, config, PrimeOrderHeuristicDevice.Gpu); + public static PrimeOrderResult Calculate(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config) + => Calculate(prime, previousOrder, divisorData, config, PrimeOrderHeuristicDevice.Gpu); public static PrimeOrderResult Calculate( ulong prime, ulong? previousOrder, + in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config, PrimeOrderHeuristicDevice device) { using var scope = UsePow2Mode(device); - return CalculateInternal(prime, previousOrder, config); + return CalculateInternal(prime, previousOrder, divisorData, config); } - private static PrimeOrderResult CalculateInternal(ulong prime, ulong? previousOrder, PrimeOrderSearchConfig config) + private static PrimeOrderResult CalculateInternal(ulong prime, ulong? previousOrder, MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config) { if (prime <= 3UL) { @@ -167,7 +168,6 @@ internal static partial class PrimeOrderCalculator } ulong phi = prime - 1UL; - MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(prime); if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out PrimeOrderResult gpuResult)) { @@ -448,8 +448,8 @@ internal static partial class PrimeOrderCalculator DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)"); int index = 0; - const int MaxGpuBatchSize = 256; - const int StackGpuBatchSize = 32; + const int MaxGpuBatchSize = 1024; + const int StackGpuBatchSize = 1024; Span stackGpuRemainders = stackalloc ulong[StackGpuBatchSize]; while (index < candidateCount && powUsed < powBudget) { @@ -473,7 +473,7 @@ internal static partial class PrimeOrderCalculator if (batchSize <= StackGpuBatchSize) { Span localRemainders = stackGpuRemainders.Slice(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders); + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData); if (status == GpuPow2ModStatus.Success) { gpuSuccess = true; @@ -484,7 +484,7 @@ internal static partial class PrimeOrderCalculator { gpuPool = ArrayPool.Shared.Rent(batchSize); Span pooledRemainders = gpuPool.AsSpan(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders); + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData); if (status == GpuPow2ModStatus.Success) { pooledGpuRemainders = pooledRemainders; @@ -775,7 +775,7 @@ internal static partial class PrimeOrderCalculator { if (IsGpuPow2Allowed) { - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder); + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData); if (status == GpuPow2ModStatus.Success) { return remainder == 1UL; diff --git a/PerfectNumbers.Core/UInt128Numbers.cs b/PerfectNumbers.Core/UInt128Numbers.cs index bccaa5b..3ed08be 100644 --- a/PerfectNumbers.Core/UInt128Numbers.cs +++ b/PerfectNumbers.Core/UInt128Numbers.cs @@ -1,3 +1,5 @@ +using PerfectNumbers.Core.Gpu; + namespace PerfectNumbers.Core { public static class UInt128Numbers @@ -9,6 +11,8 @@ namespace PerfectNumbers.Core public static readonly UInt128 OneShiftedLeft64 = UInt128.One << 64; public static readonly UInt128 OneShiftedLeft64x2 = Two * OneShiftedLeft64; + public static readonly GpuUInt128 TwoGpu = new(2UL); + // TODO: Populate this table from the UInt128 residue benchmarks so GPU kernels can pull precomputed constants // without re-deriving them at runtime. // TODO: Extend the table with cached Mod3/Mod5 folding constants so both CPU and GPU residue helpers can reuse diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index c374bba..f46b61c 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -12,680 +12,6 @@ public static class ULongExtensions private const ulong Pow2WindowFallbackThreshold = 32UL; private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1); - private struct Pow2OddPowerTable - { - public ulong Element0; - public ulong Element1; - public ulong Element2; - public ulong Element3; - public ulong Element4; - public ulong Element5; - public ulong Element6; - public ulong Element7; - public ulong Element8; - public ulong Element9; - public ulong Element10; - public ulong Element11; - public ulong Element12; - public ulong Element13; - public ulong Element14; - public ulong Element15; - public ulong Element16; - public ulong Element17; - public ulong Element18; - public ulong Element19; - public ulong Element20; - public ulong Element21; - public ulong Element22; - public ulong Element23; - public ulong Element24; - public ulong Element25; - public ulong Element26; - public ulong Element27; - public ulong Element28; - public ulong Element29; - public ulong Element30; - public ulong Element31; - public ulong Element32; - public ulong Element33; - public ulong Element34; - public ulong Element35; - public ulong Element36; - public ulong Element37; - public ulong Element38; - public ulong Element39; - public ulong Element40; - public ulong Element41; - public ulong Element42; - public ulong Element43; - public ulong Element44; - public ulong Element45; - public ulong Element46; - public ulong Element47; - public ulong Element48; - public ulong Element49; - public ulong Element50; - public ulong Element51; - public ulong Element52; - public ulong Element53; - public ulong Element54; - public ulong Element55; - public ulong Element56; - public ulong Element57; - public ulong Element58; - public ulong Element59; - public ulong Element60; - public ulong Element61; - public ulong Element62; - public ulong Element63; - public ulong Element64; - public ulong Element65; - public ulong Element66; - public ulong Element67; - public ulong Element68; - public ulong Element69; - public ulong Element70; - public ulong Element71; - public ulong Element72; - public ulong Element73; - public ulong Element74; - public ulong Element75; - public ulong Element76; - public ulong Element77; - public ulong Element78; - public ulong Element79; - public ulong Element80; - public ulong Element81; - public ulong Element82; - public ulong Element83; - public ulong Element84; - public ulong Element85; - public ulong Element86; - public ulong Element87; - public ulong Element88; - public ulong Element89; - public ulong Element90; - public ulong Element91; - public ulong Element92; - public ulong Element93; - public ulong Element94; - public ulong Element95; - public ulong Element96; - public ulong Element97; - public ulong Element98; - public ulong Element99; - public ulong Element100; - public ulong Element101; - public ulong Element102; - public ulong Element103; - public ulong Element104; - public ulong Element105; - public ulong Element106; - public ulong Element107; - public ulong Element108; - public ulong Element109; - public ulong Element110; - public ulong Element111; - public ulong Element112; - public ulong Element113; - public ulong Element114; - public ulong Element115; - public ulong Element116; - public ulong Element117; - public ulong Element118; - public ulong Element119; - public ulong Element120; - public ulong Element121; - public ulong Element122; - public ulong Element123; - public ulong Element124; - public ulong Element125; - public ulong Element126; - public ulong Element127; - - public ulong this[int index] - { - get - { - return GetElement(index); - } - - set - { - SetElement(index, value); - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private ulong GetElement(int index) => index switch - { - 0 => Element0, - 1 => Element1, - 2 => Element2, - 3 => Element3, - 4 => Element4, - 5 => Element5, - 6 => Element6, - 7 => Element7, - 8 => Element8, - 9 => Element9, - 10 => Element10, - 11 => Element11, - 12 => Element12, - 13 => Element13, - 14 => Element14, - 15 => Element15, - 16 => Element16, - 17 => Element17, - 18 => Element18, - 19 => Element19, - 20 => Element20, - 21 => Element21, - 22 => Element22, - 23 => Element23, - 24 => Element24, - 25 => Element25, - 26 => Element26, - 27 => Element27, - 28 => Element28, - 29 => Element29, - 30 => Element30, - 31 => Element31, - 32 => Element32, - 33 => Element33, - 34 => Element34, - 35 => Element35, - 36 => Element36, - 37 => Element37, - 38 => Element38, - 39 => Element39, - 40 => Element40, - 41 => Element41, - 42 => Element42, - 43 => Element43, - 44 => Element44, - 45 => Element45, - 46 => Element46, - 47 => Element47, - 48 => Element48, - 49 => Element49, - 50 => Element50, - 51 => Element51, - 52 => Element52, - 53 => Element53, - 54 => Element54, - 55 => Element55, - 56 => Element56, - 57 => Element57, - 58 => Element58, - 59 => Element59, - 60 => Element60, - 61 => Element61, - 62 => Element62, - 63 => Element63, - 64 => Element64, - 65 => Element65, - 66 => Element66, - 67 => Element67, - 68 => Element68, - 69 => Element69, - 70 => Element70, - 71 => Element71, - 72 => Element72, - 73 => Element73, - 74 => Element74, - 75 => Element75, - 76 => Element76, - 77 => Element77, - 78 => Element78, - 79 => Element79, - 80 => Element80, - 81 => Element81, - 82 => Element82, - 83 => Element83, - 84 => Element84, - 85 => Element85, - 86 => Element86, - 87 => Element87, - 88 => Element88, - 89 => Element89, - 90 => Element90, - 91 => Element91, - 92 => Element92, - 93 => Element93, - 94 => Element94, - 95 => Element95, - 96 => Element96, - 97 => Element97, - 98 => Element98, - 99 => Element99, - 100 => Element100, - 101 => Element101, - 102 => Element102, - 103 => Element103, - 104 => Element104, - 105 => Element105, - 106 => Element106, - 107 => Element107, - 108 => Element108, - 109 => Element109, - 110 => Element110, - 111 => Element111, - 112 => Element112, - 113 => Element113, - 114 => Element114, - 115 => Element115, - 116 => Element116, - 117 => Element117, - 118 => Element118, - 119 => Element119, - 120 => Element120, - 121 => Element121, - 122 => Element122, - 123 => Element123, - 124 => Element124, - 125 => Element125, - 126 => Element126, - 127 => Element127, - _ => Element0, - }; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private void SetElement(int index, ulong value) - { - switch (index) - { - case 0: - Element0 = value; - return; - case 1: - Element1 = value; - return; - case 2: - Element2 = value; - return; - case 3: - Element3 = value; - return; - case 4: - Element4 = value; - return; - case 5: - Element5 = value; - return; - case 6: - Element6 = value; - return; - case 7: - Element7 = value; - return; - case 8: - Element8 = value; - return; - case 9: - Element9 = value; - return; - case 10: - Element10 = value; - return; - case 11: - Element11 = value; - return; - case 12: - Element12 = value; - return; - case 13: - Element13 = value; - return; - case 14: - Element14 = value; - return; - case 15: - Element15 = value; - return; - case 16: - Element16 = value; - return; - case 17: - Element17 = value; - return; - case 18: - Element18 = value; - return; - case 19: - Element19 = value; - return; - case 20: - Element20 = value; - return; - case 21: - Element21 = value; - return; - case 22: - Element22 = value; - return; - case 23: - Element23 = value; - return; - case 24: - Element24 = value; - return; - case 25: - Element25 = value; - return; - case 26: - Element26 = value; - return; - case 27: - Element27 = value; - return; - case 28: - Element28 = value; - return; - case 29: - Element29 = value; - return; - case 30: - Element30 = value; - return; - case 31: - Element31 = value; - return; - case 32: - Element32 = value; - return; - case 33: - Element33 = value; - return; - case 34: - Element34 = value; - return; - case 35: - Element35 = value; - return; - case 36: - Element36 = value; - return; - case 37: - Element37 = value; - return; - case 38: - Element38 = value; - return; - case 39: - Element39 = value; - return; - case 40: - Element40 = value; - return; - case 41: - Element41 = value; - return; - case 42: - Element42 = value; - return; - case 43: - Element43 = value; - return; - case 44: - Element44 = value; - return; - case 45: - Element45 = value; - return; - case 46: - Element46 = value; - return; - case 47: - Element47 = value; - return; - case 48: - Element48 = value; - return; - case 49: - Element49 = value; - return; - case 50: - Element50 = value; - return; - case 51: - Element51 = value; - return; - case 52: - Element52 = value; - return; - case 53: - Element53 = value; - return; - case 54: - Element54 = value; - return; - case 55: - Element55 = value; - return; - case 56: - Element56 = value; - return; - case 57: - Element57 = value; - return; - case 58: - Element58 = value; - return; - case 59: - Element59 = value; - return; - case 60: - Element60 = value; - return; - case 61: - Element61 = value; - return; - case 62: - Element62 = value; - return; - case 63: - Element63 = value; - return; - case 64: - Element64 = value; - return; - case 65: - Element65 = value; - return; - case 66: - Element66 = value; - return; - case 67: - Element67 = value; - return; - case 68: - Element68 = value; - return; - case 69: - Element69 = value; - return; - case 70: - Element70 = value; - return; - case 71: - Element71 = value; - return; - case 72: - Element72 = value; - return; - case 73: - Element73 = value; - return; - case 74: - Element74 = value; - return; - case 75: - Element75 = value; - return; - case 76: - Element76 = value; - return; - case 77: - Element77 = value; - return; - case 78: - Element78 = value; - return; - case 79: - Element79 = value; - return; - case 80: - Element80 = value; - return; - case 81: - Element81 = value; - return; - case 82: - Element82 = value; - return; - case 83: - Element83 = value; - return; - case 84: - Element84 = value; - return; - case 85: - Element85 = value; - return; - case 86: - Element86 = value; - return; - case 87: - Element87 = value; - return; - case 88: - Element88 = value; - return; - case 89: - Element89 = value; - return; - case 90: - Element90 = value; - return; - case 91: - Element91 = value; - return; - case 92: - Element92 = value; - return; - case 93: - Element93 = value; - return; - case 94: - Element94 = value; - return; - case 95: - Element95 = value; - return; - case 96: - Element96 = value; - return; - case 97: - Element97 = value; - return; - case 98: - Element98 = value; - return; - case 99: - Element99 = value; - return; - case 100: - Element100 = value; - return; - case 101: - Element101 = value; - return; - case 102: - Element102 = value; - return; - case 103: - Element103 = value; - return; - case 104: - Element104 = value; - return; - case 105: - Element105 = value; - return; - case 106: - Element106 = value; - return; - case 107: - Element107 = value; - return; - case 108: - Element108 = value; - return; - case 109: - Element109 = value; - return; - case 110: - Element110 = value; - return; - case 111: - Element111 = value; - return; - case 112: - Element112 = value; - return; - case 113: - Element113 = value; - return; - case 114: - Element114 = value; - return; - case 115: - Element115 = value; - return; - case 116: - Element116 = value; - return; - case 117: - Element117 = value; - return; - case 118: - Element118 = value; - return; - case 119: - Element119 = value; - return; - case 120: - Element120 = value; - return; - case 121: - Element121 = value; - return; - case 122: - Element122 = value; - return; - case 123: - Element123 = value; - return; - case 124: - Element124 = value; - return; - case 125: - Element125 = value; - return; - case 126: - Element126 = value; - return; - case 127: - Element127 = value; - return; - default: - // ILGPU kernels cannot throw exceptions, and callers guarantee the index range. - return; - } - } - } - public static ulong CalculateOrder(this ulong q) { if (q <= 2UL) @@ -1019,8 +345,7 @@ public static class ULongExtensions ulong result = divisor.MontgomeryOne; ulong nPrime = divisor.NPrime; - Pow2OddPowerTable oddPowers = default; - InitializeMontgomeryOddPowers(divisor, modulus, nPrime, ref oddPowers, oddPowerCount); + ulong[] oddPowers = InitializeMontgomeryOddPowers(divisor, modulus, nPrime, oddPowerCount); int index = bitLength - 1; while (index >= 0) @@ -1136,12 +461,13 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void InitializeMontgomeryOddPowers(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, ref Pow2OddPowerTable oddPowers, int oddPowerCount) + private static ulong[] InitializeMontgomeryOddPowers(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount) { + ulong[] oddPowers = new ulong[PerfectNumberConstants.MaxOddPowersCount]; oddPowers[0] = divisor.MontgomeryTwo; if (oddPowerCount == 1) { - return; + return oddPowers; } ulong square = divisor.MontgomeryTwoSquared; @@ -1150,6 +476,8 @@ public static class ULongExtensions ulong previous = oddPowers[i - 1]; oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime); } + + return oddPowers; } [MethodImpl(MethodImplOptions.AggressiveInlining)]