// given ptP and ptQ, return addition result in ptP
void point_add(struct Point *p, struct Point *q){
	// P + Q = R
	// (xp, yp) + (xq, yq) = (xr, yr)
	// Y = (yq - yp)/(xq - xp)		(isn't that some kind of distance)?
	// xr = Y^2 + Y - xp - xq
	// yr = Y(xp - xr) - yp
	// assert: neither P or Q are point at infinity, and Px != Qx
	// if P or Q is point at infinity, R = other point
	// if P = Q, double instead
	// if Px == Qx, set P to point at infinity

	if(!point_iszero(q)) {
		if(point_iszero(p)) memcpy(p, q, sizeof(struct Point));
		else {
			if(point_isequal(p, q)) point_double(p);
			else if(bigint_isequal(p->x, q->x)) memset(p, 0, sizeof(struct Point));
			else {

				BIGINT t1,t2,t3,t4;

				// Yp + Yq
				bigint_add(t1, p->y, q->y);
				// Xp + Xq
				bigint_add(t2, p->x, q->x);
				// inv(t2)
				bigint_invert(t3, t2);
				// (Py + Yq) / (Xp + Xq)
				bigint_mul(t1, t1, t3);		// t1 = slope

				// slope^2
				bigint_mul(t4, t1, t1);
				// + slope
				bigint_add(t4, t4, t1);
				// + Xp + Xq
				bigint_add(t4, t4, t2);	// we have Xres

				// Xp - Xres
				bigint_add(t3, p->x, t4);
				// * slope
				bigint_mul(t3, t3, t1);
				// + Yq
				bigint_add(t3, t3, p->y);

				// copy Xres and Y res to P
				memcpy(p->x, t4, sizeof t4);
				bigint_add(p->y, t3, p->x);

			}
		}
	}
}