Computing (a^x)mod n using addition chaining. Algorithm in C++

Question

unsigned long qe2(unsigned long x, unsigned long y , unsigned long n)
{
    unsigned long s,t,u;

    s=1; t=x; u=y;

    while(u)
    {
        if(u&1)
            s = (s*t)%n;
        u>>=1;
        t= ( t*t)%n;
    }
    return s;
}

I was readin up on cryptography and , in the book Applied Cryptography by Bruice Schneier , I found the above algorithm to calculate (x^y)mod n at a low computational cost. It basically uses an algorithm called addition chaining to decrease the amount of multiplications in the calculation. I am able to use the process on pen and paper , but despite reading the above code again and again , I have not been able to understand how it works. So please point me in the right direction ( some kind of link to an article ) where I can analyse the above algorithm, or it will be really helpful if you could give the explanation here.

P.S. : the code is not explained in the book.

Answer 1

The exponent y is written as a sum of powers of two, i.e. in binary.

Consider a practical example: (x**11) % M. Mathematically,

(x**11) % M == ((((((x**1) % M) * x**2) % M) * x**8) % M)

This is useful, because a simple loop is sufficient to calculate the powers of x that are powers of two. For example, if you want to calculate x**(2**i):

for (j = 0; j < i; j++)
    x = x * x;

We can examine the logic in detail, if we look at a function that calculates basis**exponent:

unsigned long power(unsigned long basis,
                    unsigned long exponent)
{
    unsigned long result = 1u;
    unsigned long factor = basis;
    unsigned long i;

    for (i = 1u; i < exponent; i = i * 2u) {
        /* i is some power of two,
           and factor = basis**i.
        */

        /* If i is in the sum (of powers of two) for exponent,
           we include the corresponding power of basis
           in the product. */
        if (i & exponent)
            result = result * factor;

        /* Update factor for the next i. */
        factor = factor * factor;
    }

    return result;
}

Note that the test (i & exponent) is a binary AND operation, which is true if the result is nonzero, false if zero. Because i is a power of two, in binary it has a single 1 with all other binary digits zero, so it essentially tests if the exponent, written in binary, has a 1 in that position.

OP's code is simpler, because instead of using separate variables i and factor, it shifts the exponent right (dropping the rightmost binary digit), and uses the basis itself. That is,

unsigned long power(unsigned long basis,
                    unsigned long exponent)
{
    unsigned long result = 1u;

    while (exponent > 0u) {
        if (exponent & 1u)
            result = result * basis;
        basis = basis * basis;
        exponent = exponent >> 1;
    }

    return result;
}

The modulo operation is the final wrinkle, but it too is trivial. If you compute a product modulo some positive integer, you can apply the modulo operator to each term, and to each temporary result, without affecting the result:

(a * b * c * d * ... * z) % m
= ( (a % m) * (b % m) * (c % m) * ... * (z % m) ) % m
= ( ... (((((a % m) * (b % m)) % m) * (c % m)) % m) * ... * (z % m)) % m

This is obviously useful, because you can calculate the product of any number of terms, with all terms smaller than m, and all intermediate terms smaller than m*m.

Time complexity of this kind of exponentiation is O(log N), where N is the exponent.