How to implement c=m^e mod n for enormous numbers?

Question

I'm trying to figure out how to implement RSA crypto from scratch (just for the intellectual exercise), and i'm stuck on this point:

For encryption, c = m^e mod n

Now, e is normally 65537. m and n are 1024-bit integers (eg 128-byte arrays). This is obviously too big for standard methods. How would you implement this?

I've been reading a bit about exponentiation here but it just isn't clicking for me:

Wikipedia-Exponentiation by squaring

This Chapter (see section 14.85)

Thanks.

edit: Also found this - is this more what i should be looking at? Wikipedia- Modular Exponentiation

Answer 1

Exponentiation by squaring:

Let's take an example. You want to find 17²³. Note that 23 is 10111 in binary. Let's try to build it up from left to right.

           // a      exponent in binary

a = 17     //17^1          1

a = a * a  //17^2         10

a = a * a  //17^4        100
a = a * 17 //17^5        101

a = a * a  //17^10      1010
a = a * 17 //17^11      1011

a = a * a  //17^22     10110
a = a * 17 //17^23     10111

When you square, you double the exponent (shift left by 1 bit). When you multiply by m, you add 1 to the exponent.

If you want to reduce modulo n, you can do it after each multiplication (rather than leaving it to the end, which would make the numbers get very large).

65537 is 10000000000000001 in binary which makes all of this pretty easy. It's basically

a = m
repeat 16 times:
    a = a * a
    a = a mod n
a = a * m
a = a mod n

where of course a, n and m are "big integers". a needs to be at least 2048 bits as it can get as large as (n-1)².

Answer 2

For an efficient algorithm you need to combine the exponentiation by squaring with repeated application of mod after each step.

For odd e this holds:

m^e mod n = m ⋅ m^e-1 mod n

For even e:

m^e mod n = (m^e/2 mod n)² mod n

With m¹ = m as a base case this defines a recursive way to do efficient modular exponentiation.

But even with an algorithm like this, because m and n will be very large, you will still need to use a type/library that can handle integers of such sizes.

Answer 3

result = 1
while e>0:
  if (e & 1) != 0:
    result = result * m
    result = result mod n
  m = m*m
  m = m mod n
  e = e>>1
return result

This checks bits in the exponent starting with the least significant bit. Each time we move up a bit it corresponds to doubling the power of m - hence we shift e and square m. The result only gets the power of m multiplied in if the exponent has a 1 bit in that position. All multiplications need to be reduced mod n.

As an example, consider m^13. 11 = 1101 in binary. so this is the same as m^8 * m^4 * m. Notice the powers 8,4,(not 2),1 which is the same as the bits 1101. And then recall that m^8 = (m^4)^2 and m^4 = (m^2)^2.