2158 lines
62 KiB
C#
2158 lines
62 KiB
C#
using System;
|
||
using System.Numerics;
|
||
using System.Runtime.CompilerServices;
|
||
using ILGPU.Algorithms;
|
||
|
||
namespace PerfectNumbers.Core.Gpu;
|
||
|
||
public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
|
||
{
|
||
public ulong High;
|
||
|
||
public ulong Low;
|
||
|
||
private const int NativeModuloChunkBits = 8;
|
||
private const int NativeModuloChunkBitsMinusOne = NativeModuloChunkBits - 1;
|
||
private const int NativeModuloBitMaskTableSize = 1024;
|
||
private const ulong NativeModuloChunkMask = (1UL << NativeModuloChunkBits) - 1UL;
|
||
private static readonly ulong[] NativeModuloBitMasks = CreateNativeModuloBitMasks();
|
||
|
||
private const int Pow2WindowSizeBits = 8;
|
||
private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1);
|
||
private const int Pow2WindowMaxExponent = (1 << Pow2WindowSizeBits) - 1;
|
||
|
||
private struct Pow2OddPowerTable
|
||
{
|
||
public GpuUInt128 Element0;
|
||
public GpuUInt128 Element1;
|
||
public GpuUInt128 Element2;
|
||
public GpuUInt128 Element3;
|
||
public GpuUInt128 Element4;
|
||
public GpuUInt128 Element5;
|
||
public GpuUInt128 Element6;
|
||
public GpuUInt128 Element7;
|
||
public GpuUInt128 Element8;
|
||
public GpuUInt128 Element9;
|
||
public GpuUInt128 Element10;
|
||
public GpuUInt128 Element11;
|
||
public GpuUInt128 Element12;
|
||
public GpuUInt128 Element13;
|
||
public GpuUInt128 Element14;
|
||
public GpuUInt128 Element15;
|
||
public GpuUInt128 Element16;
|
||
public GpuUInt128 Element17;
|
||
public GpuUInt128 Element18;
|
||
public GpuUInt128 Element19;
|
||
public GpuUInt128 Element20;
|
||
public GpuUInt128 Element21;
|
||
public GpuUInt128 Element22;
|
||
public GpuUInt128 Element23;
|
||
public GpuUInt128 Element24;
|
||
public GpuUInt128 Element25;
|
||
public GpuUInt128 Element26;
|
||
public GpuUInt128 Element27;
|
||
public GpuUInt128 Element28;
|
||
public GpuUInt128 Element29;
|
||
public GpuUInt128 Element30;
|
||
public GpuUInt128 Element31;
|
||
public GpuUInt128 Element32;
|
||
public GpuUInt128 Element33;
|
||
public GpuUInt128 Element34;
|
||
public GpuUInt128 Element35;
|
||
public GpuUInt128 Element36;
|
||
public GpuUInt128 Element37;
|
||
public GpuUInt128 Element38;
|
||
public GpuUInt128 Element39;
|
||
public GpuUInt128 Element40;
|
||
public GpuUInt128 Element41;
|
||
public GpuUInt128 Element42;
|
||
public GpuUInt128 Element43;
|
||
public GpuUInt128 Element44;
|
||
public GpuUInt128 Element45;
|
||
public GpuUInt128 Element46;
|
||
public GpuUInt128 Element47;
|
||
public GpuUInt128 Element48;
|
||
public GpuUInt128 Element49;
|
||
public GpuUInt128 Element50;
|
||
public GpuUInt128 Element51;
|
||
public GpuUInt128 Element52;
|
||
public GpuUInt128 Element53;
|
||
public GpuUInt128 Element54;
|
||
public GpuUInt128 Element55;
|
||
public GpuUInt128 Element56;
|
||
public GpuUInt128 Element57;
|
||
public GpuUInt128 Element58;
|
||
public GpuUInt128 Element59;
|
||
public GpuUInt128 Element60;
|
||
public GpuUInt128 Element61;
|
||
public GpuUInt128 Element62;
|
||
public GpuUInt128 Element63;
|
||
public GpuUInt128 Element64;
|
||
public GpuUInt128 Element65;
|
||
public GpuUInt128 Element66;
|
||
public GpuUInt128 Element67;
|
||
public GpuUInt128 Element68;
|
||
public GpuUInt128 Element69;
|
||
public GpuUInt128 Element70;
|
||
public GpuUInt128 Element71;
|
||
public GpuUInt128 Element72;
|
||
public GpuUInt128 Element73;
|
||
public GpuUInt128 Element74;
|
||
public GpuUInt128 Element75;
|
||
public GpuUInt128 Element76;
|
||
public GpuUInt128 Element77;
|
||
public GpuUInt128 Element78;
|
||
public GpuUInt128 Element79;
|
||
public GpuUInt128 Element80;
|
||
public GpuUInt128 Element81;
|
||
public GpuUInt128 Element82;
|
||
public GpuUInt128 Element83;
|
||
public GpuUInt128 Element84;
|
||
public GpuUInt128 Element85;
|
||
public GpuUInt128 Element86;
|
||
public GpuUInt128 Element87;
|
||
public GpuUInt128 Element88;
|
||
public GpuUInt128 Element89;
|
||
public GpuUInt128 Element90;
|
||
public GpuUInt128 Element91;
|
||
public GpuUInt128 Element92;
|
||
public GpuUInt128 Element93;
|
||
public GpuUInt128 Element94;
|
||
public GpuUInt128 Element95;
|
||
public GpuUInt128 Element96;
|
||
public GpuUInt128 Element97;
|
||
public GpuUInt128 Element98;
|
||
public GpuUInt128 Element99;
|
||
public GpuUInt128 Element100;
|
||
public GpuUInt128 Element101;
|
||
public GpuUInt128 Element102;
|
||
public GpuUInt128 Element103;
|
||
public GpuUInt128 Element104;
|
||
public GpuUInt128 Element105;
|
||
public GpuUInt128 Element106;
|
||
public GpuUInt128 Element107;
|
||
public GpuUInt128 Element108;
|
||
public GpuUInt128 Element109;
|
||
public GpuUInt128 Element110;
|
||
public GpuUInt128 Element111;
|
||
public GpuUInt128 Element112;
|
||
public GpuUInt128 Element113;
|
||
public GpuUInt128 Element114;
|
||
public GpuUInt128 Element115;
|
||
public GpuUInt128 Element116;
|
||
public GpuUInt128 Element117;
|
||
public GpuUInt128 Element118;
|
||
public GpuUInt128 Element119;
|
||
public GpuUInt128 Element120;
|
||
public GpuUInt128 Element121;
|
||
public GpuUInt128 Element122;
|
||
public GpuUInt128 Element123;
|
||
public GpuUInt128 Element124;
|
||
public GpuUInt128 Element125;
|
||
public GpuUInt128 Element126;
|
||
public GpuUInt128 Element127;
|
||
|
||
public GpuUInt128 this[int index]
|
||
{
|
||
readonly get
|
||
{
|
||
return GetElement(index);
|
||
}
|
||
|
||
set
|
||
{
|
||
SetElement(index, value);
|
||
}
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private readonly GpuUInt128 GetElement(int index) => index switch
|
||
{
|
||
0 => Element0,
|
||
1 => Element1,
|
||
2 => Element2,
|
||
3 => Element3,
|
||
4 => Element4,
|
||
5 => Element5,
|
||
6 => Element6,
|
||
7 => Element7,
|
||
8 => Element8,
|
||
9 => Element9,
|
||
10 => Element10,
|
||
11 => Element11,
|
||
12 => Element12,
|
||
13 => Element13,
|
||
14 => Element14,
|
||
15 => Element15,
|
||
16 => Element16,
|
||
17 => Element17,
|
||
18 => Element18,
|
||
19 => Element19,
|
||
20 => Element20,
|
||
21 => Element21,
|
||
22 => Element22,
|
||
23 => Element23,
|
||
24 => Element24,
|
||
25 => Element25,
|
||
26 => Element26,
|
||
27 => Element27,
|
||
28 => Element28,
|
||
29 => Element29,
|
||
30 => Element30,
|
||
31 => Element31,
|
||
32 => Element32,
|
||
33 => Element33,
|
||
34 => Element34,
|
||
35 => Element35,
|
||
36 => Element36,
|
||
37 => Element37,
|
||
38 => Element38,
|
||
39 => Element39,
|
||
40 => Element40,
|
||
41 => Element41,
|
||
42 => Element42,
|
||
43 => Element43,
|
||
44 => Element44,
|
||
45 => Element45,
|
||
46 => Element46,
|
||
47 => Element47,
|
||
48 => Element48,
|
||
49 => Element49,
|
||
50 => Element50,
|
||
51 => Element51,
|
||
52 => Element52,
|
||
53 => Element53,
|
||
54 => Element54,
|
||
55 => Element55,
|
||
56 => Element56,
|
||
57 => Element57,
|
||
58 => Element58,
|
||
59 => Element59,
|
||
60 => Element60,
|
||
61 => Element61,
|
||
62 => Element62,
|
||
63 => Element63,
|
||
64 => Element64,
|
||
65 => Element65,
|
||
66 => Element66,
|
||
67 => Element67,
|
||
68 => Element68,
|
||
69 => Element69,
|
||
70 => Element70,
|
||
71 => Element71,
|
||
72 => Element72,
|
||
73 => Element73,
|
||
74 => Element74,
|
||
75 => Element75,
|
||
76 => Element76,
|
||
77 => Element77,
|
||
78 => Element78,
|
||
79 => Element79,
|
||
80 => Element80,
|
||
81 => Element81,
|
||
82 => Element82,
|
||
83 => Element83,
|
||
84 => Element84,
|
||
85 => Element85,
|
||
86 => Element86,
|
||
87 => Element87,
|
||
88 => Element88,
|
||
89 => Element89,
|
||
90 => Element90,
|
||
91 => Element91,
|
||
92 => Element92,
|
||
93 => Element93,
|
||
94 => Element94,
|
||
95 => Element95,
|
||
96 => Element96,
|
||
97 => Element97,
|
||
98 => Element98,
|
||
99 => Element99,
|
||
100 => Element100,
|
||
101 => Element101,
|
||
102 => Element102,
|
||
103 => Element103,
|
||
104 => Element104,
|
||
105 => Element105,
|
||
106 => Element106,
|
||
107 => Element107,
|
||
108 => Element108,
|
||
109 => Element109,
|
||
110 => Element110,
|
||
111 => Element111,
|
||
112 => Element112,
|
||
113 => Element113,
|
||
114 => Element114,
|
||
115 => Element115,
|
||
116 => Element116,
|
||
117 => Element117,
|
||
118 => Element118,
|
||
119 => Element119,
|
||
120 => Element120,
|
||
121 => Element121,
|
||
122 => Element122,
|
||
123 => Element123,
|
||
124 => Element124,
|
||
125 => Element125,
|
||
126 => Element126,
|
||
127 => Element127,
|
||
_ => Element0,
|
||
};
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private void SetElement(int index, GpuUInt128 value)
|
||
{
|
||
switch (index)
|
||
{
|
||
case 0:
|
||
Element0 = value;
|
||
return;
|
||
case 1:
|
||
Element1 = value;
|
||
return;
|
||
case 2:
|
||
Element2 = value;
|
||
return;
|
||
case 3:
|
||
Element3 = value;
|
||
return;
|
||
case 4:
|
||
Element4 = value;
|
||
return;
|
||
case 5:
|
||
Element5 = value;
|
||
return;
|
||
case 6:
|
||
Element6 = value;
|
||
return;
|
||
case 7:
|
||
Element7 = value;
|
||
return;
|
||
case 8:
|
||
Element8 = value;
|
||
return;
|
||
case 9:
|
||
Element9 = value;
|
||
return;
|
||
case 10:
|
||
Element10 = value;
|
||
return;
|
||
case 11:
|
||
Element11 = value;
|
||
return;
|
||
case 12:
|
||
Element12 = value;
|
||
return;
|
||
case 13:
|
||
Element13 = value;
|
||
return;
|
||
case 14:
|
||
Element14 = value;
|
||
return;
|
||
case 15:
|
||
Element15 = value;
|
||
return;
|
||
case 16:
|
||
Element16 = value;
|
||
return;
|
||
case 17:
|
||
Element17 = value;
|
||
return;
|
||
case 18:
|
||
Element18 = value;
|
||
return;
|
||
case 19:
|
||
Element19 = value;
|
||
return;
|
||
case 20:
|
||
Element20 = value;
|
||
return;
|
||
case 21:
|
||
Element21 = value;
|
||
return;
|
||
case 22:
|
||
Element22 = value;
|
||
return;
|
||
case 23:
|
||
Element23 = value;
|
||
return;
|
||
case 24:
|
||
Element24 = value;
|
||
return;
|
||
case 25:
|
||
Element25 = value;
|
||
return;
|
||
case 26:
|
||
Element26 = value;
|
||
return;
|
||
case 27:
|
||
Element27 = value;
|
||
return;
|
||
case 28:
|
||
Element28 = value;
|
||
return;
|
||
case 29:
|
||
Element29 = value;
|
||
return;
|
||
case 30:
|
||
Element30 = value;
|
||
return;
|
||
case 31:
|
||
Element31 = value;
|
||
return;
|
||
case 32:
|
||
Element32 = value;
|
||
return;
|
||
case 33:
|
||
Element33 = value;
|
||
return;
|
||
case 34:
|
||
Element34 = value;
|
||
return;
|
||
case 35:
|
||
Element35 = value;
|
||
return;
|
||
case 36:
|
||
Element36 = value;
|
||
return;
|
||
case 37:
|
||
Element37 = value;
|
||
return;
|
||
case 38:
|
||
Element38 = value;
|
||
return;
|
||
case 39:
|
||
Element39 = value;
|
||
return;
|
||
case 40:
|
||
Element40 = value;
|
||
return;
|
||
case 41:
|
||
Element41 = value;
|
||
return;
|
||
case 42:
|
||
Element42 = value;
|
||
return;
|
||
case 43:
|
||
Element43 = value;
|
||
return;
|
||
case 44:
|
||
Element44 = value;
|
||
return;
|
||
case 45:
|
||
Element45 = value;
|
||
return;
|
||
case 46:
|
||
Element46 = value;
|
||
return;
|
||
case 47:
|
||
Element47 = value;
|
||
return;
|
||
case 48:
|
||
Element48 = value;
|
||
return;
|
||
case 49:
|
||
Element49 = value;
|
||
return;
|
||
case 50:
|
||
Element50 = value;
|
||
return;
|
||
case 51:
|
||
Element51 = value;
|
||
return;
|
||
case 52:
|
||
Element52 = value;
|
||
return;
|
||
case 53:
|
||
Element53 = value;
|
||
return;
|
||
case 54:
|
||
Element54 = value;
|
||
return;
|
||
case 55:
|
||
Element55 = value;
|
||
return;
|
||
case 56:
|
||
Element56 = value;
|
||
return;
|
||
case 57:
|
||
Element57 = value;
|
||
return;
|
||
case 58:
|
||
Element58 = value;
|
||
return;
|
||
case 59:
|
||
Element59 = value;
|
||
return;
|
||
case 60:
|
||
Element60 = value;
|
||
return;
|
||
case 61:
|
||
Element61 = value;
|
||
return;
|
||
case 62:
|
||
Element62 = value;
|
||
return;
|
||
case 63:
|
||
Element63 = value;
|
||
return;
|
||
case 64:
|
||
Element64 = value;
|
||
return;
|
||
case 65:
|
||
Element65 = value;
|
||
return;
|
||
case 66:
|
||
Element66 = value;
|
||
return;
|
||
case 67:
|
||
Element67 = value;
|
||
return;
|
||
case 68:
|
||
Element68 = value;
|
||
return;
|
||
case 69:
|
||
Element69 = value;
|
||
return;
|
||
case 70:
|
||
Element70 = value;
|
||
return;
|
||
case 71:
|
||
Element71 = value;
|
||
return;
|
||
case 72:
|
||
Element72 = value;
|
||
return;
|
||
case 73:
|
||
Element73 = value;
|
||
return;
|
||
case 74:
|
||
Element74 = value;
|
||
return;
|
||
case 75:
|
||
Element75 = value;
|
||
return;
|
||
case 76:
|
||
Element76 = value;
|
||
return;
|
||
case 77:
|
||
Element77 = value;
|
||
return;
|
||
case 78:
|
||
Element78 = value;
|
||
return;
|
||
case 79:
|
||
Element79 = value;
|
||
return;
|
||
case 80:
|
||
Element80 = value;
|
||
return;
|
||
case 81:
|
||
Element81 = value;
|
||
return;
|
||
case 82:
|
||
Element82 = value;
|
||
return;
|
||
case 83:
|
||
Element83 = value;
|
||
return;
|
||
case 84:
|
||
Element84 = value;
|
||
return;
|
||
case 85:
|
||
Element85 = value;
|
||
return;
|
||
case 86:
|
||
Element86 = value;
|
||
return;
|
||
case 87:
|
||
Element87 = value;
|
||
return;
|
||
case 88:
|
||
Element88 = value;
|
||
return;
|
||
case 89:
|
||
Element89 = value;
|
||
return;
|
||
case 90:
|
||
Element90 = value;
|
||
return;
|
||
case 91:
|
||
Element91 = value;
|
||
return;
|
||
case 92:
|
||
Element92 = value;
|
||
return;
|
||
case 93:
|
||
Element93 = value;
|
||
return;
|
||
case 94:
|
||
Element94 = value;
|
||
return;
|
||
case 95:
|
||
Element95 = value;
|
||
return;
|
||
case 96:
|
||
Element96 = value;
|
||
return;
|
||
case 97:
|
||
Element97 = value;
|
||
return;
|
||
case 98:
|
||
Element98 = value;
|
||
return;
|
||
case 99:
|
||
Element99 = value;
|
||
return;
|
||
case 100:
|
||
Element100 = value;
|
||
return;
|
||
case 101:
|
||
Element101 = value;
|
||
return;
|
||
case 102:
|
||
Element102 = value;
|
||
return;
|
||
case 103:
|
||
Element103 = value;
|
||
return;
|
||
case 104:
|
||
Element104 = value;
|
||
return;
|
||
case 105:
|
||
Element105 = value;
|
||
return;
|
||
case 106:
|
||
Element106 = value;
|
||
return;
|
||
case 107:
|
||
Element107 = value;
|
||
return;
|
||
case 108:
|
||
Element108 = value;
|
||
return;
|
||
case 109:
|
||
Element109 = value;
|
||
return;
|
||
case 110:
|
||
Element110 = value;
|
||
return;
|
||
case 111:
|
||
Element111 = value;
|
||
return;
|
||
case 112:
|
||
Element112 = value;
|
||
return;
|
||
case 113:
|
||
Element113 = value;
|
||
return;
|
||
case 114:
|
||
Element114 = value;
|
||
return;
|
||
case 115:
|
||
Element115 = value;
|
||
return;
|
||
case 116:
|
||
Element116 = value;
|
||
return;
|
||
case 117:
|
||
Element117 = value;
|
||
return;
|
||
case 118:
|
||
Element118 = value;
|
||
return;
|
||
case 119:
|
||
Element119 = value;
|
||
return;
|
||
case 120:
|
||
Element120 = value;
|
||
return;
|
||
case 121:
|
||
Element121 = value;
|
||
return;
|
||
case 122:
|
||
Element122 = value;
|
||
return;
|
||
case 123:
|
||
Element123 = value;
|
||
return;
|
||
case 124:
|
||
Element124 = value;
|
||
return;
|
||
case 125:
|
||
Element125 = value;
|
||
return;
|
||
case 126:
|
||
Element126 = value;
|
||
return;
|
||
case 127:
|
||
Element127 = value;
|
||
return;
|
||
default:
|
||
// ILGPU kernels cannot throw exceptions, and callers guarantee the index range.
|
||
return;
|
||
}
|
||
}
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static int GetScalarBitLength(ulong value)
|
||
{
|
||
if (value == 0UL)
|
||
{
|
||
return 0;
|
||
}
|
||
|
||
return 64 - XMath.LeadingZeroCount(value);
|
||
}
|
||
|
||
public static readonly GpuUInt128 Zero = new();
|
||
public static readonly GpuUInt128 One = new(1UL);
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public GpuUInt128()
|
||
{
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public GpuUInt128(ulong high, ulong low)
|
||
{
|
||
High = high;
|
||
Low = low;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public GpuUInt128(ulong low)
|
||
{
|
||
Low = low;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public GpuUInt128(UInt128 value)
|
||
{
|
||
High = (ulong)(value >> 64);
|
||
Low = (ulong)value;
|
||
}
|
||
|
||
public readonly bool IsZero
|
||
{
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
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);
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static implicit operator GpuUInt128(ulong value) => new(value);
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static implicit operator UInt128(GpuUInt128 value) =>
|
||
((UInt128)value.High << 64) | value.Low;
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static GpuUInt128 operator +(GpuUInt128 left, GpuUInt128 right)
|
||
{
|
||
left.Add(right);
|
||
return left;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static GpuUInt128 operator -(GpuUInt128 left, GpuUInt128 right)
|
||
{
|
||
left.Sub(right);
|
||
return left;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static GpuUInt128 operator ^(GpuUInt128 left, GpuUInt128 right) =>
|
||
new(left.High ^ right.High, left.Low ^ right.Low);
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static GpuUInt128 operator *(GpuUInt128 left, GpuUInt128 right)
|
||
{
|
||
Multiply(left, right, out var high, out var low);
|
||
left.High = high;
|
||
left.Low = low;
|
||
return left;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static GpuUInt128 operator <<(GpuUInt128 value, int shift)
|
||
{
|
||
shift &= 127;
|
||
if (shift == 0)
|
||
{
|
||
return value;
|
||
}
|
||
|
||
ulong high = value.High;
|
||
ulong low = value.Low;
|
||
|
||
if (shift >= 64)
|
||
{
|
||
int longShift = shift - 64;
|
||
value.High = low << longShift;
|
||
value.Low = 0UL;
|
||
return value;
|
||
}
|
||
|
||
int inverseShift = 64 - shift;
|
||
value.High = (high << shift) | (low >> inverseShift);
|
||
value.Low = low << shift;
|
||
return value;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static GpuUInt128 operator >>(GpuUInt128 value, int shift)
|
||
{
|
||
shift &= 127;
|
||
if (shift == 0)
|
||
{
|
||
return value;
|
||
}
|
||
|
||
ulong high = value.High;
|
||
ulong low = value.Low;
|
||
|
||
if (shift >= 64)
|
||
{
|
||
int longShift = shift - 64;
|
||
value.Low = longShift == 0 ? high : high >> longShift;
|
||
value.High = 0UL;
|
||
return value;
|
||
}
|
||
|
||
int inverseShift = 64 - shift;
|
||
value.Low = (low >> shift) | (high << inverseShift);
|
||
value.High = high >> shift;
|
||
return value;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static bool operator ==(GpuUInt128 left, GpuUInt128 right) => left.Equals(right);
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static bool operator !=(GpuUInt128 left, GpuUInt128 right) => !left.Equals(right);
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static bool operator <(GpuUInt128 left, GpuUInt128 right) => left.CompareTo(right) < 0;
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static bool operator >(GpuUInt128 left, GpuUInt128 right) => left.CompareTo(right) > 0;
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static bool operator <=(GpuUInt128 left, GpuUInt128 right) => left.CompareTo(right) <= 0;
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static bool operator >=(GpuUInt128 left, GpuUInt128 right) => left.CompareTo(right) >= 0;
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static GpuUInt128 operator ++(GpuUInt128 value)
|
||
{
|
||
value.Add(1UL);
|
||
return value;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public readonly int CompareTo(GpuUInt128 other)
|
||
{
|
||
if (High < other.High)
|
||
{
|
||
return -1;
|
||
}
|
||
|
||
if (High > other.High)
|
||
{
|
||
return 1;
|
||
}
|
||
|
||
if (Low < other.Low)
|
||
{
|
||
return -1;
|
||
}
|
||
|
||
if (Low > other.Low)
|
||
{
|
||
return 1;
|
||
}
|
||
|
||
return 0;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public readonly bool Equals(GpuUInt128 other) => High == other.High && Low == other.Low;
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public override readonly bool Equals(object? obj) => obj is GpuUInt128 value && Equals(value);
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public override readonly int GetHashCode() => HashCode.Combine(High, Low);
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void Add(GpuUInt128 other)
|
||
{
|
||
ulong originalLow = Low;
|
||
ulong low = originalLow + other.Low;
|
||
ulong carry = low < originalLow ? 1UL : 0UL;
|
||
High = High + other.High + carry;
|
||
Low = low;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void Add(ulong value)
|
||
{
|
||
ulong low = Low + value;
|
||
High += (low < Low ? 1UL : 0UL);
|
||
Low = low;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void AddMod(GpuUInt128 value, GpuUInt128 modulus)
|
||
{
|
||
if (modulus.High == 0UL)
|
||
{
|
||
AddMod(value, modulus.Low);
|
||
return;
|
||
}
|
||
|
||
Add(value);
|
||
if (CompareTo(modulus) >= 0)
|
||
{
|
||
Sub(modulus);
|
||
}
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void AddMod(ulong value, GpuUInt128 modulus)
|
||
{
|
||
if (modulus.High == 0UL)
|
||
{
|
||
AddMod(value, modulus.Low);
|
||
return;
|
||
}
|
||
|
||
Add(value);
|
||
if (CompareTo(modulus) >= 0)
|
||
{
|
||
Sub(modulus);
|
||
}
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static GpuUInt128 Pow2Minus1Mod(ulong exponent, GpuUInt128 modulus)
|
||
{
|
||
if (modulus.IsZero)
|
||
{
|
||
return Zero;
|
||
}
|
||
|
||
GpuUInt128 pow = Pow2Mod(exponent, modulus);
|
||
if (pow.IsZero)
|
||
{
|
||
pow = modulus;
|
||
}
|
||
|
||
pow.Sub(One);
|
||
return pow;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static void Pow2Minus1ModBatch(GpuUInt128 modulus, ReadOnlySpan<ulong> exponents, Span<GpuUInt128> results)
|
||
{
|
||
int count = exponents.Length;
|
||
for (int i = 0; i < count; i++)
|
||
{
|
||
results[i] = Pow2Minus1Mod(exponents[i], modulus);
|
||
}
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void AddMod(GpuUInt128 value, ulong modulus)
|
||
{
|
||
// TODO: Pre-reduce the operands via the Montgomery ladder used in MulMod64Benchmarks so the GPU
|
||
// compatible shim stops paying for `%` on every call; the InlineUInt128 helper ran 6–82× faster on
|
||
// large 64-bit workloads while preserving compatibility with the CPU scanner.
|
||
ulong a = Low % modulus;
|
||
ulong b = value.Low % modulus;
|
||
ulong sum = a + b;
|
||
if (sum >= modulus || sum < a)
|
||
{
|
||
sum -= modulus;
|
||
}
|
||
|
||
High = 0UL;
|
||
Low = sum;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void AddMod(ulong value, ulong modulus)
|
||
{
|
||
// TODO: Fold these operands with the ImmediateModulo helper once the GPU shim exposes it, avoiding
|
||
// repeated `%` reductions that the benchmarks showed are far slower than the Montgomery-based path
|
||
// for dense operands.
|
||
ulong a = Low % modulus;
|
||
ulong b = value % modulus;
|
||
ulong sum = a + b;
|
||
if (sum >= modulus || sum < a)
|
||
{
|
||
sum -= modulus;
|
||
}
|
||
|
||
High = 0UL;
|
||
Low = sum;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void Sub(GpuUInt128 other)
|
||
{
|
||
ulong borrow = Low < other.Low ? 1UL : 0UL;
|
||
High = High - other.High - borrow;
|
||
Low -= other.Low;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static ulong MulHigh(ulong x, ulong y)
|
||
{
|
||
ulong xLow = (uint)x;
|
||
ulong yLow = (uint)y;
|
||
ulong lowProduct = xLow * yLow;
|
||
|
||
ulong yHigh = y >> 32;
|
||
ulong cross = xLow * yHigh;
|
||
|
||
ulong xHigh = x >> 32;
|
||
ulong result = xHigh * yHigh;
|
||
ulong temp = xHigh * yLow;
|
||
|
||
// Keeping the wide partial sum in a dedicated local prevents the JIT from
|
||
// materialising it on the stack before the final carry propagation. The
|
||
// additional store looks redundant in C#, but it shortens the generated
|
||
// instruction sequence by avoiding an extra temporary and results in a
|
||
// measurable throughput win in the MulHigh benchmarks.
|
||
result += cross >> 32;
|
||
result += temp >> 32;
|
||
|
||
lowProduct >>= 32;
|
||
lowProduct += (uint)cross;
|
||
lowProduct += (uint)temp;
|
||
result += lowProduct >> 32;
|
||
return result;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void Mul64(GpuUInt128 other)
|
||
{
|
||
// Multiply this.Low (assumed 64-bit value) by full 128-bit other
|
||
ulong operand = Low;
|
||
Low = operand * other.Low;
|
||
ulong highProduct = operand * other.High;
|
||
High = highProduct + MulHigh(operand, other.Low);
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void MulMod(in GpuUInt128 value, in GpuUInt128 modulus)
|
||
{
|
||
GpuUInt128 multiplicand = this;
|
||
GpuUInt128 multiplier = value;
|
||
|
||
High = 0UL;
|
||
Low = 0UL;
|
||
|
||
while (!multiplier.IsZero)
|
||
{
|
||
if ((multiplier.Low & 1UL) != 0UL)
|
||
{
|
||
AddMod(multiplicand, modulus);
|
||
}
|
||
|
||
multiplicand.ShiftLeft(1);
|
||
if (multiplicand.CompareTo(modulus) >= 0)
|
||
{
|
||
multiplicand.Sub(modulus);
|
||
}
|
||
|
||
multiplier.ShiftRight(1);
|
||
}
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static GpuUInt128 Pow2Mod(ulong exponent, GpuUInt128 modulus)
|
||
{
|
||
if (modulus.IsZero || modulus == One)
|
||
{
|
||
return new();
|
||
}
|
||
|
||
if (exponent == 0UL)
|
||
{
|
||
return new GpuUInt128(1UL);
|
||
}
|
||
|
||
Pow2OddPowerTable oddPowers = default;
|
||
GpuUInt128 power = new GpuUInt128(1UL);
|
||
int oddIndex = 0;
|
||
|
||
for (int bit = 1; bit <= Pow2WindowMaxExponent; bit++)
|
||
{
|
||
power += power;
|
||
if (power.CompareTo(modulus) >= 0)
|
||
{
|
||
power.Sub(modulus);
|
||
}
|
||
|
||
if ((bit & 1) != 0)
|
||
{
|
||
oddPowers[oddIndex] = power;
|
||
oddIndex++;
|
||
}
|
||
}
|
||
|
||
GpuUInt128 result = new GpuUInt128(1UL);
|
||
int bitLength = GetScalarBitLength(exponent);
|
||
int index = bitLength - 1;
|
||
|
||
while (index >= 0)
|
||
{
|
||
if (((exponent >> index) & 1UL) == 0UL)
|
||
{
|
||
result.SquareMod(modulus);
|
||
index--;
|
||
continue;
|
||
}
|
||
|
||
int windowStart = index - Pow2WindowSizeBits + 1;
|
||
if (windowStart < 0)
|
||
{
|
||
windowStart = 0;
|
||
}
|
||
|
||
while (((exponent >> windowStart) & 1UL) == 0UL)
|
||
{
|
||
windowStart++;
|
||
}
|
||
|
||
int windowBitCount = index - windowStart + 1;
|
||
for (int square = 0; square < windowBitCount; square++)
|
||
{
|
||
result.SquareMod(modulus);
|
||
}
|
||
|
||
ulong mask = (1UL << windowBitCount) - 1UL;
|
||
ulong windowValue = (exponent >> windowStart) & mask;
|
||
int tableIndex = (int)((windowValue - 1UL) >> 1);
|
||
|
||
GpuUInt128 factor = oddPowers[tableIndex];
|
||
result.MulMod(factor, modulus);
|
||
|
||
index = windowStart - 1;
|
||
}
|
||
|
||
return result;
|
||
}
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static int TrailingZeroCount(GpuUInt128 value)
|
||
{
|
||
if (value.Low != 0UL)
|
||
{
|
||
return BitOperations.TrailingZeroCount(value.Low);
|
||
}
|
||
|
||
return 64 + BitOperations.TrailingZeroCount(value.High);
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public static GpuUInt128 BinaryGcd(GpuUInt128 u, GpuUInt128 v)
|
||
{
|
||
// TODO: Replace this scalar binary GCD with the branchless reduction from
|
||
// GpuUInt128BinaryGcdBenchmarks so CPU fallbacks stay aligned with the GPU kernel
|
||
// performance when resolving large divisor residues.
|
||
if (u.IsZero)
|
||
{
|
||
return v;
|
||
}
|
||
|
||
if (v.IsZero)
|
||
{
|
||
return u;
|
||
}
|
||
|
||
int shift = TrailingZeroCount(new GpuUInt128(u.High | v.High, u.Low | v.Low));
|
||
int zu = TrailingZeroCount(u);
|
||
u >>= zu;
|
||
|
||
do
|
||
{
|
||
int zv = TrailingZeroCount(v);
|
||
v >>= zv;
|
||
if (u > v)
|
||
{
|
||
(u, v) = (v, u);
|
||
}
|
||
|
||
v -= u;
|
||
}
|
||
while (!v.IsZero);
|
||
|
||
return u << shift;
|
||
}
|
||
|
||
|
||
public void Sub(ulong value)
|
||
{
|
||
ulong borrow = Low < value ? 1UL : 0UL;
|
||
High -= borrow;
|
||
Low -= value;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void SubMod(GpuUInt128 value, GpuUInt128 modulus)
|
||
{
|
||
if (CompareTo(value) >= 0)
|
||
{
|
||
Sub(value);
|
||
}
|
||
else
|
||
{
|
||
Add(modulus);
|
||
Sub(value);
|
||
}
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void SubMod(ulong value, GpuUInt128 modulus)
|
||
{
|
||
ulong low = Low;
|
||
ulong high = High;
|
||
if (high == 0UL && low < value)
|
||
{
|
||
low += modulus.Low;
|
||
high = modulus.High + (low < modulus.Low ? 1UL : 0UL);
|
||
}
|
||
|
||
ulong borrow = low < value ? 1UL : 0UL;
|
||
High = high - borrow;
|
||
Low = low - value;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void Mul(GpuUInt128 other)
|
||
{
|
||
Multiply(this, other, out var high, out var low);
|
||
High = high;
|
||
Low = low;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void Mul(ulong value)
|
||
{
|
||
var (highPart, lowPart) = Mul64(Low, value);
|
||
High = highPart + High * value;
|
||
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)
|
||
{
|
||
High ^= other.High;
|
||
Low ^= other.Low;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void Xor(ulong value)
|
||
{
|
||
Low ^= value;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void ShiftLeft(int shift)
|
||
{
|
||
shift &= 127;
|
||
if (shift == 0)
|
||
{
|
||
return;
|
||
}
|
||
|
||
ulong high = High;
|
||
ulong low = Low;
|
||
|
||
if (shift >= 64)
|
||
{
|
||
High = low << (shift - 64);
|
||
Low = 0UL;
|
||
return;
|
||
}
|
||
|
||
int inverseShift = 64 - shift;
|
||
High = (high << shift) | (low >> inverseShift);
|
||
Low = low << shift;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void ShiftRight(int shift)
|
||
{
|
||
shift &= 127;
|
||
if (shift == 0)
|
||
{
|
||
return;
|
||
}
|
||
|
||
ulong high = High;
|
||
ulong low = Low;
|
||
|
||
if (shift >= 64)
|
||
{
|
||
Low = shift == 64 ? high : high >> (shift - 64);
|
||
High = 0UL;
|
||
return;
|
||
}
|
||
|
||
int inverseShift = 64 - shift;
|
||
Low = (low >> shift) | (high << inverseShift);
|
||
High = high >> shift;
|
||
}
|
||
|
||
public void MulMod(GpuUInt128 other, GpuUInt128 modulus)
|
||
{
|
||
if (modulus.High == 0UL)
|
||
{
|
||
Low = MulMod(other, modulus.Low);
|
||
High = 0UL;
|
||
return;
|
||
}
|
||
|
||
MultiplyFull(this, other, out var p3, out var p2, out var p1, out var p0);
|
||
ReduceProductBitwise(p3, p2, p1, p0, modulus);
|
||
// TODO(MOD-OPT): Replace this bitwise long-division style reduction
|
||
// with a faster algorithm suitable for GPU kernels without '%' support.
|
||
// Options:
|
||
// - Montgomery reduction for 128-bit moduli (R=2^128), requires pre-
|
||
// computed modulus-dependent constants (n' and R^2 mod n). This is
|
||
// ideal for NTT primes (odd modulus). Cache per-modulus constants.
|
||
// - Barrett reduction with a 256/128 quotient approximation using only
|
||
// multiplies and shifts. Cache mu = floor(2^k / n) for k=256.
|
||
// Implement a fast path for common 64-bit NTT moduli and a separate
|
||
// path for 128-bit moduli. This will significantly reduce the cost of
|
||
// each butterfly in NTT and LL steps.
|
||
// TODO(MOD-OPT): Plumb constants through caches in NttGpuMath.SquareCacheEntry
|
||
// (e.g., Montgomery n', R2) and provide device-friendly accessors.
|
||
}
|
||
|
||
/// <summary>
|
||
/// Modular multiplication using <see cref="BigInteger"/> reduction. This
|
||
/// method is intended for validation only and should not be used inside
|
||
/// GPU kernels.
|
||
/// </summary>
|
||
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);
|
||
High = (ulong)(reduced >> 64);
|
||
Low = (ulong)reduced;
|
||
}
|
||
|
||
/// <summary>
|
||
/// Experimental limb-based reduction. The current implementation performs
|
||
/// repeated subtractions and becomes extremely slow for large remainders.
|
||
/// Kept for future optimization work.
|
||
/// </summary>
|
||
internal void MulModByLimb(GpuUInt128 other, GpuUInt128 modulus)
|
||
{
|
||
// TODO: Relocate this limb-based reducer to the benchmark project once the production
|
||
// pipeline switches to the faster allocating legacy path demonstrated in the benchmarks.
|
||
MultiplyFull(this, other, out var p3, out var p2, out var p1, out var p0);
|
||
|
||
High = p3;
|
||
Low = p2;
|
||
while (CompareTo(modulus) >= 0)
|
||
{
|
||
Sub(modulus);
|
||
}
|
||
|
||
ulong limb = p1;
|
||
for (int i = 0; i < 2; i++)
|
||
{
|
||
ShiftLeft(64);
|
||
Low = limb;
|
||
while (CompareTo(modulus) >= 0)
|
||
{
|
||
Sub(modulus);
|
||
}
|
||
|
||
limb = p0;
|
||
}
|
||
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void MulMod(ulong value, GpuUInt128 modulus)
|
||
{
|
||
if (modulus.High == 0UL)
|
||
{
|
||
Low = MulMod(value, modulus.Low);
|
||
High = 0UL;
|
||
return;
|
||
}
|
||
|
||
MultiplyFull(this, value, out var p3, out var p2, out var p1, out var p0);
|
||
ReduceProductBitwise(p3, p2, p1, p0, modulus);
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public readonly ulong MulMod(GpuUInt128 other, ulong modulus)
|
||
{
|
||
ulong modulusLocal = modulus;
|
||
ulong a = Low % modulusLocal;
|
||
ulong b = other.Low % modulusLocal;
|
||
|
||
if (a == 0UL || b == 0UL)
|
||
{
|
||
return 0UL;
|
||
}
|
||
|
||
return b <= ulong.MaxValue / a
|
||
? (a * b) % modulusLocal
|
||
: MulMod64(a, b, modulusLocal);
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public readonly ulong MulMod(ulong value, ulong modulus)
|
||
{
|
||
ulong a = Low % modulus;
|
||
ulong b = value % modulus;
|
||
|
||
if (a == 0UL || b == 0UL)
|
||
{
|
||
return 0UL;
|
||
}
|
||
|
||
ulong ulongRange = ulong.MaxValue / a;
|
||
return b <= ulongRange
|
||
? (a * b) % modulus
|
||
: MulMod64(a, b, modulus);
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public readonly ulong MulModSimplified(ulong value, ulong modulus)
|
||
{
|
||
ulong modulusLocal = modulus;
|
||
ulong a = Low % modulusLocal;
|
||
ulong b = value % modulusLocal;
|
||
|
||
if (a == 0UL || b == 0UL)
|
||
{
|
||
return 0UL;
|
||
}
|
||
|
||
return MulMod64(a, b, modulusLocal);
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public readonly ulong MulModWithNativeModulo(ulong value, ulong modulus)
|
||
{
|
||
// TODO: Drop this native-modulo path from production after migrating callers to MulMod,
|
||
// which benchmarked 4-8× faster on dense operands and still wins on mixed workloads.
|
||
ulong multiplicand = Low % modulus;
|
||
var remainder = value % modulus;
|
||
|
||
if (multiplicand == 0UL || remainder == 0UL)
|
||
{
|
||
return 0UL;
|
||
}
|
||
|
||
ulong result = 0UL;
|
||
while (true)
|
||
{
|
||
ulong chunk = remainder & NativeModuloChunkMask;
|
||
if (chunk != 0UL)
|
||
{
|
||
chunk = MultiplyChunkModulo(multiplicand, chunk, modulus);
|
||
result = (result + chunk) % modulus;
|
||
}
|
||
|
||
remainder >>= NativeModuloChunkBits;
|
||
if (remainder == 0UL)
|
||
{
|
||
break;
|
||
}
|
||
|
||
multiplicand = ShiftLeftByNativeChunk(multiplicand, modulus);
|
||
}
|
||
|
||
return result;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private void ReduceProductBitwise(ulong p3, ulong p2, ulong p1, ulong p0, GpuUInt128 modulus)
|
||
{
|
||
GpuUInt128 remainder = new();
|
||
int bit;
|
||
// TODO: Can we modify these loops to process multiple bits at a time? E.g. 64-bit chunks.
|
||
for (bit = 63; bit >= 0; bit--)
|
||
{
|
||
remainder.ShiftLeft(1);
|
||
remainder.Low |= (p3 >> bit) & 1UL;
|
||
if (remainder.CompareTo(modulus) >= 0)
|
||
{
|
||
remainder.Sub(modulus);
|
||
}
|
||
}
|
||
|
||
for (bit = 63; bit >= 0; bit--)
|
||
{
|
||
remainder.ShiftLeft(1);
|
||
remainder.Low |= (p2 >> bit) & 1UL;
|
||
if (remainder.CompareTo(modulus) >= 0)
|
||
{
|
||
remainder.Sub(modulus);
|
||
}
|
||
}
|
||
|
||
for (bit = 63; bit >= 0; bit--)
|
||
{
|
||
remainder.ShiftLeft(1);
|
||
remainder.Low |= (p1 >> bit) & 1UL;
|
||
if (remainder.CompareTo(modulus) >= 0)
|
||
{
|
||
remainder.Sub(modulus);
|
||
}
|
||
}
|
||
|
||
for (bit = 63; bit >= 0; bit--)
|
||
{
|
||
remainder.ShiftLeft(1);
|
||
remainder.Low |= (p0 >> bit) & 1UL;
|
||
if (remainder.CompareTo(modulus) >= 0)
|
||
{
|
||
remainder.Sub(modulus);
|
||
}
|
||
}
|
||
|
||
High = remainder.High;
|
||
Low = remainder.Low;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static ulong ShiftLeftByNativeChunk(ulong value, ulong modulus)
|
||
{
|
||
// TODO: Collapse this eight-step shift ladder into the ProcessEightBitWindows helper once it lands so
|
||
// we reuse the precomputed window residues instead of emitting `% modulus` after every shift.
|
||
value = (value << 1) % modulus;
|
||
value = (value << 1) % modulus;
|
||
value = (value << 1) % modulus;
|
||
value = (value << 1) % modulus;
|
||
value = (value << 1) % modulus;
|
||
value = (value << 1) % modulus;
|
||
value = (value << 1) % modulus;
|
||
value = (value << 1) % modulus;
|
||
|
||
return value;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static ulong MultiplyChunkModulo(ulong value, ulong chunk, ulong modulus)
|
||
{
|
||
ulong result = 0UL;
|
||
ulong[] nativeModuloBitMasks = NativeModuloBitMasks;
|
||
for (int bit = 0; bit < NativeModuloChunkBits; bit++)
|
||
{
|
||
if ((chunk & nativeModuloBitMasks[bit]) != 0UL)
|
||
{
|
||
result = (result + value) % modulus;
|
||
}
|
||
|
||
if (bit == NativeModuloChunkBitsMinusOne)
|
||
{
|
||
break;
|
||
}
|
||
|
||
value = (value << 1) % modulus;
|
||
}
|
||
|
||
return result;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public readonly ulong MulModMontgomery64(GpuUInt128 other, ulong modulus, ulong nPrime, ulong r2)
|
||
{
|
||
// TODO: Retire this struct-based Montgomery path from production after adopting the extension
|
||
// helper, which benchmarks 6-7× faster across dense and near-modulus operands.
|
||
ulong modulusLocal = modulus;
|
||
ulong a = Low % modulusLocal;
|
||
ulong b = other.Low % modulusLocal;
|
||
|
||
if (a == 0UL || b == 0UL)
|
||
{
|
||
return 0UL;
|
||
}
|
||
|
||
ulong aR = MontMul64(a, r2, modulusLocal, nPrime);
|
||
ulong bR = MontMul64(b, r2, modulusLocal, nPrime);
|
||
ulong cR = MontMul64(aR, bR, modulusLocal, nPrime);
|
||
return MontMul64(cR, 1UL, modulusLocal, nPrime);
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public readonly ulong MulModMontgomery64(ulong other, ulong modulus, ulong nPrime, ulong r2)
|
||
{
|
||
// TODO: Same as above—migrate callers to the scalar extension to avoid this 6-7× slowdown.
|
||
ulong modulusLocal = modulus;
|
||
ulong a = Low % modulusLocal;
|
||
ulong b = other % modulusLocal;
|
||
|
||
if (a == 0UL || b == 0UL)
|
||
{
|
||
return 0UL;
|
||
}
|
||
|
||
ulong aR = MontMul64(a, r2, modulusLocal, nPrime);
|
||
ulong bR = MontMul64(b, r2, modulusLocal, nPrime);
|
||
ulong cR = MontMul64(aR, bR, modulusLocal, nPrime);
|
||
return MontMul64(cR, 1UL, modulusLocal, nPrime);
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void SquareMod(GpuUInt128 modulus)
|
||
{
|
||
if (modulus.High == 0UL)
|
||
{
|
||
Low = MulMod(Low, modulus.Low);
|
||
High = 0UL;
|
||
return;
|
||
}
|
||
|
||
SquareFull(this, out var p3, out var p2, out var p1, out var p0);
|
||
|
||
// TODO: This should operate on the instance itself, not on a copy. Avoid creating new instances anywhere.
|
||
GpuUInt128 remainder = new();
|
||
int bit;
|
||
// TODO: Can we modify these loops to process multiple bits at a time? E.g. 64-bit chunks.
|
||
for (bit = 63; bit >= 0; bit--)
|
||
{
|
||
remainder <<= 1;
|
||
remainder = new(remainder.High, remainder.Low | ((p3 >> bit) & 1UL));
|
||
if (remainder.CompareTo(modulus) >= 0)
|
||
{
|
||
remainder.Sub(modulus);
|
||
}
|
||
}
|
||
|
||
for (bit = 63; bit >= 0; bit--)
|
||
{
|
||
remainder <<= 1;
|
||
remainder = new(remainder.High, remainder.Low | ((p2 >> bit) & 1UL));
|
||
if (remainder.CompareTo(modulus) >= 0)
|
||
{
|
||
remainder.Sub(modulus);
|
||
}
|
||
}
|
||
|
||
for (bit = 63; bit >= 0; bit--)
|
||
{
|
||
remainder <<= 1;
|
||
remainder = new(remainder.High, remainder.Low | ((p1 >> bit) & 1UL));
|
||
if (remainder.CompareTo(modulus) >= 0)
|
||
{
|
||
remainder.Sub(modulus);
|
||
}
|
||
}
|
||
|
||
for (bit = 63; bit >= 0; bit--)
|
||
{
|
||
remainder <<= 1;
|
||
remainder = new(remainder.High, remainder.Low | ((p0 >> bit) & 1UL));
|
||
if (remainder.CompareTo(modulus) >= 0)
|
||
{
|
||
remainder.Sub(modulus);
|
||
}
|
||
}
|
||
|
||
High = remainder.High;
|
||
Low = remainder.Low;
|
||
}
|
||
|
||
public void ModPow(GpuUInt128 exponent, GpuUInt128 modulus)
|
||
{
|
||
// TODO: Swap this copy-heavy path for the pooled base/exponent ladder used in
|
||
// GpuUInt128MulModBenchmarks so the GPU shim keeps reusing buffers instead of
|
||
// allocating temporary structs during residue scans.
|
||
GpuUInt128 result = new(1UL);
|
||
GpuUInt128 baseValue = this;
|
||
|
||
// TODO: Replace the single-bit square-and-multiply loop with the 64-bit windowed
|
||
// ladder measured fastest in GpuUInt128MulModBenchmarks to align with the GPU kernel
|
||
// implementation once the shared helper lands.
|
||
while (!exponent.IsZero)
|
||
{
|
||
if ((exponent.Low & 1UL) != 0UL)
|
||
{
|
||
result.MulMod(baseValue, modulus);
|
||
}
|
||
|
||
exponent >>= 1;
|
||
baseValue.MulMod(baseValue, modulus);
|
||
}
|
||
|
||
High = result.High;
|
||
Low = result.Low;
|
||
}
|
||
|
||
public void ModPow(ulong exponent, GpuUInt128 modulus)
|
||
{
|
||
// TODO: Share the pooled ladder state from GpuUInt128MulModBenchmarks here as well so
|
||
// Lucas–Lehmer batches avoid constructing throwaway temporaries on every exponent.
|
||
GpuUInt128 result = new(1UL);
|
||
GpuUInt128 baseValue = this;
|
||
|
||
// TODO: Upgrade this loop to the same 64-bit windowed ladder proven fastest in
|
||
// GpuUInt128MulModBenchmarks so the scalar helper matches the GPU-optimized path.
|
||
while (exponent != 0UL)
|
||
{
|
||
if ((exponent & 1UL) != 0UL)
|
||
{
|
||
result.MulMod(baseValue, modulus);
|
||
}
|
||
|
||
exponent >>= 1;
|
||
baseValue.MulMod(baseValue, modulus);
|
||
}
|
||
|
||
High = result.High;
|
||
Low = result.Low;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void ModInv(ulong modulus)
|
||
{
|
||
// TODO: Replace this Fermat inversion with the Montgomery ladder highlighted in
|
||
// GpuUInt128Montgomery64Benchmarks so we avoid instantiating a temporary modulus and
|
||
// reuse the pooled reduction constants.
|
||
ModPow(modulus - 2UL, new GpuUInt128(modulus));
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public void ModInv(GpuUInt128 modulus)
|
||
{
|
||
if (modulus.High == 0UL)
|
||
{
|
||
ModInv(modulus.Low);
|
||
return;
|
||
}
|
||
|
||
var exponent = new GpuUInt128(modulus);
|
||
exponent.Sub(2UL);
|
||
ModPow(exponent, modulus);
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
public override readonly string ToString() =>
|
||
((UInt128)this).ToString();
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static void MultiplyFull(GpuUInt128 left, GpuUInt128 right, out ulong p3, out ulong p2, out ulong p1, out ulong p0)
|
||
{
|
||
// TODO: Reuse variables to reduce register pressure following the fused-limb layout
|
||
// from GpuUInt128MulModByLimbBenchmarks so the multiply helper matches the fastest
|
||
// GPU-compatible scalar routine.
|
||
var (h0, l0) = Mul64(left.Low, right.Low);
|
||
var (h1, l1) = Mul64(left.Low, right.High);
|
||
var (h2, l2) = Mul64(left.High, right.Low);
|
||
var (h3, l3) = Mul64(left.High, right.High);
|
||
|
||
p0 = l0;
|
||
ulong carry = 0UL;
|
||
ulong sum = h0;
|
||
sum += l1;
|
||
if (sum < l1)
|
||
{
|
||
carry++;
|
||
}
|
||
|
||
sum += l2;
|
||
if (sum < l2)
|
||
{
|
||
carry++;
|
||
}
|
||
|
||
p1 = sum;
|
||
sum = h1;
|
||
sum += h2;
|
||
ulong carry2 = sum < h2 ? 1UL : 0UL;
|
||
sum += l3;
|
||
if (sum < l3)
|
||
{
|
||
carry2++;
|
||
}
|
||
|
||
sum += carry;
|
||
if (sum < carry)
|
||
{
|
||
carry2++;
|
||
}
|
||
|
||
p2 = sum;
|
||
p3 = h3 + carry2;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static void MultiplyFull(GpuUInt128 left, ulong right, out ulong p3, out ulong p2, out ulong p1, out ulong p0)
|
||
{
|
||
var (highLow, lowLow) = Mul64(left.Low, right);
|
||
var (highHigh, lowHigh) = Mul64(left.High, right);
|
||
|
||
p0 = lowLow;
|
||
|
||
ulong mid = highLow + lowHigh;
|
||
ulong carry = mid < highLow ? 1UL : 0UL;
|
||
p1 = mid;
|
||
|
||
ulong upper = highHigh + carry;
|
||
p3 = upper < highHigh ? 1UL : 0UL;
|
||
p2 = upper;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static void Multiply(GpuUInt128 left, GpuUInt128 right, out ulong high, out ulong low)
|
||
{
|
||
// TODO: Reuse variables to reduce register pressure.
|
||
var (h0, l0) = Mul64(left.Low, right.Low);
|
||
low = l0;
|
||
// TODO: Since we ignore the first result element, can we create a version of the function which calculates and returns only the second element?
|
||
(_, l0) = Mul64(left.Low, right.High);
|
||
var (_, l2) = Mul64(left.High, right.Low);
|
||
|
||
// TODO: Why not just modify left instance directly instead using out parameters?
|
||
high = h0 + l0 + l2;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static (ulong High, ulong Low) Mul64(ulong left, ulong right)
|
||
{
|
||
// TODO: Reuse variables to reduce register pressure.
|
||
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);
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static void SquareFull(GpuUInt128 value, out ulong p3, out ulong p2, out ulong p1, out ulong p0)
|
||
{
|
||
// TODO: Reuse variables to reduce register pressure.
|
||
// Compute (H*2^64 + L)^2 = H^2*2^128 + 2*H*L*2^64 + L^2
|
||
ulong L = value.Low;
|
||
ulong H = value.High;
|
||
var (hLL, lLL) = Mul64(L, L); // L^2 -> (hLL,lLL)
|
||
var (hHH, lHH) = Mul64(H, H); // H^2 -> (hHH,lHH)
|
||
var (hLH, lLH) = Mul64(L, H); // L*H -> (hLH,lLH)
|
||
|
||
// double LH: (hLH,lLH) << 1
|
||
ulong dLH_low = lLH << 1;
|
||
ulong carry = (lLH >> 63) & 1UL;
|
||
ulong dLH_high = (hLH << 1) | carry;
|
||
|
||
// Assemble 256-bit result parts
|
||
p0 = lLL;
|
||
|
||
// p1 = hLL + dLH_low (with carry)
|
||
ulong sum1 = hLL + dLH_low;
|
||
ulong c1 = sum1 < hLL ? 1UL : 0UL;
|
||
p1 = sum1;
|
||
|
||
// p2 = lHH + dLH_high + c1
|
||
ulong sum2 = lHH + dLH_high;
|
||
ulong c2 = sum2 < lHH ? 1UL : 0UL;
|
||
sum2 += c1;
|
||
if (sum2 < c1)
|
||
{
|
||
c2++;
|
||
}
|
||
p2 = sum2;
|
||
|
||
// p3 = hHH + c2
|
||
p3 = hHH + c2;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static ulong MulMod64(ulong a, ulong b, ulong modulus)
|
||
{
|
||
ulong aReduced = a % modulus;
|
||
ulong bReduced = b % modulus;
|
||
if (aReduced == 0UL || bReduced == 0UL)
|
||
{
|
||
return 0UL;
|
||
}
|
||
|
||
// TODO: Can we identify where we deal with Mersenne exponent in EvenPerfectBitScanner and directly use MulModMersenne there, removing this branch?
|
||
if (TryGetMersenneExponent(modulus, out int exponent))
|
||
{
|
||
return MulModMersenne(aReduced, bReduced, modulus, exponent);
|
||
}
|
||
|
||
// Fallback for general moduli using the branch-free shift-add reducer.
|
||
ulong result = 0UL;
|
||
ulong x = aReduced;
|
||
ulong y = bReduced;
|
||
// TODO: Can we modify this loop to process multiple bits at a time? E.g. 64-bit chunks.
|
||
while (y != 0UL)
|
||
{
|
||
if ((y & 1UL) != 0UL)
|
||
{
|
||
result += x;
|
||
if (result >= modulus)
|
||
{
|
||
result -= modulus;
|
||
}
|
||
}
|
||
|
||
x <<= 1;
|
||
if (x >= modulus)
|
||
{
|
||
x -= modulus;
|
||
}
|
||
|
||
y >>= 1;
|
||
}
|
||
|
||
return result;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static ulong MulModMersenne(ulong a, ulong b, ulong modulus, int exponent)
|
||
{
|
||
// TODO: Reuse variables to reduce register pressure.
|
||
Mul64Parts(a, b, out ulong productHigh, out ulong productLow);
|
||
ulong mask = exponent == 64 ? ulong.MaxValue : modulus;
|
||
|
||
ulong currentHigh = productHigh;
|
||
ulong currentLow = productLow;
|
||
|
||
do
|
||
{
|
||
ulong lower = currentLow & mask;
|
||
ShiftRight128(currentHigh, currentLow, exponent, out ulong shiftedHigh, out ulong shiftedLow);
|
||
currentLow = lower + shiftedLow;
|
||
currentHigh = shiftedHigh;
|
||
if (currentLow < lower)
|
||
{
|
||
currentHigh++;
|
||
}
|
||
}
|
||
while (currentHigh != 0UL);
|
||
|
||
ulong result = currentLow;
|
||
// TODO: Can we modify this loop to process multiple bits at a time? E.g. 64-bit chunks.
|
||
while (result >= modulus)
|
||
{
|
||
result -= modulus;
|
||
}
|
||
|
||
return result;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static bool TryGetMersenneExponent(ulong modulus, out int exponent)
|
||
{
|
||
if (modulus <= 1UL)
|
||
{
|
||
exponent = 0;
|
||
return false;
|
||
}
|
||
|
||
modulus++;
|
||
if (modulus == 0UL)
|
||
{
|
||
exponent = 64;
|
||
return true;
|
||
}
|
||
|
||
if ((modulus & (modulus - 1UL)) != 0UL)
|
||
{
|
||
exponent = 0;
|
||
return false;
|
||
}
|
||
|
||
int bits = 0;
|
||
// TODO: Can we modify this loop to process multiple bits at a time? E.g. 64-bit chunks.
|
||
while ((modulus & 1UL) == 0UL)
|
||
{
|
||
modulus >>= 1;
|
||
bits++;
|
||
}
|
||
|
||
exponent = bits;
|
||
return true;
|
||
}
|
||
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static void ShiftRight128(ulong high, ulong low, int shift, out ulong newHigh, out ulong newLow)
|
||
{
|
||
if (shift == 0)
|
||
{
|
||
newHigh = high;
|
||
newLow = low;
|
||
return;
|
||
}
|
||
|
||
if (shift < 64)
|
||
{
|
||
newLow = (low >> shift) | (high << (64 - shift));
|
||
newHigh = high >> shift;
|
||
return;
|
||
}
|
||
|
||
if (shift == 64)
|
||
{
|
||
newLow = high;
|
||
newHigh = 0UL;
|
||
return;
|
||
}
|
||
|
||
int extra = shift - 64;
|
||
newLow = high >> extra;
|
||
newHigh = 0UL;
|
||
}
|
||
|
||
// Montgomery core for 64-bit operands. Returns a*b*R^{-1} mod modulus.
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static ulong MontMul64(ulong aR, ulong bR, ulong modulus, ulong nPrime)
|
||
{
|
||
// TODO: Reuse variables to reduce register pressure.
|
||
// t = aR * bR (128-bit)
|
||
Mul64Parts(aR, bR, out ulong tHigh, out ulong tLow);
|
||
// m = (tLow * nPrime) mod 2^64 (low 64 bits only)
|
||
Mul64Parts(tLow, nPrime, out _, out ulong mLow);
|
||
// u = (t + m * modulus) >> 64
|
||
Mul64Parts(mLow, modulus, out ulong mmHigh, out ulong mmLow);
|
||
ulong carry = 0UL;
|
||
ulong low = tLow + mmLow;
|
||
if (low < tLow)
|
||
{
|
||
carry = 1UL;
|
||
}
|
||
|
||
ulong high = tHigh + mmHigh + carry;
|
||
ulong u = high; // (t + m*n) >> 64
|
||
if (u >= modulus)
|
||
{
|
||
u -= modulus;
|
||
}
|
||
|
||
return u;
|
||
}
|
||
|
||
// 64x64 -> 128 multiply into (high, low)
|
||
[MethodImpl(MethodImplOptions.AggressiveInlining)]
|
||
private static void Mul64Parts(ulong a, ulong b, out ulong high, out ulong low)
|
||
{
|
||
// TODO: Reuse variables to reduce register pressure.
|
||
ulong a0 = (uint)a;
|
||
ulong a1 = a >> 32;
|
||
ulong b0 = (uint)b;
|
||
ulong b1 = b >> 32;
|
||
|
||
ulong lo = a0 * b0;
|
||
ulong mid1 = a1 * b0;
|
||
ulong mid2 = a0 * b1;
|
||
ulong hi = a1 * b1;
|
||
|
||
ulong carry = (lo >> 32) + (uint)mid1 + (uint)mid2;
|
||
low = (lo & 0xFFFFFFFFUL) | (carry << 32);
|
||
hi += (mid1 >> 32) + (mid2 >> 32) + (carry >> 32);
|
||
high = hi;
|
||
}
|
||
|
||
private static ulong[] CreateNativeModuloBitMasks()
|
||
{
|
||
ulong[] masks = new ulong[NativeModuloBitMaskTableSize];
|
||
for (int bit = 0; bit < NativeModuloBitMaskTableSize; bit++)
|
||
{
|
||
masks[bit] = bit < 64 ? 1UL << bit : 0UL;
|
||
}
|
||
|
||
return masks;
|
||
}
|
||
}
|
||
|
||
// TODO: Check if the TODO below is still relevant.
|
||
|
||
// TODO(MOD-OPT): Montgomery/Barrett integration plan
|
||
// - Introduce caches of modulus-dependent constants:
|
||
// * Montgomery: n' (-(n^{-1}) mod 2^64 or 2^128), R2 = (R^2 mod n)
|
||
// * Barrett: mu = floor(2^k / n) for k ∈ {128, 192, 256}
|
||
// - Add fast-path for 64-bit NTT primes (modulus.High == 0UL) using pure 64-bit ops.
|
||
// - Expose helpers to retrieve/calc constants once per modulus and reuse in kernels.
|
||
// - Wire these into MulMod and SquareMod hot paths under feature toggles.
|
||
// - Ensure ILGPU compatibility (no BigInteger, no % inside kernels).
|
||
|