From c1dc2e14d5a133f5eed356ab5952f5178f7dab58 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 11 Sep 2025 02:36:16 +0200 Subject: [PATCH 1/2] Allocate unique p ranges per thread --- EvenPerfectBitScanner/Program.cs | 124 ++++++++++++++++--------------- 1 file changed, 65 insertions(+), 59 deletions(-) diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 4d482d1..526c71b 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -589,46 +589,40 @@ internal static class Program } Console.WriteLine("Starting scan..."); - _state = ((long)currentP << 3) | (long)remainder; - Task[] tasks = new Task[threadCount]; + _state = ((long)currentP << 3) | (long)remainder; + Task[] tasks = new Task[threadCount]; - for (int i = 0; i < threadCount; i++) - { - tasks[i] = Task.Run(() => - { - ulong[] buffer = new ulong[blockSize]; - while (!Volatile.Read(ref _limitReached)) - { - int count = 0; - while (count < blockSize && !Volatile.Read(ref _limitReached)) - { - buffer[count] = ReserveNextP(); - count++; - } + for (int i = 0; i < threadCount; i++) + { + tasks[i] = Task.Run(() => + { + ulong[] buffer = new ulong[blockSize]; + while (!Volatile.Read(ref _limitReached)) + { + int count = ReserveBlock(buffer, blockSize); + if (count == 0) + { + break; + } - if (count == 0) - { - break; - } + for (int j = 0; j < count && !Volatile.Read(ref _limitReached); j++) + { + ulong p = buffer[j]; + if (useFilter && !filter.Contains(p)) + { + continue; + } + else if (useFilter) + { + Console.WriteLine($"Testing {p}"); + } - for (int j = 0; j < count && !Volatile.Read(ref _limitReached); j++) - { - ulong p = buffer[j]; - if (useFilter && !filter.Contains(p)) - { - continue; - } - else if (useFilter) - { - Console.WriteLine($"Testing {p}"); - } - - bool isPerfect = IsEvenPerfectCandidate(p, out bool searchedMersenne, out bool detailedCheck); - PrintResult(p, searchedMersenne, detailedCheck, isPerfect); - } - } - }); - } + bool isPerfect = IsEvenPerfectCandidate(p, divisorCyclesSearchLimit, out bool searchedMersenne, out bool detailedCheck); + PrintResult(p, searchedMersenne, detailedCheck, isPerfect); + } + } + }); + } Task.WaitAll(tasks); FlushBuffer(); @@ -742,24 +736,36 @@ internal static class Program return $"even_perfect_bit_scan_inc-{inc}_thr-{threads}_blk-{block}_mers-{mers}_mersdev-{mersenneDevice}_ntt-{ntt}_red-{red}_primesdev-{primesDevice}_{order}_orderdev-{orderDevice}_{gcd}_{work}_gputh-{gpuPrimeThreads}_llslice-{llSlice}_scanb-{gpuScanBatch}_warm-{warmupLimit}.csv"; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static unsafe ulong ReserveNextP() - { - while (true) - { - long state = Volatile.Read(ref _state); - ulong current = (ulong)state >> 3; - ulong remainder = ((ulong)state) & 7UL; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static unsafe int ReserveBlock(ulong[] buffer, int blockSize) + { + while (true) + { + long state = Volatile.Read(ref _state); + ulong current = (ulong)state >> 3; + ulong remainder = ((ulong)state) & 7UL; - ulong next = _transformP(current, ref remainder); - long newState = ((long)next << 3) | (long)remainder; - long original = Interlocked.CompareExchange(ref _state, newState, state); - if (original == state) - { - return current; - } - } - } + ulong p = current; + int count = 0; + while (count < blockSize && !Volatile.Read(ref _limitReached)) + { + buffer[count++] = p; + p = _transformP(p, ref remainder); + } + + if (count == 0) + { + return 0; + } + + long newState = ((long)p << 3) | (long)remainder; + long original = Interlocked.CompareExchange(ref _state, newState, state); + if (original == state) + { + return count; + } + } + } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void PrintResult(ulong currentP, bool searchedMersenne, bool detailedCheck, bool isPerfect) @@ -948,11 +954,11 @@ internal static class Program return remainder == 1UL; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static bool IsEvenPerfectCandidate(ulong p) - { - return IsEvenPerfectCandidate(p, out _, out _); - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static bool IsEvenPerfectCandidate(ulong p) + { + return IsEvenPerfectCandidate(p, 0, out _, out _); + } internal static bool IsEvenPerfectCandidate(ulong p, int divisorCyclesSearchLimit, out bool searchedMersenne, out bool detailedCheck) { From ab856c49950e577641f2d5eb041442b07c77b8b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 11 Sep 2025 18:38:42 +0200 Subject: [PATCH 2/2] Make GpuUInt128 mutable and update callers --- .../Gpu/GpuUInt128Pow2Minus1ModTests.cs | 12 +- .../Gpu/GpuUInt128Tests.cs | 34 +-- PerfectNumbers.Core/Gpu/GpuUInt128.cs | 276 ++++++++++-------- .../Gpu/MersenneNumberDivisorGpuTester.cs | 11 +- .../Gpu/MersenneNumberLucasLehmerGpuTester.cs | 19 +- PerfectNumbers.Core/Gpu/NttGpuMath.cs | 118 ++++---- 6 files changed, 267 insertions(+), 203 deletions(-) diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Pow2Minus1ModTests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Pow2Minus1ModTests.cs index f8d9c28..b3dc8a1 100644 --- a/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Pow2Minus1ModTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Pow2Minus1ModTests.cs @@ -13,7 +13,17 @@ public class GpuUInt128Pow2Minus1ModTests ulong exponent = 23UL; GpuUInt128 modulus = new(0UL, 97UL); var pow = GpuUInt128.Pow2Mod(exponent, modulus); - var expected = pow.IsZero ? modulus.Sub(GpuUInt128.One) : pow.Sub(GpuUInt128.One); + GpuUInt128 expected; + if (pow.IsZero) + { + expected = modulus; + expected.Sub(GpuUInt128.One); + } + else + { + expected = pow; + expected.Sub(GpuUInt128.One); + } var actual = GpuUInt128.Pow2Minus1Mod(exponent, modulus); actual.Should().Be(expected); } diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs index 04b68a1..fb2636f 100644 --- a/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs @@ -18,7 +18,7 @@ public class GpuUInt128Tests var a = new GpuUInt128(high1, low1); var b = new GpuUInt128(high2, low2); var expected = (UInt128)new GpuUInt128(high1, low1) * (UInt128)new GpuUInt128(high2, low2); - a = a.Mul(b); + a.Mul(b); ((UInt128)a).Should().Be(expected); } } @@ -39,7 +39,7 @@ public class GpuUInt128Tests var a = new GpuUInt128((UInt128)new GpuUInt128(high1, low1) % (UInt128)modulus); var b = new GpuUInt128((UInt128)new GpuUInt128(high2, low2) % (UInt128)modulus); var expected = MulModReference((UInt128)a, (UInt128)b, (UInt128)modulus); - a = a.MulMod(b, modulus); + a.MulMod(b, modulus); ((UInt128)a).Should().Be(expected); } } @@ -60,7 +60,7 @@ public class GpuUInt128Tests UInt128 bBig = (UInt128)bRaw % (UInt128)modulus; ulong b = (ulong)bBig; var expected = MulModReference((UInt128)a, bBig, (UInt128)modulus); - a = a.MulMod(b, modulus); + a.MulMod(b, modulus); ((UInt128)a).Should().Be(expected); } } @@ -79,7 +79,7 @@ public class GpuUInt128Tests var a = new GpuUInt128(aVal); var b = new GpuUInt128(bVal); UInt128 expected = ((UInt128)aVal * bVal) % modulus; - a = a.MulMod(b, modulus); + a.MulMod(b, modulus); ((UInt128)a).Should().Be(expected); } } @@ -101,7 +101,7 @@ public class GpuUInt128Tests var b = new GpuUInt128((UInt128)new GpuUInt128(high2, low2) % (UInt128)modulus); var expected = MulModReference((UInt128)a, (UInt128)b, (UInt128)modulus); var result = new GpuUInt128(a.High, a.Low); - result = result.MulModBigInteger(b, modulus); + result.MulModBigInteger(b, modulus); ((UInt128)result).Should().Be(expected); } } @@ -119,7 +119,7 @@ public class GpuUInt128Tests var modulus = new GpuUInt128(modHigh, modLow); var a = new GpuUInt128((UInt128)new GpuUInt128(valueHigh, valueLow) % (UInt128)modulus); var expected = MulModReference((UInt128)a, (UInt128)a, (UInt128)modulus); - a = a.SquareMod(modulus); + a.SquareMod(modulus); ((UInt128)a).Should().Be(expected); } } @@ -142,7 +142,7 @@ public class GpuUInt128Tests var exponent = new GpuUInt128(expHigh, expLow); var expected = (UInt128)System.Numerics.BigInteger.ModPow((System.Numerics.BigInteger)(UInt128)value, (System.Numerics.BigInteger)(UInt128)exponent, (System.Numerics.BigInteger)(UInt128)modulus); - value = value.ModPow(exponent, modulus); + value.ModPow(exponent, modulus); ((UInt128)value).Should().Be(expected); } } @@ -163,7 +163,7 @@ public class GpuUInt128Tests var value = new GpuUInt128((UInt128)new GpuUInt128(baseHigh, baseLow) % (UInt128)modulus); var expected = (UInt128)System.Numerics.BigInteger.ModPow((System.Numerics.BigInteger)(UInt128)value, exponent, (System.Numerics.BigInteger)(UInt128)modulus); - value = value.ModPow(exponent, modulus); + value.ModPow(exponent, modulus); ((UInt128)value).Should().Be(expected); } } @@ -184,7 +184,7 @@ public class GpuUInt128Tests var a = new GpuUInt128((UInt128)new GpuUInt128(valueHigh, valueLow) % (UInt128)modulus); var b = new GpuUInt128((UInt128)new GpuUInt128(subHigh, subLow) % (UInt128)modulus); UInt128 expected = ((UInt128)a + (UInt128)modulus - (UInt128)b) % (UInt128)modulus; - a = a.SubMod(b, modulus); + a.SubMod(b, modulus); ((UInt128)a).Should().Be(expected); } } @@ -205,7 +205,7 @@ public class GpuUInt128Tests UInt128 bBig = (UInt128)subRaw % (UInt128)modulus; ulong b = (ulong)bBig; UInt128 expected = ((UInt128)a + (UInt128)modulus - bBig) % (UInt128)modulus; - a = a.SubMod(b, modulus); + a.SubMod(b, modulus); ((UInt128)a).Should().Be(expected); } } @@ -226,7 +226,7 @@ public class GpuUInt128Tests var a = new GpuUInt128((UInt128)new GpuUInt128(valueHigh, valueLow) % (UInt128)modulus); var b = new GpuUInt128((UInt128)new GpuUInt128(addHigh, addLow) % (UInt128)modulus); UInt128 expected = ((UInt128)a + (UInt128)b) % (UInt128)modulus; - a = a.AddMod(b, modulus); + a.AddMod(b, modulus); ((UInt128)a).Should().Be(expected); } } @@ -245,7 +245,7 @@ public class GpuUInt128Tests var a = new GpuUInt128(aVal); var b = new GpuUInt128(bVal); UInt128 expected = ((UInt128)aVal + bVal) % modulus; - a = a.AddMod(b, modulus); + a.AddMod(b, modulus); ((UInt128)a).Should().Be(expected); } } @@ -263,7 +263,7 @@ public class GpuUInt128Tests ulong bVal = bRaw % modulus; var a = new GpuUInt128(aVal); UInt128 expected = ((UInt128)aVal + bVal) % modulus; - a = a.AddMod(bVal, modulus); + a.AddMod(bVal, modulus); ((UInt128)a).Should().Be(expected); } } @@ -282,11 +282,11 @@ public class GpuUInt128Tests var value = new GpuUInt128((UInt128)new GpuUInt128(valueHigh, valueLow) % (UInt128)modulus); if (value.IsZero) { - value = value.Add(1UL); + value.Add(1UL); } var expected = (UInt128)System.Numerics.BigInteger.ModPow((System.Numerics.BigInteger)(UInt128)value, (System.Numerics.BigInteger)(UInt128)modulus - 2, (System.Numerics.BigInteger)(UInt128)modulus); - value = value.ModInv(modulus); + value.ModInv(modulus); ((UInt128)value).Should().Be(expected); } } @@ -302,11 +302,11 @@ public class GpuUInt128Tests var value = new GpuUInt128((UInt128)(raw % modulus)); if (value.IsZero) { - value = value.Add(1UL); + value.Add(1UL); } var expected = (UInt128)System.Numerics.BigInteger.ModPow((System.Numerics.BigInteger)(UInt128)value, (System.Numerics.BigInteger)modulus - 2, (System.Numerics.BigInteger)modulus); - value = value.ModInv(modulus); + value.ModInv(modulus); ((UInt128)value).Should().Be(expected); } } diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs index 5ccd4b1..725780f 100644 --- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs +++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs @@ -3,11 +3,11 @@ using System.Runtime.CompilerServices; namespace PerfectNumbers.Core.Gpu; -public readonly struct GpuUInt128 : IComparable, IEquatable +public struct GpuUInt128 : IComparable, IEquatable { - public readonly ulong High; + public ulong High; - public readonly ulong Low; + public ulong Low; public static readonly GpuUInt128 Zero = new(0UL, 0UL); public static readonly GpuUInt128 One = new(0UL, 1UL); @@ -129,7 +129,11 @@ public readonly struct GpuUInt128 : IComparable, IEquatable=(GpuUInt128 left, GpuUInt128 right) => left.CompareTo(right) >= 0; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static GpuUInt128 operator ++(GpuUInt128 value) => value.Add(1UL); + public static GpuUInt128 operator ++(GpuUInt128 value) + { + value.Add(1UL); + return value; + } [MethodImpl(MethodImplOptions.AggressiveInlining)] public readonly int CompareTo(GpuUInt128 other) @@ -167,51 +171,51 @@ public readonly struct GpuUInt128 : IComparable, IEquatable HashCode.Combine(High, Low); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 Add(GpuUInt128 other) + public void Add(GpuUInt128 other) { ulong low = Low + other.Low; - return new(High + other.High + (low < Low ? 1UL : 0UL), low); + High = High + other.High + (low < Low ? 1UL : 0UL); + Low = low; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 Add(ulong value) + public void Add(ulong value) { - var low = Low + value; - return new(High + (low < Low ? 1UL : 0UL), low); + ulong low = Low + value; + High = High + (low < Low ? 1UL : 0UL); + Low = low; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 AddMod(GpuUInt128 value, GpuUInt128 modulus) + public void AddMod(GpuUInt128 value, GpuUInt128 modulus) { if (modulus.High == 0UL) { - return AddMod(value, modulus.Low); + AddMod(value, modulus.Low); + return; } - value = Add(value); - if (value.CompareTo(modulus) >= 0) + Add(value); + if (CompareTo(modulus) >= 0) { - value = value.Sub(modulus); + Sub(modulus); } - - return value; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 AddMod(ulong value, GpuUInt128 modulus) + public void AddMod(ulong value, GpuUInt128 modulus) { if (modulus.High == 0UL) { - return AddMod(value, modulus.Low); + AddMod(value, modulus.Low); + return; } - var result = Add(value); - if (result.CompareTo(modulus) >= 0) + Add(value); + if (CompareTo(modulus) >= 0) { - result = result.Sub(modulus); + Sub(modulus); } - - return result; } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -336,7 +340,7 @@ public readonly struct GpuUInt128 : IComparable, IEquatable, IEquatable, IEquatable, IEquatable= modulus) { - a = a.Sub(modulus); + a.Sub(modulus); } return a; @@ -424,7 +431,7 @@ public readonly struct GpuUInt128 : IComparable, IEquatable= modulus) { - a = a.Sub(modulus); + a.Sub(modulus); } // b >>= 1 @@ -516,76 +523,95 @@ public readonly struct GpuUInt128 : IComparable, IEquatable= 0) { - return Sub(value); + Sub(value); } else { - return Add(modulus).Sub(value); + Add(modulus); + Sub(value); } } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 SubMod(ulong value, GpuUInt128 modulus) + public void SubMod(ulong value, GpuUInt128 modulus) { - ulong low, 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; - } + ulong low, 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; + } ulong borrow = low < value ? 1UL : 0UL; - return new (high - borrow, low - value); + High = high - borrow; + Low = low - value; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 Mul(GpuUInt128 other) + public void Mul(GpuUInt128 other) { Multiply(this, other, out var high, out var low); - return new(high, low); + High = high; + Low = low; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 Mul(ulong value) + public void Mul(ulong value) { var (highPart, lowPart) = Mul64(Low, value); - return new(highPart + High * value, lowPart); + High = highPart + High * value; + Low = lowPart; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 Xor(GpuUInt128 other) => new(High ^ other.High, Low ^ other.Low); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Xor(GpuUInt128 other) + { + High ^= other.High; + Low ^= other.Low; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 Xor(ulong value) => new(High, Low ^ value); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Xor(ulong value) + { + Low ^= value; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 ShiftLeft(int shift) => this << shift; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void ShiftLeft(int shift) + { + this = this << shift; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 ShiftRight(int shift) => this >> shift; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void ShiftRight(int shift) + { + this = this >> shift; + } - public GpuUInt128 MulMod(GpuUInt128 other, GpuUInt128 modulus) + public void MulMod(GpuUInt128 other, GpuUInt128 modulus) { if (modulus.High == 0UL) { - return MulMod(other, modulus.Low); + MulMod(other, modulus.Low); + return; } MultiplyFull(this, other, out var p3, out var p2, out var p1, out var p0); @@ -598,7 +624,7 @@ public readonly struct GpuUInt128 : IComparable, IEquatable> bit) & 1UL)); if (remainder.CompareTo(modulus) >= 0) { - remainder = remainder.Sub(modulus); + remainder.Sub(modulus); } } @@ -608,7 +634,7 @@ public readonly struct GpuUInt128 : IComparable, IEquatable> bit) & 1UL)); if (remainder.CompareTo(modulus) >= 0) { - remainder = remainder.Sub(modulus); + remainder.Sub(modulus); } } @@ -618,7 +644,7 @@ public readonly struct GpuUInt128 : IComparable, IEquatable> bit) & 1UL)); if (remainder.CompareTo(modulus) >= 0) { - remainder = remainder.Sub(modulus); + remainder.Sub(modulus); } } @@ -628,11 +654,12 @@ public readonly struct GpuUInt128 : IComparable, IEquatable> bit) & 1UL)); if (remainder.CompareTo(modulus) >= 0) { - remainder = remainder.Sub(modulus); + remainder.Sub(modulus); } } - return remainder; + High = remainder.High; + Low = remainder.Low; // TODO(MOD-OPT): Replace this bitwise long-division style reduction // with a faster algorithm suitable for GPU kernels without '%' support. // Options: @@ -653,13 +680,14 @@ public readonly struct GpuUInt128 : IComparable, IEquatable - internal GpuUInt128 MulModBigInteger(GpuUInt128 other, GpuUInt128 modulus) + internal void MulModBigInteger(GpuUInt128 other, GpuUInt128 modulus) { var left = (BigInteger)(UInt128)this; var right = (BigInteger)(UInt128)other; var mod = (BigInteger)(UInt128)modulus; var reduced = (UInt128)((left * right) % mod); - return new((ulong)(reduced >> 64), (ulong)reduced); + High = (ulong)(reduced >> 64); + Low = (ulong)reduced; } /// @@ -667,14 +695,14 @@ public readonly struct GpuUInt128 : IComparable, IEquatable - internal GpuUInt128 MulModByLimb(GpuUInt128 other, GpuUInt128 modulus) + internal void MulModByLimb(GpuUInt128 other, GpuUInt128 modulus) { MultiplyFull(this, other, out var p3, out var p2, out var p1, out var p0); GpuUInt128 remainder = new(p3, p2); while (remainder.CompareTo(modulus) >= 0) { - remainder = remainder.Sub(modulus); + remainder.Sub(modulus); } ulong limb = p1; @@ -684,33 +712,32 @@ public readonly struct GpuUInt128 : IComparable, IEquatable= 0) { - remainder = remainder.Sub(modulus); + remainder.Sub(modulus); } limb = p0; } - return remainder; + High = remainder.High; + Low = remainder.Low; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 MulMod(ulong value, GpuUInt128 modulus) + public void MulMod(ulong value, GpuUInt128 modulus) { if (modulus.High == 0UL) { - return MulMod(value, modulus.Low); + MulMod(value, modulus.Low); + return; } GpuUInt128 other = new(0UL, value); - return MulMod(other, modulus); + MulMod(other, modulus); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 MulMod(GpuUInt128 other, ulong modulus) + public void MulMod(GpuUInt128 other, ulong modulus) { - // Fast-path: use single 64-bit multiply when it cannot overflow - // (enabled by ILGPU support for 64-bit % in kernels). Otherwise fall - // back to the generic branchless double-and-add implementation. ulong a = Low % modulus; ulong b = other.Low % modulus; ulong result; @@ -720,7 +747,6 @@ public readonly struct GpuUInt128 : IComparable, IEquatable, IEquatable, IEquatable, IEquatable> bit) & 1UL)); if (remainder.CompareTo(modulus) >= 0) { - remainder = remainder.Sub(modulus); + remainder.Sub(modulus); } } @@ -822,7 +850,7 @@ public readonly struct GpuUInt128 : IComparable, IEquatable> bit) & 1UL)); if (remainder.CompareTo(modulus) >= 0) { - remainder = remainder.Sub(modulus); + remainder.Sub(modulus); } } @@ -832,7 +860,7 @@ public readonly struct GpuUInt128 : IComparable, IEquatable> bit) & 1UL)); if (remainder.CompareTo(modulus) >= 0) { - remainder = remainder.Sub(modulus); + remainder.Sub(modulus); } } @@ -842,14 +870,15 @@ public readonly struct GpuUInt128 : IComparable, IEquatable> bit) & 1UL)); if (remainder.CompareTo(modulus) >= 0) { - remainder = remainder.Sub(modulus); + remainder.Sub(modulus); } } - return remainder; + High = remainder.High; + Low = remainder.Low; } - public GpuUInt128 ModPow(GpuUInt128 exponent, GpuUInt128 modulus) + public void ModPow(GpuUInt128 exponent, GpuUInt128 modulus) { GpuUInt128 result = new(0UL, 1UL); GpuUInt128 baseValue = this; @@ -858,17 +887,18 @@ public readonly struct GpuUInt128 : IComparable, IEquatable>= 1; - baseValue = baseValue.MulMod(baseValue, modulus); + baseValue.MulMod(baseValue, modulus); } - return result; + High = result.High; + Low = result.Low; } - public GpuUInt128 ModPow(ulong exponent, GpuUInt128 modulus) + public void ModPow(ulong exponent, GpuUInt128 modulus) { GpuUInt128 result = new(0UL, 1UL); GpuUInt128 baseValue = this; @@ -877,29 +907,35 @@ public readonly struct GpuUInt128 : IComparable, IEquatable>= 1; - baseValue = baseValue.MulMod(baseValue, modulus); + baseValue.MulMod(baseValue, modulus); } - return result; + High = result.High; + Low = result.Low; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 ModInv(ulong modulus) => ModPow(modulus - 2UL, new GpuUInt128(0UL, modulus)); + public void ModInv(ulong modulus) + { + ModPow(modulus - 2UL, new GpuUInt128(0UL, modulus)); + } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public GpuUInt128 ModInv(GpuUInt128 modulus) + public void ModInv(GpuUInt128 modulus) { if (modulus.High == 0UL) { - return ModInv(modulus.Low); + ModInv(modulus.Low); + return; } - var exponent = new GpuUInt128(modulus).Sub(2UL); - return ModPow(exponent, modulus); + var exponent = new GpuUInt128(modulus); + exponent.Sub(2UL); + ModPow(exponent, modulus); } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs index a6a38d0..59b5d59 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs @@ -48,10 +48,11 @@ public sealed class MersenneNumberDivisorGpuTester ulong pow2x = 1UL << x; ulong y = divisor - pow2x; var modVal = new GpuUInt128(0UL, divisor); - var baseVal = new GpuUInt128(0UL, pow2x); // 2^x ≡ -y (mod divisor) - var part1 = baseVal.ModPow(q, modVal); - var part2 = GpuUInt128.Pow2Mod(r, modVal); - var rem = part1.MulMod(part2, modVal); - result[0] = rem.High == 0UL && rem.Low == 1UL ? (byte)1 : (byte)0; + var baseVal = new GpuUInt128(0UL, pow2x); // 2^x ≡ -y (mod divisor) + var part1 = baseVal; + part1.ModPow(q, modVal); + var part2 = GpuUInt128.Pow2Mod(r, modVal); + part1.MulMod(part2, modVal); + result[0] = part1.High == 0UL && part1.Low == 1UL ? (byte)1 : (byte)0; } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs index 6aa29f4..f169ce7 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs @@ -526,7 +526,7 @@ public class MersenneNumberLucasLehmerGpuTester int target = i - limbCount; var sum = value[target]; var original = sum; - sum = sum.Add(limb); + sum.Add(limb); value[target] = sum; value[i] = new GpuUInt128(0UL, 0UL); if (sum.CompareTo(original) < 0) @@ -536,7 +536,7 @@ public class MersenneNumberLucasLehmerGpuTester { var next = value[j]; var nextOriginal = next; - next = next.Add(1UL); + next.Add(1UL); value[j] = next; if (next.CompareTo(nextOriginal) >= 0) { @@ -569,7 +569,7 @@ public class MersenneNumberLucasLehmerGpuTester { var cur = value[j]; var curOriginal = cur; - cur = cur.Add(carryBits); + cur.Add(carryBits); value[j] = cur; carryBits = cur.CompareTo(curOriginal) < 0 ? new GpuUInt128(0UL, 1UL) : new GpuUInt128(0UL, 0UL); j++; @@ -756,7 +756,8 @@ public class MersenneNumberLucasLehmerGpuTester ulong limit = exponent - 2UL; for (ulong i = 0UL; i < limit; i++) { - s = SquareModMersenne128(s, exponent).SubMod(2UL, modulus); + s = SquareModMersenne128(s, exponent); + s.SubMod(2UL, modulus); } states[idx] = s; @@ -773,14 +774,16 @@ public class MersenneNumberLucasLehmerGpuTester { for (; i < limit; i++) { - s = SquareModMersenne128(s, exponent).SubMod(2UL, modulus); + s = SquareModMersenne128(s, exponent); + s.SubMod(2UL, modulus); } } else { for (; i < limit; i++) { - s = s.SquareMod(modulus).SubMod(2UL, modulus); + s.SquareMod(modulus); + s.SubMod(2UL, modulus); } } @@ -797,7 +800,7 @@ public class MersenneNumberLucasLehmerGpuTester // v = low128 + high128 (fold by 128 bits once) var vHigh = new GpuUInt128(q3, q2); var vLow = new GpuUInt128(q1, q0); - vLow = vLow.Add(vHigh); + vLow.Add(vHigh); // Mask to keep only p low bits and fold the carry bits (v >> p) int topBits = (int)(p % 128UL); @@ -817,7 +820,7 @@ public class MersenneNumberLucasLehmerGpuTester while (!carryBits.IsZero) { before = vLow; - vLow = vLow.Add(carryBits); + vLow.Add(carryBits); carryBits = vLow.CompareTo(before) < 0 ? GpuUInt128.One : GpuUInt128.Zero; } diff --git a/PerfectNumbers.Core/Gpu/NttGpuMath.cs b/PerfectNumbers.Core/Gpu/NttGpuMath.cs index 5466748..86f217d 100644 --- a/PerfectNumbers.Core/Gpu/NttGpuMath.cs +++ b/PerfectNumbers.Core/Gpu/NttGpuMath.cs @@ -77,14 +77,24 @@ public static class NttGpuMath Buffer = accelerator.Allocate1D(length); Root = new GpuUInt128(primitiveRoot); RootInv = new GpuUInt128(primitiveRoot); - RootInv = modulus.High == 0UL - ? RootInv.ModInv(modulus.Low) - : RootInv.ModInv(modulus); + if (modulus.High == 0UL) + { + RootInv.ModInv(modulus.Low); + } + else + { + RootInv.ModInv(modulus); + } NInv = new GpuUInt128((UInt128)length); - NInv = modulus.High == 0UL - ? NInv.ModInv(modulus.Low) - : NInv.ModInv(modulus); + if (modulus.High == 0UL) + { + NInv.ModInv(modulus.Low); + } + else + { + NInv.ModInv(modulus); + } // Precompute stage twiddle factors (forward and inverse) and upload to GPU. Bits = (int)Math.Log2(length); @@ -102,20 +112,20 @@ public static class NttGpuMath StageOffsets[stage - 1] = offset; UInt128 exp = (modValue - 1UL) / (ulong)len; var wLen = new GpuUInt128(Root); - wLen = wLen.ModPow((ulong)exp, modulus); + wLen.ModPow((ulong)exp, modulus); var w = new GpuUInt128(0UL, 1UL); for (int j = 0; j < half; j++) { forward[offset + j] = w; - w = w.MulMod(wLen, modulus); + w.MulMod(wLen, modulus); } var wLenInv = new GpuUInt128(RootInv); - wLenInv = wLenInv.ModPow((ulong)exp, modulus); + wLenInv.ModPow((ulong)exp, modulus); var wi = new GpuUInt128(0UL, 1UL); for (int j = 0; j < half; j++) { inverse[offset + j] = wi; - wi = wi.MulMod(wLenInv, modulus); + wi.MulMod(wLenInv, modulus); } offset += half; } @@ -393,7 +403,9 @@ public static class NttGpuMath private static void MulKernel(Index1D index, ArrayView a, ArrayView b, GpuUInt128 modulus) { - a[index] = a[index].MulMod(b[index], modulus); + var val = a[index]; + val.MulMod(b[index], modulus); + a[index] = val; } private static void StageKernel(Index1D index, ArrayView data, int len, int half, int stageOffset, ArrayView twiddles, GpuUInt128 modulus) @@ -410,10 +422,10 @@ public static class NttGpuMath var u = data[i1]; var v = data[i2]; var w = twiddles[stageOffset + j]; - v = v.MulMod(w, modulus); + v.MulMod(w, modulus); var sum = new GpuUInt128(u); - sum = sum.AddMod(v, modulus); - u = u.SubMod(v, modulus); + sum.AddMod(v, modulus); + u.SubMod(v, modulus); data[i1] = sum; data[i2] = u; } @@ -421,7 +433,7 @@ public static class NttGpuMath private static void ScaleKernel(Index1D index, ArrayView data, GpuUInt128 scale, GpuUInt128 modulus) { var v = data[index]; - v = v.MulMod(scale, modulus); + v.MulMod(scale, modulus); data[index] = v; } @@ -515,8 +527,8 @@ public static class NttGpuMath var sum = new GpuUInt128(u); var mulred = new GpuUInt128(resHigh, resLow); - sum = sum.AddMod(mulred, new GpuUInt128(modHigh, modLow)); - u = u.SubMod(mulred, new GpuUInt128(modHigh, modLow)); + sum.AddMod(mulred, new GpuUInt128(modHigh, modLow)); + u.SubMod(mulred, new GpuUInt128(modHigh, modLow)); data[i1] = sum; data[i2] = u; } @@ -735,9 +747,9 @@ public static class NttGpuMath { var term = input[i]; var power = new GpuUInt128(root); - power = power.ModPow((ulong)index.X * (ulong)i, modulus); - term = term.MulMod(power, modulus); - sum = sum.AddMod(term, modulus); + power.ModPow((ulong)index.X * (ulong)i, modulus); + term.MulMod(power, modulus); + sum.AddMod(term, modulus); } output[index] = sum; @@ -753,12 +765,12 @@ public static class NttGpuMath { var term = input[i]; var power = new GpuUInt128(rootInv); - power = power.ModPow((ulong)index.X * (ulong)i, modulus); - term = term.MulMod(power, modulus); - sum = sum.AddMod(term, modulus); + power.ModPow((ulong)index.X * (ulong)i, modulus); + term.MulMod(power, modulus); + sum.AddMod(term, modulus); } - sum = sum.MulMod(nInv, modulus); + sum.MulMod(nInv, modulus); output[index] = sum; } @@ -831,7 +843,7 @@ public static class NttGpuMath { UInt128 expBase = (modValue - 1UL) / (ulong)len; var wLen = new GpuUInt128(primitiveRoot); - wLen = wLen.ModPow((ulong)expBase, modulus); + wLen.ModPow((ulong)expBase, modulus); for (int i = 0; i < n; i += len) { var w = new GpuUInt128(0UL, 1UL); @@ -840,13 +852,13 @@ public static class NttGpuMath { var u = values[i + j]; var v = values[i + j + half]; - v = v.MulMod(w, modulus); + v.MulMod(w, modulus); var sum = new GpuUInt128(u); - sum = sum.AddMod(v, modulus); - u = u.SubMod(v, modulus); + sum.AddMod(v, modulus); + u.SubMod(v, modulus); values[i + j] = sum; values[i + j + half] = u; - w = w.MulMod(wLen, modulus); + w.MulMod(wLen, modulus); } } } @@ -879,7 +891,7 @@ public static class NttGpuMath UInt128 modValue = modulus; UInt128 expBase = (modValue - 1UL) / (ulong)n; var root = new GpuUInt128(primitiveRoot); - root = root.ModPow((ulong)expBase, modulus); + root.ModPow((ulong)expBase, modulus); kernel(n, inputBuffer.View, outputBuffer.View, n, modulus, root); accelerator.Synchronize(); outputBuffer.View.CopyToCPU(ref outputArray[0], n); @@ -900,33 +912,33 @@ public static class NttGpuMath var root = new GpuUInt128(primitiveRoot); if (modulus.High == 0UL) { - root = root.ModInv(modulus.Low); + root.ModInv(modulus.Low); } else { - root = root.ModInv(modulus); + root.ModInv(modulus); } for (int len = 2; len <= n; len <<= 1) { - UInt128 expBase = (modValue - 1UL) / (ulong)len; - var wLen = new GpuUInt128(root); - wLen = wLen.ModPow((ulong)expBase, modulus); + UInt128 expBase = (modValue - 1UL) / (ulong)len; + var wLen = new GpuUInt128(root); + wLen.ModPow((ulong)expBase, modulus); for (int i = 0; i < n; i += len) { var w = new GpuUInt128(0UL, 1UL); int half = len >> 1; for (int j = 0; j < half; j++) { - var u = values[i + j]; - var v = values[i + j + half]; - v = v.MulMod(w, modulus); - var sum = new GpuUInt128(u); - sum = sum.AddMod(v, modulus); - u = u.SubMod(v, modulus); - values[i + j] = sum; - values[i + j + half] = u; - w = w.MulMod(wLen, modulus); + var u = values[i + j]; + var v = values[i + j + half]; + v.MulMod(w, modulus); + var sum = new GpuUInt128(u); + sum.AddMod(v, modulus); + u.SubMod(v, modulus); + values[i + j] = sum; + values[i + j + half] = u; + w.MulMod(wLen, modulus); } } } @@ -934,16 +946,18 @@ public static class NttGpuMath var nInv = new GpuUInt128((UInt128)n); if (modulus.High == 0UL) { - nInv = nInv.ModInv(modulus.Low); + nInv.ModInv(modulus.Low); } else { - nInv = nInv.ModInv(modulus); + nInv.ModInv(modulus); } for (int i = 0; i < n; i++) { - values[i] = values[i].MulMod(nInv, modulus); + var val = values[i]; + val.MulMod(nInv, modulus); + values[i] = val; } } @@ -974,25 +988,25 @@ public static class NttGpuMath UInt128 modValue = modulus; UInt128 expBase = (modValue - 1UL) / (ulong)n; var root = new GpuUInt128(primitiveRoot); - root = root.ModPow((ulong)expBase, modulus); + root.ModPow((ulong)expBase, modulus); var rootInv = new GpuUInt128(root); if (modulus.High == 0UL) { - rootInv = rootInv.ModInv(modulus.Low); + rootInv.ModInv(modulus.Low); } else { - rootInv = rootInv.ModInv(modulus); + rootInv.ModInv(modulus); } var nInv = new GpuUInt128((UInt128)n); if (modulus.High == 0UL) { - nInv = nInv.ModInv(modulus.Low); + nInv.ModInv(modulus.Low); } else { - nInv = nInv.ModInv(modulus); + nInv.ModInv(modulus); } kernel(n, inputBuffer.View, outputBuffer.View, n, modulus, rootInv, nInv); accelerator.Synchronize();