void pointZZ_pAdd(PointZZ_p * rop, const PointZZ_p * op1, const PointZZ_p * op2, const CurveZZ_p * curve) {
    // handle identity element cases
    if(pointZZ_pIsIdentityElement(op1) && pointZZ_pIsIdentityElement(op2)) {
        return pointZZ_pSetToIdentityElement(rop);
    } else if(pointZZ_pIsIdentityElement(op1)) {
        mpz_set(rop->x, op2->x);
        mpz_set(rop->y, op2->y);
        return;
    } else if(pointZZ_pIsIdentityElement(op2)) {
        mpz_set(rop->x, op1->x);
        mpz_set(rop->y, op1->y);
        return;
    }

    // use doubling algorithm if points are equal
    if(pointZZ_pEqual(op1, op2)) {
        pointZZ_pDouble(rop, op1, curve);
        return;
    }

    // check if points sum to identity element
    mpz_t negy;
    mpz_init(negy);
    mpz_sub(negy, curve->p, op2->y);
    if(mpz_cmp(op1->x, op2->x) == 0 && mpz_cmp(op1->y, negy) == 0) {
        mpz_clear(negy);
        return pointZZ_pSetToIdentityElement(rop);
    }


    mpz_t xdiff, ydiff, lambda;
    mpz_inits(xdiff, ydiff, lambda, NULL);

    // calculate lambda
    mpz_sub(ydiff, op2->y, op1->y);
    mpz_sub(xdiff, op2->x, op1->x);
    mpz_invert(xdiff, xdiff, curve->p);  // TODO check status
    mpz_mul(lambda, ydiff, xdiff);
    mpz_mod(lambda, lambda, curve->p);

    // calculate resulting x coord
    mpz_mul(rop->x, lambda, lambda);
    mpz_sub(rop->x, rop->x, op1->x);
    mpz_sub(rop->x, rop->x, op2->x);
    mpz_mod(rop->x, rop->x, curve->p);

    //calculate resulting y coord
    mpz_sub(rop->y, op1->x, rop->x);
    mpz_mul(rop->y, lambda, rop->y);
    mpz_sub(rop->y, rop->y, op1->y);
    mpz_mod(rop->y, rop->y, curve->p);

    mpz_clears(negy, xdiff, ydiff, lambda, NULL);
}
void pointZZ_pAdd(PointZZ_p * rop, const PointZZ_p * op1, const PointZZ_p * op2, const CurveZZ_p * curve) {
     
    if(pointZZ_pIsIdentityElement(op1) && pointZZ_pIsIdentityElement(op2)) {
        return pointZZ_pSetToIdentityElement(rop);
    } else if(pointZZ_pIsIdentityElement(op1)) {
        mpz_set(rop->x, op2->x);
        mpz_set(rop->y, op2->y);
        return;
    } else if(pointZZ_pIsIdentityElement(op2)) {
        mpz_set(rop->x, op1->x);
        mpz_set(rop->y, op1->y);
        return;
    }

     
    if(pointZZ_pEqual(op1, op2)) {
        pointZZ_pDouble(rop, op1, curve);
        return;
    }

     
    mpz_t negy;
    mpz_init(negy);
    mpz_sub(negy, curve->p, op2->y);
    if(mpz_cmp(op1->x, op2->x) == 0 && mpz_cmp(op1->y, negy) == 0) {
        mpz_clear(negy);
        return pointZZ_pSetToIdentityElement(rop);
    }


    mpz_t xdiff, ydiff, lambda;
    mpz_inits(xdiff, ydiff, lambda, NULL);

     
    mpz_sub(ydiff, op2->y, op1->y);
    mpz_sub(xdiff, op2->x, op1->x);
    mpz_invert(xdiff, xdiff, curve->p);   
    mpz_mul(lambda, ydiff, xdiff);
    mpz_mod(lambda, lambda, curve->p);

     
    mpz_mul(rop->x, lambda, lambda);
    mpz_sub(rop->x, rop->x, op1->x);
    mpz_sub(rop->x, rop->x, op2->x);
    mpz_mod(rop->x, rop->x, curve->p);

     
    mpz_sub(rop->y, op1->x, rop->x);
    mpz_mul(rop->y, lambda, rop->y);
    mpz_sub(rop->y, rop->y, op1->y);
    mpz_mod(rop->y, rop->y, curve->p);

    mpz_clears(negy, xdiff, ydiff, lambda, NULL);
}