even-perfect-numbers-scanner/PerfectNumbers.Core/Gpu/Kernels/DivisorCycleKernels.cs

using System.Runtime.CompilerServices;
using ILGPU;
using ILGPU.Runtime;

namespace PerfectNumbers.Core.Gpu;

internal static class DivisorCycleKernels
{
    private const byte ByteZero = 0;
    private const byte ByteOne = 1;

    // GPU kernel for divisor cycle calculation
    public static void GpuDivisorCycleKernel(
        Index1D index,
        ArrayView1D<ulong, Stride1D.Dense> divisors,
        ArrayView1D<ulong, Stride1D.Dense> outCycles)
    {
        int i = index.X;
        outCycles[i] = CalculateCycleLengthGpu(divisors[i]);
    }

    // GPU-friendly version of cycle length calculation
    /// <summary>
    /// GPU-friendly cycle calculator that unrolls sixteen doubling steps; it wins the 8,388,607 benchmark and stays within ~2%
    /// of the octo loop at divisor 131,071.
    /// </summary>
    public static ulong CalculateCycleLengthGpu(ulong divisor)
    {
        if ((divisor & (divisor - 1UL)) == 0UL)
            return 1UL;

        ulong order = 1UL;
        ulong pow = 2UL;

        while (true)
        {
            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;

            if (GpuStep(ref pow, divisor, ref order))
                return order;
        }
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static bool GpuStep(ref ulong pow, ulong divisor, ref ulong order)
    {
        pow += pow;
        if (pow >= divisor)
            pow -= divisor;

        order++;
        return pow == 1UL;
    }

    public static void GpuAdvanceDivisorCyclesKernel(
        Index1D index,
        int steps,
        ArrayView1D<ulong, Stride1D.Dense> divisors,
        ArrayView1D<ulong, Stride1D.Dense> pow,
        ArrayView1D<ulong, Stride1D.Dense> order,
        ArrayView1D<ulong, Stride1D.Dense> cycles,
        ArrayView1D<byte, Stride1D.Dense> status)
    {
        if (status[index] != ByteZero)
        {
            return;
        }

        ulong divisor = divisors[index];

        if ((divisor & (divisor - 1UL)) == 0UL)
        {
            cycles[index] = 1UL;
            status[index] = ByteOne;
            pow[index] = 1UL;
            order[index] = 1UL;
            return;
        }

        ulong currentPow = pow[index];
        ulong currentOrder = order[index];

        if (divisor <= 3UL)
        {
            while (currentPow != 1UL)
            {
                currentPow += currentPow;
                if (currentPow >= divisor)
                {
                    currentPow -= divisor;
                }

                currentOrder++;
            }

            cycles[index] = currentOrder;
            status[index] = ByteOne;
            pow[index] = currentPow;
            order[index] = currentOrder;
            return;
        }

        do
        {
            currentPow += currentPow;
            if (currentPow >= divisor)
            {
                currentPow -= divisor;
            }

            currentOrder++;

            if (currentPow == 1UL)
            {
                cycles[index] = currentOrder;
                status[index] = ByteOne;
                pow[index] = currentPow;
                order[index] = currentOrder;
                return;
            }
        }
        while (--steps != 0);

        pow[index] = currentPow;
        order[index] = currentOrder;
    }
}