I was trying to implement the double and addition formulas for elliptic points on a Montgomery elliptic curve. I came across this weird thing which should definitely not be happening.
I took a point $P_1$, and doubled it twice. This should be $4P_1$. Let's call this point $P_4'$.
Next I doubled $P_1$ to get $P_2$. Then carried out $P_1 + P_2$ to get $P_3$, and then again added $P_3$ and $P_1$ to get $P_4$.
Now I was expecting both these points to be the same, since both are simply $4$ times $P_1$.
For example, let us take a curve $y^2 = x^3 + 6x^2 + x$ (Montgomery Curve), over $F_{2425967623052370772757633156976982469681}$
Let a point $P_1$ be $[2 :: 1]$ (projective coordinates, $y$ dropped)
Doubling this twice, gives the point $[339112225 :: 126909216]$ ($P_4'$).
Adding it stepwise to reach $4P_1$, gives the point $[106258781030400 :: 39766241378304]$
Question: Shouldn't $P_4$ be equal to $P_4'$ ? Since both are equal to $4P_1$?
I am using the following formulas for addition and doubling.
$$X_{m+n} = (z_{m-n}\cdot((x_m - z_m)\cdot(x_n + z_n) + (x_m + z_m)\cdot(x_n - z_n))^2 \mod{n})\\ Z_{m+n} = (x_{m-n}\cdot((x_m - z_m)\cdot(x_n + z_n) - (x_m + z_m)\cdot(x_n - z_n))^2 \mod{n})\\ X_{2n} = (((x_n + z_n)^2)\cdot((x_n - z_n)^2) \mod{n})\\ Z_{2n} = ((4\cdot x_n\cdot z_n)\cdot((x_n - z_n)^2 + {a+2\over4}\cdot 4\cdot x_n\cdot z_n) \mod{n})$$
These formulas were mentioned here http://en.wikipedia.org/wiki/Montgomery_curve
I have implemented a python version of the addition and double, and carried out the above arithmetic in this piece of code.