Issues with Elliptic Curve Point Addition 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 &gt; 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 &gt; 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 &gt; 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(&amp;x, &amp;y, &amp;z);

    extendedEuclid(ans, pow(z, 2), p);
    x = fmod(x*ans[1], p);
    if(x &lt; 0)
    {
        x = p + x;
    }

    extendedEuclid(ans, pow(z, 3), p);
    y = fmod(y*ans[1], p);
    if(y &lt; 0)
    {
        y = p + y;
    }

    printf("x = %d\n", x);
    printf("y = %d\n", y);

    return 0;
}
All arithmetic operations in jacobianDouble should be in the finite field $\mathbb{F}_p$, i.e. all results should mod $p$, not plain integer arithmetics. Note division Y/2 is defined as $Y\cdot 2^{-1}$ where $2^{-1}$ is the multiplicative inverse of 2 modulo $p$.