logo

Square-and-Multiply Algorithm Proof 📂Number Theory

Square-and-Multiply Algorithm Proof

Algorithm

Given a natural number a,k,ma,k,m, bak(modm)b \equiv a^{k} \pmod{m} can be computed as follows.


Step 1. Binary Expansion of kk

Represent ui=0u_{i} = 0 or ui=1u_{i} = 1 as follows. k=i=0rui2i=u0+2u1++2rur k = \sum_{i=0}^{r} u_{i} 2^{i} = u_{0} + 2 u_{1} + \cdots + 2^r u_{r}


Step 2.

Calculate a2r(a2r1)2Ar12Ar(modm)a^{2^{r}} \equiv ( a^{2^{r-1}} )^2 \equiv A_{r-1}^2 \equiv A_{r} \pmod{m} as follows: aA0(modm)a2(a1)2A02A1(modm)a4(a2)2A12A2(modm)a8(a4)2A22A3(modm) \begin{align*} a & & & \equiv A_{0} \pmod{m} \\ a^{2} & \equiv ( a^1 )^2 & \equiv A_{0}^2 & \equiv A_{1} \pmod{m} \\ a^{4} & \equiv ( a^2 )^2 & \equiv A_{1}^2 & \equiv A_{2} \pmod{m} \\ a^{8} & \equiv ( a^4 )^2 & \equiv A_{2}^2 & \equiv A_{3} \pmod{m} \\ & & & \vdots \end{align*}


Step 3.

Calculating A0u0A1u1ArurA_{0}^{u_{0}} \cdot A_{1}^{u_{1}} \cdot \cdots A_{r}^{u_{r}} yields: akA0u0A1u1Arur(modm) a^{k} \equiv A_{0}^{u_{0}} \cdot A_{1}^{u_{1}} \cdot \cdots A_{r}^{u_{r}} \pmod{m}

Explanation

Although the advancement of computers has made it possible to compute powers quickly, humans are never satisfied. The method of repeated squaring is a way to get the answer without directly computing immensely large numbers in modular operations. By dividing by (modm)\pmod{m} after each square, the number falls below mm, thus reducing the amount of computation.

Example

For instance, let’s calculate 7327(mod853)7^{327} \pmod{853}. Taking the common logarithm of 73277^{327} gives 327log73270.8450=276.315327 \log 7 \approx 327 \cdot 0.8450 = 276.315, which is a large number with a digit count of 277277. It’s too big to multiply honestly and then divide, so let’s use the above algorithm. If you do the binary expansion of 327327 as described in the algorithm, it looks like this: 327=256+64+4+2+1=28+26+22+21+20 327 = 256 + 64 + 4 +2 +1 = 2^{8} + 2^{6} + 2^2 + 2^1 + 2^0 If you calculate according to Step 2, 777(mod853)72(71)27249(mod853)74(72)24922401695(mod853)78(74)26952227(mod853)716(78)22272349(mod853)732(716)23492675(mod853)764(732)26752123(mod853)7128(764)21232628(mod853)7256(7128)26282298(mod853) \begin{align*} 7 & & \equiv 7 & \equiv 7 \pmod{853} \\ 7^{2} & \equiv ( 7^1 )^2 & \equiv 7^2 & \equiv 49 \pmod{853} \\ 7^{4} & \equiv ( 7^2 )^2 & \equiv 49^2 & \equiv 2401 \equiv 695 \pmod{853} \\ 7^{8} & \equiv ( 7^4 )^2 & \equiv 695^2 & \equiv 227 \pmod{853} \\ 7^{16} & \equiv ( 7^8 )^2 & \equiv 227^2 & \equiv 349 \pmod{853} \\ 7^{32} & \equiv ( 7^{16} )^2 & \equiv 349^2 & \equiv 675 \pmod{853} \\ 7^{64} & \equiv ( 7^{32} )^2 & \equiv 675^2 & \equiv 123 \pmod{853} \\ 7^{128} & \equiv ( 7^{64} )^2 & \equiv 123^2 & \equiv 628 \pmod{853} \\ 7^{256} & \equiv ( 7^{128} )^2 & \equiv 628^2 & \equiv 298 \pmod{853} \end{align*} Thus, 7327=7256764747271=298123695497286(mod853) 7^{327} = 7^{256} \cdot 7^{64} \cdot 7^{4} \cdot 7^{2} \cdot 7^{1} = 298 \cdot 123 \cdot 695 \cdot 49 \cdot 7 \equiv 286 \pmod{853} This method might still feel hard, but the amount of computation is quite reasonable compared to multiplying 327327 times individually.

Proof

ak=au0+2u1++2rur=au0a2u1a2rurA0u0A1u1Arur(modm) \begin{align*} a^{k} =& a^{u_{0} + 2u_{1} + \cdots + 2^r u_{r} } \\ =& a^{u_{0}} a^{ 2u_{1} } \cdots a^{ 2^r u_{r} } \\ \equiv & A_{0}^{u_{0}} A_{1}^{u_{1}} \cdots A_{r}^{u_{r}} \pmod{m} \end{align*}

Code

Below is an example code for calculating 7327(mod853)7^{327} \pmod{853} using the repeated squaring method, written in R language. Note that the power option accommodates negative integers, and if mod is a prime number, entering power=-1 allows you to find the multiplicative inverse of the given base.

FPM<-function(base,power,mod) #It is equal to (base^power)%%mod
{
  i<-0
  if (power<0) {
    while((base*i)%%mod != 1) {i=i+1}
    base<-i
    power<-(-power)}
  
  if (power==0) {return(1)}
  if (power==1) {return(base%%mod)}
  
  n<-0
  while(power>=2^n) {n=n+1}
  
  A<-rep(1,n)
  A[1]=base
  
  for(i in 1:(n-1)) {A[i+1]=(A[i]^2)%%mod}
  
  for(i in n:1) {
    if(power>=2^(i-1)) {power=power-2^(i-1)}
    else {A[i]=1} }
  
  for(i in 2:n) {A[1]=(A[1]*A[i])%%mod}
  
  return(A[1])
}
 
FPM(7,327,853)
FPM(7,-1,11)