Efficient calculation of large exponents using the Pingala's algorithm

I was reading about how to efficiently calculate large exponents $$\pmod m$$ in the context of primality testing. A way is binary expansion. So if we need to calculate if $$a$$ is a witness of the non-primality of $$91$$ we would have to compute: $$a^{90}$$ and an efficient way to do that would be to notice that the binary expansion of $$90 = 64 + 16 + 8 + 2$$ and hence we can just create a table and using the successive squaring method we could calculate $$x^{90}$$ with just $$7$$ multiplications since (e.g. for $$x = 3$$):

$$x^{2^e}$$ $$\pmod {91}$$
$$3^1$$ 1
$$3^2$$ $$9$$
$$3^4$$ $$81$$
$$3^8$$ $$9$$
$$3^{16}$$ $$81$$
$$3^{32}$$ $$9$$
$$3^{64}$$ $$81$$

So we have: $$3^{90} = 3^{64}\cdot 3^{16}\cdot 3^{8}\cdot 3^{2} \equiv 1 \pmod {91}$$

So far, it is clear but then I came about what seems to be a more generic (?) and efficient approach for calculating large exponent's which is the Pingala's algorithm.
From what I understand, the algorithm (working on the exponent) proceeds to halve it (same process as when we calculate the binary version of a decimal) but if in the process of halving the number if odd, it first subtracts $$1$$ and then halves.
Then it follows the reverse process (bottom up) for the multiplication of a number where if the number was halved, we square the current number, otherwise we square and then multiply by the number.
Here is an example for simplifying $$3^{666} \pmod {667}$$. The left column goes top-down where if a number is even is halved, if a number is odd, first we subtract by $$1$$ and then half.
The right column goes bottom up, if we have divided a number we square, if we also subtracted we multiply by $$3$$ after squaring

Exponent $$e$$ $$3^e \pmod {667}$$
$$666$$ $$660$$
$$333$$ $$188$$
$$332$$ $$285$$
$$166$$ $$187$$
$$83$$ $$39$$
$$82$$ $$13$$
$$41$$ $$512$$
$$40$$ $$393$$
$$20$$ $$547$$
$$10$$ $$353$$
$$5$$ $$243$$
$$4$$ $$81$$
$$2$$ $$9$$
$$1$$ $$3$$

Apparently this process simplifies $$3^{666} \equiv 660 \pmod {667}$$

But I am lost understanding what is the idea/intuition here. I thought it does binary expansion, but can't understand why $$1$$ is subtracted. Additionally following the idea bottom up I re-wrote as follows:

$$(((((((((3^2)^2\cdot 3)^2)^2)^2\cdot 3)^2\cdot 3)^2)^2)\cdot 3)^2 \Leftrightarrow (((3^{108}\cdot 3^{32} \cdot 3^{16} \cdot 3^2) \cdot 3)^4 \cdot 3)^2 \Leftrightarrow 3^{197} \cdot 3^{256} \cdot 3^{128} \cdot 3^{16} \cdot 3^{8} \cdot 3^2 \Leftrightarrow 3^{607}$$

So I don't really understand the process, and how come my calculations are different.

• This is easiest to follow using the binary representation of the exponent, and it is known as binary exponentiation or square-and-multiply. See for example Exponentiation by squaring and Fast exponentiation algorithm - How to arrive at it?.
– dxiv
Commented Jul 18, 2021 at 0:22
• @dxiv: Yes I did understand the binary representation as mentioned in the example with $90$ in my post. I don't understand the subtraction with $1$ and why that table bottom up gives $\equiv 660 \pmod 667$, why my following of the multiplications leds to $3^{607}$ and how do we know that $\equiv 660$ is the correct answer
– Jim
Commented Jul 18, 2021 at 9:38
• @dxiv: did I misunderstood your comment? My question/confusion is about the pingala's algorithm and why we subtract $1$ and the final result.
– Jim
Commented Jul 19, 2021 at 7:45
• Subtracting $1$ is part of the conversion to binary, think at how it works for $6=1010_2$ for example. As for the final result, I am not sure I follow your table. The algorithm is worked out step by step in the previous links, and $666=2^1 + 2^3 + 2^4 + 2^7 + 2^9=2+8+16+128+512$ so $3^{666}=3^2\cdot 3^8 \cdot 3^{16}\cdot 3^{128}\cdot 3^{512}\,$. Maybe you should choose a smaller exponent to check that you got the steps right.
– dxiv
Commented Jul 19, 2021 at 7:55
• @dxiv: Ok, I think I got it. Thank you for your help
– Jim
Commented Jul 19, 2021 at 14:56

Suppose we need to compute $$a^{22}$$. We can do it by first computing $$a^{11}$$ and then squaring it to give $$(a^{11})^2 = a^{22}$$. We can compute $$a^{11}$$ similarly but since 11 is odd, we can express $$a^{11} = a^{10}a$$, and compute $$a^{10}$$ by first computing $$a^5$$ and squaring it. This recursive process can be repeated till the trivial subproblem of computing $$a^1$$. The subproblems created in this example have exponents as (in order of their creation): 11, 5, 2, 1.
Did you note how we handled the odd exponents $$e$$? By "subtracting 1" and halving the remaining even number $$(e-1)$$. And this action of subtracting is handled by the multiplication with $$a$$ after squaring.
The algorithm may be easy to understand and reason about if we consider the binary representation of $$b$$. It effectively processes the bits of $$b$$ from MSB to LSB: squaring for each bit, and additionally multiplying with $$a$$ for every 1 bit. For elaboration, you may refer to this article (written by me), section "Left-to-Right Binary Algorithm".