even-perfect-numbers-scanner/EvenPerfectBitScanner.Benchmarks/MontgomeryMultiplyBenchmarks.cs

using System;
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs;

namespace EvenPerfectBitScanner.Benchmarks;

[MemoryDiagnoser]
[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 10)]
public class MontgomeryMultiplyBenchmarks
{
    private const int MaxBatchSize = 1024;

    private ulong[] _operandsA = Array.Empty<ulong>();
    private ulong[] _operandsB = Array.Empty<ulong>();
    private ulong[] _moduli = Array.Empty<ulong>();
    private ulong[] _nPrimes = Array.Empty<ulong>();

    [Params(64, 256, 1024)]
    public int BatchSize { get; set; }

    [GlobalSetup]
    public void Setup()
    {
        MontgomerySampleData data = MontgomerySampleDataCache.Instance;
        _operandsA = data.OperandsA;
        _operandsB = data.OperandsB;
        _moduli = data.Moduli;
        _nPrimes = data.NPrimes;
    }

    private sealed class MontgomerySampleData
    {
        public MontgomerySampleData(ulong[] operandsA, ulong[] operandsB, ulong[] moduli, ulong[] nPrimes)
        {
            OperandsA = operandsA;
            OperandsB = operandsB;
            Moduli = moduli;
            NPrimes = nPrimes;
        }

        public ulong[] OperandsA { get; }
        public ulong[] OperandsB { get; }
        public ulong[] Moduli { get; }
        public ulong[] NPrimes { get; }
    }

    private static class MontgomerySampleDataCache
    {
        private static readonly Lazy<MontgomerySampleData> Cache = new(Create);

        public static MontgomerySampleData Instance => Cache.Value;

        private static MontgomerySampleData Create()
        {
            Console.WriteLine("Preparing Montgomery multiply benchmark sample data...");

            Random random = new(1979);
            var operandsA = new ulong[MaxBatchSize];
            var operandsB = new ulong[MaxBatchSize];
            var moduli = new ulong[MaxBatchSize];
            var nPrimes = new ulong[MaxBatchSize];

            for (int i = 0; i < MaxBatchSize; i++)
            {
                ulong modulus = NextOddModulus(random);
                moduli[i] = modulus;
                nPrimes[i] = ComputeMontgomeryNPrime(modulus);
                operandsA[i] = (ulong)random.NextInt64(0, (long)modulus);
                operandsB[i] = (ulong)random.NextInt64(0, (long)modulus);
            }

            MontgomerySampleData data = new MontgomerySampleData(operandsA, operandsB, moduli, nPrimes);

            Console.WriteLine("Finished preparing Montgomery multiply benchmark sample data.");

            return data;
        }

        private static ulong NextOddModulus(Random random)
        {
            ulong modulus = (ulong)random.NextInt64(3, long.MaxValue);
            if ((modulus & 1UL) == 0UL)
            {
                modulus++;
            }

            return modulus;
        }
    }

    /// <summary>
    /// Baseline Montgomery multiply used in earlier revisions; measured 247 ns (64 batch), 974 ns (256), and 3.94 μs (1024).
    /// </summary>
    /// <remarks>
    /// Observed means: BatchSize 64 → 247.0 ns (1.00×), 256 → 974.2 ns, 1024 → 3,940.6 ns.
    /// </remarks>
    [Benchmark(Baseline = true)]
    public ulong OriginalImplementation()
    {
        ulong checksum = 0UL;
        for (int i = 0; i < BatchSize; i++)
        {
            checksum ^= OriginalMontgomeryMultiply(_operandsA[i], _operandsB[i], _moduli[i], _nPrimes[i]);
        }

        return checksum;
    }

    /// <summary>
    /// Optimized Montgomery multiply; consistently 3–5% faster with 234.9 ns (64), 942.6 ns (256), and 3.75 μs (1024).
    /// </summary>
    /// <remarks>
    /// Observed means: BatchSize 64 → 234.9 ns (0.95×), 256 → 942.6 ns, 1024 → 3,751.6 ns.
    /// </remarks>
    [Benchmark]
    public ulong OptimizedImplementation()
    {
        ulong checksum = 0UL;
        for (int i = 0; i < BatchSize; i++)
        {
            checksum ^= OptimizedMontgomeryMultiply(_operandsA[i], _operandsB[i], _moduli[i], _nPrimes[i]);
        }

        return checksum;
    }

    private static ulong OriginalMontgomeryMultiply(ulong a, ulong b, ulong modulus, ulong nPrime)
    {
        ulong tLow = unchecked(a * b);
        ulong tHigh = MultiplyHigh(a, b);
        ulong m = unchecked(tLow * nPrime);
        ulong mnLow = unchecked(m * modulus);
        ulong mnHigh = MultiplyHigh(m, modulus);

        ulong sumLow = unchecked(tLow + mnLow);
        ulong carry = sumLow < tLow ? 1UL : 0UL;
        ulong sumHigh = unchecked(tHigh + mnHigh + carry);

        ulong result = sumHigh;
        if (result >= modulus)
        {
            result -= modulus;
        }

        return result;
    }

    private static ulong OptimizedMontgomeryMultiply(ulong a, ulong b, ulong modulus, ulong nPrime)
    {
        ulong tLow = unchecked(a * b);
        ulong m = unchecked(tLow * nPrime);
        ulong sumLow = unchecked(tLow + m * modulus);
        ulong carry = sumLow < tLow ? 1UL : 0UL;
        ulong sumHigh = unchecked(MultiplyHigh(a, b) + MultiplyHigh(m, modulus) + carry);

        ulong result = sumHigh;
        if (result >= modulus)
        {
            result -= modulus;
        }

        return result;
    }

    private static ulong ComputeMontgomeryNPrime(ulong modulus)
    {
        ulong inv = modulus;
        inv *= unchecked(2UL - modulus * inv);
        inv *= unchecked(2UL - modulus * inv);
        inv *= unchecked(2UL - modulus * inv);
        inv *= unchecked(2UL - modulus * inv);
        inv *= unchecked(2UL - modulus * inv);
        inv *= unchecked(2UL - modulus * inv);
        return unchecked(0UL - inv);
    }

    private static ulong MultiplyHigh(ulong x, ulong y)
    {
        ulong xLow = (uint)x;
        ulong xHigh = x >> 32;
        ulong yLow = (uint)y;
        ulong yHigh = y >> 32;

        ulong lowProduct = xLow * yLow;
        ulong mid = xHigh * yLow + (lowProduct >> 32);
        ulong midLow = mid & 0xFFFFFFFFUL;
        ulong midHigh = mid >> 32;
        midLow += xLow * yHigh;

        ulong resultHigh = xHigh * yHigh + midHigh + (midLow >> 32);
        return resultHigh;
    }
}