From 604e91b3e3879ececd86c651570c33047a404019 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 1 Oct 2025 00:41:27 +0200 Subject: [PATCH] Return scalar MulMod results and benchmark native reducers --- .../GpuPow2ModBenchmarks.cs | 69 +++ .../GpuUInt128AddBenchmarks.cs | 57 ++ .../GpuUInt128BinaryGcdBenchmarks.cs | 81 +++ .../GpuUInt128Montgomery64Benchmarks.cs | 76 +++ .../GpuUInt128Mul64Benchmarks.cs | 57 ++ .../GpuUInt128MulModBenchmarks.cs | 102 ++++ .../GpuUInt128MulModByLimbBenchmarks.cs | 144 +++++ .../GpuUInt128NativeModuloBenchmarks.cs | 43 ++ .../GpuUInt128ScalarMulModBenchmarks.cs | 57 ++ .../GpuUInt128SubModScalarBenchmarks.cs | 68 +++ .../MulMod64Benchmarks.cs | 6 +- .../Gpu/GpuUInt128Tests.cs | 4 +- PerfectNumbers.Core/Gpu/GpuUInt128.cs | 533 +++++++++--------- .../PerfectNumbers.Core.csproj | 1 + PerfectNumbers.Core/ULongExtensions.cs | 6 +- 15 files changed, 1043 insertions(+), 261 deletions(-) create mode 100644 EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs create mode 100644 EvenPerfectBitScanner.Benchmarks/GpuUInt128AddBenchmarks.cs create mode 100644 EvenPerfectBitScanner.Benchmarks/GpuUInt128BinaryGcdBenchmarks.cs create mode 100644 EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs create mode 100644 EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs create mode 100644 EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs create mode 100644 EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModByLimbBenchmarks.cs create mode 100644 EvenPerfectBitScanner.Benchmarks/GpuUInt128NativeModuloBenchmarks.cs create mode 100644 EvenPerfectBitScanner.Benchmarks/GpuUInt128ScalarMulModBenchmarks.cs create mode 100644 EvenPerfectBitScanner.Benchmarks/GpuUInt128SubModScalarBenchmarks.cs diff --git a/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs new file mode 100644 index 0000000..d618610 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs @@ -0,0 +1,69 @@ +using System.Collections.Generic; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +[SimpleJob(RuntimeMoniker.Net80)] +[MemoryDiagnoser] +public class GpuPow2ModBenchmarks +{ + private static readonly Pow2ModInput[] Inputs = new Pow2ModInput[] + { + new(63UL, new GpuUInt128(0UL, 97UL), "SmallExponentSmallModulus"), + new(1_048_575UL, new GpuUInt128(0UL, 0xFFFF_FFFBUL), "MediumExponentPrimeModulus"), + new(6_710_886_400UL, new GpuUInt128(0x0000_0000_0000_0001UL, 0xFFFF_FFFF_FFFF_FFC3UL), "LargeExponentHighWordModulus"), + new(ulong.MaxValue, new GpuUInt128(0xFFFF_FFFF_FFFF_FFFBUL, 0xFFFF_FFFF_FFFF_FFC5UL), "FullWidthExponentFullWidthModulus"), + }; + + [ParamsSource(nameof(GetInputs))] + public Pow2ModInput Input { get; set; } + + public static IEnumerable GetInputs() => Inputs; + + [Benchmark(Baseline = true)] + public GpuUInt128 ProcessEightBitWindows() + { + return GpuUInt128.Pow2Mod(Input.Exponent, Input.Modulus); + } + + [Benchmark] + public GpuUInt128 ProcessSingleBits() + { + return Pow2ModSingleBit(Input.Exponent, Input.Modulus); + } + + private static GpuUInt128 Pow2ModSingleBit(ulong exponent, in GpuUInt128 modulus) + { + if (modulus.IsZero || modulus == GpuUInt128.One) + { + return new GpuUInt128(); + } + + GpuUInt128 result = new(1UL); + GpuUInt128 baseVal = new(2UL); + + ulong e = exponent; + while (e != 0UL) + { + if ((e & 1UL) != 0UL) + { + result.MulMod(baseVal, modulus); + } + + baseVal.MulMod(baseVal, modulus); + e >>= 1; + } + + return result; + } + + public readonly record struct Pow2ModInput(ulong Exponent, GpuUInt128 Modulus, string Name) + { + public override string ToString() + { + return Name; + } + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128AddBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128AddBenchmarks.cs new file mode 100644 index 0000000..5274cb9 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128AddBenchmarks.cs @@ -0,0 +1,57 @@ +using System.Collections.Generic; +using System.Runtime.CompilerServices; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +[SimpleJob(RuntimeMoniker.Net80)] +[MemoryDiagnoser] +public class GpuUInt128AddBenchmarks +{ + private static readonly AddInput[] Inputs = new AddInput[] + { + new(new GpuUInt128(0UL, 0UL), new GpuUInt128(0UL, 1UL), "IncrementLow"), + new(new GpuUInt128(0UL, ulong.MaxValue), new GpuUInt128(0UL, 1UL), "CarryIntoHigh"), + new(new GpuUInt128(1UL, ulong.MaxValue), new GpuUInt128(3UL, 7UL), "MixedOperands"), + new(new GpuUInt128(ulong.MaxValue, ulong.MaxValue), new GpuUInt128(ulong.MaxValue, ulong.MaxValue), "AllBitsSet"), + }; + + [ParamsSource(nameof(GetInputs))] + public AddInput Input { get; set; } + + public static IEnumerable GetInputs() => Inputs; + + [Benchmark(Baseline = true)] + public GpuUInt128 CarryMaterialisedInLocal() + { + GpuUInt128 value = Input.Left; + value.Add(Input.Right); + return value; + } + + [Benchmark] + public GpuUInt128 CarryInExpression() + { + GpuUInt128 value = Input.Left; + AddInline(ref value, Input.Right); + return value; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void AddInline(ref GpuUInt128 left, GpuUInt128 right) + { + ulong low = left.Low + right.Low; + left.High = left.High + right.High + (low < left.Low ? 1UL : 0UL); + left.Low = low; + } + + public readonly record struct AddInput(GpuUInt128 Left, GpuUInt128 Right, string Name) + { + public override string ToString() + { + return Name; + } + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128BinaryGcdBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128BinaryGcdBenchmarks.cs new file mode 100644 index 0000000..062af47 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128BinaryGcdBenchmarks.cs @@ -0,0 +1,81 @@ +using System.Collections.Generic; +using System.Runtime.CompilerServices; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +[SimpleJob(RuntimeMoniker.Net80)] +[MemoryDiagnoser] +public class GpuUInt128BinaryGcdBenchmarks +{ + private static readonly BinaryGcdInput[] Inputs = new BinaryGcdInput[] + { + new(new GpuUInt128(0UL, 48UL), new GpuUInt128(0UL, 18UL), "SmallOperands"), + new( + new GpuUInt128(0x0123_4567_89AB_CDEFUL, 0xFEDC_BA98_7654_3210UL), + new GpuUInt128(0x0FED_CBA9_8765_4321UL, 0x0123_4567_89AB_CDEFUL), + "HighEntropy"), + new(new GpuUInt128(0UL, ulong.MaxValue), new GpuUInt128(0UL, ulong.MaxValue - 1UL), "LowWordHeavy"), + new(new GpuUInt128(0x8000_0000_0000_0000UL, 0UL), new GpuUInt128(0x7FFF_FFFF_FFFF_FFFFUL, 0UL), "HighWordOnly"), + }; + + [ParamsSource(nameof(GetInputs))] + public BinaryGcdInput Input { get; set; } + + public static IEnumerable GetInputs() => Inputs; + + [Benchmark(Baseline = true)] + public GpuUInt128 ReusedVariables() + { + return GpuUInt128.BinaryGcd(Input.Left, Input.Right); + } + + [Benchmark] + public GpuUInt128 TemporaryStructPerIteration() + { + return BinaryGcdWithTemporary(Input.Left, Input.Right); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static GpuUInt128 BinaryGcdWithTemporary(GpuUInt128 u, GpuUInt128 v) + { + if (u.IsZero) + { + return v; + } + + if (v.IsZero) + { + return u; + } + + int shift = GpuUInt128.TrailingZeroCount(new GpuUInt128(u.High | v.High, u.Low | v.Low)); + int zu = GpuUInt128.TrailingZeroCount(u); + u >>= zu; + + do + { + int zv = GpuUInt128.TrailingZeroCount(v); + v >>= zv; + if (u > v) + { + (u, v) = (v, u); + } + + v -= u; + } + while (!v.IsZero); + + return u << shift; + } + + public readonly record struct BinaryGcdInput(GpuUInt128 Left, GpuUInt128 Right, string Name) + { + public override string ToString() + { + return Name; + } + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs new file mode 100644 index 0000000..d24d3ab --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs @@ -0,0 +1,76 @@ +using System; +using System.Collections.Generic; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core.Gpu; +using PerfectNumbers.Core; + +namespace EvenPerfectBitScanner.Benchmarks; + +[SimpleJob(RuntimeMoniker.Net80)] +[MemoryDiagnoser] +public class GpuUInt128Montgomery64Benchmarks +{ + private static readonly MontgomeryInput[] Inputs = new MontgomeryInput[] + { + Create(0x0000_0000_0000_0065UL, 0x0000_0000_0000_0037UL, 0x0000_0000_0000_007FUL, "TinyOperands"), + Create(0x7FFF_FFFF_FFFF_FFABUL, 0x7FFF_FFFF_FFFF_FFC5UL, 0xFFFF_FFFF_FFFF_FFC5UL, "NearModulus"), + Create(0x0123_4567_89AB_CDEFUL, 0x0FED_CBA9_8765_4321UL, 0xFFFF_FFFF_FFFF_FF4FUL, "MixedMagnitude"), + Create(0xFFFF_FFFF_FFFF_FF01UL, 0xFFFF_FFFF_FFFF_FDCBUL, 0xFFFF_FFFF_FFFF_FFCFUL, "DenseOperands"), + }; + + [ParamsSource(nameof(GetInputs))] + public MontgomeryInput Input { get; set; } + + public static IEnumerable GetInputs() => Inputs; + + [Benchmark(Baseline = true)] + public ulong ExtensionMontgomeryMultiply() + { + return Input.Left.MontgomeryMultiply(Input.Right, Input.Modulus, Input.NPrime); + } + + [Benchmark] + public ulong GpuStructMontgomeryMultiply() + { + GpuUInt128 state = new(Input.Left); + return state.MulModMontgomery64(Input.Right, Input.Modulus, Input.NPrime, Input.R2); + } + + private static MontgomeryInput Create(ulong left, ulong right, ulong modulus, string name) + { + ulong adjustedModulus = (modulus & 1UL) == 0UL ? modulus | 1UL : modulus; + ulong nPrime = ComputeMontgomeryNPrime(adjustedModulus); + ulong r2 = ComputeMontgomeryR2(adjustedModulus); + return new MontgomeryInput(left, right, adjustedModulus, nPrime, r2, name); + } + + private static ulong ComputeMontgomeryR2(ulong modulus) + { + UInt128 r = UInt128.One << 64; + return (ulong)((r * r) % modulus); + } + + private static ulong ComputeMontgomeryNPrime(ulong modulus) + { + ulong inv = modulus; + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + return unchecked(0UL - inv); + } + + public readonly record struct MontgomeryInput( + ulong Left, + ulong Right, + ulong Modulus, + ulong NPrime, + ulong R2, + string Name) + { + public override string ToString() => Name; + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs new file mode 100644 index 0000000..e822560 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs @@ -0,0 +1,57 @@ +using System.Collections.Generic; +using System.Runtime.CompilerServices; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +[SimpleJob(RuntimeMoniker.Net80)] +[MemoryDiagnoser] +public class GpuUInt128Mul64Benchmarks +{ + private static readonly Mul64Input[] Inputs = new Mul64Input[] + { + new(new GpuUInt128(0UL, 1UL), new GpuUInt128(0UL, 1UL), "TinyOperands"), + new(new GpuUInt128(0UL, ulong.MaxValue), new GpuUInt128(0UL, ulong.MaxValue), "LowWordMax"), + new(new GpuUInt128(1UL, 0UL), new GpuUInt128(0UL, ulong.MaxValue), "HighByLow"), + new(new GpuUInt128(ulong.MaxValue, ulong.MaxValue), new GpuUInt128(ulong.MaxValue, ulong.MaxValue), "AllBitsSet"), + }; + + [ParamsSource(nameof(GetInputs))] + public Mul64Input Input { get; set; } + + public static IEnumerable GetInputs() => Inputs; + + [Benchmark(Baseline = true)] + public GpuUInt128 HighProductMaterialisedInLocal() + { + GpuUInt128 value = Input.Multiplicand; + value.Mul64(Input.Multiplier); + return value; + } + + [Benchmark] + public GpuUInt128 HighProductInline() + { + GpuUInt128 value = Input.Multiplicand; + Mul64Inline(ref value, Input.Multiplier); + return value; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void Mul64Inline(ref GpuUInt128 left, GpuUInt128 right) + { + ulong operand = left.Low; + left.Low = operand * right.Low; + left.High = operand * right.High + GpuUInt128.MulHigh(operand, right.Low); + } + + public readonly record struct Mul64Input(GpuUInt128 Multiplicand, GpuUInt128 Multiplier, string Name) + { + public override string ToString() + { + return Name; + } + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs new file mode 100644 index 0000000..6dc0547 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs @@ -0,0 +1,102 @@ +using System.Collections.Generic; +using System.Runtime.CompilerServices; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +[SimpleJob(RuntimeMoniker.Net80)] +[MemoryDiagnoser] +public class GpuUInt128MulModBenchmarks +{ + private static readonly MulModInput[] Inputs = new MulModInput[] + { + new( + new GpuUInt128(0UL, 5UL), + new GpuUInt128(0UL, 7UL), + new GpuUInt128(0UL, 97UL), + "TinyOperands"), + new( + new GpuUInt128(0UL, ulong.MaxValue - 3UL), + new GpuUInt128(0UL, ulong.MaxValue - 7UL), + new GpuUInt128(0UL, ulong.MaxValue - 11UL), + "LowWordHeavy"), + new( + new GpuUInt128(0x1234_5678_9ABC_DEF0UL, 0x0FED_CBA9_8765_4321UL), + new GpuUInt128(0x0ACE_BDF0_1357_9BDFUL, 0x0246_8ACE_ECA8_6421UL), + new GpuUInt128(0x1FFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFF7UL), + "HighWordModulus"), + new( + new GpuUInt128(0x8000_0000_0000_0000UL, 0x0000_0000_0000_0001UL), + new GpuUInt128(0x7FFF_FFFF_FFFF_FFFFUL, 0x8000_0000_0000_0000UL), + new GpuUInt128(0x7FFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFF1UL), + "MixedMagnitude"), + }; + + [ParamsSource(nameof(GetInputs))] + public MulModInput Input { get; set; } + + public static IEnumerable GetInputs() => Inputs; + + [Benchmark(Baseline = true)] + public GpuUInt128 InPlaceMulMod() + { + GpuUInt128 value = Input.Left; + value.MulMod(Input.Right, Input.Modulus); + return value; + } + + [Benchmark] + public GpuUInt128 AllocatePerIteration() + { + return MulModInline(Input.Left, Input.Right, Input.Modulus); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static GpuUInt128 MulModInline(GpuUInt128 left, GpuUInt128 right, GpuUInt128 modulus) + { + GpuUInt128 a = left; + GpuUInt128 b = right; + GpuUInt128 result = new(); + + while (!b.IsZero) + { + if ((b.Low & 1UL) != 0UL) + { + result.Add(a); + if (result.CompareTo(modulus) >= 0) + { + result.Sub(modulus); + } + } + + a <<= 1; + if (a.CompareTo(modulus) >= 0) + { + a.Sub(modulus); + } + + if (b.High == 0UL) + { + b = new GpuUInt128(b.Low >> 1); + } + else + { + ulong low = (b.Low >> 1) | (b.High << 63); + ulong high = b.High >> 1; + b = new GpuUInt128(high, low); + } + } + + return result; + } + + public readonly record struct MulModInput(GpuUInt128 Left, GpuUInt128 Right, GpuUInt128 Modulus, string Name) + { + public override string ToString() + { + return Name; + } + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModByLimbBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModByLimbBenchmarks.cs new file mode 100644 index 0000000..99745a8 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModByLimbBenchmarks.cs @@ -0,0 +1,144 @@ +using System.Collections.Generic; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +[SimpleJob(RuntimeMoniker.Net80)] +[MemoryDiagnoser] +public class GpuUInt128MulModByLimbBenchmarks +{ + private static readonly MulModByLimbInput[] Inputs = new MulModByLimbInput[] + { + new( + new GpuUInt128(0UL, 17UL), + new GpuUInt128(0UL, 19UL), + new GpuUInt128(0x0000_0000_0000_0001UL, 0x0000_0000_0000_00FFUL), + "TinyOperands"), + new( + new GpuUInt128(0x0123_4567_89AB_CDEFUL, 0xFEDC_BA98_7654_3210UL), + new GpuUInt128(0x0F1E_2D3C_4B5A_6978UL, 0x8776_6554_4332_2110UL), + new GpuUInt128(0x0FFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFC3UL), + "MixedMagnitude"), + new( + new GpuUInt128(0x8000_0000_0000_0000UL, 0x0000_0000_0000_0001UL), + new GpuUInt128(0x7FFF_FFFF_FFFF_FFFFUL, 0x8000_0000_0000_0000UL), + new GpuUInt128(0x7FFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFF1UL), + "HighWordDominant"), + }; + + [ParamsSource(nameof(GetInputs))] + public MulModByLimbInput Input { get; set; } + + public static IEnumerable GetInputs() => Inputs; + + [Benchmark(Baseline = true)] + public GpuUInt128 LegacyAllocating() + { + return MulModByLimbLegacy(Input.Left, Input.Right, Input.Modulus); + } + + [Benchmark] + public GpuUInt128 InPlaceReduction() + { + GpuUInt128 value = Input.Left; + value.MulModByLimb(Input.Right, Input.Modulus); + return value; + } + + private static GpuUInt128 MulModByLimbLegacy(GpuUInt128 left, GpuUInt128 right, GpuUInt128 modulus) + { + var (p3, p2, p1, p0) = MultiplyFullLegacy(left, right); + + GpuUInt128 remainder = new(p3, p2); + while (remainder.CompareTo(modulus) >= 0) + { + remainder.Sub(modulus); + } + + ulong limb = p1; + for (int i = 0; i < 2; i++) + { + remainder <<= 64; + remainder = new GpuUInt128(remainder.High, limb); + while (remainder.CompareTo(modulus) >= 0) + { + remainder.Sub(modulus); + } + + limb = p0; + } + + return remainder; + } + + private static (ulong P3, ulong P2, ulong P1, ulong P0) MultiplyFullLegacy(GpuUInt128 left, GpuUInt128 right) + { + var (h0, l0) = Mul64Legacy(left.Low, right.Low); + var (h1, l1) = Mul64Legacy(left.Low, right.High); + var (h2, l2) = Mul64Legacy(left.High, right.Low); + var (h3, l3) = Mul64Legacy(left.High, right.High); + + ulong carry = 0UL; + ulong sum0 = l0; + ulong sum1 = h0; + + sum1 += l1; + if (sum1 < l1) + { + carry++; + } + + sum1 += l2; + if (sum1 < l2) + { + carry++; + } + + ulong sum2 = h1 + h2; + ulong carry2 = sum2 < h2 ? 1UL : 0UL; + sum2 += l3; + if (sum2 < l3) + { + carry2++; + } + + sum2 += carry; + if (sum2 < carry) + { + carry2++; + } + + ulong p1 = sum1; + ulong p2 = sum2; + ulong p3 = h3 + carry2; + + return (p3, p2, p1, sum0); + } + + private static (ulong High, ulong Low) Mul64Legacy(ulong left, ulong right) + { + ulong a0 = (uint)left; + ulong a1 = left >> 32; + ulong b0 = (uint)right; + ulong b1 = right >> 32; + + ulong lo = a0 * b0; + ulong mid1 = a1 * b0; + b0 = a0 * b1; + b1 *= a1; + + a0 = (lo >> 32) + (uint)mid1 + (uint)b0; + a1 = (lo & 0xFFFF_FFFFUL) | (a0 << 32); + b1 += (mid1 >> 32) + (b0 >> 32) + (a0 >> 32); + + return (b1, a1); + } + + public readonly record struct MulModByLimbInput(GpuUInt128 Left, GpuUInt128 Right, GpuUInt128 Modulus, string Name) + { + public override string ToString() => Name; + } +} + diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128NativeModuloBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128NativeModuloBenchmarks.cs new file mode 100644 index 0000000..f15cab9 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128NativeModuloBenchmarks.cs @@ -0,0 +1,43 @@ +using System.Collections.Generic; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +[SimpleJob(RuntimeMoniker.Net80)] +[MemoryDiagnoser] +public class GpuUInt128NativeModuloBenchmarks +{ + private static readonly NativeModuloInput[] Inputs = new NativeModuloInput[] + { + new(0x0000_0000_0000_0003UL, 0x0000_0000_0000_0005UL, 0x0000_0000_0000_007FUL, "TinyOperands"), + new(0x7FFF_FFFF_FFFF_FFC3UL, 0x7FFF_FFFF_FFFF_FF21UL, 0x7FFF_FFFF_FFFF_FFC5UL, "NearModulus"), + new(0xFFFF_FFFF_FFFF_FFFBUL, 0xFFFF_FFFF_FFFF_FFCFUL, 0xFFFF_FFFF_FFFF_FFF5UL, "DenseOperands"), + new(0x0123_4567_89AB_CDEFUL, 0x0FED_CBA9_8765_4321UL, 0x1FFF_FFFF_FFFF_FFFBUL, "MixedMagnitude"), + }; + + [ParamsSource(nameof(GetInputs))] + public NativeModuloInput Input { get; set; } + + public static IEnumerable GetInputs() => Inputs; + + [Benchmark(Baseline = true)] + public ulong ImmediateModulo() + { + GpuUInt128 state = new(Input.Left % Input.Modulus); + return state.MulMod(Input.Right, Input.Modulus); + } + + [Benchmark] + public ulong DeferredNativeModulo() + { + GpuUInt128 state = new(Input.Left); + return state.MulModWithNativeModulo(Input.Right, Input.Modulus); + } + + public readonly record struct NativeModuloInput(ulong Left, ulong Right, ulong Modulus, string Name) + { + public override string ToString() => Name; + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128ScalarMulModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128ScalarMulModBenchmarks.cs new file mode 100644 index 0000000..15a1d23 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128ScalarMulModBenchmarks.cs @@ -0,0 +1,57 @@ +using System.Collections.Generic; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +[SimpleJob(RuntimeMoniker.Net80)] +[MemoryDiagnoser] +public class GpuUInt128ScalarMulModBenchmarks +{ + private static readonly ScalarMulModInput[] Inputs = new ScalarMulModInput[] + { + new( + new GpuUInt128(0x0000_0000_0000_0000UL, 0x0000_0000_0000_0011UL), + 0x0000_0000_0000_0013UL, + new GpuUInt128(0x0000_0000_0000_0001UL, 0x0000_0000_0000_00FFUL), + "TinyOperands"), + new( + new GpuUInt128(0x0123_4567_89AB_CDEFUL, 0xFEDC_BA98_7654_3210UL), + 0x0F1E_2D3C_4B5A_6978UL, + new GpuUInt128(0x0FFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFC3UL), + "MixedMagnitude"), + new( + new GpuUInt128(0x8000_0000_0000_0000UL, 0x0000_0000_0000_0001UL), + 0x7FFF_FFFF_FFFF_FFFBUL, + new GpuUInt128(0x7FFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFF1UL), + "HighWordDominant"), + }; + + [ParamsSource(nameof(GetInputs))] + public ScalarMulModInput Input { get; set; } + + public static IEnumerable GetInputs() => Inputs; + + [Benchmark(Baseline = true)] + public GpuUInt128 StructAllocating() + { + GpuUInt128 value = Input.Left; + value.MulMod(new GpuUInt128(Input.Right), Input.Modulus); + return value; + } + + [Benchmark] + public GpuUInt128 SpecializedScalar() + { + GpuUInt128 value = Input.Left; + value.MulMod(Input.Right, Input.Modulus); + return value; + } + + public readonly record struct ScalarMulModInput(GpuUInt128 Left, ulong Right, GpuUInt128 Modulus, string Name) + { + public override string ToString() => Name; + } +} + diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128SubModScalarBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128SubModScalarBenchmarks.cs new file mode 100644 index 0000000..6258958 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128SubModScalarBenchmarks.cs @@ -0,0 +1,68 @@ +using System.Collections.Generic; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +[SimpleJob(RuntimeMoniker.Net80)] +[MemoryDiagnoser] +public class GpuUInt128SubModScalarBenchmarks +{ + private static readonly SubModInput[] Inputs = new SubModInput[] + { + new(new GpuUInt128(0UL, 3UL), 5UL, new GpuUInt128(0UL, 11UL), "UnderflowAddsModulus"), + new(new GpuUInt128(0UL, ulong.MaxValue - 3UL), 7UL, new GpuUInt128(0UL, ulong.MaxValue - 11UL), "LargeLowWord"), + new(new GpuUInt128(17UL, 2UL), 1UL, new GpuUInt128(3UL, 7UL), "HighWordBorrow"), + }; + + [ParamsSource(nameof(GetInputs))] + public SubModInput Input { get; set; } + + public static IEnumerable GetInputs() => Inputs; + + [Benchmark(Baseline = true)] + public GpuUInt128 CurrentInPlace() + { + GpuUInt128 value = Input.Left; + value.SubMod(Input.Scalar, Input.Modulus); + return value; + } + + [Benchmark] + public GpuUInt128 LegacyTemporaries() + { + GpuUInt128 value = Input.Left; + LegacySubMod(ref value, Input.Scalar, Input.Modulus); + return value; + } + + private static void LegacySubMod(ref GpuUInt128 value, ulong scalar, GpuUInt128 modulus) + { + ulong low; + ulong high; + if (value.High == 0UL && value.Low < scalar) + { + low = value.Low + modulus.Low; + ulong carry = low < value.Low ? 1UL : 0UL; + high = value.High + modulus.High + carry; + } + else + { + high = value.High; + low = value.Low; + } + + ulong borrow = low < scalar ? 1UL : 0UL; + value.High = high - borrow; + value.Low = low - scalar; + } + + public readonly record struct SubModInput(GpuUInt128 Left, ulong Scalar, GpuUInt128 Modulus, string Name) + { + public override string ToString() + { + return Name; + } + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs index 0ac0722..4d181c9 100644 --- a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs @@ -57,8 +57,7 @@ public class MulMod64Benchmarks public ulong GpuCompatibleMulModExtension() { GpuUInt128 gpuUInt128 = new(Input.Left); - gpuUInt128.MulMod(Input.Right, Input.Modulus); - return gpuUInt128.Low; + return gpuUInt128.MulMod(Input.Right, Input.Modulus); } /// @@ -74,8 +73,7 @@ public class MulMod64Benchmarks public ulong GpuCompatibleMulModSimplifiedExtension() { GpuUInt128 gpuUInt128 = new(Input.Left); - gpuUInt128.MulModSimplified(Input.Right, Input.Modulus); - return gpuUInt128.Low; + return gpuUInt128.MulModSimplified(Input.Right, Input.Modulus); } /// diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs index fb2636f..d6778d9 100644 --- a/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs @@ -79,8 +79,8 @@ public class GpuUInt128Tests var a = new GpuUInt128(aVal); var b = new GpuUInt128(bVal); UInt128 expected = ((UInt128)aVal * bVal) % modulus; - a.MulMod(b, modulus); - ((UInt128)a).Should().Be(expected); + ulong result = a.MulMod(b, modulus); + ((UInt128)result).Should().Be(expected); } } diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs index f614eba..ff9cbaa 100644 --- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs +++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs @@ -62,20 +62,15 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public static GpuUInt128 operator +(GpuUInt128 left, GpuUInt128 right) { - var low = left.Low + right.Low; - var high = left.High + right.High + (low < left.Low ? 1UL : 0UL); - // TODO: This should operate and return the instance itself, not a copy. - return new(high, low); + left.Add(right); + return left; } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static GpuUInt128 operator -(GpuUInt128 left, GpuUInt128 right) { - var borrow = left.Low < right.Low ? 1UL : 0UL; - var low = left.Low - right.Low; - var high = left.High - right.High - borrow; - // TODO: This should operate and return the instance itself, not a copy. - return new(high, low); + left.Sub(right); + return left; } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -86,52 +81,61 @@ public struct GpuUInt128 : IComparable, IEquatable public static GpuUInt128 operator *(GpuUInt128 left, GpuUInt128 right) { Multiply(left, right, out var high, out var low); - // TODO: This should operate and return the instance itself, not a copy. - return new(high, low); + left.High = high; + left.Low = low; + return left; } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static GpuUInt128 operator <<(GpuUInt128 value, int shift) { - // TODO: Reuse variables to reduce register pressure. shift &= 127; if (shift == 0) { return value; } + ulong high = value.High; + ulong low = value.Low; + if (shift >= 64) { - ulong high = value.Low << (shift - 64); - return new(high, 0UL); + int longShift = shift - 64; + value.High = low << longShift; + value.Low = 0UL; + return value; } - ulong highPart = (value.High << shift) | (value.Low >> (64 - shift)); - ulong lowPart = value.Low << shift; - // TODO: This should operate and return the instance itself, not a copy. - return new(highPart, lowPart); + int inverseShift = 64 - shift; + value.High = (high << shift) | (low >> inverseShift); + value.Low = low << shift; + return value; } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static GpuUInt128 operator >>(GpuUInt128 value, int shift) { - // TODO: Reuse variables to reduce register pressure. shift &= 127; if (shift == 0) { return value; } + ulong high = value.High; + ulong low = value.Low; + if (shift >= 64) { - ulong low = value.High >> (shift - 64); - return new(low); + int longShift = shift - 64; + value.Low = longShift == 0 ? high : high >> longShift; + value.High = 0UL; + return value; } - ulong highPart = value.High >> shift; - ulong lowPart = (value.Low >> shift) | (value.High << (64 - shift)); - // TODO: This should operate and return the instance itself, not a copy. - return new(highPart, lowPart); + int inverseShift = 64 - shift; + value.Low = (low >> shift) | (high << inverseShift); + value.High = high >> shift; + return value; } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -197,9 +201,10 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public void Add(GpuUInt128 other) { - ulong low = Low + other.Low; - // TODO: Let's add a version and a benchmark where we assign (low < Low ? 1UL : 0UL) to a local, first. - High = High + other.High + (low < Low ? 1UL : 0UL); + ulong originalLow = Low; + ulong low = originalLow + other.Low; + ulong carry = low < originalLow ? 1UL : 0UL; + High = High + other.High + carry; Low = low; } @@ -405,23 +410,29 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MulHigh(ulong x, ulong y) { - // TODO: Reuse variables to reduce register pressure. ulong xLow = (uint)x; - ulong xHigh = x >> 32; ulong yLow = (uint)y; - ulong yHigh = y >> 32; + ulong lowProduct = xLow * yLow; - ulong w1 = xLow * yHigh; - ulong w2 = xHigh * yLow; - ulong w3 = xLow * yLow; + ulong yHigh = y >> 32; + ulong cross = xLow * yHigh; + + ulong xHigh = x >> 32; + ulong result = xHigh * yHigh; + ulong temp = xHigh * yLow; // 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. - ulong result = (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32); - result += ((w3 >> 32) + (uint)w1 + (uint)w2) >> 32; + result += cross >> 32; + result += temp >> 32; + + lowProduct >>= 32; + lowProduct += (uint)cross; + lowProduct += (uint)temp; + result += lowProduct >> 32; return result; } @@ -430,59 +441,35 @@ public struct GpuUInt128 : IComparable, IEquatable { // Multiply this.Low (assumed 64-bit value) by full 128-bit other ulong operand = Low; - // otherLow = other.Low; - Low = operand * other.Low; - // Let's try a version and add a benchmark for it where operand * other.High is assigned to a local, first. - High = operand * other.High + MulHigh(operand, other.Low); + ulong highProduct = operand * other.High; + High = highProduct + MulHigh(operand, other.Low); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static GpuUInt128 AddMod(GpuUInt128 a, GpuUInt128 b, GpuUInt128 modulus) + public void MulMod(in GpuUInt128 value, in GpuUInt128 modulus) { - a.Add(b); - if (a >= modulus) + GpuUInt128 multiplicand = this; + GpuUInt128 multiplier = value; + + High = 0UL; + Low = 0UL; + + while (!multiplier.IsZero) { - a.Sub(modulus); + if ((multiplier.Low & 1UL) != 0UL) + { + AddMod(multiplicand, modulus); + } + + multiplicand.ShiftLeft(1); + if (multiplicand.CompareTo(modulus) >= 0) + { + multiplicand.Sub(modulus); + } + + multiplier.ShiftRight(1); } - - // TODO: Don't return anything. Return void. - return a; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static GpuUInt128 MulMod(GpuUInt128 a, GpuUInt128 b, GpuUInt128 modulus) - { - // TODO: This should operate on the instance itself, not on a copy. Avoid creating new instances anywhere. We should return void. - GpuUInt128 result = new(); - while (!b.IsZero) - { - if ((b.Low & 1UL) != 0UL) - { - result = AddMod(result, a, modulus); - } - - // a = (a << 1) % modulus - a <<= 1; - if (a >= modulus) - { - a.Sub(modulus); - } - - // b >>= 1 - if (b.High == 0UL) - { - b = new GpuUInt128(b.Low >> 1); - } - else - { - ulong low = (b.Low >> 1) | (b.High << 63); - ulong high = b.High >> 1; - b = new GpuUInt128(high, low); - } - } - - return result; } @@ -494,20 +481,31 @@ public struct GpuUInt128 : IComparable, IEquatable return new(); } + const int WindowSizeBits = 8; + const ulong WindowMask = (1UL << WindowSizeBits) - 1UL; + GpuUInt128 result = new(1UL); - GpuUInt128 baseVal = new(2UL); + GpuUInt128 baseVal = new(2UL); ulong e = exponent; - // TODO: Can we modify this loop to process multiple bits at a time? E.g. 64-bit chunks. while (e != 0UL) { - if ((e & 1UL) != 0UL) + int remainingBits = BitOperations.Log2(e) + 1; + int bitsToProcess = remainingBits >= WindowSizeBits ? WindowSizeBits : remainingBits; + ulong chunk = e & (WindowMask >> (WindowSizeBits - bitsToProcess)); + + for (int bitIndex = 0; bitIndex < bitsToProcess; bitIndex++) { - result = MulMod(result, baseVal, modulus); + if ((chunk & 1UL) != 0UL) + { + result.MulMod(baseVal, modulus); + } + + baseVal.MulMod(baseVal, modulus); + chunk >>= 1; } - baseVal = MulMod(baseVal, baseVal, modulus); - e >>= 1; + e >>= bitsToProcess; } return result; @@ -527,7 +525,6 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public static GpuUInt128 BinaryGcd(GpuUInt128 u, GpuUInt128 v) { - // TODO: Reuse variables to reduce register pressure. if (u.IsZero) { return v; @@ -538,8 +535,12 @@ public struct GpuUInt128 : IComparable, IEquatable return u; } - int shift = TrailingZeroCount(new GpuUInt128(u.High | v.High, u.Low | v.Low)); - // remove factors of 2 + ulong combinedLow = u.Low | v.Low; + ulong combinedHigh = u.High | v.High; + int shift = combinedLow != 0UL + ? BitOperations.TrailingZeroCount(combinedLow) + : 64 + BitOperations.TrailingZeroCount(combinedHigh); + int zu = TrailingZeroCount(u); u >>= zu; @@ -584,18 +585,12 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public void SubMod(ulong value, GpuUInt128 modulus) { - // TODO: Reuse variables to reduce register pressure. - ulong low, high; - if (High == 0UL && Low < value) + ulong low = Low; + ulong high = High; + if (high == 0UL && low < value) { - low = Low + modulus.Low; - ulong carry = low < Low ? 1UL : 0UL; - high = High + modulus.High + carry; - } - else - { - high = High; - low = Low; + low += modulus.Low; + high = modulus.High + (low < modulus.Low ? 1UL : 0UL); } ulong borrow = low < value ? 1UL : 0UL; @@ -635,73 +630,62 @@ public struct GpuUInt128 : IComparable, IEquatable [MethodImpl(MethodImplOptions.AggressiveInlining)] public void ShiftLeft(int shift) { - // TODO: This should operate on the instance itself, not on a copy. Don't use operator. Implement your own logic. - this <<= shift; + shift &= 127; + if (shift == 0) + { + return; + } + + ulong high = High; + ulong low = Low; + + if (shift >= 64) + { + High = low << (shift - 64); + Low = 0UL; + return; + } + + int inverseShift = 64 - shift; + High = (high << shift) | (low >> inverseShift); + Low = low << shift; } [MethodImpl(MethodImplOptions.AggressiveInlining)] public void ShiftRight(int shift) { - // TODO: This should operate on the instance itself, not on a copy. Don't use operator. Implement your own logic. - this >>= shift; + shift &= 127; + if (shift == 0) + { + return; + } + + ulong high = High; + ulong low = Low; + + if (shift >= 64) + { + Low = shift == 64 ? high : high >> (shift - 64); + High = 0UL; + return; + } + + int inverseShift = 64 - shift; + Low = (low >> shift) | (high << inverseShift); + High = high >> shift; } public void MulMod(GpuUInt128 other, GpuUInt128 modulus) { if (modulus.High == 0UL) { - MulMod(other, modulus.Low); + Low = MulMod(other, modulus.Low); + High = 0UL; return; } MultiplyFull(this, other, out var p3, out var p2, out var p1, out var p0); - - // TODO: This should operate on the instance itself, not on a copy. Avoid creating new instances anywhere. - GpuUInt128 remainder = new(); - int bit; - // TODO: Can we modify these loops to process multiple bits at a time? E.g. 64-bit chunks. - for (bit = 63; bit >= 0; bit--) - { - remainder <<= 1; - remainder = new(remainder.High, remainder.Low | ((p3 >> bit) & 1UL)); - if (remainder.CompareTo(modulus) >= 0) - { - remainder.Sub(modulus); - } - } - - for (bit = 63; bit >= 0; bit--) - { - remainder <<= 1; - remainder = new(remainder.High, remainder.Low | ((p2 >> bit) & 1UL)); - if (remainder.CompareTo(modulus) >= 0) - { - remainder.Sub(modulus); - } - } - - for (bit = 63; bit >= 0; bit--) - { - remainder <<= 1; - remainder = new(remainder.High, remainder.Low | ((p1 >> bit) & 1UL)); - if (remainder.CompareTo(modulus) >= 0) - { - remainder.Sub(modulus); - } - } - - for (bit = 63; bit >= 0; bit--) - { - remainder <<= 1; - remainder = new(remainder.High, remainder.Low | ((p0 >> bit) & 1UL)); - if (remainder.CompareTo(modulus) >= 0) - { - remainder.Sub(modulus); - } - } - - High = remainder.High; - Low = remainder.Low; + ReduceProductBitwise(p3, p2, p1, p0, modulus); // TODO(MOD-OPT): Replace this bitwise long-division style reduction // with a faster algorithm suitable for GPU kernels without '%' support. // Options: @@ -741,28 +725,26 @@ public struct GpuUInt128 : IComparable, IEquatable { MultiplyFull(this, other, out var p3, out var p2, out var p1, out var p0); - // TODO: This should operate on the instance itself, not on a copy. Avoid creating new instances anywhere. - GpuUInt128 remainder = new(p3, p2); - while (remainder.CompareTo(modulus) >= 0) + High = p3; + Low = p2; + while (CompareTo(modulus) >= 0) { - remainder.Sub(modulus); + Sub(modulus); } ulong limb = p1; for (int i = 0; i < 2; i++) { - remainder <<= 64; - remainder = new(remainder.High, limb); - while (remainder.CompareTo(modulus) >= 0) + ShiftLeft(64); + Low = limb; + while (CompareTo(modulus) >= 0) { - remainder.Sub(modulus); + Sub(modulus); } limb = p0; } - High = remainder.High; - Low = remainder.Low; } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -770,110 +752,149 @@ public struct GpuUInt128 : IComparable, IEquatable { if (modulus.High == 0UL) { - MulMod(value, modulus.Low); + Low = MulMod(value, modulus.Low); + High = 0UL; return; } - // TODO: We shouldn't need a new instance to multiply by value. - GpuUInt128 other = new(value); - MulMod(other, modulus); + MultiplyFull(this, value, out var p3, out var p2, out var p1, out var p0); + ReduceProductBitwise(p3, p2, p1, p0, modulus); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void MulMod(GpuUInt128 other, ulong modulus) + public ulong MulMod(GpuUInt128 other, ulong modulus) { - // TODO: Since we always set High = 0, we should return ulong for low instead of modifying the instance. - High = 0UL; - ulong a = Low % modulus; - ulong b = other.Low % modulus; + ulong modulusLocal = modulus; + ulong a = Low % modulusLocal; + ulong b = other.Low % modulusLocal; + if (a == 0UL || b == 0UL) { - Low = 0UL; - } - else if (b <= ulong.MaxValue / a) - { - Low = (a * b) % modulus; - } - else - { - Low = MulMod64(a, b, modulus); + return 0UL; } + + return b <= ulong.MaxValue / a + ? (a * b) % modulusLocal + : MulMod64(a, b, modulusLocal); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void MulMod(ulong value, ulong modulus) + public ulong MulMod(ulong value, ulong modulus) { - // TODO: Since we always set High = 0, we should return ulong for low instead of modifying the instance. - High = 0UL; + ulong modulusLocal = modulus; + ulong a = Low % modulusLocal; + ulong b = value % modulusLocal; - ulong a = Low % modulus; - ulong b = value % modulus; if (a == 0UL || b == 0UL) { - Low = 0UL; - } - else if (b <= ulong.MaxValue / a) - { - Low = (a * b) % modulus; - } - else - { - Low = MulMod64(a, b, modulus); + return 0UL; } + + return b <= ulong.MaxValue / a + ? (a * b) % modulusLocal + : MulMod64(a, b, modulusLocal); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void MulModSimplified(ulong value, ulong modulus) + public ulong MulModSimplified(ulong value, ulong modulus) { - // TODO: Since we always set High = 0, we should return ulong for low instead of modifying the instance. - High = 0UL; - ulong a = Low % modulus; - ulong b = value % modulus; + ulong modulusLocal = modulus; + ulong a = Low % modulusLocal; + ulong b = value % modulusLocal; + if (a == 0UL || b == 0UL) { - Low = 0UL; - } - else - { - Low = MulMod64(a, b, modulus); + return 0UL; } + + return MulMod64(a, b, modulusLocal); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void MulModWithNativeModulo(ulong value, ulong modulus) + public ulong MulModWithNativeModulo(ulong value, ulong modulus) { - // TODO: Since we always set High = 0, we should return ulong for low instead of modifying the instance. - // TODO: Reuse variables to reduce register pressure. - High = 0UL; + ulong modulusLocal = modulus; + ulong multiplicand = Low % modulusLocal; + ulong multiplier = value % modulusLocal; - ulong a = Low % modulus; - ulong b = value % modulus; - if (a == 0UL || b == 0UL) + if (multiplicand == 0UL || multiplier == 0UL) { - Low = 0UL; - return; + return 0UL; } ulong result = 0UL; + ulong current = multiplicand; + ulong remaining = multiplier; + while (true) { - ulong chunk = b & NativeModuloChunkMask; + ulong chunk = remaining & NativeModuloChunkMask; if (chunk != 0UL) { - ulong chunkContribution = MultiplyChunkModulo(a, chunk, modulus); - result = (result + chunkContribution) % modulus; + ulong chunkContribution = MultiplyChunkModulo(current, chunk, modulusLocal); + result = (result + chunkContribution) % modulusLocal; } - b >>= NativeModuloChunkBits; - if (b == 0UL) + remaining >>= NativeModuloChunkBits; + if (remaining == 0UL) { break; } - a = ShiftLeftByNativeChunk(a, modulus); + current = ShiftLeftByNativeChunk(current, modulusLocal); } - Low = result; + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private void ReduceProductBitwise(ulong p3, ulong p2, ulong p1, ulong p0, GpuUInt128 modulus) + { + GpuUInt128 remainder = new(); + int bit; + // TODO: Can we modify these loops to process multiple bits at a time? E.g. 64-bit chunks. + for (bit = 63; bit >= 0; bit--) + { + remainder.ShiftLeft(1); + remainder.Low |= (p3 >> bit) & 1UL; + if (remainder.CompareTo(modulus) >= 0) + { + remainder.Sub(modulus); + } + } + + for (bit = 63; bit >= 0; bit--) + { + remainder.ShiftLeft(1); + remainder.Low |= (p2 >> bit) & 1UL; + if (remainder.CompareTo(modulus) >= 0) + { + remainder.Sub(modulus); + } + } + + for (bit = 63; bit >= 0; bit--) + { + remainder.ShiftLeft(1); + remainder.Low |= (p1 >> bit) & 1UL; + if (remainder.CompareTo(modulus) >= 0) + { + remainder.Sub(modulus); + } + } + + for (bit = 63; bit >= 0; bit--) + { + remainder.ShiftLeft(1); + remainder.Low |= (p0 >> bit) & 1UL; + if (remainder.CompareTo(modulus) >= 0) + { + remainder.Sub(modulus); + } + } + + High = remainder.High; + Low = remainder.Low; } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -915,47 +936,39 @@ public struct GpuUInt128 : IComparable, IEquatable } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void MulModMontgomery64(GpuUInt128 other, ulong modulus, ulong nPrime, ulong r2) + public ulong MulModMontgomery64(GpuUInt128 other, ulong modulus, ulong nPrime, ulong r2) { - // TODO: Since we always set High = 0, we should return ulong for low instead of modifying the instance. - // TODO: Reuse variables to reduce register pressure. - ulong a = Low % modulus; - ulong b = other.Low % modulus; + ulong modulusLocal = modulus; + ulong a = Low % modulusLocal; + ulong b = other.Low % modulusLocal; + if (a == 0UL || b == 0UL) { - High = 0UL; - Low = 0UL; - return; + return 0UL; } - ulong aR = MontMul64(a, r2, modulus, nPrime); - ulong bR = MontMul64(b, r2, modulus, nPrime); - ulong cR = MontMul64(aR, bR, modulus, nPrime); - ulong c = MontMul64(cR, 1UL, modulus, nPrime); - High = 0UL; - Low = c; + ulong aR = MontMul64(a, r2, modulusLocal, nPrime); + ulong bR = MontMul64(b, r2, modulusLocal, nPrime); + ulong cR = MontMul64(aR, bR, modulusLocal, nPrime); + return MontMul64(cR, 1UL, modulusLocal, nPrime); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void MulModMontgomery64(ulong other, ulong modulus, ulong nPrime, ulong r2) + public ulong MulModMontgomery64(ulong other, ulong modulus, ulong nPrime, ulong r2) { - // TODO: Since we always set High = 0, we should return ulong for low instead of modifying the instance. - // TODO: Reuse variables to reduce register pressure. - ulong a = Low % modulus; - ulong b = other % modulus; + ulong modulusLocal = modulus; + ulong a = Low % modulusLocal; + ulong b = other % modulusLocal; + if (a == 0UL || b == 0UL) { - High = 0UL; - Low = 0UL; - return; + return 0UL; } - ulong aR = MontMul64(a, r2, modulus, nPrime); - ulong bR = MontMul64(b, r2, modulus, nPrime); - ulong cR = MontMul64(aR, bR, modulus, nPrime); - ulong c = MontMul64(cR, 1UL, modulus, nPrime); - High = 0UL; - Low = c; + ulong aR = MontMul64(a, r2, modulusLocal, nPrime); + ulong bR = MontMul64(b, r2, modulusLocal, nPrime); + ulong cR = MontMul64(aR, bR, modulusLocal, nPrime); + return MontMul64(cR, 1UL, modulusLocal, nPrime); } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -963,7 +976,8 @@ public struct GpuUInt128 : IComparable, IEquatable { if (modulus.High == 0UL) { - MulMod(Low, modulus.Low); + Low = MulMod(Low, modulus.Low); + High = 0UL; return; } @@ -1130,6 +1144,23 @@ public struct GpuUInt128 : IComparable, IEquatable p3 = h3 + carry2; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void MultiplyFull(GpuUInt128 left, ulong right, out ulong p3, out ulong p2, out ulong p1, out ulong p0) + { + var (highLow, lowLow) = Mul64(left.Low, right); + var (highHigh, lowHigh) = Mul64(left.High, right); + + p0 = lowLow; + + ulong mid = highLow + lowHigh; + ulong carry = mid < highLow ? 1UL : 0UL; + p1 = mid; + + ulong upper = highHigh + carry; + p3 = upper < highHigh ? 1UL : 0UL; + p2 = upper; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void Multiply(GpuUInt128 left, GpuUInt128 right, out ulong high, out ulong low) { diff --git a/PerfectNumbers.Core/PerfectNumbers.Core.csproj b/PerfectNumbers.Core/PerfectNumbers.Core.csproj index f4ef795..fa27778 100644 --- a/PerfectNumbers.Core/PerfectNumbers.Core.csproj +++ b/PerfectNumbers.Core/PerfectNumbers.Core.csproj @@ -17,6 +17,7 @@ + diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 291d4b3..a2eee9b 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -215,16 +215,14 @@ public static class ULongExtensions public static ulong MulMod64GpuCompatible(this ulong a, ulong b, ulong modulus) { GpuUInt128 state = new(a % modulus); - state.MulMod(b, modulus); - return state.Low; + return state.MulMod(b, modulus); } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MulMod64GpuCompatibleDeferred(this ulong a, ulong b, ulong modulus) { GpuUInt128 state = new(a); - state.MulModWithNativeModulo(b, modulus); - return state.Low; + return state.MulModWithNativeModulo(b, modulus); } [MethodImpl(MethodImplOptions.AggressiveInlining)]