Extension methods benchmarked. Slower alternatives replaced.

This commit is contained in:
Marek Lemańczyk committed 2025-09-23 22:01:09 +02:00
1 parent 9becb8df67
commit 6e99ee39e3
40 files changed
+929 -1025

No files matched your search

@@ -1,5 +1,4 @@
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks; namespace EvenPerfectBitScanner.Benchmarks;
@@ -15,38 +14,45 @@ public class Mod10UInt128Benchmarks
return Value % 10UL; return Value % 10UL;
} }
/// <summary>
/// Fastest
/// </summary>
[Benchmark] [Benchmark]
public ulong ModWithHighLowModulo() public ulong ModWithHighLowModulo()
{ {
// Split and fold under mod 10: 2^64 ≡ 6 (mod 10) return ModWithHighLowModuloHelper(Value);
ulong value = (ulong)Value;
ulong high = value >> 64;
ulong low = value;
ulong highRem = high % 10UL;
ulong lowRem = low % 10UL;
ulong combined = lowRem + highRem * 6UL;
return combined % 10UL;
} }
[Benchmark] private static ulong ModWithHighLowModuloHelper(UInt128 value128)
public ulong ModWithLoopModulo() {
{ // Split and fold under mod 10: 2^64 ≡ 6 (mod 10)
UInt128 value = Value; ulong value = (ulong)value128;
return ((value % 10UL) + ((value >> 64) % 10UL) * 6UL) % 10UL;
}
[Benchmark]
public ulong ModWithLoopModulo()
{
return ModWithLoopModuloHelper(Value);
}
private static ulong ModWithLoopModuloHelper(UInt128 value)
{
UInt128 zero = UInt128.Zero; UInt128 zero = UInt128.Zero;
if (value == zero) if (value == zero)
return 0UL; return 0UL;
ulong result = (ulong)value; ulong result = (ulong)value;
UInt128 high = value >> 64; value >>= 64;
while (high != zero) while (value != zero)
{ {
// 2^64 ≡ 6 (mod 10) // 2^64 ≡ 6 (mod 10)
result = (result + (ulong)high * 6UL) % 10UL; result = (result + (ulong)value * 6UL) % 10UL;
high >>= 64; value >>= 64;
} }
return result % 10UL; return result % 10UL;
} }
} }
@@ -1,3 +1,4 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -8,10 +9,16 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod10UIntBenchmarks public class Mod10UIntBenchmarks
{ {
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod10(uint value) => value - (uint)((value * UIntExtensions.Mod5Mask) >> 35) * 10U;
[Params(10U, 2047U, 65535U, 2147483647U)] [Params(10U, 2047U, 65535U, 2147483647U)]
public uint Value { get; set; } public uint Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public uint ModuloOperator() public uint ModuloOperator()
{ {
return Value % 10U; return Value % 10U;
@@ -20,7 +27,7 @@ public class Mod10UIntBenchmarks
[Benchmark] [Benchmark]
public uint ExtensionMethod() public uint ExtensionMethod()
{ {
return Value.Mod10(); return Mod10(Value);
} }
} }
@@ -1,4 +1,3 @@
using System.Numerics;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -16,7 +15,10 @@ public class Mod10ULongBenchmarks
return Value % 10UL; return Value % 10UL;
} }
[Benchmark] /// <summary>
/// Fastest
/// </summary>
[Benchmark]
public ulong ModMethodModulo() public ulong ModMethodModulo()
{ {
return Value.Mod10(); return Value.Mod10();
@@ -1,4 +1,3 @@
using System.Numerics;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -10,7 +9,10 @@ public class Mod10_8_5_3Benchmarks
[Params(3UL, 8191UL, 131071UL, 2147483647UL)] [Params(3UL, 8191UL, 131071UL, 2147483647UL)]
public ulong Value { get; set; } public ulong Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public (ulong Mod10, ulong Mod8, ulong Mod3, ulong Mod5) ModuloOperator() public (ulong Mod10, ulong Mod8, ulong Mod3, ulong Mod5) ModuloOperator()
{ {
ulong value = Value; ulong value = Value;
@@ -1,4 +1,3 @@
using System.Numerics;
using System.Runtime.CompilerServices; using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -36,10 +35,13 @@ public class Mod10_8_5_3PropertiesBenchmarks
return rem >= 5UL ? rem - 5UL : rem; return rem >= 5UL ? rem - 5UL : rem;
} }
[Params(3UL, 8191UL, 131071UL, 2147483647UL)] [Params(3UL, 6UL, 128UL, 256UL, 2_047UL, 8_191UL, 65_535UL, 131_071UL, 2_147_483_647UL, ulong.MaxValue - 1_024UL, ulong.MaxValue - 511UL, ulong.MaxValue - 255UL, ulong.MaxValue - 127UL, ulong.MaxValue - 63UL, ulong.MaxValue - 31UL, ulong.MaxValue - 17UL)]
public ulong Value { get; set; } public ulong Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public (ulong Mod10, ulong Mod8, ulong Mod5, ulong Mod3) LegacyModulo() public (ulong Mod10, ulong Mod8, ulong Mod5, ulong Mod3) LegacyModulo()
{ {
UInt128 value = Value; UInt128 value = Value;
@@ -53,12 +55,7 @@ public class Mod10_8_5_3PropertiesBenchmarks
[Benchmark] [Benchmark]
public (ulong Mod10, ulong Mod8, ulong Mod3, ulong Mod5) ModMethodModulo() public (ulong Mod10, ulong Mod8, ulong Mod3, ulong Mod5) ModMethodModulo()
{ {
Value.Mod10_8_5_3(out var mod10, out ulong mod8, out ulong mod5, out ulong mod3); Value.Mod10_8_5_3(out Mod10R, out Mod8R, out Mod5R, out Mod3R);
Mod10R = mod10;
Mod8R = mod8;
Mod5R = mod5;
Mod3R = mod3;
return (Mod10R, Mod8R, Mod3R, Mod5R); return (Mod10R, Mod8R, Mod3R, Mod5R);
} }
} }
@@ -1,4 +1,3 @@
using System.Numerics;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -14,19 +13,21 @@ public class Mod10_8_5_3StepsBenchmarks
public (ulong Mod10, ulong Mod8, ulong Mod3, ulong Mod5) LegacyModulo() public (ulong Mod10, ulong Mod8, ulong Mod3, ulong Mod5) LegacyModulo()
{ {
UInt128 value = Value; UInt128 value = Value;
ulong mod10 = (ulong)(value % 10UL) % 10UL; ulong mod10 = (((ulong)(value % 10UL)) << 1) % 10UL;
ulong mod8 = (ulong)(value % 8UL); ulong mod8 = (((ulong)(value & 7UL)) << 1) & 7UL;
ulong mod3 = (ulong)(value % 3UL); ulong mod5 = (((ulong)(value % 5UL)) << 1) % 5UL;
ulong mod5 = (ulong)(value % 5UL); ulong mod3 = (((ulong)(value % 3UL)) << 1) % 3UL;
return (mod10, mod8, mod3, mod5); return (mod10, mod8, mod3, mod5);
} }
[Benchmark] /// <summary>
/// Fastest
/// </summary>
[Benchmark]
public (ulong Mod10, ulong Mod8, ulong Mod3, ulong Mod5) ModMethodModulo() public (ulong Mod10, ulong Mod8, ulong Mod3, ulong Mod5) ModMethodModulo()
{ {
Value.Mod10_8_5_3(out ulong mod10, out ulong mod8, out ulong mod5, out ulong mod3); Value.Mod10_8_5_3Steps(out ulong mod10, out ulong mod8, out ulong mod5, out ulong mod3);
// Step residues for q += 2*p // Step residues for q += 2*p
mod10 %= 10UL;
return (mod10, mod8, mod3, mod5); return (mod10, mod8, mod3, mod5);
} }
} }
@@ -1,6 +1,6 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks; namespace EvenPerfectBitScanner.Benchmarks;
@@ -8,10 +8,27 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod11UIntBenchmarks public class Mod11UIntBenchmarks
{ {
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod11(uint value)
{
uint remainder = 0U;
while (value != 0U)
{
remainder += value & 1023U;
value >>= 10;
remainder -= 11U * (remainder / 11U);
}
return remainder;
}
[Params(11U, 2047U, 65535U, 2147483647U)] [Params(11U, 2047U, 65535U, 2147483647U)]
public uint Value { get; set; } public uint Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public uint ModuloOperator() public uint ModuloOperator()
{ {
return Value % 11U; return Value % 11U;
@@ -20,7 +37,7 @@ public class Mod11UIntBenchmarks
[Benchmark] [Benchmark]
public uint ExtensionMethod() public uint ExtensionMethod()
{ {
return Value.Mod11(); return Mod11(Value);
} }
} }
@@ -1,6 +1,6 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks; namespace EvenPerfectBitScanner.Benchmarks;
@@ -8,10 +8,25 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod11ULongBenchmarks public class Mod11ULongBenchmarks
{ {
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod11(ulong value)
{
ulong remainder = ((uint)value % 11U) + ((uint)(value >> 32) % 11U) << 2;
while (remainder >= 11UL)
{
remainder -= 11UL;
}
return remainder;
}
[Params(11UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 127UL)] [Params(11UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 127UL)]
public ulong Value { get; set; } public ulong Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public ulong ModuloOperator() public ulong ModuloOperator()
{ {
return Value % 11UL; return Value % 11UL;
@@ -20,7 +35,7 @@ public class Mod11ULongBenchmarks
[Benchmark] [Benchmark]
public ulong ExtensionMethod() public ulong ExtensionMethod()
{ {
return Value.Mod11(); return Mod11(Value);
} }
} }
@@ -1,32 +1,92 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks; namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod11_7_5_3Benchmarks public class Mod11_7_5_3Benchmarks
{ {
private static readonly byte[] Mod7ByteCoefficients = [1, 4, 2, 1, 4, 2, 1, 4];
private static readonly byte[] Mod11ByteCoefficients = [1, 3, 9, 5, 4, 1, 3, 9];
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong ReduceMod7(uint value)
{
while (value >= 7U)
{
value = (value >> 3) + (value & 7U);
if (value >= 7U && value < 14U)
{
value -= 7U;
}
}
return value;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong ReduceMod11(uint value)
{
while (value >= 11U)
{
value = (value & 15U) + ((value >> 4) * 5U);
if (value >= 11U && value < 22U)
{
value -= 11U;
}
}
return value;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void Mod11_7_5_3(ulong value, out ulong mod11, out ulong mod7, out ulong mod5, out ulong mod3)
{
uint mod3And5Accumulator = 0U;
uint mod7Accumulator = 0U;
uint mod11Accumulator = 0U;
int index = 0;
byte current;
ulong temp = value;
byte[] mod7ByteCoefficients = Mod7ByteCoefficients;
byte[] mod11ByteCoefficients = Mod11ByteCoefficients;
do
{
current = (byte)temp;
mod3And5Accumulator += current;
mod7Accumulator += (uint)(current * mod7ByteCoefficients[index]);
mod11Accumulator += (uint)(current * mod11ByteCoefficients[index]);
temp >>= 8;
index++;
}
while (temp != 0UL);
mod3 = mod3And5Accumulator % 3;
mod5 = mod3And5Accumulator % 5;
mod7 = (uint)ReduceMod7(mod7Accumulator);
mod11 = (uint)ReduceMod11(mod11Accumulator);
}
[Params(3UL, 8191UL, 131071UL, 2147483647UL)] [Params(3UL, 8191UL, 131071UL, 2147483647UL)]
public ulong Value { get; set; } public ulong Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public (ulong Mod11, ulong Mod7, ulong Mod5, ulong Mod3) ModuloOperator() public (ulong Mod11, ulong Mod7, ulong Mod5, ulong Mod3) ModuloOperator()
{ {
ulong value = Value; ulong value = Value;
return (value % 11UL, value % 7UL, value % 5UL, value % 3UL); return (value % 11UL, value % 7UL, value % 5UL, value % 3UL);
} }
[Benchmark]
public (ulong Mod11, ulong Mod7, ulong Mod5, ulong Mod3) ExtensionMethods()
{
ulong value = Value;
return (value.Mod11(), value.Mod7(), value.Mod5(), value.Mod3());
}
[Benchmark] [Benchmark]
public (ulong Mod11, ulong Mod7, ulong Mod5, ulong Mod3) CombinedMethod() public (ulong Mod11, ulong Mod7, ulong Mod5, ulong Mod3) CombinedMethod()
{ {
Value.Mod11_7_5_3(out ulong mod11, out ulong mod7, out ulong mod5, out ulong mod3); Mod11_7_5_3(Value, out ulong mod11, out ulong mod7, out ulong mod5, out ulong mod3);
return (mod11, mod7, mod5, mod3); return (mod11, mod7, mod5, mod3);
} }
} }
@@ -1,4 +1,3 @@
using System.Numerics;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -12,12 +11,16 @@ public class Mod128UInt128Benchmarks
[Params(128UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 511UL)] [Params(128UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 511UL)]
public UInt128 Value { get; set; } public UInt128 Value { get; set; }
[Benchmark(Baseline = true)] [Benchmark(Baseline = true)]
public ulong ModuloOperator() public ulong ModuloOperator()
{ {
return (ulong)(Value % 128UL); return (ulong)(Value % 128UL);
} }
/// <summary>
/// Fastest
/// </summary>
/// <returns></returns>
[Benchmark] [Benchmark]
public ulong ExtensionMethod() public ulong ExtensionMethod()
{ {
@@ -1,6 +1,6 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks; namespace EvenPerfectBitScanner.Benchmarks;
@@ -8,10 +8,17 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod128UIntBenchmarks public class Mod128UIntBenchmarks
{ {
[Params(128U, 2047U, 65535U, 2147483647U)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod128(uint value) => value & 127U;
[Params(3UL, 5UL, 6UL, 8UL, 10UL, 11UL, 32UL, 64UL, 128UL, 256UL, 2_047UL, 8_191UL, 65_535UL, 131_071UL, 2_147_483_647UL, uint.MaxValue - 1_024UL, uint.MaxValue - 511UL, uint.MaxValue - 255UL, uint.MaxValue - 127UL, uint.MaxValue - 63UL, uint.MaxValue - 31UL, uint.MaxValue - 17UL, uint.MaxValue)]
public uint Value { get; set; } public uint Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest most of the times. It's usually much faster than it's slower.
/// </summary>
/// <returns></returns>
[Benchmark(Baseline = true)]
public uint ModuloOperator() public uint ModuloOperator()
{ {
return Value % 128U; return Value % 128U;
@@ -20,7 +27,7 @@ public class Mod128UIntBenchmarks
[Benchmark] [Benchmark]
public uint ExtensionMethod() public uint ExtensionMethod()
{ {
return Value.Mod128(); return Mod128(Value);
} }
} }
@@ -8,15 +8,19 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod128ULongBenchmarks public class Mod128ULongBenchmarks
{ {
[Params(128UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 255UL)] [Params(3UL, 5UL, 6UL, 8UL, 10UL, 11UL, 32UL, 64UL, 128UL, 256UL, 2_047UL, 8_191UL, 65_535UL, 131_071UL, 2_147_483_647UL, uint.MaxValue - 1_024UL, uint.MaxValue - 511UL, uint.MaxValue - 255UL, uint.MaxValue - 127UL, uint.MaxValue - 63UL, uint.MaxValue - 31UL, uint.MaxValue - 17UL, uint.MaxValue)]
public ulong Value { get; set; } public ulong Value { get; set; }
[Benchmark(Baseline = true)] [Benchmark(Baseline = true)]
public ulong ModuloOperator() public ulong ModuloOperator()
{ {
return Value % 128UL; return Value % 128UL;
} }
/// <summary>
/// Usually fastest for high values. It's slower for some values and faster for other.
/// </summary>
/// <returns></returns>
[Benchmark] [Benchmark]
public ulong ExtensionMethod() public ulong ExtensionMethod()
{ {
@@ -20,7 +20,7 @@ public class Mod3UIntBenchmarks
[Benchmark] [Benchmark]
public uint ExtensionMethod() public uint ExtensionMethod()
{ {
return Value.Mod3(); return UIntExtensions.Mod3(Value);
} }
} }
@@ -1,3 +1,4 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -8,10 +9,22 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod3ULongBenchmarks public class Mod3ULongBenchmarks
{ {
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong Mod3(ulong value) =>
(
(uint)(value & ULongExtensions.WordBitMask) +
(uint)((value >> 16) & ULongExtensions.WordBitMask) +
(uint)((value >> 32) & ULongExtensions.WordBitMask) +
(uint)((value >> 48) & ULongExtensions.WordBitMask)
) % 3;
[Params(3UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 17UL)] [Params(3UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 17UL)]
public ulong Value { get; set; } public ulong Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public ulong ModuloOperator() public ulong ModuloOperator()
{ {
return Value % 3UL; return Value % 3UL;
@@ -20,7 +33,7 @@ public class Mod3ULongBenchmarks
[Benchmark] [Benchmark]
public ulong ExtensionMethod() public ulong ExtensionMethod()
{ {
return Value.Mod3(); return Mod3(Value);
} }
} }
@@ -8,7 +8,7 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod5UIntBenchmarks public class Mod5UIntBenchmarks
{ {
[Params(5U, 2047U, 65535U, 2147483647U)] [Params(3UL, 5UL, 6UL, 8UL, 10UL, 11UL, 32UL, 64UL, 128UL, 256UL, 2_047UL, 8_191UL, 65_535UL, 131_071UL, 2_147_483_647UL, uint.MaxValue - 1_024UL, uint.MaxValue - 511UL, uint.MaxValue - 255UL, uint.MaxValue - 127UL, uint.MaxValue - 63UL, uint.MaxValue - 31UL, uint.MaxValue - 17UL, uint.MaxValue)]
public uint Value { get; set; } public uint Value { get; set; }
[Benchmark(Baseline = true)] [Benchmark(Baseline = true)]
@@ -17,7 +17,10 @@ public class Mod5UIntBenchmarks
return Value % 5U; return Value % 5U;
} }
[Benchmark] /// <summary>
/// Mixed results. Sometimes it's faster, in higher ranges, sometimes it's slower.
/// </summary>
[Benchmark]
public uint ExtensionMethod() public uint ExtensionMethod()
{ {
return Value.Mod5(); return Value.Mod5();
@@ -1,3 +1,4 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -8,10 +9,22 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod5ULongBenchmarks public class Mod5ULongBenchmarks
{ {
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod5(ulong value) =>
(
(uint)(value & ULongExtensions.WordBitMask) +
(uint)((value >> 16) & ULongExtensions.WordBitMask) +
(uint)((value >> 32) & ULongExtensions.WordBitMask) +
(uint)((value >> 48) & ULongExtensions.WordBitMask)
) % 5;
[Params(5UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 31UL)] [Params(5UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 31UL)]
public ulong Value { get; set; } public ulong Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public ulong ModuloOperator() public ulong ModuloOperator()
{ {
return Value % 5UL; return Value % 5UL;
@@ -20,7 +33,7 @@ public class Mod5ULongBenchmarks
[Benchmark] [Benchmark]
public ulong ExtensionMethod() public ulong ExtensionMethod()
{ {
return Value.Mod5(); return Mod5(Value);
} }
} }
@@ -1,4 +1,4 @@
using System.Numerics; using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -9,6 +9,11 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod6ComparisonBenchmarks public class Mod6ComparisonBenchmarks
{ {
private static readonly byte[] Mod6Lookup = [0, 3, 4, 1, 2, 5];
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod6(ulong value) => Mod6Lookup[(int)(((value % 3UL) << 1) | (value & 1UL))];
private const int Iterations = 10_000_000; private const int Iterations = 10_000_000;
private const ulong Divisor = 6UL; private const ulong Divisor = 6UL;
private static readonly ulong FastDivMul = (ulong)(((UInt128)1 << 64) / Divisor); private static readonly ulong FastDivMul = (ulong)(((UInt128)1 << 64) / Divisor);
@@ -71,7 +76,7 @@ public class Mod6ComparisonBenchmarks
for (int i = 0; i < Iterations; i++) for (int i = 0; i < Iterations; i++)
{ {
checksum += value.Mod6(); checksum += Mod6(value);
value += (ulong)((i & 15) + 1); value += (ulong)((i & 15) + 1);
} }
@@ -1,6 +1,6 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks; namespace EvenPerfectBitScanner.Benchmarks;
@@ -8,10 +8,18 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod6UIntBenchmarks public class Mod6UIntBenchmarks
{ {
private static readonly byte[] Mod6Lookup = [0, 3, 4, 1, 2, 5];
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod6(uint value) => Mod6Lookup[(int)(((value % 3U) << 1) | (value & 1U))];
[Params(6U, 2047U, 65535U, 2147483647U)] [Params(6U, 2047U, 65535U, 2147483647U)]
public uint Value { get; set; } public uint Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public uint ModuloOperator() public uint ModuloOperator()
{ {
return Value % 6U; return Value % 6U;
@@ -20,7 +28,7 @@ public class Mod6UIntBenchmarks
[Benchmark] [Benchmark]
public uint ExtensionMethod() public uint ExtensionMethod()
{ {
return Value.Mod6(); return Mod6(Value);
} }
} }
@@ -1,4 +1,4 @@
using System.Numerics; using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -9,13 +9,20 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod6ULongBenchmarks public class Mod6ULongBenchmarks
{ {
private static readonly byte[] Mod6Lookup = [0, 3, 4, 1, 2, 5];
private static readonly ulong Divisor = 6UL; private static readonly ulong Divisor = 6UL;
private static readonly ulong FastDivMul = (ulong)(((UInt128)1 << 64) / Divisor); private static readonly ulong FastDivMul = (ulong)(((UInt128)1 << 64) / Divisor);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod6(ulong value) => Mod6Lookup[(int)(((value % 3UL) << 1) | (value & 1UL))];
[Params(3UL, 131071UL, ulong.MaxValue)] [Params(3UL, 131071UL, ulong.MaxValue)]
public ulong Value { get; set; } public ulong Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public ulong ModuloOperator() public ulong ModuloOperator()
{ {
return Value % Divisor; return Value % Divisor;
@@ -38,6 +45,6 @@ public class Mod6ULongBenchmarks
[Benchmark] [Benchmark]
public ulong ExtensionMethod() public ulong ExtensionMethod()
{ {
return Value.Mod6(); return Mod6(Value);
} }
} }
@@ -17,7 +17,11 @@ public class Mod7UInt128Benchmarks
return (ulong)(Value % 7UL); return (ulong)(Value % 7UL);
} }
[Benchmark] /// <summary>
/// Fastest
/// </summary>
/// <returns></returns>
[Benchmark]
public ulong ExtensionMethod() public ulong ExtensionMethod()
{ {
return Value.Mod7(); return Value.Mod7();
@@ -1,6 +1,6 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks; namespace EvenPerfectBitScanner.Benchmarks;
@@ -8,10 +8,30 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod7UIntBenchmarks public class Mod7UIntBenchmarks
{ {
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod7(uint value)
{
uint remainder = 0U;
while (value != 0U)
{
remainder += value & 7U;
value >>= 3;
if (remainder >= 7U)
{
remainder -= 7U;
}
}
return remainder >= 7U ? remainder - 7U : remainder;
}
[Params(7U, 2047U, 65535U, 2147483647U)] [Params(7U, 2047U, 65535U, 2147483647U)]
public uint Value { get; set; } public uint Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public uint ModuloOperator() public uint ModuloOperator()
{ {
return Value % 7U; return Value % 7U;
@@ -20,7 +40,7 @@ public class Mod7UIntBenchmarks
[Benchmark] [Benchmark]
public uint ExtensionMethod() public uint ExtensionMethod()
{ {
return Value.Mod7(); return Mod7(Value);
} }
} }
@@ -1,6 +1,6 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks; namespace EvenPerfectBitScanner.Benchmarks;
@@ -8,10 +8,25 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod7ULongBenchmarks public class Mod7ULongBenchmarks
{ {
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod7(ulong value)
{
ulong remainder = ((uint)value % 7U) + (((uint)(value >> 32) % 7U) << 2);
while (remainder >= 7UL)
{
remainder -= 7UL;
}
return remainder;
}
[Params(7UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 63UL)] [Params(7UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 63UL)]
public ulong Value { get; set; } public ulong Value { get; set; }
[Benchmark(Baseline = true)] /// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)]
public ulong ModuloOperator() public ulong ModuloOperator()
{ {
return Value % 7UL; return Value % 7UL;
@@ -20,7 +35,7 @@ public class Mod7ULongBenchmarks
[Benchmark] [Benchmark]
public ulong ExtensionMethod() public ulong ExtensionMethod()
{ {
return Value.Mod7(); return Mod7(Value);
} }
} }
@@ -1,4 +1,3 @@
using System.Numerics;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -18,7 +17,10 @@ public class Mod8UInt128Benchmarks
return (ulong)(Value % 8UL); return (ulong)(Value % 8UL);
} }
[Benchmark] /// <summary>
/// Fastest
/// </summary>
[Benchmark]
public ulong ExtensionMethod() public ulong ExtensionMethod()
{ {
return Value.Mod8(); return Value.Mod8();
@@ -17,6 +17,9 @@ public class Mod8UIntBenchmarks
return Value % 8U; return Value % 8U;
} }
/// <summary>
/// Fastest
/// </summary>
[Benchmark] [Benchmark]
public uint ExtensionMethod() public uint ExtensionMethod()
{ {
@@ -1,6 +1,6 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs; using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks; namespace EvenPerfectBitScanner.Benchmarks;
@@ -8,19 +8,25 @@ namespace EvenPerfectBitScanner.Benchmarks;
[MemoryDiagnoser] [MemoryDiagnoser]
public class Mod8ULongBenchmarks public class Mod8ULongBenchmarks
{ {
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod8(ulong value) => value & 7UL;
[Params(8UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 127UL)] [Params(8UL, 8191UL, 131071UL, 2147483647UL, ulong.MaxValue - 127UL)]
public ulong Value { get; set; } public ulong Value { get; set; }
/// <summary>
/// Fastest
/// </summary>
[Benchmark(Baseline = true)] [Benchmark(Baseline = true)]
public ulong ModuloOperator() public ulong ModuloOperator()
{ {
return Value % 8UL; return Value % 8UL;
} }
[Benchmark] [Benchmark]
public ulong ExtensionMethod() public ulong ExtensionMethod()
{ {
return Value.Mod8(); return Mod8(Value);
} }
} }
+1 -2
View File
@@ -3,7 +3,6 @@ using System.Globalization;
using System.Numerics; using System.Numerics;
using System.Runtime.CompilerServices; using System.Runtime.CompilerServices;
using System.Text; using System.Text;
using System.Threading;
using Open.Collections; using Open.Collections;
using Open.Numeric.Primes; using Open.Numeric.Primes;
using PerfectNumbers.Core; using PerfectNumbers.Core;
@@ -1674,7 +1673,7 @@ internal static class Program
// Optional: RLE blacklist and binary-threshold filters on p (safe only when configured) // Optional: RLE blacklist and binary-threshold filters on p (safe only when configured)
if (p <= _rleHardMaxP) if (p <= _rleHardMaxP)
{ {
if (!_rleOnlyLast7 || p % 10UL == 7UL) if (!_rleOnlyLast7 || p.Mod10() == 7UL)
{ {
if (RleBlacklist.IsLoaded() && RleBlacklist.Matches(p)) if (RleBlacklist.IsLoaded() && RleBlacklist.Matches(p))
{ {
@@ -13,7 +13,7 @@ public class UIntExtensionsTests
[InlineData(uint.MaxValue)] [InlineData(uint.MaxValue)]
public void Mod3_matches_operator(uint value) public void Mod3_matches_operator(uint value)
{ {
value.Mod3().Should().Be(value % 3U); UIntExtensions.Mod3(value).Should().Be(value % 3U);
} }
[Theory] [Theory]
@@ -27,28 +27,6 @@ public class UIntExtensionsTests
value.Mod5().Should().Be(value % 5U); value.Mod5().Should().Be(value % 5U);
} }
[Theory]
[Trait("Category", "Fast")]
[InlineData(0U)]
[InlineData(1U)]
[InlineData(1234567890U)]
[InlineData(uint.MaxValue)]
public void Mod6_matches_operator(uint value)
{
value.Mod6().Should().Be(value % 6U);
}
[Theory]
[Trait("Category", "Fast")]
[InlineData(0U)]
[InlineData(1U)]
[InlineData(1234567890U)]
[InlineData(uint.MaxValue)]
public void Mod7_matches_operator(uint value)
{
value.Mod7().Should().Be(value % 7U);
}
[Theory] [Theory]
[Trait("Category", "Fast")] [Trait("Category", "Fast")]
[InlineData(0U)] [InlineData(0U)]
@@ -59,38 +37,5 @@ public class UIntExtensionsTests
{ {
value.Mod8().Should().Be(value % 8U); value.Mod8().Should().Be(value % 8U);
} }
[Theory]
[Trait("Category", "Fast")]
[InlineData(0U)]
[InlineData(1U)]
[InlineData(1234567890U)]
[InlineData(uint.MaxValue)]
public void Mod10_matches_operator(uint value)
{
value.Mod10().Should().Be(value % 10U);
}
[Theory]
[Trait("Category", "Fast")]
[InlineData(0U)]
[InlineData(1U)]
[InlineData(1234567890U)]
[InlineData(uint.MaxValue)]
public void Mod11_matches_operator(uint value)
{
value.Mod11().Should().Be(value % 11U);
}
[Theory]
[Trait("Category", "Fast")]
[InlineData(0U)]
[InlineData(1U)]
[InlineData(1234567890U)]
[InlineData(uint.MaxValue)]
public void Mod128_matches_operator(uint value)
{
value.Mod128().Should().Be(value % 128U);
}
} }
@@ -21,66 +21,6 @@ public class ULongExtensionsTests
ULongExtensions.Mod10(value).Should().Be(expected); ULongExtensions.Mod10(value).Should().Be(expected);
} }
[Theory]
[Trait("Category", "Fast")]
[InlineData(0UL)]
[InlineData(ulong.MaxValue)]
[InlineData(1234567890123456789UL)]
public void Mod3_matches_operator(ulong value)
{
value.Mod3().Should().Be(value % 3UL);
}
[Theory]
[Trait("Category", "Fast")]
[InlineData(0UL)]
[InlineData(ulong.MaxValue)]
[InlineData(1234567890123456789UL)]
public void Mod5_matches_operator(ulong value)
{
value.Mod5().Should().Be(value % 5UL);
}
[Theory]
[Trait("Category", "Fast")]
[InlineData(0UL)]
[InlineData(ulong.MaxValue)]
[InlineData(1234567890123456789UL)]
public void Mod6_matches_operator(ulong value)
{
value.Mod6().Should().Be(value % 6UL);
}
[Theory]
[Trait("Category", "Fast")]
[InlineData(0UL)]
[InlineData(ulong.MaxValue)]
[InlineData(1234567890123456789UL)]
public void Mod7_matches_operator(ulong value)
{
value.Mod7().Should().Be(value % 7UL);
}
[Theory]
[Trait("Category", "Fast")]
[InlineData(0UL)]
[InlineData(ulong.MaxValue)]
[InlineData(1234567890123456789UL)]
public void Mod8_matches_operator(ulong value)
{
value.Mod8().Should().Be(value % 8UL);
}
[Theory]
[Trait("Category", "Fast")]
[InlineData(0UL)]
[InlineData(ulong.MaxValue)]
[InlineData(1234567890123456789UL)]
public void Mod11_matches_operator(ulong value)
{
value.Mod11().Should().Be(value % 11UL);
}
[Theory] [Theory]
[Trait("Category", "Fast")] [Trait("Category", "Fast")]
[InlineData(0UL)] [InlineData(0UL)]
@@ -111,22 +51,6 @@ public class ULongExtensionsTests
step3.Should().Be(expected3); step3.Should().Be(expected3);
} }
[Theory]
[Trait("Category", "Fast")]
[InlineData(0UL)]
[InlineData(1UL)]
[InlineData(1234567890123456789UL)]
[InlineData(ulong.MaxValue)]
public void Mod11_7_5_3_matches_operator(ulong value)
{
value.Mod11_7_5_3(out ulong mod11, out ulong mod7, out ulong mod5, out ulong mod3);
mod11.Should().Be(value % 11UL);
mod7.Should().Be(value % 7UL);
mod5.Should().Be(value % 5UL);
mod3.Should().Be(value % 3UL);
}
[Theory] [Theory]
[Trait("Category", "Fast")] [Trait("Category", "Fast")]
[InlineData(2UL, true)] [InlineData(2UL, true)]
@@ -9,10 +9,10 @@ public class MersenneNumberLucasLehmerCpuTester
// - If p ≡ 1 (mod 4) and p shares a factor with (p-1), reject fast. // - If p ≡ 1 (mod 4) and p shares a factor with (p-1), reject fast.
// - If p is divisible by divisors specific to numbers ending with 1 or 7, reject fast. // - If p is divisible by divisors specific to numbers ending with 1 or 7, reject fast.
if ( if (
(exponent.Mod3() == 0UL && exponent != 3UL) || ((exponent % 3UL) == 0UL && exponent != 3UL) ||
(exponent.Mod5() == 0UL && exponent != 5UL) || ((exponent % 5UL) == 0UL && exponent != 5UL) ||
(exponent.Mod7() == 0UL && exponent != 7UL) || ((exponent % 7UL) == 0UL && exponent != 7UL) ||
(exponent.Mod11() == 0UL && exponent != 11UL) ((exponent % 11UL) == 0UL && exponent != 11UL)
) )
{ {
return false; return false;
+10 -11
View File
@@ -1,6 +1,5 @@
using System.Collections.Concurrent; using System.Collections.Concurrent;
using System.Runtime.CompilerServices; using System.Runtime.CompilerServices;
using System.Threading;
using ILGPU; using ILGPU;
using ILGPU.Runtime; using ILGPU.Runtime;
using static PerfectNumbers.Core.Gpu.GpuContextPool; using static PerfectNumbers.Core.Gpu.GpuContextPool;
@@ -200,8 +199,8 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
ArrayView1D<uint, Stride1D.Dense> smallCycles) ArrayView1D<uint, Stride1D.Dense> smallCycles)
{ {
ulong idx = (ulong)index.X; ulong idx = (ulong)index.X;
ulong idxMod3 = idx.Mod3(); ulong idxMod3 = idx % 3UL;
ulong idxMod5 = idx.Mod5(); ulong idxMod5 = idx % 5UL;
// residue automaton // residue automaton
ulong step10 = (exponent.Mod10() << 1).Mod10(); ulong step10 = (exponent.Mod10() << 1).Mod10();
ulong r10 = q0m10 + (step10 * idx).Mod10(); ulong r10 = q0m10 + (step10 * idx).Mod10();
@@ -217,7 +216,7 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
} }
else else
{ {
ulong step3 = (exponent.Mod3() << 1).Mod3(); ulong step3 = ((exponent % 3UL) << 1) % 3UL;
ulong r3 = q0m3 + (step3 * idxMod3); r3 -= (r3 >= 3UL) ? 3UL : 0UL; ulong r3 = q0m3 + (step3 * idxMod3); r3 -= (r3 >= 3UL) ? 3UL : 0UL;
if (r3 == 0UL) if (r3 == 0UL)
{ {
@@ -225,7 +224,7 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
} }
else else
{ {
ulong step5Loc = (exponent.Mod5() << 1).Mod5(); ulong step5Loc = ((exponent % 5UL) << 1) % 5UL;
ulong r5 = q0m5 + (step5Loc * idxMod5); if (r5 >= 5UL) r5 -= 5UL; ulong r5 = q0m5 + (step5Loc * idxMod5); if (r5 >= 5UL) r5 -= 5UL;
if (r5 == 0UL) if (r5 == 0UL)
{ {
@@ -360,8 +359,8 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
ArrayView1D<ulong, Stride1D.Dense> smallPrimesPow2LastSeven) ArrayView1D<ulong, Stride1D.Dense> smallPrimesPow2LastSeven)
{ {
ulong idx = (ulong)index.X; ulong idx = (ulong)index.X;
ulong idxMod3 = idx.Mod3(); ulong idxMod3 = idx % 3UL;
ulong idxMod5 = idx.Mod5(); ulong idxMod5 = idx % 5UL;
ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL;
bool shouldCheck = r10 != 5UL; bool shouldCheck = r10 != 5UL;
if (shouldCheck) if (shouldCheck)
@@ -493,8 +492,8 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
ResidueAutomatonArgs ra, ArrayView<int> found, ArrayView1D<uint, Stride1D.Dense> smallCycles) ResidueAutomatonArgs ra, ArrayView<int> found, ArrayView1D<uint, Stride1D.Dense> smallCycles)
{ {
ulong idx = (ulong)index.X; ulong idx = (ulong)index.X;
ulong idxMod3 = idx.Mod3(); ulong idxMod3 = idx % 3UL;
ulong idxMod5 = idx.Mod5(); ulong idxMod5 = idx % 5UL;
ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL;
bool shouldCheck = r10 != 5UL; bool shouldCheck = r10 != 5UL;
if (shouldCheck) if (shouldCheck)
@@ -579,8 +578,8 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
ResidueAutomatonArgs ra, ArrayView<int> found, ArrayView1D<uint, Stride1D.Dense> smallCycles) ResidueAutomatonArgs ra, ArrayView<int> found, ArrayView1D<uint, Stride1D.Dense> smallCycles)
{ {
ulong idx = (ulong)index.X; ulong idx = (ulong)index.X;
ulong idxMod3 = idx.Mod3(); ulong idxMod3 = idx % 3UL;
ulong idxMod5 = idx.Mod5(); ulong idxMod5 = idx % 5UL;
ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL; ulong r10 = ra.Q0M10 + (ra.Step10 * idx).Mod10(); r10 -= (r10 >= 10UL) ? 10UL : 0UL;
bool shouldCheck = r10 != 5UL; bool shouldCheck = r10 != 5UL;
if (shouldCheck) if (shouldCheck)
@@ -1,6 +1,5 @@
using System.Buffers; using System.Buffers;
using System.Runtime.InteropServices; using System.Runtime.InteropServices;
using ILGPU;
using ILGPU.Runtime; using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core.Gpu;
@@ -10,98 +9,95 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u
private readonly GpuKernelType _kernelType = kernelType; private readonly GpuKernelType _kernelType = kernelType;
private readonly bool _useGpuOrder = useGpuOrder; private readonly bool _useGpuOrder = useGpuOrder;
public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime) public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime)
{ {
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder); var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var execution = gpuLease.EnterExecutionScope(); var execution = gpuLease.EnterExecutionScope();
var accelerator = gpuLease.Accelerator; var accelerator = gpuLease.Accelerator;
var stream = gpuLease.Stream; var stream = gpuLease.Stream;
int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy
UInt128 kStart = 1UL; UInt128 kStart = 1UL;
ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL; ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL;
byte last = lastIsSeven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters byte last = lastIsSeven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters
var pow2Kernel = gpuLease.Pow2ModKernel; var pow2Kernel = gpuLease.Pow2ModKernel;
var incKernel = gpuLease.IncrementalKernel; var incKernel = gpuLease.IncrementalKernel;
ulong step10 = (exponent.Mod10() << 1).Mod10(); exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
ulong step8 = ((exponent & 7UL) << 1) & 7UL; GpuUInt128 twoPGpu = (GpuUInt128)twoP;
ulong step3 = (exponent.Mod3() << 1).Mod3(); var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
ulong step5 = (exponent.Mod5() << 1).Mod5(); ResiduePrimeViews primeViews = default;
GpuUInt128 twoPGpu = (GpuUInt128)twoP; if (_kernelType == GpuKernelType.Pow2Mod)
var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator); {
ResiduePrimeViews primeViews = default; primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
if (_kernelType == GpuKernelType.Pow2Mod) }
{
primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
}
var orderBuffer = accelerator.Allocate1D<ulong>(batchSize); var orderBuffer = accelerator.Allocate1D<ulong>(batchSize);
// Avoid giant stack allocations that can trigger StackOverflow when batchSize is large. // Avoid giant stack allocations that can trigger StackOverflow when batchSize is large.
// Rent a reusable array from the shared pool instead. // Rent a reusable array from the shared pool instead.
ulong[] orderArray = ArrayPool<ulong>.Shared.Rent(batchSize); ulong[] orderArray = ArrayPool<ulong>.Shared.Rent(batchSize);
UInt128 remaining; UInt128 remaining;
int currentSize; int currentSize;
int i; int i;
UInt128 q = UInt128.Zero; UInt128 q = UInt128.Zero;
try try
{ {
while (kStart <= maxK && Volatile.Read(ref isPrime)) while (kStart <= maxK && Volatile.Read(ref isPrime))
{ {
remaining = maxK - kStart + UInt128.One; remaining = maxK - kStart + UInt128.One;
currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining; currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining;
Span<ulong> orders = orderArray.AsSpan(0, currentSize); Span<ulong> orders = orderArray.AsSpan(0, currentSize);
// Precompute residue automaton bases for this batch // Precompute residue automaton bases for this batch
UInt128 q0 = twoP * kStart + UInt128.One; UInt128 q0 = twoP * kStart + UInt128.One;
if (_kernelType == GpuKernelType.Pow2Mod) if (_kernelType == GpuKernelType.Pow2Mod)
{ {
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3); q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5); var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
pow2Kernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, divMul, pow2Kernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, divMul,
kernelArgs, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven, kernelArgs, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven,
primeViews.LastOnePow2, primeViews.LastSevenPow2); primeViews.LastOnePow2, primeViews.LastSevenPow2);
} }
else else
{ {
q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3); q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
incKernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, divMul, incKernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, divMul,
q0m10, q0m8, q0m3, q0m5, orderBuffer.View, smallCyclesView); q0m10, q0m8, q0m3, q0m5, orderBuffer.View, smallCyclesView);
} }
stream.Synchronize(); stream.Synchronize();
orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize); orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize);
if (_kernelType == GpuKernelType.Pow2Mod) if (_kernelType == GpuKernelType.Pow2Mod)
{ {
for (i = 0; i < currentSize; i++) for (i = 0; i < currentSize; i++)
{ {
if (orders[i] != 0UL) if (orders[i] != 0UL)
{ {
Volatile.Write(ref isPrime, false); Volatile.Write(ref isPrime, false);
break; break;
} }
} }
} }
else else
{ {
q = twoP * kStart + UInt128.One; q = twoP * kStart + UInt128.One;
for (i = 0; i < currentSize && Volatile.Read(ref isPrime); i++, q += twoP) for (i = 0; i < currentSize && Volatile.Read(ref isPrime); i++, q += twoP)
{ {
if (orders[i] != 0UL && q.IsPrimeCandidate()) if (orders[i] != 0UL && q.IsPrimeCandidate())
{ {
Volatile.Write(ref isPrime, false); Volatile.Write(ref isPrime, false);
} }
} }
} }
kStart += (UInt128)currentSize; kStart += (UInt128)currentSize;
} }
} }
finally finally
{ {
ArrayPool<ulong>.Shared.Return(orderArray); ArrayPool<ulong>.Shared.Return(orderArray);
orderBuffer.Dispose(); orderBuffer.Dispose();
execution.Dispose(); execution.Dispose();
gpuLease.Dispose(); gpuLease.Dispose();
} }
} }
} }
@@ -51,7 +51,7 @@ public class MersenneNumberLucasLehmerGpuTester
// Early rejections aligned with incremental/order sieves, but safe for small p: // Early rejections aligned with incremental/order sieves, but safe for small p:
// - If 3 | p and p != 3, then 7 | M_p -> composite. // - If 3 | p and p != 3, then 7 | M_p -> composite.
// - If p ≡ 1 (mod 4) and p shares a factor with (p-1), reject fast. // - If p ≡ 1 (mod 4) and p shares a factor with (p-1), reject fast.
if (exponent.Mod3() == 0UL && exponent != 3UL) if ((exponent % 3UL) == 0UL && exponent != 3UL)
{ {
return false; return false;
} }
@@ -25,12 +25,7 @@ public class MersenneNumberOrderGpuTester(GpuKernelType kernelType, bool useGpuO
}; };
var foundBuffer = accelerator.Allocate1D<int>(1); var foundBuffer = accelerator.Allocate1D<int>(1);
UInt128 exponent128 = exponent; exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
exponent128.Mod10_8_5_3(out ulong exponentMod10, out ulong exponentMod8, out ulong exponentMod5, out ulong exponentMod3);
ulong step10 = (exponentMod10 << 1).Mod10();
ulong step8 = (exponentMod8 << 1) & 7UL;
ulong step3 = (exponentMod3 << 1).Mod3();
ulong step5 = (exponentMod5 << 1).Mod5();
try try
{ {
UInt128 remaining; UInt128 remaining;
@@ -24,7 +24,11 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
UInt128 kStart = 1UL; UInt128 kStart = 1UL;
byte last = lastIsSeven ? (byte)1 : (byte)0; byte last = lastIsSeven ? (byte)1 : (byte)0;
var kernel = gpuLease.Pow2ModKernel; var kernel = gpuLease.Pow2ModKernel;
exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3); ulong step3 = exponent % 3UL,
step5 = exponent % 5UL,
step8 = exponent % 8UL,
step10 = exponent.Mod10();
GpuUInt128 twoPGpu = (GpuUInt128)twoP; GpuUInt128 twoPGpu = (GpuUInt128)twoP;
var orderBuffer = accelerator.Allocate1D<ulong>(batchSize); var orderBuffer = accelerator.Allocate1D<ulong>(batchSize);
+24 -24
View File
@@ -189,7 +189,7 @@ public class MersenneDivisorCycles
for (int i = 0; i < localCycleIndex; i++) for (int i = 0; i < localCycleIndex; i++)
{ {
(d, cycle) = localCycles[i]; (d, cycle) = localCycles[i];
if (d.Mod3() == 0UL || d.Mod5() == 0UL || d.Mod7() == 0UL || d.Mod11() == 0UL) if ((d % 3UL) == 0UL || (d % 5UL) == 0UL || (d % 7UL) == 0UL || (d % 11UL) == 0UL)
continue; continue;
writer.Write(d); writer.Write(d);
@@ -222,15 +222,15 @@ public class MersenneDivisorCycles
int count = (int)Math.Min(batchSizeUL, maxDivisor), i, idx; int count = (int)Math.Min(batchSizeUL, maxDivisor), i, idx;
ulong[] divisors, outCycles, validDivisors; ulong[] divisors, outCycles, validDivisors;
divisors = pool.Rent(count); divisors = pool.Rent(count);
outCycles = pool.Rent(count); outCycles = pool.Rent(count);
var lease = GpuKernelPool.GetKernel(useGpuOrder: true); var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var accelerator = lease.Accelerator; var accelerator = lease.Accelerator;
var stream = accelerator.CreateStream(); var stream = accelerator.CreateStream();
var kernel = accelerator.LoadAutoGroupedStreamKernel< var kernel = accelerator.LoadAutoGroupedStreamKernel<
Index1D, Index1D,
ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>,
ArrayView1D<ulong, Stride1D.Dense>>(GpuDivisorCycleKernel); ArrayView1D<ulong, Stride1D.Dense>>(GpuDivisorCycleKernel);
using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: true); using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: true);
@@ -249,7 +249,7 @@ public class MersenneDivisorCycles
{ {
while (skipCount > 0L && start <= maxDivisor) while (skipCount > 0L && start <= maxDivisor)
{ {
if ((start & 1UL) == 0UL || start.Mod3() == 0UL || start.Mod5() == 0UL || start.Mod7() == 0UL || start.Mod11() == 0UL) if ((start & 1UL) == 0UL || (start % 3UL) == 0UL || (start % 5UL) == 0UL || (start % 7UL) == 0UL || (start % 11UL) == 0UL)
{ {
start++; start++;
continue; continue;
@@ -280,10 +280,10 @@ public class MersenneDivisorCycles
var bufferDiv = accelerator.Allocate1D(validDivisors); var bufferDiv = accelerator.Allocate1D(validDivisors);
var bufferCycle = accelerator.Allocate1D<ulong>(idx); var bufferCycle = accelerator.Allocate1D<ulong>(idx);
kernel( kernel(
idx, idx,
bufferDiv.View, bufferDiv.View,
bufferCycle.View); bufferCycle.View);
accelerator.Synchronize(); accelerator.Synchronize();
@@ -296,7 +296,7 @@ public class MersenneDivisorCycles
{ {
d = divisors[i]; d = divisors[i];
if (d.Mod3() == 0UL || d.Mod5() == 0UL || d.Mod7() == 0UL || d.Mod11() == 0UL) if ((d % 3UL) == 0UL || (d % 5UL) == 0UL || (d % 7UL) == 0UL || (d % 11UL) == 0UL)
continue; continue;
writer.Write(d); writer.Write(d);
@@ -308,16 +308,16 @@ public class MersenneDivisorCycles
start = end + 1UL; start = end + 1UL;
} }
} }
finally finally
{ {
pool.Return(divisors, clearArray: false); pool.Return(divisors, clearArray: false);
pool.Return(outCycles, clearArray: false); pool.Return(outCycles, clearArray: false);
} }
stream.Dispose(); stream.Dispose();
lease.Dispose(); lease.Dispose();
} }
const byte ByteZero = 0; const byte ByteZero = 0;
const byte ByteOne = 1; const byte ByteOne = 1;
+1 -1
View File
@@ -203,7 +203,7 @@ public sealed class MersenneNumberTester(
{ {
// Safe early rejections using known orders for tiny primes: // Safe early rejections using known orders for tiny primes:
// 7 | M_p iff 3 | p. Avoid rejecting p == 3 where M_3 == 7 is prime. // 7 | M_p iff 3 | p. Avoid rejecting p == 3 where M_3 == 7 is prime.
if (exponent.Mod3() == 0UL && exponent != 3UL) if ((exponent % 3UL) == 0UL && exponent != 3UL)
{ {
return false; return false;
} }
+2 -2
View File
@@ -53,7 +53,7 @@ public sealed class PrimeTester(bool useInternal = false)
{ {
result = false; result = false;
} }
else if (n > 5UL && n.Mod5() == 0UL) else if (n > 5UL && (n % 5UL) == 0UL)
{ {
result = false; result = false;
} }
@@ -260,7 +260,7 @@ public sealed class PrimeTester(bool useInternal = false)
return; return;
} }
if ((n & 1UL) == 0UL || (n > 5UL && n.Mod5() == 0UL)) if ((n & 1UL) == 0UL || (n > 5UL && (n % 5UL) == 0UL))
{ {
results[index] = 0; results[index] = 0;
return; return;
+149 -171
View File
@@ -96,109 +96,147 @@ public static class UInt128Extensions
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsPrimeCandidate(this UInt128 n) public static bool IsPrimeCandidate(this UInt128 n)
{ {
UInt128 p; UInt128 p, zero = UInt128.Zero;
uint[] smallPrimes = PrimesGenerator.SmallPrimes; uint[] smallPrimes = PrimesGenerator.SmallPrimes;
ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2; ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2;
ulong i, smallPrimesCount = (ulong)smallPrimes.Length; ulong i, smallPrimesCount = (ulong)smallPrimes.Length;
for (i = 0UL; i < smallPrimesCount; i++) for (i = 0UL; i < smallPrimesCount; i++)
{ {
p = smallPrimes[i]; if (smallPrimesPow2[i] > n)
if (smallPrimesPow2[i] > n) {
{ break;
break; }
}
if (n % p == UInt128.Zero) p = smallPrimes[i];
{ if (n % p == zero)
return n == p; {
} return n == p;
} }
}
return true; return true;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void Mod10_8_5_3(this UInt128 value, out ulong mod10, out ulong mod8, out ulong mod5, out ulong mod3) public static ulong Mod3(this UInt128 value)
{ {
UInt128 temp = UInt128.Zero; ulong remainder = ((ulong)value) % 3UL + ((ulong)(value >> 64)) % 3UL;
uint byteSum = 0U; return remainder >= 3UL ? remainder - 3UL : remainder;
byte current = 0; }
uint accumulator = 0U;
temp = value;
do
{
current = (byte)temp;
byteSum += current;
temp >>= 8;
}
while (temp != UInt128.Zero);
mod8 = (ulong)value & 7UL;
accumulator = byteSum.Mod5();
mod5 = accumulator;
mod10 = (mod8 & 1UL) == 0UL
? accumulator switch
{
0U => 0UL,
1U => 6UL,
2U => 2UL,
3U => 8UL,
_ => 4UL,
}
: accumulator switch
{
0U => 5UL,
1U => 1UL,
2U => 7UL,
3U => 3UL,
_ => 9UL,
};
// Reusing accumulator to hold the mod 3 reduction now that the mod 5 value has been stored.
accumulator = byteSum.Mod3();
mod3 = accumulator;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 Pow(this UInt128 value, ulong exponent) public static ulong Mod5(this UInt128 value)
{ {
UInt128 result = UInt128.One; ulong remainder = (((ulong)value) % 5UL) + (((ulong)(value >> 64)) % 5UL);
UInt128 baseValue = value; return remainder >= 5UL ? remainder - 5UL : remainder;
}
while (exponent != 0UL) [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod6(this UInt128 value) => ((value.Mod3() << 1) | ((ulong)value & 1UL)) switch
{
0UL => 0UL,
1UL => 3UL,
2UL => 4UL,
3UL => 1UL,
4UL => 2UL,
_ => 5UL,
};
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod7(this UInt128 value)
{
ulong remainder = (((ulong)value) % 7UL) + (((ulong)(value >> 64)) % 7UL) << 1;
while (remainder >= 7UL)
{ {
if ((exponent & 1UL) != 0UL) remainder -= 7UL;
{
result *= baseValue;
}
baseValue *= baseValue;
exponent >>= 1;
} }
return result; return remainder;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod8(this UInt128 value) => (ulong)value & 7UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod10(this UInt128 value128)
{
// Split and fold under mod 10: 2^64 ≡ 6 (mod 10)
ulong value = (ulong)value128;
return ((value % 10UL) + ((value >> 64) % 10UL) * 6UL) % 10UL;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void Mod10_8_5_3(this UInt128 value, out ulong mod10, out ulong mod8, out ulong mod5, out ulong mod3)
{
byte current;
uint byteSum = 0U;
mod8 = (ulong)value & 7UL;
UInt128 zero = UInt128.Zero;
do
{
current = (byte)value;
byteSum += current;
value >>= 8;
}
while (value != zero);
mod5 = byteSum.Mod5();
mod10 = (mod8 & 1UL) == 0UL
? mod5 switch
{
0U => 0UL,
1U => 6UL,
2U => 2UL,
3U => 8UL,
_ => 4UL,
}
: mod5 switch
{
0U => 5UL,
1U => 1UL,
2U => 7UL,
3U => 3UL,
_ => 9UL,
};
mod3 = byteSum % 3U;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod11(this UInt128 value)
{
ulong remainder = (((ulong)value) % 11UL) + (((ulong)(value >> 64)) % 11UL) * 5UL;
while (remainder >= 11UL)
{
remainder -= 11UL;
}
return remainder;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod128(this UInt128 value) => (ulong)value & 127UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 ModPow(this UInt128 value, UInt128 exponent, UInt128 modulus) public static UInt128 ModPow(this UInt128 value, UInt128 exponent, UInt128 modulus)
{ {
UInt128 zero = UInt128.Zero; UInt128 zero = UInt128.Zero;
UInt128 result = UInt128.One; UInt128 one = UInt128.One;
UInt128 baseValue = value % modulus; UInt128 result = UInt128.One;
value %= modulus;
while (exponent != zero) while (exponent != zero)
{ {
if ((exponent & UInt128.One) != zero) if ((exponent & one) != zero)
{ {
result = MulMod(result, baseValue, modulus); result = MulMod(result, value, modulus);
} }
baseValue = MulMod(baseValue, baseValue, modulus); value = MulMod(value, value, modulus);
exponent >>= 1; exponent >>= 1;
} }
@@ -206,87 +244,8 @@ public static class UInt128Extensions
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod10(this UInt128 value) public static UInt128 Mul64(this UInt128 a, UInt128 b)
{ {
UInt128 zero = UInt128.Zero;
if (value == zero)
return 0UL;
ulong result = (ulong)value;
UInt128 high = value >> 64;
while (high != zero)
{
// 2^64 ≡ 6 (mod 10)
result = (result + (ulong)high * 6UL) % 10UL;
high >>= 64;
}
return result % 10UL;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod8(this UInt128 value) => (ulong)value & 7UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod3(this UInt128 value)
{
ulong remainder = ((ulong)value % 3UL) + ((ulong)(value >> 64) % 3UL);
return remainder >= 3UL ? remainder - 3UL : remainder;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod5(this UInt128 value)
{
ulong remainder = ((ulong)value % 5UL) + ((ulong)(value >> 64) % 5UL);
return remainder >= 5UL ? remainder - 5UL : remainder;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod6(this UInt128 value)
{
return ((value.Mod3() << 1) | ((ulong)value & 1UL)) switch
{
0UL => 0UL,
1UL => 3UL,
2UL => 4UL,
3UL => 1UL,
4UL => 2UL,
_ => 5UL,
};
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod7(this UInt128 value)
{
ulong remainder = ((ulong)value % 7UL) + ((ulong)(value >> 64) % 7UL) * 2UL;
while (remainder >= 7UL)
{
remainder -= 7UL;
}
return remainder;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod11(this UInt128 value)
{
ulong remainder = ((ulong)value % 11UL) + ((ulong)(value >> 64) % 11UL) * 5UL;
while (remainder >= 11UL)
{
remainder -= 11UL;
}
return remainder;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod128(this UInt128 value) => (ulong)value & 127UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 Mul64(this UInt128 a, UInt128 b)
{
ulong aLow = (ulong)a; ulong aLow = (ulong)a;
ulong bLow = (ulong)b; ulong bLow = (ulong)b;
return ((UInt128)(aLow * (ulong)(b >> 64) + aLow.MulHigh(bLow)) << 64) | (aLow * bLow); return ((UInt128)(aLow * (ulong)(b >> 64) + aLow.MulHigh(bLow)) << 64) | (aLow * bLow);
@@ -295,15 +254,16 @@ public static class UInt128Extensions
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 MulMod(this UInt128 a, UInt128 b, UInt128 modulus) public static UInt128 MulMod(this UInt128 a, UInt128 b, UInt128 modulus)
{ {
UInt128 zero = UInt128.Zero; UInt128 one = UInt128.One,
UInt128 result = UInt128.Zero; result = UInt128.Zero,
zero = UInt128.Zero;
while (b != zero) while (b != zero)
{ {
if ((b & UInt128.One) != zero) if ((b & one) != zero)
{ {
result += a; result += a;
if (result >= modulus) if (result >= modulus)
{ {
result -= modulus; result -= modulus;
} }
@@ -320,4 +280,22 @@ public static class UInt128Extensions
return result; return result;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 Pow(this UInt128 value, ulong exponent)
{
UInt128 result = UInt128.One;
while (exponent != 0UL)
{
if ((exponent & 1UL) != 0UL)
{
result *= value;
}
value *= value;
exponent >>= 1;
}
return result;
}
} }
+8 -20
View File
@@ -4,28 +4,16 @@ namespace PerfectNumbers.Core;
public static class UIntExtensions public static class UIntExtensions
{ {
[MethodImpl(MethodImplOptions.AggressiveInlining)] private const ulong Mod3Mask = 0xAAAAAAABUL;
public static uint Mod3(this uint value) => value % 3U; public const ulong Mod5Mask = 0xCCCCCCCDUL;
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod5(this uint value) => value % 5U; public static uint Mod3(uint value) => value - (uint)((value * Mod3Mask) >> 33) * 3U;
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod6(this uint value) => value % 6U; public static uint Mod5(this uint value) => value - (uint)((value * Mod5Mask) >> 34) * 5U;
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod7(this uint value) => value % 7U; public static uint Mod8(this uint value) => value & 7U;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod8(this uint value) => value & 7U;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod10(this uint value) => value % 10U;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod11(this uint value) => value % 11U;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static uint Mod128(this uint value) => value & 127U;
} }
+322 -466
View File
@@ -48,196 +48,110 @@ public static class ULongExtensions
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong FastDiv64(this ulong value, ulong divisor, ulong mul) public static ulong FastDiv64(this ulong value, ulong divisor, ulong mul)
{ {
ulong q = value.MulHigh(mul); mul = value.MulHigh(mul);
if ((UInt128)value - q.Mul64(divisor) >= divisor) if ((UInt128)value - mul.Mul64(divisor) >= divisor)
{ {
q++; mul++;
} }
return q; return mul;
} }
private static readonly byte[] Mod7ByteCoefficients = { 1, 4, 2, 1, 4, 2, 1, 4 }; public const ulong WordBitMask = 0xFFFFUL;
private static readonly byte[] Mod11ByteCoefficients = { 1, 3, 9, 5, 4, 1, 3, 9 };
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void Mod10_8_5_3(this ulong value, out ulong mod10, out ulong mod8, out ulong mod5, out ulong mod3) public static bool IsPrimeCandidate(this ulong n)
{ {
ulong temp = 0UL; int i = 0;
uint byteSum = 0U; uint[] smallPrimes = PrimesGenerator.SmallPrimes;
byte current = 0; ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2;
uint accumulator = 0U; int len = smallPrimes.Length;
ulong p;
for (; i < len; i++)
{
if (smallPrimesPow2[i] > n)
{
break;
}
temp = value; p = smallPrimes[i];
do if ((n % p) == 0UL)
{ {
current = (byte)temp; return n == p;
byteSum += current; }
temp >>= 8; }
}
while (temp != 0UL);
mod8 = value & 7UL; return true;
}
accumulator = byteSum.Mod5(); [MethodImpl(MethodImplOptions.AggressiveInlining)]
mod5 = accumulator; public static ulong Mod10(this ulong value) => (value & 1UL) == 0UL
mod10 = (mod8 & 1UL) == 0UL ? value % 5UL switch
? accumulator switch {
{ 0UL => 0UL,
0U => 0UL, 1UL => 6UL,
1U => 6UL, 2UL => 2UL,
2U => 2UL, 3UL => 8UL,
3U => 8UL, _ => 4UL,
_ => 4UL, }
} : value % 5UL switch
: accumulator switch {
{ 0UL => 5UL,
0U => 5UL, 1UL => 1UL,
1U => 1UL, 2UL => 7UL,
2U => 7UL, 3UL => 3UL,
3U => 3UL, _ => 9UL,
_ => 9UL, };
};
// Reusing accumulator to hold the mod 3 reduction now that the mod 5 value has been stored. [MethodImpl(MethodImplOptions.AggressiveInlining)]
accumulator = byteSum.Mod3(); public static ulong Mod128(this ulong value) => value & 127UL;
mod3 = accumulator;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void Mod11_7_5_3(this ulong value, out ulong mod11, out ulong mod7, out ulong mod5, out ulong mod3) public static void Mod10_8_5_3(this ulong value, out ulong mod10, out ulong mod8, out ulong mod5, out ulong mod3)
{ {
uint mod3Accumulator = 0U; mod8 = value & 7UL;
uint mod5Accumulator = 0U;
uint mod7Accumulator = 0U;
uint mod11Accumulator = 0U;
ulong temp = 0UL;
int index = 0;
byte current = 0;
temp = value; byte current;
do uint byteSum = 0U;
{ do
current = (byte)temp; {
mod3Accumulator += current; current = (byte)value;
mod5Accumulator += current; byteSum += current;
mod7Accumulator += (uint)(current * Mod7ByteCoefficients[index]); value >>= 8;
mod11Accumulator += (uint)(current * Mod11ByteCoefficients[index]); }
while (value != 0UL);
temp >>= 8; mod3 = byteSum % 3U;
index++; mod5 = byteSum % 5U;
} mod10 = (mod8 & 1UL) == 0UL
while (temp != 0UL); ? mod5 switch
{
0UL => 0UL,
1UL => 6UL,
2UL => 2UL,
3UL => 8UL,
_ => 4UL,
}
: mod5 switch
{
0UL => 5UL,
1UL => 1UL,
2UL => 7UL,
3UL => 3UL,
_ => 9UL,
};
}
// Reusing mod3Accumulator to store the reduced mod 3 remainder before exposing it. [MethodImpl(MethodImplOptions.AggressiveInlining)]
mod3Accumulator = mod3Accumulator.Mod3(); public static void Mod10_8_5_3Steps(this ulong value, out ulong step10, out ulong step8, out ulong step5, out ulong step3)
// Reusing mod5Accumulator to store the reduced mod 5 remainder before exposing it. {
mod5Accumulator = mod5Accumulator.Mod5(); value.Mod10_8_5_3(out step10, out step8, out step5, out step3);
// Reusing mod7Accumulator to store the reduced mod 7 remainder before exposing it.
mod7Accumulator = (uint)ReduceMod7(mod7Accumulator);
// Reusing mod11Accumulator to store the reduced mod 11 remainder before exposing it.
mod11Accumulator = (uint)ReduceMod11(mod11Accumulator);
mod3 = mod3Accumulator; step10 = (step10 << 1).Mod10();
mod5 = mod5Accumulator; step8 = (step8 << 1) & 7UL;
mod7 = mod7Accumulator; step5 = (step5 << 1) % 5UL;
mod11 = mod11Accumulator; step3 = (step3 << 1) % 3UL;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void Mod10_8_5_3Steps(this ulong value, out ulong step10, out ulong step8, out ulong step5, out ulong step3)
{
ulong mod10 = 0UL;
ulong mod8 = 0UL;
ulong mod5 = 0UL;
ulong mod3 = 0UL;
value.Mod10_8_5_3(out mod10, out mod8, out mod5, out mod3);
step10 = (mod10 << 1).Mod10();
step8 = ((mod8 << 1) & 7UL);
step5 = (mod5 << 1).Mod5();
step3 = (mod3 << 1).Mod3();
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsPrimeCandidate(this ulong n)
{
int i = 0;
ulong p = 0UL;
uint[] smallPrimes = PrimesGenerator.SmallPrimes;
ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2;
int len = smallPrimes.Length;
for (; i < len; i++)
{
if (smallPrimesPow2[i] > n)
{
break;
}
p = smallPrimes[i];
if ((n % p) == 0UL)
{
return n == p;
}
}
return true;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod10(this ulong value) => value % 10UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod8(this ulong value) => value & 7UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod3(this ulong value) => value % 3UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod5(this ulong value) => value % 5UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod6(this ulong value) => value % 6UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod7(this ulong value) => value % 7UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod11(this ulong value) => value % 11UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mod128(this ulong value) => value & 127UL;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong ReduceMod7(uint value)
{
while (value >= 7U)
{
value = (value >> 3) + (value & 7U);
if (value >= 7U && value < 14U)
{
value -= 7U;
}
}
return value;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong ReduceMod11(uint value)
{
while (value >= 11U)
{
value = (value & 15U) + ((value >> 4) * 5U);
if (value >= 11U && value < 22U)
{
value -= 11U;
}
}
return value;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 Mul64(this ulong a, ulong b) => ((UInt128)a.MulHigh(b) << 64) | (UInt128)(a * b); public static UInt128 Mul64(this ulong a, ulong b) => ((UInt128)a.MulHigh(b) << 64) | (UInt128)(a * b);
@@ -252,384 +166,326 @@ public static class ULongExtensions
ulong w1 = xLow * yHigh; ulong w1 = xLow * yHigh;
ulong w2 = xHigh * yLow; ulong w2 = xHigh * yLow;
ulong carry =
(((xLow * yLow) >> 32) +
(uint)w1 +
(uint)w2) >> 32;
return (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32) + carry; return (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32) +
(
((xLow * yLow) >> 32) +
(uint)w1 +
(uint)w2
) >> 32;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong ModPow64(this ulong value, ulong exponent, ulong modulus) public static ulong ModPow64(this ulong value, ulong exponent, ulong modulus)
{ {
ulong result = 1UL; ulong result = 1UL;
ulong baseValue = value % modulus; value %= modulus;
while (exponent != 0UL) while (exponent != 0UL)
{ {
if ((exponent & 1UL) != 0UL) if ((exponent & 1UL) != 0UL)
{ {
result = MulMod64(result, baseValue, modulus); result = MulMod64(result, value, modulus);
} }
baseValue = MulMod64(baseValue, baseValue, modulus); value = MulMod64(value, value, modulus);
exponent >>= 1; exponent >>= 1;
} }
return result; return result;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong MulMod64(this ulong a, ulong b, ulong modulus) public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus);
{
if (modulus <= 1UL)
{
return 0UL;
}
if (TryGetMersenneExponent(modulus, out int exponent))
{
return MulModMersenne(a, b, modulus, exponent);
}
return (ulong)(((UInt128)a * b) % modulus);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool TryGetMersenneExponent(ulong modulus, out int exponent)
{
if (modulus == 0UL)
{
exponent = 0;
return false;
}
ulong plusOne = modulus + 1UL;
if (plusOne == 0UL)
{
exponent = 64;
return true;
}
if ((plusOne & (plusOne - 1UL)) != 0UL)
{
exponent = 0;
return false;
}
exponent = BitOperations.TrailingZeroCount(plusOne);
return true;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong MulModMersenne(ulong a, ulong b, ulong modulus, int exponent)
{
UInt128 mask = (UInt128)modulus;
UInt128 product = (UInt128)a * b;
UInt128 folded = (product & mask) + (product >> exponent);
while (folded > mask)
{
folded = (folded & mask) + (folded >> exponent);
}
ulong result = (ulong)folded;
if (result >= modulus)
{
result -= modulus;
}
return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 PowMod(this ulong exponent, UInt128 modulus) public static UInt128 PowMod(this ulong exponent, UInt128 modulus)
{ {
UInt128 result = UInt128.One; UInt128 result = UInt128.One;
ulong exponentLoopIndex = 0UL; ulong exponentLoopIndex = 0UL;
// Return 1 because 2^0 = 1 // Return 1 because 2^0 = 1
if (exponent == 0UL) if (exponent == 0UL)
return result; return result;
// Any number mod 1 is 0 // Any number mod 1 is 0
if (modulus == UInt128.One) if (modulus == UInt128.One)
return UInt128.Zero; return UInt128.Zero;
// For small exponents, do classic method // For small exponents, do classic method
if (exponent < 64 || modulus < 4) if (exponent < 64 || modulus < 4)
{ {
for (; exponentLoopIndex < exponent; exponentLoopIndex++) for (; exponentLoopIndex < exponent; exponentLoopIndex++)
{ {
result <<= 1; result <<= 1;
if (result >= modulus) if (result >= modulus)
result -= modulus; result -= modulus;
} }
return result; return result;
} }
// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness // Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
if ((modulus & (modulus - 1)) == 0) if ((modulus & (modulus - 1)) == 0)
{ {
result = UInt128.One << (int)(exponent & 127); result = UInt128.One << (int)(exponent & 127);
return result & (modulus - 1); return result & (modulus - 1);
} }
// Reusing exponentLoopIndex to iterate again over the exponent range for the general-case accumulation. // Reusing exponentLoopIndex to iterate again over the exponent range for the general-case accumulation.
exponentLoopIndex = 0UL; exponentLoopIndex = 0UL;
// Reusing result after resetting it for the general modulus accumulation phase. // Reusing result after resetting it for the general modulus accumulation phase.
result = UInt128.One; result = UInt128.One;
for (; exponentLoopIndex < exponent; exponentLoopIndex++) for (; exponentLoopIndex < exponent; exponentLoopIndex++)
{ {
result <<= 1; result <<= 1;
if (result >= modulus) if (result >= modulus)
result -= modulus; result -= modulus;
} }
return result; return result;
} }
/// <summary> /// <summary>
/// Computes 2^exponent mod modulus using a known cycle length. /// Computes 2^exponent mod modulus using a known cycle length.
/// </summary> /// </summary>
public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, ulong cycleLength) public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, ulong cycleLength)
{ {
UInt128 result = UInt128.One; UInt128 result = UInt128.One;
ulong exponentLoopIndex = 0UL; ulong exponentLoopIndex = 0UL;
ulong rotationCount = 0UL;
// Return 1 because 2^0 = 1 // Return 1 because 2^0 = 1
if (exponent == 0UL) if (exponent == 0UL)
return result; return result;
// Any number mod 1 is 0 // Any number mod 1 is 0
if (modulus == UInt128.One) if (modulus == UInt128.One)
return UInt128.Zero; return UInt128.Zero;
// For small exponents, do classic method // For small exponents, do classic method
if (exponent < 64 || modulus < 4) if (exponent < 64 || modulus < 4)
{ {
for (; exponentLoopIndex < exponent; exponentLoopIndex++) for (; exponentLoopIndex < exponent; exponentLoopIndex++)
{ {
result <<= 1; result <<= 1;
if (result >= modulus) if (result >= modulus)
result -= modulus; result -= modulus;
} }
return result; return result;
} }
// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness // Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
if ((modulus & (modulus - 1)) == 0) if ((modulus & (modulus - 1)) == 0)
{ {
result = UInt128.One << (int)(exponent & 127); result = UInt128.One << (int)(exponent & 127);
return result & (modulus - 1); return result & (modulus - 1);
} }
// Reusing exponentLoopIndex to iterate over the rotation count for the cycle-aware path. // Reusing exponentLoopIndex to iterate over the rotation count for the cycle-aware path.
exponentLoopIndex = 0UL; exponentLoopIndex = 0UL;
// Reusing result after resetting it for the rotation-accumulation pass. // Reusing result after resetting it for the rotation-accumulation pass.
result = UInt128.One; result = UInt128.One;
rotationCount = exponent % cycleLength; ulong rotationCount = exponent % cycleLength;
for (; exponentLoopIndex < rotationCount; exponentLoopIndex++) for (; exponentLoopIndex < rotationCount; exponentLoopIndex++)
{ {
result <<= 1; result <<= 1;
if (result >= modulus) if (result >= modulus)
result -= modulus; result -= modulus;
} }
return result; return result;
} }
/// <summary> /// <summary>
/// Computes 2^exponent mod modulus using a known cycle length. /// Computes 2^exponent mod modulus using a known cycle length.
/// </summary> /// </summary>
public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, UInt128 cycleLength) public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, UInt128 cycleLength)
{ {
UInt128 result = UInt128.One; UInt128 result = UInt128.One;
ulong exponentLoopIndex = 0UL; ulong exponentLoopIndex = 0UL;
UInt128 rotationIndex = UInt128.Zero;
UInt128 rotationCount = UInt128.Zero;
// Return 1 because 2^0 = 1 // Return 1 because 2^0 = 1
if (exponent == UInt128.Zero) if (exponent == UInt128.Zero)
return result; return result;
// Any number mod 1 is 0 // Any number mod 1 is 0
if (modulus == UInt128.One) if (modulus == UInt128.One)
return UInt128.Zero; return UInt128.Zero;
// For small exponents, do classic method // For small exponents, do classic method
if (exponent < UInt128Numbers.SixtyFour || modulus < UInt128Numbers.Four) if (exponent < UInt128Numbers.SixtyFour || modulus < UInt128Numbers.Four)
{ {
for (; exponentLoopIndex < exponent; exponentLoopIndex++) for (; exponentLoopIndex < exponent; exponentLoopIndex++)
{ {
result <<= 1; result <<= 1;
if (result >= modulus) if (result >= modulus)
result -= modulus; result -= modulus;
} }
return result; return result;
} }
// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness // Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
if ((modulus & (modulus - UInt128.One)) == UInt128.Zero) if ((modulus & (modulus - UInt128.One)) == UInt128.Zero)
{ {
result = UInt128.One << (int)(exponent & UInt128Numbers.OneHundredTwentySeven); result = UInt128.One << (int)(exponent & UInt128Numbers.OneHundredTwentySeven);
return result & (modulus - 1); return result & (modulus - 1);
} }
// Reusing result after resetting it for the rotation-driven accumulation phase. // Reusing result after resetting it for the rotation-driven accumulation phase.
result = UInt128.One; result = UInt128.One;
rotationCount = exponent % cycleLength; UInt128 rotationCount = exponent % cycleLength;
rotationIndex = UInt128.Zero; UInt128 rotationIndex = UInt128.Zero;
while (rotationIndex < rotationCount) while (rotationIndex < rotationCount)
{ {
result <<= 1; result <<= 1;
if (result >= modulus) if (result >= modulus)
result -= modulus; result -= modulus;
rotationIndex += UInt128.One; rotationIndex += UInt128.One;
} }
return result; return result;
} }
/// <summary> /// <summary>
/// Computes 2^exponent mod modulus using a known cycle length. /// Computes 2^exponent mod modulus using a known cycle length.
/// </summary> /// </summary>
public static UInt128 PowModWithCycle(this UInt128 exponent, UInt128 modulus, UInt128 cycleLength) public static UInt128 PowModWithCycle(this UInt128 exponent, UInt128 modulus, UInt128 cycleLength)
{ {
UInt128 result = UInt128.One; UInt128 one = UInt128.One,
ulong exponentLoopIndex = 0UL; result = one,
UInt128 rotationIndex = UInt128.Zero; zero = UInt128.Zero;
UInt128 rotationCount = UInt128.Zero; ulong exponentLoopIndex = 0UL;
// Return 1 because 2^0 = 1 // Return 1 because 2^0 = 1
if (exponent == UInt128.Zero) if (exponent == zero)
return result; return result;
// Any number mod 1 is 0 // Any number mod 1 is 0
if (modulus == UInt128.One) if (modulus == one)
return UInt128.Zero; return zero;
// For small exponents, do classic method // For small exponents, do classic method
if (exponent < UInt128Numbers.SixtyFour || modulus < UInt128Numbers.Four) if (exponent < UInt128Numbers.SixtyFour || modulus < UInt128Numbers.Four)
{ {
for (; exponentLoopIndex < exponent; exponentLoopIndex++) for (; exponentLoopIndex < exponent; exponentLoopIndex++)
{ {
result <<= 1; result <<= 1;
if (result >= modulus) if (result >= modulus)
result -= modulus; result -= modulus;
} }
return result; return result;
} }
// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness // Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
if ((modulus & (modulus - UInt128.One)) == UInt128.Zero) if ((modulus & (modulus - one)) == zero)
{ {
result = UInt128.One << (int)(exponent & UInt128Numbers.OneHundredTwentySeven); result = one << (int)(exponent & UInt128Numbers.OneHundredTwentySeven);
return result & (modulus - 1); return result & (modulus - 1);
} }
// Reusing result after resetting it for the rotation-driven accumulation phase. // Reusing result after resetting it for the rotation-driven accumulation phase.
result = UInt128.One; result = one;
rotationCount = exponent % cycleLength; UInt128 rotationCount = exponent % cycleLength;
rotationIndex = UInt128.Zero;
while (rotationIndex < rotationCount)
{
result <<= 1;
if (result >= modulus)
result -= modulus;
rotationIndex += UInt128.One; // We're reusing "zero" as rotation index for just a little better performance
} while (zero < rotationCount)
{
result <<= 1;
if (result >= modulus)
result -= modulus;
return result; zero += one;
} }
return result;
}
/// <summary> /// <summary>
/// Computes 2^exponent mod modulus using iterative CRT composition from mod 10 up to modulus. /// Computes 2^exponent mod modulus using iterative CRT composition from mod 10 up to modulus.
/// Only for modulus >= 10 and reasonable size. /// Only for modulus >= 10 and reasonable size.
/// </summary> /// </summary>
public static UInt128 PowModCrt(this ulong exponent, UInt128 modulus, MersenneDivisorCycles cycles) public static UInt128 PowModCrt(this ulong exponent, UInt128 modulus, MersenneDivisorCycles cycles)
{ {
if (modulus < 10) if (modulus < 10)
return PowMod(exponent, modulus); // fallback to classic return PowMod(exponent, modulus); // fallback to classic
// Use cycle length 4 for mod 10 // Use cycle length 4 for mod 10
UInt128 result = PowModWithCycle(exponent, 10, 4); UInt128 currentModulus = 10,
UInt128 currentModulus = 10; cycle,
UInt128 modulusCandidate = 11; modulusCandidate = 11,
UInt128 cycle = UInt128.Zero; remainderForCandidate,
UInt128 remainderForCandidate = UInt128.Zero; result = PowModWithCycle(exponent, 10, 4),
zero = UInt128.Zero;
for (; modulusCandidate <= modulus; modulusCandidate++) for (; modulusCandidate <= modulus; modulusCandidate++)
{ {
cycle = MersenneDivisorCycles.GetCycle(modulusCandidate); cycle = MersenneDivisorCycles.GetCycle(modulusCandidate);
remainderForCandidate = cycle > UInt128.Zero remainderForCandidate = cycle > zero
? PowModWithCycle(exponent, modulusCandidate, cycle) ? PowModWithCycle(exponent, modulusCandidate, cycle)
: PowMod(exponent, modulusCandidate); : PowMod(exponent, modulusCandidate);
// Solve x ≡ result mod currentModulus // Solve x ≡ result mod currentModulus
// x ≡ remM mod m // x ≡ remM mod m
// Find x mod (currentModulus * m) // Find x mod (currentModulus * m)
// Since currentModulus and m are coprime, use CRT: // Since currentModulus and m are coprime, use CRT:
// x = result + currentModulus * t, where t ≡ (remM - result) * inv(currentModulus, m) mod m // x = result + currentModulus * t, where t ≡ (remM - result) * inv(currentModulus, m) mod m
result += currentModulus * ((remainderForCandidate + modulusCandidate - (result % modulusCandidate)) * ModInverse(currentModulus, modulusCandidate) % modulusCandidate); result += currentModulus * ((remainderForCandidate + modulusCandidate - (result % modulusCandidate)) * ModInverse(currentModulus, modulusCandidate) % modulusCandidate);
currentModulus *= modulusCandidate; currentModulus *= modulusCandidate;
if (currentModulus >= modulus) if (currentModulus >= modulus)
break; break;
} }
return result % modulus; return result % modulus;
} }
// Helper: modular inverse (extended Euclidean algorithm) // Helper: modular inverse (extended Euclidean algorithm)
private static UInt128 ModInverse(UInt128 a, UInt128 m) private static UInt128 ModInverse(UInt128 a, UInt128 m)
{ {
UInt128 m0 = m, temp; UInt128 m0 = m,
UInt128 x0 = 0, x1 = 1; originalA,
UInt128 originalA = UInt128.Zero; originalM,
UInt128 originalM = UInt128.Zero; temp,
if (m == 1) x0 = 0,
{ x1 = 1;
return 0;
}
while (a > 1) if (m == 1)
{ {
originalA = a; return 0;
originalM = m; }
m = originalA % originalM;
a = originalM;
temp = x0;
x0 = x1 - (originalA / originalM) * x0;
x1 = temp;
}
if (x1 < 0) while (a > 1)
{ {
x1 += m0; originalA = a;
} originalM = m;
m = originalA % originalM;
a = originalM;
temp = x0;
x0 = x1 - (originalA / originalM) * x0;
x1 = temp;
}
return x1; if (x1 < 0)
} {
x1 += m0;
}
return x1;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool SharesFactorWithExponentMinusOne(this ulong exponent) public static bool SharesFactorWithExponentMinusOne(this ulong exponent)
{ {
ulong value = exponent - 1UL; ulong prime, value = exponent - 1UL;
value >>= BitOperations.TrailingZeroCount(value); value >>= BitOperations.TrailingZeroCount(value);
uint[] smallPrimes = PrimesGenerator.SmallPrimes; uint[] smallPrimes = PrimesGenerator.SmallPrimes;
ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2; ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2;
int i = 1, smallPrimesLength = smallPrimes.Length; int i = 1, smallPrimesLength = smallPrimes.Length;
ulong prime;
for (; i < smallPrimesLength && smallPrimesPow2[i] <= value; i++) for (; i < smallPrimesLength && smallPrimesPow2[i] <= value; i++)
{ {