ModPow in LSL

// modular reduction: compute ab mod p
// assumes 0 ≤ a < p, 0 ≤ b < p, p > 2
integer _reduce(integer a, integer b, integer p) {
    /* calculate ab using long multiplication of 16bit limbs */

    // split a, b into 16bit limbs
    integer al = a & 0xFFFF;
    integer ah = (a>>16) & 0xFFFF;
    integer bl = b & 0xFFFF;
    integer bh = (b>>16) & 0xFFFF;

    // intermediate limb products
    integer p0 = al * bl;
    integer p1 = ah * bl;
    integer p2 = al * bh;

    // calculate upper limb of the 64bit product's lower 32 bits.
    // any overflow will be carry to add to the upper 32bits of he 64bit product
    a = ((p0 >> 16) & 0xFFFF) + (p1 & 0xFFFF) + (p2 & 0xFFFF);

    // 64bit product ab, split across 2 32bit integers (store in ah, al)
    al = (p0 & 0xFFFF) | (a << 16);
    ah = (a >> 16) + ((p1 >> 16) & 0xFFFF) + ((p2 >> 16) & 0xFFFF) + ah * bh;

    // if we just have a 31 bit product, can just use builtin operator
    if(!ah && (al > 0)) return al % p;

    /* now divide by p to calculate remainder */

    // determine number of bits set (store in a)
    if(ah) {
        // bitcount is 32 + number of bits in ah
        a = 32;
        b = ah;
    } else if(al & 0x80000000) a = 32; // msb set, so bitcount is 32
    else b = al; // bitcount is number of bits in al
    if(b) do ++a; while(b = b >> 1);

    integer r; // remainder R to calculate

    // 64 bit mask M (store in p0,p1)
    // Used to select the dividend bit to bring down
    p0 = 0; p1 = 0;
    if(a > 32)
        p0 = 1 << (a - 33);
    else
        p1 = 1 << (a - 1);

    // long binary division, ignoring quotient
    // starting with the MSB of the dividend, pull
    // it down to build the remainder. Once the
    // remainder exceeds the divisor, subtract the
    // divisor to get the next remainder to continue
    // building with pulled down dividend bits.
    // Once all the dividend bits have ben pulled
    // down you have the resulting remainder
    for(; a; --a) {
        // pull down the dividend bit and build the
        // remainder
        if(p0)
            r = (r << 1) | ((ah & p0) != 0);
        else
            r = (r << 1) | ((al & p1) != 0);

        if((r >= p) || (r < 0))
            // remainder exceeded divisor,
            // so subtract divisor
            r -= p;

        // next dividend bit selector (M = M >> 1)
        if(p1) p1 = (p1 >> 1) & 0x7FFFFFFF;
        else if(!(p0 = p0 >> 1)) p1 = 0x80000000;
    }

    return r;
}

// compute xⁿ mod p
// uses exponentiation by squaring method and
// utilises fact that xy mod n = (x mod n)(y mod n) mod n
integer ModPow(integer x, integer n, integer p){
    // don't support -ve powers
    if(n < 0) return 0;
    // modulo less than 2 makes no sense
    if(p < 2) return 0;

    // if modulo fits within 16bits, we can trivially
    // use * and % operators, no need to use the
    // expensive _reduce UDF
    integer t = p < 0x10000;

    integer r = 1; // result

    // trivial reduction
    if(x >= p) x = x % p;

    if(n) do {
        if(n&1)
        if(r>1) {
            // bit set, so current square mod p
            //  is a factor of the remainder
            if(t) r = (r * x) % p;
            else r = _reduce(r, x, p);
        } else r = x;

        // calculate next square mod p
        if(n = n >> 1)
        if(t) x = (x * x) % p;
        else x = _reduce(x, x, p);
    } while(n);

    return r;
}

default
{
    state_entry()
    {
        integer x; integer n; integer p;
        integer e; integer r;
        list d = [
            // X. N, P, Expected
            // Expected taken from https://www.dcode.fr/modular-exponentiation
            2386073, 2, 2147483647, 365213132,
            54023227, 188543765, 128767565, 72911872,
            105840436, 187237174, 131390926, 37550416,
            41160897, 165236876, 146705624, 57948809,
            61936308, 110630338 , 65008019, 11358659,
            80557, 2111897, 7177, 2959
            ];
        integer l = llGetListLength(d) - 1;
        float t;

        integer i = -1;
        llResetTime();
        do {
            x = llList2Integer(d, ++i);
            n = llList2Integer(d, ++i);
            p = llList2Integer(d, ++i);
            e = llList2Integer(d, ++i);
            r = ModPow(x, n, p);
            llOwnerSay((string)[
                "ModPow(x=", x, "; n=", n, "; p=", p,
                ")\n=> expected=", e, "; actual=", r, " ",
                llGetSubString("✘✔", r==e, r==e)
                ]);
        } while (i < l);
        t = llGetTime();
        llOwnerSay((string)[ "time taken: ", t, " seconds" ]);

        i = -1;
        llResetTime();
        do {
            x = llList2Integer(d, ++i);
            n = llList2Integer(d, ++i);
            p = llList2Integer(d, ++i);
            e = llList2Integer(d, ++i);
            r = llModPow(x, n, p);
            llOwnerSay((string)[
                "llModPow(x=", x, "; n=", n, "; p=", p,
                ")\n=> expected=", e, "; actual=", r, " ",
                llGetSubString("✘✔", r==e, r==e)
                ]);
        } while (i < l);
        t = llGetTime();
        llOwnerSay((string)[ "time taken: ", t, " seconds" ]);
    }
}