I'm trying to implement elliptic curve point addition in hardware. I've managed to create a working module using Affine coordinates however I've been reading about how points can be added by using Jacobian coordinates more efficently due to needing less modular inversions. I then found this repeated point doubling algorithm which I've attempted to implement in C.
I take the Jacobian coordinates returned from this algorithm and convert them into Affine coordinates using the modular inversions $x = X/Z^2$ and $y = Y/Z^3$. I get the correct point on the curve when calculating $2P$ with $m = 1$ however if I set $m > 1$ to calculate $4P$, $8P$ etc. I get the wrong points.
I'm wondering if I'm doing something obviously wrong? I've included my c code below. I'm fairly new to all this so I'm sorry if there is some obvious errors.
int extendedEuclid(int *c, int phi, int e) {
int a[2] = {1,0};
int b[2] = {0,1};
int q;
while(1) {
q = phi / e;
phi = phi % e;
a[0] = a[0] - q*a[1];
b[0] = b[0] - q*b[1];
if(phi == 0) {
c[0] = e;
c[1] = a[1];
c[2] = b[1];
return;
};
q = e / phi;
e = e % phi;
a[1] = a[1] - q*a[0];
b[1] = b[1] - q*b[0];
if(e == 0) {
c[0] = phi;
c[1] = a[0];
c[2] = b[0];
return;
};
};
}
int jacobianDouble(int* x, int* y, int* z)
{
int a, b, w, m = 1;
*y = 2*(*y);
w = pow(*z, 4);
while(m > 0)
{
a = 3*(pow(*x, 2) - w);
b = *x*pow(*y, 2);
*x = pow(a, 2) - 2*b;
*z = *z*(*y);
m = m - 1;
if(m > 0)
{
w = w*pow(*y, 4);
}
*y = 2*a*(b - *x) - pow(*y, 4);
}
*y = *y/2;
}
int main() {
int x = 5, y = 4, z = 1, p = 97;
int ans[2];
jacobianDouble(&x, &y, &z);
extendedEuclid(ans, pow(z, 2), p);
x = fmod(x*ans[1], p);
if(x < 0)
{
x = p + x;
}
extendedEuclid(ans, pow(z, 3), p);
y = fmod(y*ans[1], p);
if(y < 0)
{
y = p + y;
}
printf("x = %d\n", x);
printf("y = %d\n", y);
return 0;
}