31-Bit RSA Encryption: u_modPow(b, e, m)

I have become very, VERY annoyed recently with the 'llModPow()' function. Why? Well:

1. It can only use up to 16-bit keys. This offers next to zero security, but more importantly I can't even find a utility to generate 16-bit keys (they are too small). I could have written one myself, but I'd rather take advantage of an existing one.
   
2. It (or the documentation) is BROKEN. See /54/50/42686/1.html#post1297155

So, with a little work I have written a (SLOW) 31-bit implementation. It should accept values in the range [0, 2^31) for the modulus and both exponents. Each 'u_modPow()' operation takes about 45 seconds (so a full encryption/decyption round takes about 90 seconds), so it is far from the most speedy encryption you will find.

I have tested it about 30 times (that's 30 function calls) with random values, but I haven't yet done thorough stress testing (e.g. using values of 2^31-1). You are welcome to do your own testing (and if you do, please post the results!). In successive posts I will provide both a Java program which can be used to verify 'u_modPow()' results, and an LSL test program that uses RSA encryption/decryption with random message integers.

// NOTE
// ----
// All funcitons beginning with the 'u_' prefix assume their arguments are
// non-negative integers in the range [0, 2^31-1].  Passing values outside
// this range results in unspecified behavior.
//
// All results from these functions will either also be in this range or
// will overflow into two integer words.  In the latter case, the least
// significant integer word will be in the above range (the sign bit will be
// zero).  This means that if 'msw' is the most significant integer word and
// 'lsw' is the least significant word, the full number is actually:
//     msw * 2^31 + lsw
// (notice 2^31, NOT 2^32).

integer SIGN_BIT_MASK = 0x80000000;

// Returns the full product of two unsigned integers.  The first (index 0)
// integer is the most significant integer of the product.  The second
// (index 1) integer is the least significant integer of the product.
list u_intProd(integer i1, integer i2)
{
    integer i1Low = i1 & 0xFFFF;
    integer i1High = i1 >> 16;
    integer i2Low = i2 & 0xFFFF;
    integer i2High = i2 >> 16;
    
    //llOwnerSay("DEBUG - intProd: i1Low = "+(string)i1Low);
    //llOwnerSay("DEBUG - intProd: i1High = "+(string)i1High);
    //llOwnerSay("DEBUG - intProd: i2Low = "+(string)i2Low);
    //llOwnerSay("DEBUG - intProd: i2High = "+(string)i2High);
    
    integer lowLow = i1Low*i2Low;
    integer lowHigh = i1Low*i2High;
    integer highLow = i1High*i2Low;
    integer highHigh = i1High*i2High;
    
    integer lowLowMSW = (lowLow >> 31) & 0x1;
    integer lowLowLSW = lowLow & ~SIGN_BIT_MASK;
    integer lowHighMSW = (lowHigh >> 15) & 0xFFFF;
    integer lowHighLSW = (lowHigh << 16) & ~SIGN_BIT_MASK;
    integer highLowMSW = (highLow >> 15) & 0xFFFF;
    integer highLowLSW = (highLow << 16) & ~SIGN_BIT_MASK;
    
    //llOwnerSay("DEBUG - intProd: lowLowMSW = "+(string)lowLowMSW);
    //llOwnerSay("DEBUG - intProd: lowLowLSW = "+(string)lowLowLSW);
    //llOwnerSay("DEBUG - intProd: lowHighMSW = "+(string)lowHighMSW);
    //llOwnerSay("DEBUG - intProd: lowHighLSW = "+(string)lowHighLSW);
    //llOwnerSay("DEBUG - intProd: highLowMSW = "+(string)highLowMSW);
    //llOwnerSay("DEBUG - intProd: highLowLSW = "+(string)highLowLSW);
    
    integer msw = lowLowMSW+(highHigh<<1)+lowHighMSW+highLowMSW;
    
    integer accum0 = lowLowLSW+lowHighLSW;
    if (accum0 & SIGN_BIT_MASK)
    {
        accum0 = accum0 & ~SIGN_BIT_MASK;
        ++msw;
    }
    
    integer lsw = accum0+highLowLSW;
    if (lsw & SIGN_BIT_MASK)
    {
        lsw = lsw & ~SIGN_BIT_MASK;
        ++msw;
    }
    
    //llOwnerSay("DEBUG - i1 = "+(string)i1);
    //llOwnerSay("DEBUG - i2 = "+(string)i2);
    //llOwnerSay("DEBUG - msw = "+(string)msw);
    //llOwnerSay("DEBUG - lsw = "+(string)lsw);
    
    return [ msw, lsw ];
}

// Returns the sum of a big number (product of integers) and a little number
// (integer).  The first (index 0) integer is the most significant integer
// of the sum, and the second (index 1) integer is the least significant
// integer of the sum.
list u_bigSum(integer msw1, integer lsw1, integer i2)
{
    integer msw = msw1;
    integer lsw = lsw1+i2;
    
    if (lsw & SIGN_BIT_MASK)
    {
        lsw = lsw & ~SIGN_BIT_MASK;
        ++msw;
    }
    
    return [ msw, lsw ];
}

// Returns the modulo of a large dividend (product of two integers) and a
// small divisor (integer).
integer u_bigMod(integer dividendMSW, integer dividendLSW, integer divisor)
{
    integer msw = dividendMSW%divisor;
    integer lsw = dividendLSW%divisor;
    
    //llOwnerSay("DEBUG - bigMod: dividendMSW = "+(string)dividendMSW);
    //llOwnerSay("DEBUG - bigMod: dividendLSW = "+(string)dividendLSW);
    //llOwnerSay("DEBUG - bigMod: divisor = "+(string)divisor);
    
    // This might be slow, but it WILL finish.  Proof (assume msw > 0):
    //     x = msw * 2^31 + lsw
    //     x % y = (msw * 2^31) % y + lsw % y
    //         = (msw % y) * (2^31 % y) + lsw % y
    //     msw % y <= msw
    //     lsw % y <= lsw
    //     y <= 2^31 - 1
    //     2^31 % y < 2^31 - 1
    // thus:
    //     (msw % y) * (2^31 % y) + lsw % y < msw * 2^31 + lsw
    while (msw != 0)
    {
        //llOwnerSay("DEBUG - msw = "+(string)msw);
        
        integer wsmod = (0x7FFFFFFF%divisor)+1;
        if (wsmod == divisor)  // Added correction Feb. 8, 2008
        {
            wsmod = 0;
        }
        
        list prod = u_intProd(wsmod, msw);
        integer prodMSW = llList2Integer(prod, 0);
        integer prodLSW = llList2Integer(prod, 1);
        
        list sum = u_bigSum(prodMSW, prodLSW, lsw);
        msw = llList2Integer(sum, 0)%divisor;
        lsw = llList2Integer(sum, 1)%divisor;
        
        //llOwnerSay("DEBUG - bigMod: msw = "+(string)msw);
        //llOwnerSay("DEBUG - bigMod: lsw = "+(string)lsw);
    }
    
    return lsw;
}

// Returns (a*b)%c.
integer u_modProd(integer a, integer b, integer c)
{
    list prod = u_intProd(a, b);
    integer prodMSW = llList2Integer(prod, 0);
    integer prodLSW = llList2Integer(prod, 1);
    
    integer result =  u_bigMod(prodMSW, prodLSW, c);
    
    //llOwnerSay("DEBUG modProd - a = "+(string)a);
    //llOwnerSay("DEBUG modProd - b = "+(string)b);
    //llOwnerSay("DEBUG modProd - c = "+(string)c);
    //llOwnerSay("DEBUG modProd - prodMSW = "+(string)prodMSW);
    //llOwnerSay("DEBUG modProd - prodLSW = "+(string)prodLSW);
    //llOwnerSay("DEBUG modProd - result = "+(string)result);
    
    return result;
}

// Returns (a^b)%c.
integer u_modPow(integer a, integer b, integer c)
{
    integer result = 1;
    integer currentPow = a;
    integer exp = b;
    
    while (exp != 0)
    {
        //llOwnerSay("DEBUG - modPow: exp = "+(string)exp);
        
        if (exp & 0x1)
        {
            result = u_modProd(result, currentPow, c);
        }
        
        //llOwnerSay("DEBUG - modPow: currentPow = "+(string)currentPow);
        //llOwnerSay("DEBUG - modPow: result = "+(string)result);
        
        currentPow = u_modProd(currentPow, currentPow, c);
        exp = exp >> 1;
    }
    
    return result;
}

default
{
    touch_start(integer total_number)
    {
        integer N_TESTS = 10;
        
        float startTime = llGetTime();
        
        integer i;
        for (i = 0; i < N_TESTS; ++i)
        {
            integer a = (integer)llFloor(llFrand(0x7FFFFFFF));
            integer b = (integer)llFloor(llFrand(0x7FFFFFFF));
            integer c = (integer)llFloor(llFrand(0x7FFFFFFF));
            
            llOwnerSay("a = "+(string)a);
            llOwnerSay("b = "+(string)b);
            llOwnerSay("c = "+(string)c);
            
            integer result = u_modPow(a, b, c);
            
            llOwnerSay("(a^b)%c = "+(string)result);
        }
        
        float endTime = llGetTime();
        
        float timeDiff = endTime - startTime;
        float avgTime = timeDiff / N_TESTS;
        
        llOwnerSay("Average execution time: "+(string)avgTime+" s");
    }
}

Here's some test code (replace the contents of the above 'touch_start' event handler with it or whatever) that does RSA encryption/decryption on random 31-bit values (non-negative and less than the modulus) to verify the functionality (again not stress-tested).

NOTE: Feel free to test as is or with different RSA keys, and MAKE SURE TO CHANGE THIS KEY IF YOU USE THIS IN A REAL APPLICATION! There are several utilities out there that can generate RSA keys. I used plainrsa-gen in Linux from the ipsec-tools package as shown below.

// Generated with plainrsa-gen -e 3 -b 31
integer RSA_MOD = 0x6d801c75;
integer RSA_PUB = 0x03;
integer RSA_PRI = 0x48ff3033;

integer N_TESTS = 10;

float startTime = llGetTime();

llOwnerSay("Starting Test");

integer i;
for (i = 0; i < N_TESTS; ++i)
{
    // Feb. 8, 2008 - Note that this is actually a signing/authentication example.  For encryption, you'd do the reverse: encrypt using the public key and decrypt using the private key.  In terms of testing performance and accuracy, it is pretty symmetrical though.
    integer clear = (integer)llFloor(llFrand(RSA_MOD));
    integer encrypt = u_modPow(clear, RSA_PRI, RSA_MOD);
    integer decrypt = u_modPow(encrypt, RSA_PUB, RSA_MOD);

    if (decrypt != clear)
    {
        llOwnerSay("Failed test "+(string)(i+1)+"/"+(string)N_TESTS);
        llOwnerSay("   clear = "+(string)clear);
        llOwnerSay("   encrypt = "+(string)encrypt);
        llOwnerSay("   decrypt = "+(string)decrypt);

        return;
    }

    llOwnerSay("Passed test "+(string)(i+1)+"/"+(string)N_TESTS);
}

float endTime = llGetTime();

float timeDiff = endTime - startTime;
float avgTime = timeDiff / N_TESTS;

llOwnerSay(
    "Average execution time of full encrypt/decrypt round: "+(string)avgTime+" s");

Here's a Java application for verifying the results of 'u_modPow()'. After compiling with:

you invoke it like:

to get the results of (a^b)%c. The numbers can be arbitrarily sized.

import java.math.BigInteger;

class ModPow
{
   public static void main(String[] args)
   {
      BigInteger i1 = new BigInteger(args[0]);
      BigInteger i2 = new BigInteger(args[1]);
      BigInteger i3 = new BigInteger(args[2]);

      BigInteger res = i1.modPow(i2, i3);

      System.out.println(
         "("+i1+" ^ "+i2+") % "+i3+" = "+res);
   }
}

I tested this in Mono a few weeks ago, and it completes in about a second, unlike the 45 or so seconds it takes in the pre-Mono LSL VM.