magicdiv.cs

using System;
using System.Diagnostics;
using System.Runtime.CompilerServices;

class Program
{
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    static int MULHI(int x, int y) => (int)(((long)x * (long)y) >> 32);

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    static int RSH(int x, int c) => x >> c;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    static int RSZ(int x, int c) => (int)((uint)x >> c);

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    static int MagicDiv(int n, int m, int s)
    {
        int r = MULHI(n, m);

        if (m < 0)
            r += n;

        return RSH(r, s) + RSZ(r, 31);
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    static int Div(int n, int d)
    {
        return n / d;
    }

    static void Test(int d)
    {
        int s;
        int m = GetSignedMagicNumberForDivide(d, out s);
        Console.WriteLine($"denominator:{d} magic:{m:X8} shift:{s}");

        int sum = 0;
        var sw = Stopwatch.StartNew();

        for (int n = 0; n < int.MaxValue; n++)
        {
            //sum += MagicDiv(n, m, s);
            //sum += Div(n, d);

            if (MagicDiv(n, m, s) != Div(n, d))
              Console.WriteLine($"bad for {n} -> {MagicDiv(n, m, s)} {Div(n, d)}");
        }

        Console.WriteLine($"checksum:{sum} time:{sw.ElapsedMilliseconds}");
    }

    static void Main()
    {
        for (int i = 2; i < int.MaxValue; i++)
            Test(i);
    }

    // From RyuJIT:
    // code to generate a magic number and shift amount for the magic number division
    // optimization.  This code is previously from UTC where it notes it was taken from
    // _The_PowerPC_Compiler_Writer's_Guide_, pages 57-58.
    // The paper it is based on is "Division by invariant integers using multiplication"
    // by Torbjorn Granlund and Peter L. Montgomery in PLDI 94
    
    static int GetSignedMagicNumberForDivide(int denom, out int shift)
    {
        const int bits = sizeof(int) * 8;
        const int bits_minus_1 = bits - 1;

        const uint two_nminus1 = 1U << bits_minus_1;

        int p;
        uint absDenom;
        uint absNc;
        uint delta;
        uint q1;
        uint r1;
        uint r2;
        uint q2;
        uint t;
        int result_magic;
        int result_shift;
        int iters = 0;

        absDenom = (uint)Math.Abs(denom);
        t = two_nminus1 + ((uint)denom >> 31);
        absNc = t - 1 - (t % absDenom);        // absolute value of nc
        p = bits_minus_1;                  // initialize p
        q1 = two_nminus1 / absNc;           // initialize q1 = 2^p / abs(nc)
        r1 = two_nminus1 - (q1 * absNc);    // initialize r1 = rem(2^p, abs(nc))
        q2 = two_nminus1 / absDenom;        // initialize q1 = 2^p / abs(denom)
        r2 = two_nminus1 - (q2 * absDenom); // initialize r1 = rem(2^p, abs(denom))

        do
        {
            iters++;
            p++;
            q1 *= 2; // update q1 = 2^p / abs(nc)
            r1 *= 2; // update r1 = rem(2^p / abs(nc))

            if (r1 >= absNc)
            { // must be unsigned comparison
                q1++;
                r1 -= absNc;
            }

            q2 *= 2; // update q2 = 2^p / abs(denom)
            r2 *= 2; // update r2 = rem(2^p / abs(denom))

            if (r2 >= absDenom)
            { // must be unsigned comparison
                q2++;
                r2 -= absDenom;
            }

            delta = absDenom - r2;
        } while (q1 < delta || (q1 == delta && r1 == 0));

        result_magic = (int)(q2 + 1); // resulting magic number
        if (denom < 0)
        {
            result_magic = -result_magic;
        }
        shift = p - bits; // resulting shift

        return result_magic;
    }
}