mlemanczyk/even-perfect-numbers-scanner
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #14 from mlemanczyk/codex/restore-mutable-gpuuint128-structure

Browse Source 
Switch GpuUInt128 back to mutable form
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-09-11 18:39:04 +02:00
2 parents 40477eca5d ab856c4995
commit 5af55131d1
6 files changed
+267 -203
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core.Tests/Gpu/GpuUInt128Pow2Minus1ModTests.cs
+11 -1
View File
@@ -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);
}
PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs
+17 -17
View File
@@ -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);
}
}
PerfectNumbers.Core/Gpu/GpuUInt128.cs
+156 -120
View File
@@ -3,11 +3,11 @@ using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core.Gpu;
public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
{
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<GpuUInt128>, IEquatable<GpuUInt1
public static bool operator >=(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<GpuUInt128>, IEquatable<GpuUInt1
public override readonly int GetHashCode() => 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<GpuUInt128>, IEquatable<GpuUInt1
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public GpuUInt128 AddMod(GpuUInt128 value, ulong modulus)
public void AddMod(GpuUInt128 value, ulong modulus)
{
ulong a = Low % modulus;
ulong b = value.Low % modulus;
@@ -346,11 +350,12 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
sum -= modulus;
}
return new(0UL, sum);
High = 0UL;
Low = sum;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public GpuUInt128 AddMod(ulong value, ulong modulus)
public void AddMod(ulong value, ulong modulus)
{
ulong a = Low % modulus;
ulong b = value % modulus;
@@ -360,14 +365,16 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
sum -= modulus;
}
return new(0UL, sum);
High = 0UL;
Low = sum;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public GpuUInt128 Sub(GpuUInt128 other)
public void Sub(GpuUInt128 other)
{
ulong borrow = Low < other.Low ? 1UL : 0UL;
return new(High - other.High - borrow, Low - other.Low);
High = High - other.High - borrow;
Low -= other.Low;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -400,10 +407,10 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GpuUInt128 AddMod(GpuUInt128 a, GpuUInt128 b, GpuUInt128 modulus)
{
a = a.Add(b);
a.Add(b);
if (a >= modulus)
{
a = a.Sub(modulus);
a.Sub(modulus);
}
return a;
@@ -424,7 +431,7 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
a <<= 1;
if (a >= modulus)
{
a = a.Sub(modulus);
a.Sub(modulus);
}
// b >>= 1
@@ -516,76 +523,95 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
}
public GpuUInt128 Sub(ulong value)
public void Sub(ulong value)
{
ulong borrow = Low < value ? 1UL : 0UL;
return new(High - borrow, Low - value);
High -= borrow;
Low -= value;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public GpuUInt128 SubMod(GpuUInt128 value, GpuUInt128 modulus)
public void SubMod(GpuUInt128 value, GpuUInt128 modulus)
{
if (CompareTo(value) >= 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<GpuUInt128>, IEquatable<GpuUInt1
remainder = new(remainder.High, remainder.Low | ((p3 >> bit) & 1UL));
if (remainder.CompareTo(modulus) >= 0)
{
remainder = remainder.Sub(modulus);
remainder.Sub(modulus);
}
}
@@ -608,7 +634,7 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
remainder = new(remainder.High, remainder.Low | ((p2 >> bit) & 1UL));
if (remainder.CompareTo(modulus) >= 0)
{
remainder = remainder.Sub(modulus);
remainder.Sub(modulus);
}
}
@@ -618,7 +644,7 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
remainder = new(remainder.High, remainder.Low | ((p1 >> bit) & 1UL));
if (remainder.CompareTo(modulus) >= 0)
{
remainder = remainder.Sub(modulus);
remainder.Sub(modulus);
}
}
@@ -628,11 +654,12 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
remainder = new(remainder.High, remainder.Low | ((p0 >> 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<GpuUInt128>, IEquatable<GpuUInt1
/// method is intended for validation only and should not be used inside
/// GPU kernels.
/// </summary>
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;
}
/// <summary>
@@ -667,14 +695,14 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
/// repeated subtractions and becomes extremely slow for large remainders.
/// Kept for future optimization work.
/// </summary>
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<GpuUInt128>, IEquatable<GpuUInt1
remainder = new(remainder.High, limb);
while (remainder.CompareTo(modulus) >= 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<GpuUInt128>, IEquatable<GpuUInt1
}
else if (b <= ulong.MaxValue / a)
{
// Safe to multiply without overflow
result = (a * b) % modulus;
}
else
@@ -728,13 +754,13 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
result = MulMod64(a, b, modulus);
}
return new(0UL, result);
High = 0UL;
Low = result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public GpuUInt128 MulMod(ulong value, ulong modulus)
public void MulMod(ulong value, ulong modulus)
{
// Same fast-path as above for 64-bit operands
ulong a = Low % modulus;
ulong b = value % modulus;
ulong result;
@@ -751,57 +777,59 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
result = MulMod64(a, b, modulus);
}
return new(0UL, result);
High = 0UL;
Low = result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public GpuUInt128 MulModMontgomery64(GpuUInt128 other, ulong modulus, ulong nPrime, ulong r2)
public void MulModMontgomery64(GpuUInt128 other, ulong modulus, ulong nPrime, ulong r2)
{
// Montgomery-based 64-bit modular multiplication without division.
// Computes (this * other) % modulus using:
// aR = MontMul(a, R^2), bR = MontMul(b, R^2), cR = MontMul(aR, bR), c = MontMul(cR, 1)
// where R = 2^64, nPrime = -mod^{-1} mod 2^64, r2 = R^2 mod mod.
ulong a = Low % modulus;
ulong b = other.Low % modulus;
if (a == 0UL || b == 0UL)
{
return Zero;
High = 0UL;
Low = 0UL;
return;
}
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);
return new(0UL, c);
High = 0UL;
Low = c;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public GpuUInt128 MulModMontgomery64(ulong other, ulong modulus, ulong nPrime, ulong r2)
public void MulModMontgomery64(ulong other, ulong modulus, ulong nPrime, ulong r2)
{
ulong a = Low % modulus;
ulong b = other % modulus;
if (a == 0UL || b == 0UL)
{
return Zero;
High = 0UL;
Low = 0UL;
return;
}
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);
return new(0UL, c);
High = 0UL;
Low = c;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public GpuUInt128 SquareMod(GpuUInt128 modulus)
public void SquareMod(GpuUInt128 modulus)
{
if (modulus.High == 0UL)
{
// Use fast 64-bit path.
return MulMod(Low, modulus.Low);
MulMod(Low, modulus.Low);
return;
}
// Specialized 128-bit squaring with bitwise reduction (faster than generic multiply).
SquareFull(this, out var p3, out var p2, out var p1, out var p0);
GpuUInt128 remainder = new(0UL, 0UL);
@@ -812,7 +840,7 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
remainder = new(remainder.High, remainder.Low | ((p3 >> bit) & 1UL));
if (remainder.CompareTo(modulus) >= 0)
{
remainder = remainder.Sub(modulus);
remainder.Sub(modulus);
}
}
@@ -822,7 +850,7 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
remainder = new(remainder.High, remainder.Low | ((p2 >> bit) & 1UL));
if (remainder.CompareTo(modulus) >= 0)
{
remainder = remainder.Sub(modulus);
remainder.Sub(modulus);
}
}
@@ -832,7 +860,7 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
remainder = new(remainder.High, remainder.Low | ((p1 >> bit) & 1UL));
if (remainder.CompareTo(modulus) >= 0)
{
remainder = remainder.Sub(modulus);
remainder.Sub(modulus);
}
}
@@ -842,14 +870,15 @@ public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt1
remainder = new(remainder.High, remainder.Low | ((p0 >> 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<GpuUInt128>, IEquatable<GpuUInt1
{
if ((exponent.Low & 1UL) != 0UL)
{
result = result.MulMod(baseValue, modulus);
result.MulMod(baseValue, modulus);
}
exponent >>= 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<GpuUInt128>, IEquatable<GpuUInt1
{
if ((exponent & 1UL) != 0UL)
{
result = result.MulMod(baseValue, modulus);
result.MulMod(baseValue, modulus);
}
exponent >>= 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)]
PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs
+6 -5
View File
@@ -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;
}
}
PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs
+11 -8
View File
@@ -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;
}
PerfectNumbers.Core/Gpu/NttGpuMath.cs
+66 -52
View File
@@ -77,14 +77,24 @@ public static class NttGpuMath
Buffer = accelerator.Allocate1D<GpuUInt128>(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<GpuUInt128> a, ArrayView<GpuUInt128> 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<GpuUInt128> data, int len, int half, int stageOffset, ArrayView<GpuUInt128> 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<GpuUInt128> 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();