Fast exponentiation algorithm - How to arrive at it?

Write $b = b_0 + b_1 2 + b_2 2^2 + ... + b_n 2^n$ where $b_k$ are the binary digits of the representation of $b$ in base $2$. Note that $n$ varies logarithmically with $b$. Then:

$$ a^b = a^{b_0} \cdot (a^2)^{b_1} \cdot (a^{2^2})^{b_2} \cdot ... \cdot (a^{2^n})^{b_n} $$

The terms where bits $b_k = 0$ evaluate to $1$ and can be skipped. The terms with $b_k = 1$ are of the form $a^{2^k}$ and can be calculated by repeatedly squaring $a$ $k$ times. This is precisely what the posted code does.

The algorithm uses the binary expansion of the exponent to reduce the number of multiplications one has to do. If you take $a$ and square it and then square it again and then square it again, you produce the numbers $a, a^2, a^4, a^8,\ldots,a^{2^k}$ until $2^{k+1}>b$ is So that takes about $\log_2 b$ multiplications. Then if you express $b$ in binary, say, $b=1101001$, then $a^b = a^{2^6+2^5+2^3+2^0} = a^{2^6}a^{2^5}a^{2^3}a^{2^0}$, and you've just computed all these powers of $a$, so multiply them together and that's about $\log_2 b$ more multiplications.