diff --git a/EvenPerfectBitScanner.Benchmarks/GpuRationalBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuRationalBenchmarks.cs new file mode 100644 index 0000000..744aed8 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/GpuRationalBenchmarks.cs @@ -0,0 +1,29 @@ +using BenchmarkDotNet.Attributes; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +[MemoryDiagnoser] +public class GpuRationalBenchmarks +{ + private static readonly GpuRational Left = new(new GpuUInt128(3UL), new GpuUInt128(4UL)); + private static readonly GpuRational Right = new(new GpuUInt128(5UL), new GpuUInt128(6UL)); + private static readonly GpuUInt128 IntegerLeft = new(3UL); + private static readonly GpuUInt128 IntegerRight = new(5UL); + + [Benchmark] + public GpuRational AddRational() + { + GpuRational value = Left; + value = value + Right; + return value; + } + + [Benchmark] + public GpuUInt128 AddInteger() + { + GpuUInt128 value = IntegerLeft; + value.Add(IntegerRight); + return value; + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs new file mode 100644 index 0000000..440cd96 --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs @@ -0,0 +1,283 @@ +using System; +using System.Buffers.Binary; +using System.Numerics; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core; +using PerfectNumbers.Core.Gpu; + +namespace EvenPerfectBitScanner.Benchmarks; + +[MemoryDiagnoser] +[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)] +public class Pow2MontgomeryModUInt128Benchmarks +{ + private const int SampleCount = 256; + private const ulong Pow2WindowFallbackThreshold = 32UL; + + private readonly UInt128[] _moduli = new UInt128[SampleCount]; + private readonly UInt128[] _exponents = new UInt128[SampleCount]; + + private readonly Random _random = new(19); + + [GlobalSetup] + public void Setup() + { + for (int i = 0; i < SampleCount; i++) + { + UInt128 modulus = NextWideOddModulus(); + _moduli[i] = modulus; + + UInt128 exponent = NextWideExponent(); + _exponents[i] = exponent; + } + } + + /// + /// Current UInt128-native windowed Montgomery ladder that avoids GpuUInt128 on CPU paths. + /// + [Benchmark(Baseline = true)] + public UInt128 UInt128Windowed() + { + UInt128 checksum = UInt128.Zero; + + for (int i = 0; i < SampleCount; i++) + { + checksum ^= _exponents[i].Pow2MontgomeryModWindowed(_moduli[i]); + } + + return checksum; + } + + /// + /// Previous CPU shim that reused the GPU-oriented GpuUInt128 implementation. + /// + [Benchmark] + public UInt128 GpuStructShim() + { + UInt128 checksum = UInt128.Zero; + + for (int i = 0; i < SampleCount; i++) + { + checksum ^= Pow2MontgomeryModWindowedGpuShim(_exponents[i], _moduli[i]); + } + + return checksum; + } + + private UInt128 NextWideOddModulus() + { + Span buffer = stackalloc byte[16]; + _random.NextBytes(buffer); + + ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer); + ulong high = BinaryPrimitives.ReadUInt64LittleEndian(buffer[8..]); + if (high == 0UL) + { + high = 1UL; + } + + UInt128 modulus = ((UInt128)high << 64) | low; + modulus |= UInt128.One; + + if (modulus <= UInt128.One) + { + modulus += 2UL; + } + + return modulus; + } + + private UInt128 NextWideExponent() + { + Span buffer = stackalloc byte[16]; + _random.NextBytes(buffer); + + ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer); + ulong high = BinaryPrimitives.ReadUInt64LittleEndian(buffer[8..]); + if (high == 0UL) + { + high = 1UL; + } + + return ((UInt128)high << 64) | low; + } + + private static UInt128 Pow2MontgomeryModWindowedGpuShim(UInt128 exponent, UInt128 modulus) + { + if (modulus == UInt128.One) + { + return UInt128.Zero; + } + + if (exponent == UInt128.Zero) + { + return UInt128.One % modulus; + } + + GpuUInt128 modulusGpu = new(modulus); + GpuUInt128 exponentGpu = new(exponent); + GpuUInt128 baseValue = new(2UL); + if (baseValue.CompareTo(modulusGpu) >= 0) + { + baseValue.Sub(modulusGpu); + } + + if (ShouldUseSingleBit(exponentGpu)) + { + GpuUInt128 singleBitResult = Pow2MontgomeryModSingleBit(exponentGpu, modulusGpu, baseValue); + return (UInt128)singleBitResult; + } + + int bitLength = exponentGpu.GetBitLength(); + int windowSize = GetWindowSize(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + Span oddPowers = stackalloc GpuUInt128[oddPowerCount]; + InitializeOddPowers(baseValue, modulusGpu, oddPowers); + + GpuUInt128 result = GpuUInt128.One; + int index = bitLength - 1; + + while (index >= 0) + { + if (!IsBitSet(exponentGpu, index)) + { + result.MulMod(result, modulusGpu); + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (!IsBitSet(exponentGpu, windowStart)) + { + windowStart++; + } + + int windowBitCount = index - windowStart + 1; + for (int square = 0; square < windowBitCount; square++) + { + result.MulMod(result, modulusGpu); + } + + ulong windowValue = ExtractWindowValue(exponentGpu, windowStart, windowBitCount); + int tableIndex = (int)((windowValue - 1UL) >> 1); + GpuUInt128 factor = oddPowers[tableIndex]; + result.MulMod(factor, modulusGpu); + + index = windowStart - 1; + } + + return (UInt128)result; + } + + private static bool ShouldUseSingleBit(in GpuUInt128 exponent) + { + if (exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold) + { + return true; + } + + return false; + } + + private static int GetWindowSize(int bitLength) + { + if (bitLength <= 8) + { + return Math.Max(bitLength, 1); + } + + if (bitLength <= 23) + { + return 4; + } + + if (bitLength <= 79) + { + return 5; + } + + if (bitLength <= 239) + { + return 6; + } + + if (bitLength <= 671) + { + return 7; + } + + return 8; + } + + private static void InitializeOddPowers(GpuUInt128 baseValue, GpuUInt128 modulus, Span oddPowers) + { + oddPowers[0] = baseValue; + if (oddPowers.Length == 1) + { + return; + } + + // Reusing baseValue to hold base^2 for the shared odd-power ladder. + baseValue.MulMod(baseValue, modulus); + for (int i = 1; i < oddPowers.Length; i++) + { + oddPowers[i] = oddPowers[i - 1]; + oddPowers[i].MulMod(baseValue, modulus); + } + } + + private static GpuUInt128 Pow2MontgomeryModSingleBit(GpuUInt128 exponent, GpuUInt128 modulus, GpuUInt128 baseValue) + { + GpuUInt128 result = GpuUInt128.One; + + while (!exponent.IsZero) + { + if ((exponent.Low & 1UL) != 0UL) + { + result.MulMod(baseValue, modulus); + } + + exponent.ShiftRight(1); + if (exponent.IsZero) + { + break; + } + + // Reusing baseValue to store the squared base for the next iteration. + baseValue.MulMod(baseValue, modulus); + } + + return result; + } + + private static bool IsBitSet(in GpuUInt128 value, int bitIndex) + { + if (bitIndex >= 64) + { + return ((value.High >> (bitIndex - 64)) & 1UL) != 0UL; + } + + return ((value.Low >> bitIndex) & 1UL) != 0UL; + } + + private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount) + { + if (windowStart != 0) + { + GpuUInt128 shifted = exponent; + shifted.ShiftRight(windowStart); + ulong mask = (1UL << windowBitCount) - 1UL; + return shifted.Low & mask; + } + + ulong directMask = (1UL << windowBitCount) - 1UL; + return exponent.Low & directMask; + } +} diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index c8da4b7..22c7781 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -1319,9 +1319,9 @@ internal static class Program int zeroCountValue = 0; int maxZeroBlockValue = 0; ComputeBitStats(p, out bitLength, out zeroCountValue, out maxZeroBlockValue); + zf = (double)zeroCountValue / bitLength; if (_zeroFracHard >= 0) { - zf = (double)zeroCountValue / bitLength; if (zf > _zeroFracHard) { return false; @@ -1330,7 +1330,6 @@ internal static class Program if (_zeroFracConj >= 0 && _maxZeroConj >= 0) { - zf = (double)zeroCountValue / bitLength; if (zf > _zeroFracConj && maxZeroBlockValue >= _maxZeroConj) { return false; diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuRationalTests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuRationalTests.cs new file mode 100644 index 0000000..9c88fb1 --- /dev/null +++ b/PerfectNumbers.Core.Tests/Gpu/GpuRationalTests.cs @@ -0,0 +1,82 @@ +using FluentAssertions; +using PerfectNumbers.Core.Gpu; +using Xunit; + +namespace PerfectNumbers.Core.Tests.Gpu; + +public class GpuRationalTests +{ + [Fact, Trait("Category", "Fast")] + public void Constructor_ReducesFraction() + { + GpuRational value = new(6UL, 8UL); + value.Numerator.Should().Be(new GpuUInt128(3UL)); + value.Denominator.Should().Be(new GpuUInt128(4UL)); + } + + [Fact, Trait("Category", "Fast")] + public void Addition_ReturnsExpectedValue() + { + GpuRational left = new(1UL, 3UL); + GpuRational right = new(1UL, 6UL); + + GpuRational sum = left + right; + + sum.Numerator.Should().Be(new GpuUInt128(1UL)); + sum.Denominator.Should().Be(new GpuUInt128(2UL)); + } + + [Fact, Trait("Category", "Fast")] + public void Subtraction_ReturnsExpectedValue() + { + GpuRational left = new(5UL, 6UL); + GpuRational right = new(1UL, 6UL); + + GpuRational difference = left - right; + + difference.Numerator.Should().Be(new GpuUInt128(2UL)); + difference.Denominator.Should().Be(new GpuUInt128(3UL)); + } + + [Fact, Trait("Category", "Fast")] + public void Multiplication_ReturnsExpectedValue() + { + GpuRational left = new(3UL, 4UL); + GpuRational right = new(2UL, 3UL); + + GpuRational product = left * right; + + product.Numerator.Should().Be(new GpuUInt128(1UL)); + product.Denominator.Should().Be(new GpuUInt128(2UL)); + } + + [Fact, Trait("Category", "Fast")] + public void Division_ReturnsExpectedValue() + { + GpuRational left = new(3UL, 5UL); + GpuRational right = new(9UL, 10UL); + + GpuRational quotient = left / right; + + quotient.Numerator.Should().Be(new GpuUInt128(2UL)); + quotient.Denominator.Should().Be(new GpuUInt128(3UL)); + } + + [Fact, Trait("Category", "Fast")] + public void CompareTo_DetectsOrdering() + { + GpuRational smaller = new(1UL, 4UL); + GpuRational larger = new(2UL, 3UL); + + smaller.CompareTo(larger).Should().BeLessThan(0); + larger.CompareTo(smaller).Should().BeGreaterThan(0); + smaller.CompareTo(new GpuRational(1UL, 4UL)).Should().Be(0); + } + + [Fact, Trait("Category", "Fast")] + public void ToDouble_ReturnsExpectedResult() + { + GpuRational value = new(1UL, 5UL); + value.ToDouble().Should().BeApproximately(0.2, 1e-12); + } +} diff --git a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs index 1bb9335..aca4483 100644 --- a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs @@ -1,4 +1,6 @@ using System; +using System.Numerics; +using System.Reflection; using FluentAssertions; using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; @@ -12,6 +14,7 @@ public class PrimeOrderGpuHeuristicsTests public PrimeOrderGpuHeuristicsTests() { PrimeOrderGpuHeuristics.OverflowRegistry.Clear(); + PrimeOrderGpuHeuristics.OverflowRegistryWide.Clear(); PrimeOrderGpuHeuristics.ResetCapabilitiesForTesting(); } @@ -165,6 +168,160 @@ public class PrimeOrderGpuHeuristicsTests } } + [Fact] + public void TryPow2ModWide_returns_overflow_for_marked_prime() + { + UInt128 prime = ((UInt128)1 << 80) + 7; + try + { + PrimeOrderGpuHeuristics.OverflowRegistryWide[prime] = 0; + + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(UInt128.One, prime, out UInt128 remainder); + + status.Should().Be(GpuPow2ModStatus.Overflow); + remainder.Should().Be(UInt128.Zero); + } + finally + { + PrimeOrderGpuHeuristics.OverflowRegistryWide.TryRemove(prime, out _); + } + } + + [Fact] + public void TryPow2ModWide_returns_success_for_supported_prime() + { + UInt128 prime = ((UInt128)1 << 96) + 11; + UInt128 exponent = ((UInt128)1 << 64) + 3; + + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out UInt128 remainder); + + status.Should().Be(GpuPow2ModStatus.Success); + UInt128 expected = (UInt128)BigInteger.ModPow(2, (BigInteger)exponent, (BigInteger)prime); + remainder.Should().Be(expected); + } + + [Theory] + [InlineData("193", "31")] + [InlineData("193", "255")] + [InlineData("79228162514264337593543950347", "36893488147419103251")] + public void Pow2ModKernelCore_matches_cpu_windowed_path(string modulusText, string exponentText) + { + UInt128 modulus = UInt128.Parse(modulusText); + UInt128 exponent = UInt128.Parse(exponentText); + + var method = typeof(PrimeOrderGpuHeuristics).GetMethod("Pow2ModKernelCore", BindingFlags.NonPublic | BindingFlags.Static); + method.Should().NotBeNull(); + + object? gpuResultObject = method!.Invoke(null, new object[] { new GpuUInt128(exponent), new GpuUInt128(modulus) }); + GpuUInt128 gpuResult = (GpuUInt128)gpuResultObject!; + UInt128 expected = exponent.Pow2MontgomeryModWindowed(modulus); + + ((UInt128)gpuResult).Should().Be(expected); + } + + [Fact] + public void TryPow2ModWide_returns_overflow_when_modulus_exceeds_capability() + { + UInt128 prime = ((UInt128)1 << 80) + 9; + var capability = new PrimeOrderGpuHeuristics.PrimeOrderGpuCapability(64, 128); + PrimeOrderGpuHeuristics.OverrideCapabilitiesForTesting(capability); + try + { + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(UInt128.One, prime, out UInt128 remainder); + + status.Should().Be(GpuPow2ModStatus.Overflow); + remainder.Should().Be(UInt128.Zero); + } + finally + { + PrimeOrderGpuHeuristics.ResetCapabilitiesForTesting(); + PrimeOrderGpuHeuristics.OverflowRegistryWide.TryRemove(prime, out _); + } + } + + [Fact] + public void TryPow2ModWide_returns_overflow_when_exponent_exceeds_capability() + { + var capability = new PrimeOrderGpuHeuristics.PrimeOrderGpuCapability(128, 4); + PrimeOrderGpuHeuristics.OverrideCapabilitiesForTesting(capability); + try + { + UInt128 exponent = UInt128.One << 8; + UInt128 prime = ((UInt128)1 << 70) + 3; + + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out UInt128 remainder); + + status.Should().Be(GpuPow2ModStatus.Overflow); + remainder.Should().Be(UInt128.Zero); + } + finally + { + PrimeOrderGpuHeuristics.ResetCapabilitiesForTesting(); + } + } + + [Fact] + public void TryPow2ModBatchWide_returns_success_for_supported_prime() + { + UInt128 prime = ((UInt128)1 << 90) + 19; + UInt128[] exponents = { UInt128.One, ((UInt128)1 << 64) + 5 }; + UInt128[] remainders = { UInt128.Zero, UInt128.Zero }; + + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders); + + status.Should().Be(GpuPow2ModStatus.Success); + for (int i = 0; i < exponents.Length; i++) + { + UInt128 expected = (UInt128)BigInteger.ModPow(2, (BigInteger)exponents[i], (BigInteger)prime); + remainders[i].Should().Be(expected); + } + } + + [Fact] + public void TryPow2ModBatchWide_returns_overflow_when_modulus_exceeds_capability() + { + UInt128 prime = ((UInt128)1 << 88) + 13; + UInt128[] exponents = { UInt128.One, (UInt128)5 }; + UInt128[] remainders = { UInt128.Zero, UInt128.Zero }; + var capability = new PrimeOrderGpuHeuristics.PrimeOrderGpuCapability(64, 128); + PrimeOrderGpuHeuristics.OverrideCapabilitiesForTesting(capability); + try + { + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders); + + status.Should().Be(GpuPow2ModStatus.Overflow); + remainders[0].Should().Be(UInt128.Zero); + remainders[1].Should().Be(UInt128.Zero); + } + finally + { + PrimeOrderGpuHeuristics.ResetCapabilitiesForTesting(); + PrimeOrderGpuHeuristics.OverflowRegistryWide.TryRemove(prime, out _); + } + } + + [Fact] + public void TryPow2ModBatchWide_returns_overflow_when_any_exponent_exceeds_capability() + { + UInt128 prime = ((UInt128)1 << 72) + 23; + UInt128[] exponents = { UInt128.One << 5, UInt128.One << 9 }; + UInt128[] remainders = { UInt128.Zero, UInt128.Zero }; + var capability = new PrimeOrderGpuHeuristics.PrimeOrderGpuCapability(128, 6); + PrimeOrderGpuHeuristics.OverrideCapabilitiesForTesting(capability); + try + { + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2ModBatch(exponents, prime, remainders); + + status.Should().Be(GpuPow2ModStatus.Overflow); + remainders[0].Should().Be(UInt128.Zero); + remainders[1].Should().Be(UInt128.Zero); + } + finally + { + PrimeOrderGpuHeuristics.ResetCapabilitiesForTesting(); + } + } + [Fact] public void TryPow2ModBatch_throws_when_remainder_span_is_too_small() { diff --git a/PerfectNumbers.Core/Gpu/GpuRational.cs b/PerfectNumbers.Core/Gpu/GpuRational.cs new file mode 100644 index 0000000..4e9c1c9 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/GpuRational.cs @@ -0,0 +1,269 @@ +using System; +using System.Numerics; +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core.Gpu; + +public struct GpuRational : IEquatable, IComparable +{ + private const double Pow64 = 18446744073709551616.0; + + public GpuUInt128 Numerator; + + public GpuUInt128 Denominator; + + public static readonly GpuRational Zero = new(GpuUInt128.Zero, GpuUInt128.One); + + public static readonly GpuRational One = new(GpuUInt128.One, GpuUInt128.One); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public GpuRational(GpuUInt128 numerator) + : this(numerator, GpuUInt128.One) + { + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public GpuRational(ulong numerator) + : this(new GpuUInt128(numerator), GpuUInt128.One) + { + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public GpuRational(GpuUInt128 numerator, GpuUInt128 denominator) + { + if (denominator.IsZero) + { + throw new DivideByZeroException("GpuRational denominator cannot be zero."); + } + + Numerator = numerator; + Denominator = denominator; + Normalize(); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public GpuRational(ulong numerator, ulong denominator) + : this(new GpuUInt128(numerator), new GpuUInt128(denominator)) + { + } + + public readonly bool IsZero + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get => Numerator.IsZero; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static implicit operator GpuRational(GpuUInt128 value) => new(value); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static implicit operator GpuRational(ulong value) => new(value); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Add(in GpuRational other) + { + GpuUInt128 leftDenominator = Denominator; + GpuUInt128 rightDenominator = other.Denominator; + + GpuUInt128 gcd = GpuUInt128.BinaryGcd(leftDenominator, rightDenominator); + if (gcd != GpuUInt128.One) + { + leftDenominator = GpuUInt128.DivideExact(leftDenominator, gcd); + rightDenominator = GpuUInt128.DivideExact(rightDenominator, gcd); + // Reusing both denominator variables to store their reduced values because the unreduced inputs are no longer needed. + } + + Numerator.Mul(rightDenominator); + GpuUInt128 scaledOtherNumerator = other.Numerator; + scaledOtherNumerator.Mul(leftDenominator); + Numerator.Add(scaledOtherNumerator); + + Denominator = leftDenominator; + // Denominator now holds the reduced left denominator; multiply by the original right denominator to form the LCM. + Denominator.Mul(other.Denominator); + Normalize(); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Sub(in GpuRational other) + { + GpuUInt128 leftDenominator = Denominator; + GpuUInt128 rightDenominator = other.Denominator; + + GpuUInt128 gcd = GpuUInt128.BinaryGcd(leftDenominator, rightDenominator); + if (gcd != GpuUInt128.One) + { + leftDenominator = GpuUInt128.DivideExact(leftDenominator, gcd); + rightDenominator = GpuUInt128.DivideExact(rightDenominator, gcd); + // Reusing the denominator variables to keep the reduced values since the raw denominators are no longer required. + } + + Numerator.Mul(rightDenominator); + GpuUInt128 scaledOtherNumerator = other.Numerator; + scaledOtherNumerator.Mul(leftDenominator); + if (Numerator.CompareTo(scaledOtherNumerator) < 0) + { + throw new InvalidOperationException("Resulting GpuRational would be negative."); + } + + Numerator.Sub(scaledOtherNumerator); + + Denominator = leftDenominator; + Denominator.Mul(other.Denominator); + Normalize(); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Mul(in GpuRational other) + { + GpuUInt128 rightDenominator = other.Denominator; + GpuUInt128 gcd = GpuUInt128.BinaryGcd(Numerator, rightDenominator); + if (gcd != GpuUInt128.One) + { + Numerator = GpuUInt128.DivideExact(Numerator, gcd); + rightDenominator = GpuUInt128.DivideExact(rightDenominator, gcd); + // Reusing rightDenominator to keep the reduced value because the original denominator is no longer required. + } + + GpuUInt128 rightNumerator = other.Numerator; + GpuUInt128 leftDenominator = Denominator; + gcd = GpuUInt128.BinaryGcd(rightNumerator, leftDenominator); + if (gcd != GpuUInt128.One) + { + rightNumerator = GpuUInt128.DivideExact(rightNumerator, gcd); + leftDenominator = GpuUInt128.DivideExact(leftDenominator, gcd); + // Reusing rightNumerator and leftDenominator after cross-cancellation keeps the intermediates compact. + } + + Numerator.Mul(rightNumerator); + Denominator = leftDenominator; + Denominator.Mul(rightDenominator); + Normalize(); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Div(in GpuRational other) + { + if (other.Numerator.IsZero) + { + throw new DivideByZeroException("Cannot divide by a zero GpuRational."); + } + + GpuUInt128 rightNumerator = other.Numerator; + GpuUInt128 gcd = GpuUInt128.BinaryGcd(Numerator, rightNumerator); + if (gcd != GpuUInt128.One) + { + Numerator = GpuUInt128.DivideExact(Numerator, gcd); + rightNumerator = GpuUInt128.DivideExact(rightNumerator, gcd); + // Reusing rightNumerator to store the reduced numerator keeps temporary allocations out of the hot path. + } + + GpuUInt128 rightDenominator = other.Denominator; + GpuUInt128 leftDenominator = Denominator; + gcd = GpuUInt128.BinaryGcd(rightDenominator, leftDenominator); + if (gcd != GpuUInt128.One) + { + rightDenominator = GpuUInt128.DivideExact(rightDenominator, gcd); + leftDenominator = GpuUInt128.DivideExact(leftDenominator, gcd); + // Reusing both denominator variables here mirrors the cross-cancellation performed on the numerators. + } + + Numerator.Mul(rightDenominator); + Denominator = leftDenominator; + Denominator.Mul(rightNumerator); + Normalize(); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static GpuRational operator +(GpuRational left, GpuRational right) + { + left.Add(right); + return left; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static GpuRational operator -(GpuRational left, GpuRational right) + { + left.Sub(right); + return left; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static GpuRational operator *(GpuRational left, GpuRational right) + { + left.Mul(right); + return left; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static GpuRational operator /(GpuRational left, GpuRational right) + { + left.Div(right); + return left; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public readonly int CompareTo(GpuRational other) + { + GpuUInt128 leftScaled = Numerator; + leftScaled.Mul(other.Denominator); + GpuUInt128 rightScaled = other.Numerator; + rightScaled.Mul(Denominator); + return leftScaled.CompareTo(rightScaled); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public readonly bool Equals(GpuRational other) => Numerator == other.Numerator && Denominator == other.Denominator; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public override readonly bool Equals(object? obj) => obj is GpuRational other && Equals(other); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public override readonly int GetHashCode() => HashCode.Combine(Numerator, Denominator); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool operator ==(GpuRational left, GpuRational right) => left.Equals(right); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool operator !=(GpuRational left, GpuRational right) => !left.Equals(right); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool operator <(GpuRational left, GpuRational right) => left.CompareTo(right) < 0; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool operator >(GpuRational left, GpuRational right) => left.CompareTo(right) > 0; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool operator <=(GpuRational left, GpuRational right) => left.CompareTo(right) <= 0; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool operator >=(GpuRational left, GpuRational right) => left.CompareTo(right) >= 0; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public readonly double ToDouble() + { + double numerator = Numerator.High == 0UL ? Numerator.Low : Numerator.High * Pow64 + Numerator.Low; + double denominator = Denominator.High == 0UL ? Denominator.Low : Denominator.High * Pow64 + Denominator.Low; + return numerator / denominator; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public override readonly string ToString() => $"{Numerator}/{Denominator}"; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private void Normalize() + { + if (Numerator.IsZero) + { + Denominator = GpuUInt128.One; + return; + } + + GpuUInt128 gcd = GpuUInt128.BinaryGcd(Numerator, Denominator); + if (gcd != GpuUInt128.One) + { + Numerator = GpuUInt128.DivideExact(Numerator, gcd); + Denominator = GpuUInt128.DivideExact(Denominator, gcd); + } + } +} diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs index e213e93..2f5bc36 100644 --- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs +++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs @@ -54,6 +54,26 @@ public struct GpuUInt128 : IComparable, IEquatable get => High == 0UL && Low == 0UL; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public readonly int GetBitLength() + { + if (High != 0UL) + { + return 64 + (64 - BitOperations.LeadingZeroCount(High)); + } + + ulong low = Low; + if (low == 0UL) + { + return 0; + } + + return 64 - BitOperations.LeadingZeroCount(low); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static int GetBitLength(GpuUInt128 value) => value.GetBitLength(); + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static implicit operator GpuUInt128(UInt128 value) => new(value); @@ -564,6 +584,43 @@ public struct GpuUInt128 : IComparable, IEquatable Low = lowPart; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static GpuUInt128 DivideExact(GpuUInt128 dividend, GpuUInt128 divisor) + { + int divisorBits = divisor.GetBitLength(); + int shift = dividend.GetBitLength() - divisorBits; + + GpuUInt128 quotient = Zero; + GpuUInt128 remainder = dividend; + GpuUInt128 shiftedDivisor = divisor << shift; + GpuUInt128 quotientBit = One << shift; + + while (true) + { + if (remainder.CompareTo(shiftedDivisor) >= 0) + { + remainder.Sub(shiftedDivisor); + quotient.Add(quotientBit); + } + + if (shift == 0) + { + break; + } + + shift--; + shiftedDivisor.ShiftRight(1); + quotientBit.ShiftRight(1); + } + + if (!remainder.IsZero) + { + throw new InvalidOperationException("DivideExact requires divisible operands."); + } + + return quotient; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public void Xor(GpuUInt128 other) { diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index ea94b27..e3da840 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -1,9 +1,12 @@ using System; +using System.Buffers; using System.Collections.Concurrent; using ILGPU; using ILGPU.Runtime; +using System.Numerics; using PerfectNumbers.Core; using System.Runtime.InteropServices; +using System.Runtime.CompilerServices; namespace PerfectNumbers.Core.Gpu; internal enum GpuPow2ModStatus @@ -16,10 +19,24 @@ internal enum GpuPow2ModStatus internal static class PrimeOrderGpuHeuristics { private static readonly ConcurrentDictionary OverflowedPrimes = new(); + private static readonly ConcurrentDictionary OverflowedPrimesWide = new(); private static readonly ConcurrentDictionary, MontgomeryDivisorData, ArrayView1D>> Pow2ModKernelCache = new(); + private static readonly ConcurrentDictionary, GpuUInt128, ArrayView1D>> Pow2ModKernelWideCache = new(); + private const int WideStackThreshold = 8; private static PrimeOrderGpuCapability s_capability = PrimeOrderGpuCapability.Default; + private const int Pow2WindowSizeBits = 8; + private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1); + private const ulong Pow2WindowFallbackThreshold = 32UL; + + [InlineArray(Pow2WindowOddPowerCount)] + private struct Pow2OddPowerTable + { + private GpuUInt128 _element0; + } + internal static ConcurrentDictionary OverflowRegistry => OverflowedPrimes; + internal static ConcurrentDictionary OverflowRegistryWide => OverflowedPrimesWide; internal static void OverrideCapabilitiesForTesting(PrimeOrderGpuCapability capability) { @@ -89,6 +106,59 @@ internal static class PrimeOrderGpuHeuristics return computed ? GpuPow2ModStatus.Success : GpuPow2ModStatus.Unavailable; } + public static GpuPow2ModStatus TryPow2Mod(UInt128 exponent, UInt128 prime, out UInt128 remainder) + { + Span exponents = stackalloc UInt128[1]; + Span remainders = stackalloc UInt128[1]; + exponents[0] = exponent; + + GpuPow2ModStatus status = TryPow2ModBatch(exponents, prime, remainders); + remainder = remainders[0]; + return status; + } + + public static GpuPow2ModStatus TryPow2ModBatch(ReadOnlySpan exponents, UInt128 prime, Span remainders) + { + 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)); + } + + Span target = remainders.Slice(0, exponents.Length); + target.Clear(); + + ConcurrentDictionary overflowRegistryWide = OverflowedPrimesWide; + + if (overflowRegistryWide.ContainsKey(prime)) + { + return GpuPow2ModStatus.Overflow; + } + + PrimeOrderGpuCapability capability = s_capability; + + if (prime.GetBitLength() > capability.ModulusBits) + { + overflowRegistryWide[prime] = 0; + return GpuPow2ModStatus.Overflow; + } + + for (int i = 0; i < exponents.Length; i++) + { + if (exponents[i].GetBitLength() > capability.ExponentBits) + { + return GpuPow2ModStatus.Overflow; + } + } + + bool computed = TryComputeOnGpuWide(exponents, prime, target); + return computed ? GpuPow2ModStatus.Success : GpuPow2ModStatus.Unavailable; + } + private static bool TryComputeOnGpu(ReadOnlySpan exponents, ulong prime, Span results) { try @@ -128,12 +198,94 @@ internal static class PrimeOrderGpuHeuristics private static void ComputePow2ModCpu(ReadOnlySpan exponents, ulong prime, Span results) { - for (int i = 0; i < exponents.Length; i++) + int length = exponents.Length; + for (int i = 0; i < length; i++) { results[i] = Pow2ModCpu(exponents[i], prime); } } + private static bool TryComputeOnGpuWide(ReadOnlySpan exponents, UInt128 prime, Span results) + { + int length = exponents.Length; + if (length == 0) + { + return true; + } + + GpuUInt128[]? rentedExponents = null; + GpuUInt128[]? rentedResults = null; + var lease = GpuKernelPool.GetKernel(useGpuOrder: true); + + try + { + using 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); + + Span exponentSpan = length <= WideStackThreshold + ? stackalloc GpuUInt128[length] + : new Span(rentedExponents = ArrayPool.Shared.Rent(length), 0, length); + + for (int i = 0; i < length; i++) + { + exponentSpan[i] = (GpuUInt128)exponents[i]; + } + + exponentBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSpan), length); + remainderBuffer.MemSetToZero(); + + GpuUInt128 modulus = (GpuUInt128)prime; + kernel(stream, length, exponentBuffer.View, modulus, remainderBuffer.View); + + stream.Synchronize(); + + Span resultSpan = length <= WideStackThreshold + ? stackalloc GpuUInt128[length] + : new Span(rentedResults = ArrayPool.Shared.Rent(length), 0, length); + + remainderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(resultSpan), length); + + for (int i = 0; i < length; i++) + { + results[i] = (UInt128)resultSpan[i]; + } + + return true; + } + catch (Exception ex) when (ex is AcceleratorException or InternalCompilerException or NotSupportedException or InvalidOperationException or AggregateException) + { + ComputePow2ModCpuWide(exponents, prime, results); + return true; + } + finally + { + if (rentedExponents is not null) + { + ArrayPool.Shared.Return(rentedExponents, clearArray: false); + } + + if (rentedResults is not null) + { + ArrayPool.Shared.Return(rentedResults, clearArray: false); + } + + lease.Dispose(); + } + } + + private static void ComputePow2ModCpuWide(ReadOnlySpan exponents, UInt128 prime, Span results) + { + int length = exponents.Length; + for (int i = 0; i < length; i++) + { + results[i] = exponents[i].Pow2MontgomeryModWindowed(prime); + } + } + private static ulong Pow2ModCpu(ulong exponent, ulong modulus) { if (modulus <= 1UL) @@ -161,8 +313,191 @@ internal static class PrimeOrderGpuHeuristics remainders[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false); } + private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) + { + return Pow2ModKernelWideCache.GetOrAdd(accelerator, static accel => + { + var loaded = accel.LoadAutoGroupedStreamKernel, GpuUInt128, ArrayView1D>(Pow2ModKernelWide); + var kernel = KernelUtil.GetKernel(loaded); + return kernel.CreateLauncherDelegate, GpuUInt128, ArrayView1D>>(); + }); + } + + private static void Pow2ModKernelWide(Index1D index, ArrayView1D exponents, GpuUInt128 modulus, ArrayView1D remainders) + { + GpuUInt128 exponent = exponents[index]; + remainders[index] = Pow2ModKernelCore(exponent, modulus); + } + + private static GpuUInt128 Pow2ModKernelCore(GpuUInt128 exponent, GpuUInt128 modulus) + { + if (modulus == GpuUInt128.One) + { + return GpuUInt128.Zero; + } + + if (exponent.IsZero) + { + return GpuUInt128.One; + } + + GpuUInt128 baseValue = new GpuUInt128(2UL); + 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 oddPowerCount = 1 << (windowSize - 1); + + Pow2OddPowerTable oddPowers = default; + InitializeOddPowers(ref oddPowers, baseValue, modulus, oddPowerCount); + + GpuUInt128 result = GpuUInt128.One; + int index = bitLength - 1; + + while (index >= 0) + { + if (!IsBitSet(exponent, index)) + { + result.MulMod(result, modulus); + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (!IsBitSet(exponent, windowStart)) + { + windowStart++; + } + + int windowBitCount = index - windowStart + 1; + for (int square = 0; square < windowBitCount; square++) + { + 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; + } + + return result; + } + + private static bool ShouldUseSingleBit(in GpuUInt128 exponent) => exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold; + + private static int GetWindowSize(int bitLength) + { + if (bitLength <= Pow2WindowSizeBits) + { + return Math.Max(bitLength, 1); + } + + if (bitLength <= 23) + { + return 4; + } + + if (bitLength <= 79) + { + return 5; + } + + if (bitLength <= 239) + { + return 6; + } + + if (bitLength <= 671) + { + return 7; + } + + 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++) + { + oddPowers[i] = oddPowers[i - 1]; + oddPowers[i].MulMod(baseValue, modulus); + } + } + + private static GpuUInt128 Pow2MontgomeryModSingleBit(GpuUInt128 exponent, GpuUInt128 modulus, GpuUInt128 baseValue) + { + GpuUInt128 result = GpuUInt128.One; + + while (!exponent.IsZero) + { + if ((exponent.Low & 1UL) != 0UL) + { + result.MulMod(baseValue, modulus); + } + + exponent.ShiftRight(1); + if (exponent.IsZero) + { + break; + } + + // Reusing baseValue to store the squared base for the next iteration. + baseValue.MulMod(baseValue, modulus); + } + + return result; + } + + private static bool IsBitSet(in GpuUInt128 value, int bitIndex) + { + if (bitIndex >= 64) + { + return ((value.High >> (bitIndex - 64)) & 1UL) != 0UL; + } + + return ((value.Low >> bitIndex) & 1UL) != 0UL; + } + + private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount) + { + if (windowStart != 0) + { + GpuUInt128 shifted = exponent; + shifted.ShiftRight(windowStart); + ulong mask = (1UL << windowBitCount) - 1UL; + return shifted.Low & mask; + } + + ulong directMask = (1UL << windowBitCount) - 1UL; + return exponent.Low & directMask; + } + internal readonly record struct PrimeOrderGpuCapability(int ModulusBits, int ExponentBits) { - public static PrimeOrderGpuCapability Default => new(128, 64); + public static PrimeOrderGpuCapability Default => new(128, 128); } } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs index 8951b4d..9e85e7b 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs @@ -959,15 +959,20 @@ internal static partial class PrimeOrderCalculator { ulong prime64 = (ulong)modulus; ulong exponent64 = (ulong)exponent; - // TODO: Extend PrimeOrderGpuHeuristics to handle 128-bit moduli so this fast path can remain on the device. GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent64, prime64, out ulong remainder); if (status == GpuPow2ModStatus.Success) { return remainder; } } - BigInteger result = BigInteger.ModPow(2, (BigInteger)exponent, (BigInteger)modulus); - return (UInt128)result; + + GpuPow2ModStatus wideStatus = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, modulus, out UInt128 wideRemainder); + if (wideStatus == GpuPow2ModStatus.Success) + { + return wideRemainder; + } + + return exponent.Pow2MontgomeryModWindowed(modulus); } private static void AddFactor(Dictionary counts, UInt128 prime, int exponent) diff --git a/PerfectNumbers.Core/UInt128Extensions.cs b/PerfectNumbers.Core/UInt128Extensions.cs index 06653dd..870a87e 100644 --- a/PerfectNumbers.Core/UInt128Extensions.cs +++ b/PerfectNumbers.Core/UInt128Extensions.cs @@ -1,10 +1,15 @@ +using System; using System.Numerics; using System.Runtime.CompilerServices; +using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; public static class UInt128Extensions { + private const int Pow2WindowSize = 8; + private const ulong Pow2WindowFallbackThreshold = 32UL; + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static UInt128 BinaryGcd(this UInt128 u, UInt128 v) { @@ -40,6 +45,315 @@ public static class UInt128Extensions return u << shift; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static int GetBitLength(this UInt128 value) + { + ulong high = (ulong)(value >> 64); + if (high != 0UL) + { + return 64 + (64 - BitOperations.LeadingZeroCount(high)); + } + + ulong low = (ulong)value; + if (low == 0UL) + { + return 0; + } + + return 64 - BitOperations.LeadingZeroCount(low); + } + + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static UInt128 Pow2MontgomeryModWindowed(this UInt128 exponent, UInt128 modulus) + { + if (modulus == UInt128.One) + { + return UInt128.Zero; + } + + if (exponent == UInt128.Zero) + { + return UInt128.One % modulus; + } + + GpuUInt128 modulusGpu = new(modulus); + GpuUInt128 baseValue = new(2UL); + if (baseValue.CompareTo(modulusGpu) >= 0) + { + baseValue.Sub(modulusGpu); + } + + if (ShouldUseSingleBit(exponent)) + { + GpuUInt128 singleBitResult = Pow2MontgomeryModSingleBit(exponent, modulusGpu, baseValue); + return (UInt128)singleBitResult; + } + + GpuUInt128 exponentGpu = new(exponent); + int bitLength = exponentGpu.GetBitLength(); + int windowSize = GetWindowSize(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + Span oddPowers = stackalloc GpuUInt128[oddPowerCount]; + InitializeOddPowers(baseValue, modulusGpu, oddPowers); + + GpuUInt128 result = GpuUInt128.One; + int index = bitLength - 1; + + while (index >= 0) + { + if (!IsBitSet(exponentGpu, index)) + { + result.MulMod(result, modulusGpu); + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (!IsBitSet(exponentGpu, windowStart)) + { + windowStart++; + } + + int windowBitCount = index - windowStart + 1; + for (int square = 0; square < windowBitCount; square++) + { + result.MulMod(result, modulusGpu); + } + + ulong windowValue = ExtractWindowValue(exponentGpu, windowStart, windowBitCount); + int tableIndex = (int)((windowValue - 1UL) >> 1); + GpuUInt128 factor = oddPowers[tableIndex]; + result.MulMod(factor, modulusGpu); + + index = windowStart - 1; + } + + return (UInt128)result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool ShouldUseSingleBit(UInt128 exponent) + { + return (exponent >> 64) == UInt128.Zero && (ulong)exponent <= Pow2WindowFallbackThreshold; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetWindowSize(int bitLength) + { + if (bitLength <= Pow2WindowSize) + { + return Math.Max(bitLength, 1); + } + + if (bitLength <= 23) + { + return 4; + } + + if (bitLength <= 79) + { + return 5; + } + + if (bitLength <= 239) + { + return 6; + } + + if (bitLength <= 671) + { + return 7; + } + + return Pow2WindowSize; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void InitializeOddPowers(GpuUInt128 baseValue, GpuUInt128 modulus, Span oddPowers) + { + oddPowers[0] = baseValue; + if (oddPowers.Length == 1) + { + return; + } + + // Reusing baseValue to hold base^2 for the shared odd-power ladder, just like the GPU helper. + baseValue.MulMod(baseValue, modulus); + for (int i = 1; i < oddPowers.Length; i++) + { + oddPowers[i] = oddPowers[i - 1]; + oddPowers[i].MulMod(baseValue, modulus); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static GpuUInt128 Pow2MontgomeryModSingleBit(UInt128 exponent, GpuUInt128 modulus, GpuUInt128 baseValue) + { + GpuUInt128 result = GpuUInt128.One; + GpuUInt128 remainingExponent = new(exponent); + + while (!remainingExponent.IsZero) + { + if ((remainingExponent.Low & 1UL) != 0UL) + { + result.MulMod(baseValue, modulus); + } + + remainingExponent.ShiftRight(1); + if (remainingExponent.IsZero) + { + break; + } + + // Reusing baseValue to hold the squared base in-place between iterations. + baseValue.MulMod(baseValue, modulus); + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool IsBitSet(in GpuUInt128 value, int bitIndex) + { + if (bitIndex >= 64) + { + return ((value.High >> (bitIndex - 64)) & 1UL) != 0UL; + } + + return ((value.Low >> bitIndex) & 1UL) != 0UL; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount) + { + if (windowStart != 0) + { + GpuUInt128 shifted = exponent; + shifted.ShiftRight(windowStart); + ulong mask = (1UL << windowBitCount) - 1UL; + return shifted.Low & mask; + } + + ulong directMask = (1UL << windowBitCount) - 1UL; + return exponent.Low & directMask; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static UInt128 ReduceProductBitwise(ulong p3, ulong p2, ulong p1, ulong p0, UInt128 modulus) + { + UInt128 remainder = UInt128.Zero; + for (int bit = 63; bit >= 0; bit--) + { + remainder <<= 1; + remainder |= (UInt128)((p3 >> bit) & 1UL); + if (remainder >= modulus) + { + remainder -= modulus; + } + } + + for (int bit = 63; bit >= 0; bit--) + { + remainder <<= 1; + remainder |= (UInt128)((p2 >> bit) & 1UL); + if (remainder >= modulus) + { + remainder -= modulus; + } + } + + for (int bit = 63; bit >= 0; bit--) + { + remainder <<= 1; + remainder |= (UInt128)((p1 >> bit) & 1UL); + if (remainder >= modulus) + { + remainder -= modulus; + } + } + + for (int bit = 63; bit >= 0; bit--) + { + remainder <<= 1; + remainder |= (UInt128)((p0 >> bit) & 1UL); + if (remainder >= modulus) + { + remainder -= modulus; + } + } + + return remainder; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void MultiplyFull(UInt128 left, UInt128 right, out ulong p3, out ulong p2, out ulong p1, out ulong p0) + { + ulong leftLow = (ulong)left; + ulong leftHigh = (ulong)(left >> 64); + ulong rightLow = (ulong)right; + ulong rightHigh = (ulong)(right >> 64); + + var partial0 = Multiply64(leftLow, rightLow); + var partial1 = Multiply64(leftLow, rightHigh); + var partial2 = Multiply64(leftHigh, rightLow); + var partial3 = Multiply64(leftHigh, rightHigh); + + p0 = partial0.Low; + ulong carry = partial0.High; + ulong sum = carry + partial1.Low; + ulong carryMid = sum < partial1.Low ? 1UL : 0UL; + sum += partial2.Low; + if (sum < partial2.Low) + { + carryMid++; + } + + p1 = sum; + sum = partial1.High + partial2.High; + ulong carryHigh = sum < partial2.High ? 1UL : 0UL; + sum += partial3.Low; + if (sum < partial3.Low) + { + carryHigh++; + } + + sum += carryMid; + if (sum < carryMid) + { + carryHigh++; + } + + p2 = sum; + p3 = partial3.High + carryHigh; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static (ulong High, ulong Low) Multiply64(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 & 0xFFFFFFFFUL) | (a0 << 32); + b1 += (mid1 >> 32) + (b0 >> 32) + (a0 >> 32); + + return (b1, a1); + } + public static ulong CalculateOrder(this UInt128 q) { if (q <= UInt128Numbers.Two) @@ -294,31 +608,31 @@ public static class UInt128Extensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public static UInt128 MulMod(this UInt128 a, UInt128 b, UInt128 modulus) { - UInt128 one = UInt128.One, - result = UInt128.Zero, - zero = UInt128.Zero; - - while (b != zero) + if (modulus <= UInt128.One) { - if ((b & one) != zero) - { - result += a; - if (result >= modulus) - { - result -= modulus; - } - } - - a <<= 1; - if (a >= modulus) - { - a -= modulus; - } - - b >>= 1; + return UInt128.Zero; } - return result; + UInt128 left = a % modulus; + UInt128 right = b % modulus; + if (left == UInt128.Zero || right == UInt128.Zero) + { + return UInt128.Zero; + } + + if ((modulus >> 64) == UInt128.Zero) + { + // Reusing left to hold the reduced product; this mirrors the inline UInt128 multiply + // that topped MulMod64Benchmarks for 64-bit moduli. + ulong modulus64 = (ulong)modulus; + left = (UInt128)(((UInt128)(ulong)left * (ulong)right) % modulus64); + return left; + } + + MultiplyFull(left, right, out var p3, out var p2, out var p1, out var p0); + // Reusing left to capture the reduced 128-bit remainder before returning. + left = ReduceProductBitwise(p3, p2, p1, p0, modulus); + return left; } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/todo.txt b/todo.txt index 489dbff..9586c58 100644 --- a/todo.txt +++ b/todo.txt @@ -1,2 +1,3 @@ -- Implement GpuRational struct to extend GpuUInt128 range. Avoid extending method calls stack. Instead prefer code duplication for a better performance. Create unit tests for it and benchmarks to compare against GpuUInt128. GpuRational should work in a similar way to ERational to avoid excessive memory usage. You have many methods available in GpuUInt128 and GpuUInt128Extensions. Many are also similar to UInt128Extensions and/or ULongExtensions. +- Prototype a UInt128-native Pow2MontgomeryModWindowed path that matches the GPU implementation once a faster modular multiplication helper is available, so we can revisit removing the GPU struct dependency without regressing benchmarks. +- Mirror the complete CPU heuristic cycle algorithm on the GPU path so both implementations stay aligned. - Extend the GPU numeric type probe to cover proposed custom types such as GpuUInt256 and GpuUInt512 once their definitions are available so we can plan the wide-order heuristics on the device.