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.

Probably the best way to see how that works, is you step through the code line by line with the debugger. Watch then how the variable values are changing at each step. — user0042, Jul 16 '17 at 13:37
"addition chain" is the general term used to describe this class of algorithm. However, this is a simplified special case known more commonly as the [right-to-left binary exponentiation algorithm](https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method), and you can find a lot of documentation using that search term. — President James K. Polk, Jul 16 '17 at 14:05
You should ask a specific question for a particular problem. Since Stack Overflow hides the Close reason from you: *"Questions asking us to recommend or find a book, tool, software library, tutorial or other off-site resource are off-topic for Stack Overflow as they tend to attract opinionated answers and spam. Instead, describe the problem and what has been done so far to solve it."* — jww, Jul 16 '17 at 16:50

Nominal Animal · Accepted Answer · 2017-07-19T18:37:40.627

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.

this helped me thanks , But i think that there's a mistake in the code : in the second example of the function power() the result is initially zero ( i think it should be 1) .. and it remains zero till the end , and so the function will return zero always. — Ratul Thakur, Jul 19 '17 at 12:38
@RatulThakur: Both `power()` functions should start with `result = 1u`, of course. Fixed. Good catch; thank you! — Nominal Animal, Jul 19 '17 at 18:41

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

1 Answers1