addlaws.c

/* addlaws.c - various addition laws for ECM
   Author: F. Morain
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <gmp.h> /* GMP header file */

#include "ecm.h" /* ecm header file */
#include "ecm-impl.h"
#include "ecm-ecm.h"
#include "mpmod.h"

#include "addlaws.h"

#if DEBUG_ADD_LAWS >= 1
void
print_mpz_from_mpres(mpres_t x, mpmod_t n)
{
    mpz_t tmp;

    mpz_init(tmp);
    mpres_get_z(tmp, x, n);
    gmp_printf("%Zd", tmp);
    mpz_clear(tmp);
}
#endif

/******************** Weierstrass section ********************/

void
pt_w_set_to_zero(ell_point_t P, mpmod_t n)
{
    mpres_set_ui(P->x, 0, n);
    mpres_set_ui(P->y, 1, n);
    mpres_set_ui(P->z, 0, n);
}

int
pt_w_is_zero(mpres_t z, mpmod_t n)
{
    return mpres_is_zero(z, n);
}

void
pt_w_set(mpres_t x0, mpres_t y0, mpres_t z0,
	    mpres_t x, mpres_t y, mpres_t z,
	    ATTRIBUTE_UNUSED mpmod_t n)
{
  mpres_set(x0, x, n);
  mpres_set(y0, y, n);
  mpres_set(z0, z, n);
}

#if DEBUG_ADD_LAWS >= 1
void
pt_w_print(mpres_t x, mpres_t y, mpres_t z, ell_curve_t E, mpmod_t n)
{
    printf("[");
    print_mpz_from_mpres(x, n);
    printf(", ");
    print_mpz_from_mpres(y, n);
    printf(", ");
    if(E->type == ECM_EC_TYPE_WEIERSTRASS && E->law == ECM_LAW_AFFINE)
	gmp_printf("%Zd", z);
    else
	print_mpz_from_mpres(z, n);
    printf("]");
}
#endif

/* [x0, y0, z0] <- [x1, y1, z1] + [x2, y2, z2] using lambda=num/den
   with buffer inv.

   (lambda*x+mu)^2+a1*x*(lambda*x+mu)+a3*(lambda*x+mu)=x^3+a2*x^2+...
   x^3+(a2-lambda^2-a1*lambda)*x^2+... = 0
   x1+x2+x3 = lambda^2+a1*lambda-a2.
   y3 = lambda*(x1-x3)-y1-a1*x3-a3
 */
static int
pt_w_common_aff(mpz_t f, mpres_t x0, mpres_t y0, mpres_t z0,
		mpres_t x1, mpres_t y1,
		mpres_t x2, mpres_t a1, mpres_t a3, mpres_t a2,
		mpmod_t n, mpres_t num, mpres_t den, mpres_t lambda)
{
    if(mpres_invert(lambda, den, n) == 0){
	mpres_gcd(f, den, n);
	return 0;
    }
    /** lambda = num/den **/
    mpres_mul(lambda, lambda, num, n);
    /** num <- (lambda+a1)*lambda **/
    mpres_add(num, lambda, a1, n);
    mpres_mul(num, num, lambda, n);
    mpres_sub(num, num, a2, n);
    /** x0 = den <- num-x1-x2 **/
    mpres_sub(den, num, x1, n);
    mpres_sub(den, den, x2, n);
    /** y0 = num <- lambda*(x1-x0)-(y1+a1*x0+a3) **/
    mpres_sub(num, x1, den, n);
    mpres_mul(num, num, lambda, n);
    mpres_sub(y0, num, y1, n);
    mpres_sub(y0, y0, a3, n);
    mpres_mul(x0, a1, den, n);
    mpres_sub(y0, y0, x0, n);
    /** finish **/
    mpres_set(x0, den, n);
    mpz_set_ui(z0, 1); /* just in case */
    return 1;
}

/* [x3, y3, z3] <- [2] * [x1, y1, z1] */
int
pt_w_duplicate(mpz_t f, mpres_t x3, mpres_t y3, mpres_t z3,
	       mpres_t x1, mpres_t y1, mpres_t z1,
	       mpmod_t n, ell_curve_t E)
{
    if(pt_w_is_zero(z1, n) == 1){
      pt_w_set(x3, y3, z3, x1, y1, z1, n);
      return 1;
    }
    if(E->type == ECM_EC_TYPE_WEIERSTRASS && E->law == ECM_LAW_AFFINE){
	/* buf[1] <- 2*y1+a1*x1+a3 */
	mpres_mul(E->buf[1], E->a1, x1, n);
	mpres_add(E->buf[1], E->buf[1], E->a3, n);
	mpres_add(E->buf[1], E->buf[1], y1, n);
	mpres_add(E->buf[1], E->buf[1], y1, n);
	if(mpres_is_zero(E->buf[1], n)){
	    /* buf1 = 0 <=> P is a [2]-torsion point */
	    mpres_set_ui(x3, 0, n);
	    mpres_set_ui(y3, 1, n);
	    mpres_set_ui(z3, 0, n);
	    return 1;
	}
	/* buf[0] <- 3*x^2+2*a2*x+a4-a1*y = (3*x+2*a2)*x+a4-a1*y */
	mpres_mul_ui(E->buf[0], x1, 3, n);
	mpres_add(E->buf[0], E->buf[0], E->a2, n);
	mpres_add(E->buf[0], E->buf[0], E->a2, n);
	mpres_mul(E->buf[0], E->buf[0], x1, n);
	mpres_add(E->buf[0], E->buf[0], E->a4, n);
	mpres_mul(E->buf[2], E->a1, y1, n);
	mpres_sub(E->buf[0], E->buf[0], E->buf[2], n);
	return pt_w_common_aff(f, x3, y3, z3, x1, y1, x1, 
			       E->a1, E->a3, E->a2, n, 
			       E->buf[0], E->buf[1], E->buf[2]);
    }
    else if(E->type == ECM_EC_TYPE_WEIERSTRASS 
	    && E->law == ECM_LAW_HOMOGENEOUS){
	/* source is dbl-2007-bl: 5M + 6S + 1*a + 7add + 3*2 + 1*3 */
	/* mapping: h = buf[0], w = buf[1], s = buf[2], RR = buf[3], B = buf[4];*/
	/*  h:=X1^2 mod p;	    # S*/
	mpres_sqr(E->buf[0], x1, n);
	/*	w:=Z1^2 mod p;*/
	mpres_sqr(E->buf[1], z1, n);
	/*	w:=a*w mod p;*/
	mpres_mul(E->buf[1], E->buf[1], E->a4, n);
	/*	s:=3*h mod p;	    # *3*/
	mpres_mul_ui(E->buf[2], E->buf[0], 3, n);
	/*	w:=w+s mod p;*/
	mpres_add(E->buf[1], E->buf[1], E->buf[2], n);
	/*	s:=Y1*Z1 mod p;*/
	mpres_mul(E->buf[2], y1, z1, n);
	/*	s:=2*s mod p;*/
	mpres_mul_ui(E->buf[2], E->buf[2], 2, n);
	/*	Z3:=s^2 mod p;*/
	mpres_sqr(z3, E->buf[2], n);
	/*	Z3:=s*Z3 mod p;*/
	mpres_mul(z3, z3, E->buf[2], n);
	/*	RR:=Y1*s mod p;	    # M*/
	mpres_mul(E->buf[3], y1, E->buf[2], n);
	/*	B:=X1+RR mod p;    # add*/
	mpres_add(E->buf[4], x1, E->buf[3], n);
	/*	B:=B^2 mod p;*/
	mpres_sqr(E->buf[4], E->buf[4], n);
	/*	RR:=RR^2 mod p;	    # S*/
	mpres_sqr(E->buf[3], E->buf[3], n);
	/*	B:=B-h mod p;*/
	mpres_sub(E->buf[4], E->buf[4], E->buf[0], n);
	/*	B:=B-RR mod p;*/
	mpres_sub(E->buf[4], E->buf[4], E->buf[3], n);
	/*	h:=w^2 mod p;*/
	mpres_sqr(E->buf[0], E->buf[1], n);
	/*	X3:=2*B mod p;*/
	mpres_mul_ui(x3, E->buf[4], 2, n);
	/*	h:=h-X3 mod p;*/
	mpres_sub(E->buf[0], E->buf[0], x3, n);
	/*	X3:=h*s mod p;	    # M*/
	mpres_mul(x3, E->buf[0], E->buf[2], n);
	/*	s:=B-h mod p;*/
	mpres_sub(E->buf[2], E->buf[4], E->buf[0], n);
	/*	s:=w*s mod p;*/
	mpres_mul(E->buf[2], E->buf[2], E->buf[1], n);
	/*	Y3:=2*RR mod p;*/
	mpres_mul_ui(y3, E->buf[3], 2, n);
	/*	Y3:=s-Y3 mod p;*/
	mpres_sub(y3, E->buf[2], y3, n);
	return 1;
    }
    return 0;
}

/* [x3, y3, z3] <- [x1, y1, z1] + [x2, y2, z2]; P3 can be either P1 or P2. */
int
pt_w_add(mpz_t f, mpres_t x3, mpres_t y3, mpres_t z3,
	 mpres_t x1, mpres_t y1, mpres_t z1,
	 mpres_t x2, mpres_t y2, mpres_t z2,
	 mpmod_t n, ell_curve_t E)
{
    if(pt_w_is_zero(z1, n)){
	pt_w_set(x3, y3, z3, x2, y2, z2, n);
	return 1;
    }
    else if(pt_w_is_zero(z2, n)){
	pt_w_set(x3, y3, z3, x1, y1, z1, n);
	return 1;
    }
    if(E->type == ECM_EC_TYPE_WEIERSTRASS && E->law == ECM_LAW_AFFINE)
	if(mpres_equal(x1, x2, n) && mpres_equal(y1, y2, n))
	    return pt_w_duplicate(f, x3, y3, z3, x1, y1, z1, n, E);
	else{
	    mpres_sub(E->buf[0], y1, y2, n);
	    mpres_sub(E->buf[1], x1, x2, n);
	    return pt_w_common_aff(f, x3, y3, z3, x1, y1, x2, 
				   E->a1, E->a3, E->a2,
				   n, E->buf[0], E->buf[1], E->buf[2]);
	}
    else if(E->type == ECM_EC_TYPE_WEIERSTRASS 
	    && E->law == ECM_LAW_HOMOGENEOUS){
	/* Cohen-Miyaji-Ono: 12M+2S+6add+1*2 */
	/* mapping: y1z2 = buf, AA = buf+1, u = buf+2, v = buf+3, R = buf+4, */
	/* vvv = buf+5; */
#if DEBUG_ADD_LAWS >= 2
	printf("y1="); print_mpz_from_mpres(y1, n); printf("\n");
	printf("y2="); print_mpz_from_mpres(y2, n); printf("\n");
	printf("z1="); print_mpz_from_mpres(z1, n); printf("\n");
	printf("z2="); print_mpz_from_mpres(z2, n); printf("\n");
#endif
	/*  Y1Z2:=Y1*Z2 mod p;	# M*/
	mpres_mul(E->buf[0], y1, z2, n);
	/*	A:=X1*Z2 mod p;	# M*/
	mpres_mul(E->buf[1], x1, z2, n);
	/*	u:=Y2*Z1 mod p;*/
	mpres_mul(E->buf[2], y2, z1, n);
	/*	u:=u-Y1Z2 mod p;*/
	mpres_sub(E->buf[2], E->buf[2], E->buf[0], n);
	/*	v:=X2*Z1 mod p;*/
	mpres_mul(E->buf[3], x2, z1, n);
	/*	v:=v-A mod p;*/
	mpres_sub(E->buf[3], E->buf[3], E->buf[1], n);
	if(mpz_sgn(E->buf[2]) == 0 && mpz_sgn(E->buf[3]) == 0){
	    /* u = 0 <=> Y2*Z1 = Y1*Z2 <=> Y2/Z2 = Y1/Z1*/
	    /* v = 0 <=> X2*Z1 = X1*Z2 <=> X2/Z2 = X1/Z1*/
	    return pt_w_duplicate(f, x3, y3, z3, x1, y1, z1, n, E);
	}
	/*	Z3:=Z1*Z2 mod p;	# M*/
	mpres_mul(z3, z1, z2, n);
	/*	X3:=u^2 mod p;*/
	mpres_sqr(x3, E->buf[2], n);
	/*	X3:=X3*Z3 mod p;*/
	mpres_mul(x3, x3, z3, n);
	/*	R:=v^2 mod p;*/
	mpres_sqr(E->buf[4], E->buf[3], n);
	/*	vvv:=v*R mod p;*/
	mpres_mul(E->buf[5], E->buf[3], E->buf[4], n);
	/*	R:=R*A mod p;*/
	mpres_mul(E->buf[4], E->buf[4], E->buf[1], n);
	/*	Y3:=2*R mod p; 		# *2*/
	mpres_mul_ui(y3, E->buf[4], 2, n);
	/*	A:=X3-vvv mod p;*/
	mpres_sub(E->buf[1], x3, E->buf[5], n);
	/*	A:=A-Y3 mod p;*/
	mpres_sub(E->buf[1], E->buf[1], y3, n);
	/*	X3:=v*A mod p;		# M*/
	mpres_mul(x3, E->buf[3], E->buf[1], n);
	/*	Y3:=R-A mod p;*/
	mpres_sub(y3, E->buf[4], E->buf[1], n);
	/*	Y3:=u*Y3 mod p;*/
	mpres_mul(y3, y3, E->buf[2], n);
	/*	A:=vvv*Y1Z2 mod p;*/
	mpres_mul(E->buf[1], E->buf[5], E->buf[0], n);
	/*	Y3:=Y3-A mod p;*/
	mpres_sub(y3, y3, E->buf[1], n);
	/*	Z3:=vvv*Z3 mod p;	# M*/
	mpres_mul(z3, z3, E->buf[5], n);
	return 1;
    }
    return 0;
}

#if USE_ADD_SUB_CHAINS > 0
/* [x3, y3, z3] <- [x1, y1, z1] - [x2, y2, z2]; P3 != P1, P3 != P2. 
   -P2 ~ -(x2/z2, y2/z2, 1) = (x2/z2, -y2/z2-a1*x/z2-a3, 1) 
                            ~ (x2, -y2-a1*x2-a3*z2, z2).
*/
int
pt_w_sub(mpz_t f, mpres_t x3, mpres_t y3, mpres_t z3,
	 mpres_t x1, mpres_t y1, mpres_t z1,
	 mpres_t x2, mpres_t y2, mpres_t z2,
	 mpmod_t n, ell_curve_t E)
{
    int res = 1;

    if(E->law == ECM_LAW_HOMOGENEOUS){
	/* FIXME: does not work for complete equation! */
	mpres_neg(y2, y2, n);
	res = pt_w_add(f, x3, y3, z3, x1, y1, z1, x2, y2, z2, n, E);
	mpres_neg(y2, y2, n);
    }
    else if(E->law == ECM_LAW_AFFINE){
	/* buf[3] not used in law, so use it */
	mpres_mul(E->buf[3], E->a1, x2, n);
	mpres_add(E->buf[3], E->buf[3], E->a3, n);
	mpres_add(E->buf[3], E->buf[3], y2, n);
	mpres_neg(E->buf[3], E->buf[3], n);
	res = pt_w_add(f, x3, y3, z3, x1, y1, z1, x2, E->buf[3], z2, n, E);
    }
    return res;
}
#endif

/******************** projective Hessian form ********************/

/* U^3+V^3+W^3 = 3*D*U*V*W, D^3 <> 1.
   O_H = [1:-1:0]
   -[u:v:w] = [v:u:w]
   Warning: there can exist two other points at infinity, namely
   [1:-omega:0] and [1:-omega^2:0] where omega^3 = 1.
*/
int
hessian_is_zero(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
{
    mpres_t tmp;
    int ret;

    if(mpz_sgn(P->z) != 0)
	return 0;
    mpres_init(tmp, n);
    mpres_add(tmp, P->x, P->y, n);
    ret = mpz_sgn(tmp) == 0;
#if 0
    if(ret)
	gmp_printf("found a third root of unity? %Zd/%Zd\n", P->x, P->y);
#endif
    mpres_clear(tmp, n);
    return ret;
}

void
hessian_set_to_zero(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
{
    mpres_set_si(P->x,  1, n);
    mpres_set_si(P->y, -1, n);
    mpres_set_si(P->z,  0, n);
}

#if DEBUG_ADD_LAWS >= 1
void
hessian_print(ell_point_t P, ell_curve_t E, mpmod_t n)
{
    pt_w_print(P->x, P->y, P->z, E, n);
}
#endif

#if USE_ADD_SUB_CHAINS > 0
/* -[u:v:w] = [v:u:w] */
void
hessian_negate(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, ATTRIBUTE_UNUSED mpmod_t n)
{
    mpz_swap(P->x, P->y); /* humf */
}
#endif

/* TODO: decrease the number of buffers? */
int
hessian_duplicate(ell_point_t R, ell_point_t P, 
		  ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
{
    /* A = buf[0], ..., G = buf[6], H = buf[7], J = buf[8] */
    /* A:=P[1]^2 mod N; */
    mpres_mul(E->buf[0], P->x, P->x, n);
    /* B:=P[2]^2 mod N; */
    mpres_mul(E->buf[1], P->y, P->y, n);
    /* C:=P[3]^2 mod N; */
    mpres_mul(E->buf[2], P->z, P->z, n);
    /* D:=(A+B) mod N; */
    mpres_add(E->buf[3], E->buf[0], E->buf[1], n);
    /* E:=(A+C) mod N; */
    mpres_add(E->buf[4], E->buf[0], E->buf[2], n);
    /* F:=(B+C) mod N; */
    mpres_add(E->buf[5], E->buf[1], E->buf[2], n);
    /* G:=((P[1]+P[2])^2-D) mod N; */
    mpres_add(E->buf[6], P->x, P->y, n);
    mpres_mul(E->buf[6], E->buf[6], E->buf[6], n);
    mpres_sub(E->buf[6], E->buf[6], E->buf[3], n);
    /* H:=((P[1]+P[3])^2-E) mod N; */
    mpres_add(E->buf[7], P->x, P->z, n);
    mpres_mul(E->buf[7], E->buf[7], E->buf[7], n);
    mpres_sub(E->buf[7], E->buf[7], E->buf[4], n);
    /* J:=((P[2]+P[3])^2-F) mod N; */
    mpres_add(E->buf[8], P->y, P->z, n);
    mpres_mul(E->buf[8], E->buf[8], E->buf[8], n);
    mpres_sub(E->buf[8], E->buf[8], E->buf[5], n);
    /* R->x = ((J-G)*(H+2*E)) mod N */
    mpres_sub(E->buf[0], E->buf[8], E->buf[6], n);
    mpres_add(E->buf[1], E->buf[7], E->buf[4], n);
    mpres_add(E->buf[1], E->buf[1], E->buf[4], n);
    mpres_mul(R->x, E->buf[0], E->buf[1], n);
    /* R->y = ((G-H)*(J+2*F)) mod N */
    mpres_sub(E->buf[0], E->buf[6], E->buf[7], n);
    mpres_add(E->buf[1], E->buf[8], E->buf[5], n);
    mpres_add(E->buf[1], E->buf[1], E->buf[5], n);
    mpres_mul(R->y, E->buf[0], E->buf[1], n);
    /* R->z = ((H-J)*(G+2*D)) mod N */
    mpres_sub(E->buf[0], E->buf[7], E->buf[8], n);
    mpres_add(E->buf[1], E->buf[6], E->buf[3], n);
    mpres_add(E->buf[1], E->buf[1], E->buf[3], n);
    mpres_mul(R->z, E->buf[0], E->buf[1], n);
    return 1;
}

/* TODO: reduce the number of buffers? */
int
hessian_plus(ell_point_t R, ell_point_t P, ell_point_t Q, 
	     ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
{
    /* P = [T1,T2,T3], Q = [T4,T5,T6] */
    /* P = Q <=> T1/T3=T4/T6 and T2/T3=T5/T6 
             <=> T1*T6=T3*T4 and T2*T6=T3*T5
     */
    /* T1 = buf[0], ..., T7 = buf[6] */
    /* T7:=(T1*T6) mod N; */
    mpres_mul(E->buf[6], P->x, Q->z, n);
    /* T1:=(T1*T5) mod N; */
    mpres_mul(E->buf[0], P->x, Q->y, n);
    /* T5:=(T3*T5) mod N; */
    mpres_mul(E->buf[4], P->z, Q->y, n);
    /* T3:=(T3*T4) mod N; */
    mpres_mul(E->buf[2], P->z, Q->x, n);
    /* T4:=(T2*T4) mod N; */
    mpres_mul(E->buf[3], P->y, Q->x, n);
    /* T2:=(T2*T6) mod N; */
    mpres_mul(E->buf[1], P->y, Q->z, n);

    if(mpres_equal(E->buf[6], E->buf[2], n)
       && mpres_equal(E->buf[4], E->buf[1], n))
	/* as a matter of that, P = Q and we need duplicate */
	return hessian_duplicate(R, P, E, n);

    /* T6:=(T2*T7) mod N; */
    mpres_mul(E->buf[5], E->buf[1], E->buf[6], n);
    /* T2:=(T2*T4) mod N; */
    mpres_mul(E->buf[1], E->buf[1], E->buf[3], n);
    /* T4:=(T3*T4) mod N; */
    mpres_mul(E->buf[3], E->buf[2], E->buf[3], n);
    /* T3:=(T3*T5) mod N; */
    mpres_mul(E->buf[2], E->buf[2], E->buf[4], n);
    /* T5:=(T1*T5) mod N; */
    mpres_mul(E->buf[4], E->buf[0], E->buf[4], n);
    /* T1:=(T1*T7) mod N; */
    mpres_mul(E->buf[0], E->buf[0], E->buf[6], n);
    /* T1:=(T1-T4) mod N; */
    mpres_sub(R->y, E->buf[0], E->buf[3], n);
    /* T2:=(T2-T5) mod N; */
    mpres_sub(R->x, E->buf[1], E->buf[4], n);
    /* T3:=(T3-T6) mod N; */
    mpres_sub(R->z, E->buf[2], E->buf[5], n);
    /* return [T2, T1, T3]; */
    return 1;
}

int
hessian_add(ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
{
    if(hessian_is_zero(P, E, n)){
	ell_point_set(R, Q, E, n);
	return 1;
    }
    else if(hessian_is_zero(Q, E, n)){
	ell_point_set(R, P, E, n);
	return 1;
    }
    else
	return hessian_plus(R, P, Q, E, n);
}

#if USE_ADD_SUB_CHAINS > 0
int
hessian_sub(ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
{
    int ret;

    hessian_negate(Q, E, n);
    ret = hessian_add(R, P, Q, E, n);
    hessian_negate(Q, E, n);
    return ret;
}
#endif

/* switch from X^3+Y^3+1=3*D*X*Y to Y^2=X^3+A*X+B
   A:=-27*D*(D^3+8);
   B:=54*(D^6-20*D^3-8);
   xi:=12*(D^3-1)/(D*u+v+1);
   x:=-9*D^2+xi*u;
   y:=3*xi*(v-1);
   OUTPUT: If a factor is found during the inversion, it is put in f and
   ECM_FACTOR_FOUND_STEP1 is returned. Otherwise, ECM_NO_FACTOR_FOUND is
   returned.
   SIDE-EFFECT: (x, y, D) <- (x_on_W, y_on_W, A_of_W)
 */
int
hessian_to_weierstrass(mpz_t f, mpres_t x, mpres_t y, mpres_t D, mpmod_t n)
{
    mpres_t D3, A, xi, tmp1, tmp2;
    int ret = ECM_NO_FACTOR_FOUND;

#if DEBUG_ADD_LAWS >= 1
    printf("P:=[");
    print_mpz_from_mpres(x, n);
    printf(", ");
    print_mpz_from_mpres(y, n);
    printf(", 1];\n");
    printf("D:=");
    print_mpz_from_mpres(D, n);
    printf(";\n");
#endif
    /* D3 <- D^3 */
    mpres_init(D3, n);
    mpres_mul(D3, D, D, n);
    mpres_mul(D3, D3, D, n);
    /* finish A */
    mpres_init(A, n);
    mpres_add_ui(A, D3, 8, n);
    mpres_mul(A, A, D, n);
    mpres_mul_ui(A, A, 27, n);
    mpres_neg(A, A, n);
    /* compute xi */
    mpres_init(xi, n);
    mpres_init(tmp1, n);
    mpres_mul(tmp1, D, x, n);
    mpres_add(tmp1, tmp1, y, n);
    mpres_add_ui(tmp1, tmp1, 1, n);
    mpres_init(tmp2, n);
    mpres_sub_ui(tmp2, D3, 1, n);
    mpres_mul_ui(tmp2, tmp2, 12, n);
    if(mpres_invert(xi, tmp1, n) == 0){
	mpres_gcd(f, tmp1, n);
	ret = ECM_FACTOR_FOUND_STEP1;
    }
    else{
	mpres_mul(xi, xi, tmp2, n);
	/* compute x */
	mpres_mul(tmp1, D, D, n);
	mpres_mul_ui(tmp1, tmp1, 9, n);
	mpres_mul(tmp2, xi, x, n);
	mpres_sub(x, tmp2, tmp1, n);
	/* compute y */
	mpres_sub_ui(tmp1, y, 1, n);
	mpres_mul(tmp1, tmp1, xi, n);
	mpres_mul_ui(y, tmp1, 3, n);
	mpres_set(D, A, n);
#if DEBUG_ADD_LAWS >= 1
	printf("WP:=[");
	print_mpz_from_mpres(x, n);
	printf(", ");
	print_mpz_from_mpres(y, n);
	printf(", 1];\n");
	printf("WA:=");
	print_mpz_from_mpres(D, n);
	printf(";\nWB:=(WP[2]^2-WP[1]^3-WA*WP[1]) mod N;WE:=[WA, WB];\n");
#endif
    }
    mpres_clear(A, n);
    mpres_clear(D3, n);
    mpres_clear(xi, n);
    mpres_clear(tmp1, n);
    mpres_clear(tmp2, n);
    return ret;
}

int
mult_by_3(mpz_t f, mpres_t x, mpres_t y, mpres_t A, mpmod_t n)
{
    ell_curve_t E;
    ell_point_t P, Q;
    int ret = ECM_NO_FACTOR_FOUND;
    mpz_t e;

    ell_curve_init_set(E, ECM_EC_TYPE_WEIERSTRASS, ECM_LAW_AFFINE, A, n);
    ell_point_init(P, E, n);
    mpres_set(P->x, x, n);
    mpres_set(P->y, y, n);
    mpres_set_ui(P->z, 1, n);
    ell_point_init(Q, E, n);
    mpz_init_set_ui(e, 3);
    if(ell_point_mul(f, Q, e, P, E, n) != 0){
	mpres_set(x, Q->x, n);
	mpres_set(y, Q->y, n);
    }
    mpz_clear(e);
    ell_point_clear(Q, E, n);
    ell_point_clear(P, E, n);
    ell_curve_clear(E, n);
    return ret;
}

/******************** projective twisted Hessian form ********************/

/* a*U^3+V^3+W^3 = d*U*V*W
   O_E = [0:-1:1]
   -[U:V:W]=[U:W:V]
*/
int
twisted_hessian_is_zero(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
{
    mpres_t tmp;
    int ret;

    if(mpz_sgn(P->x) != 0)
	return 0;
    mpres_init(tmp, n);
    mpres_add(tmp, P->y, P->z, n);
    ret = mpz_sgn(tmp) == 0;
#if 0
    if(ret)
	gmp_printf("found a third root of unity? %Zd/%Zd\n", P->x, P->y);
#endif
    mpres_clear(tmp, n);
    return ret;
}

void
twisted_hessian_set_to_zero(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
{
    mpres_set_si(P->x,  0, n);
    mpres_set_si(P->y, -1, n);
    mpres_set_si(P->z,  1, n);
}

#if DEBUG_ADD_LAWS >= 1
void
twisted_hessian_print(ell_point_t P, ell_curve_t E, mpmod_t n)
{
    pt_w_print(P->x, P->y, P->z, E, n);
}
#endif

#if USE_ADD_SUB_CHAINS > 0
/* -[u:v:w] = [u:w:v] */
void
twisted_hessian_negate(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, ATTRIBUTE_UNUSED mpmod_t n)
{
    mpz_swap(P->y, P->z); /* humf */
}
#endif

/* TODO: decrease the number of buffers? */
/* 6M+2S+1M_d: better when d is small */
int
twisted_hessian_duplicate(ell_point_t R, ell_point_t P, 
		  ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
{
    /* R = buf[0], ..., W = buf[5], C = buf[6], D = buf[7], E = buf[8] */
    /*  R:=Y1+Z1;*/
    mpres_add(E->buf[0], P->y, P->z, n);
    /*	S:=Y1-Z1;*/
    mpres_sub(E->buf[1], P->y, P->z, n);
    /*	T:=R^2 mod N;*/
    mpres_sqr(E->buf[2], E->buf[0], n);
    /*	U:=S^2 mod N;*/
    mpres_sqr(E->buf[3], E->buf[1], n);
    /*	V:=T+3*U;*/
    mpres_add(E->buf[4], E->buf[2], E->buf[3], n);
    mpres_add(E->buf[4], E->buf[4], E->buf[3], n);
    mpres_add(E->buf[4], E->buf[4], E->buf[3], n);
    /*	W:=3*T+U;*/
    mpres_add(E->buf[5], E->buf[3], E->buf[2], n);
    mpres_add(E->buf[5], E->buf[5], E->buf[2], n);
    mpres_add(E->buf[5], E->buf[5], E->buf[2], n);
    /*	C:=(R*V) mod N;*/
    mpres_mul(E->buf[6], E->buf[0], E->buf[4], n);
    /*	D:=(S*W) mod N;*/
    mpres_mul(E->buf[7], E->buf[1], E->buf[5], n);
    /*	E:=(3*C-E0[2]*X1*(W-V)) mod N;*/
    mpres_sub(E->buf[8], E->buf[5], E->buf[4], n);
    mpres_mul(E->buf[8], E->buf[8], P->x, n);
    mpres_mul(E->buf[8], E->buf[8], E->a6, n);
    mpres_sub(E->buf[8], E->buf[6], E->buf[8], n);
    mpres_add(E->buf[8], E->buf[8], E->buf[6], n);
    mpres_add(E->buf[8], E->buf[8], E->buf[6], n);
    /*	X3:=(-2*X1*D) mod N;*/
    mpres_mul(R->x, P->x, E->buf[7], n);
    mpres_add(R->x, R->x, R->x, n);
    mpres_neg(R->x, R->x, n);
    /*	Y3:=((D+E)*Z1) mod N;*/
    mpres_add(E->buf[0], E->buf[7], E->buf[8], n);
    mpres_mul(E->buf[1], E->buf[0], P->z, n);
    /*	Z3:=((D-E)*Y1) mod N;*/
    mpres_sub(E->buf[0], E->buf[7], E->buf[8], n);
    mpres_mul(R->z, E->buf[0], P->y, n);
    mpres_set(R->y, E->buf[1], n);
    return 1;
}

/* TODO: reduce the number of buffers? */
int
twisted_hessian_plus(ell_point_t R, ell_point_t P, ell_point_t Q, 
	     ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
{
    /* A = buf[0], ... F = buf[5], G = [6], H = [7], J = [8] */
    // A:=X1*Z2 mod N;
    mpres_mul(E->buf[0], P->x, Q->z, n);
    // B:=Z1*Z2 mod N;
    mpres_mul(E->buf[1], P->z, Q->z, n);
    // C:=Y1*X2 mod N;
    mpres_mul(E->buf[2], P->y, Q->x, n);
    // D:=Y1*Y2 mod N;
    mpres_mul(E->buf[3], P->y, Q->y, n);
    // E:=Z1*Y2 mod N;
    mpres_mul(E->buf[4], P->z, Q->y, n);
    // F:=E0[1]*X1*X2 mod N;
    mpres_mul(E->buf[5], P->x, Q->x, n);
    mpres_mul(E->buf[5], E->buf[5], E->a4, n);
    // Hisil
    // G := (D+B)*(A-C) mod N;
    mpres_add(E->buf[9], E->buf[3], E->buf[1], n);
    mpres_sub(E->buf[6], E->buf[0], E->buf[2], n);
    mpres_mul(E->buf[6], E->buf[6], E->buf[9], n);
    // H := (D-B)*(A+C) mod N;
    mpres_sub(E->buf[9], E->buf[3], E->buf[1], n);
    mpres_add(E->buf[7], E->buf[0], E->buf[2], n);
    mpres_mul(E->buf[7], E->buf[7], E->buf[9], n);
    // J := (D+F)*(A-E) mod N;
    mpres_add(E->buf[9], E->buf[3], E->buf[5], n);
    mpres_sub(E->buf[8], E->buf[0], E->buf[4], n);
    mpres_mul(E->buf[8], E->buf[8], E->buf[9], n);
    // K := (D-F)*(A+E) mod N;
    // this is the last use of A, so that K -> buf[0]
    mpres_sub(E->buf[9], E->buf[3], E->buf[5], n);
    mpres_add(E->buf[0], E->buf[0], E->buf[4], n);
    mpres_mul(E->buf[0], E->buf[0], E->buf[9], n);
    // X3 := G-H
    mpres_sub(R->x, E->buf[6], E->buf[7], n);
    // Y3 := K-J
    mpres_sub(R->y, E->buf[0], E->buf[8], n);
    // Z3 := (J+K-G-H-2*(B-F)*(C+E)) mod N;
    mpres_sub(E->buf[9], E->buf[1], E->buf[5], n);
    mpres_add(R->z, E->buf[2], E->buf[4], n);
    mpres_mul(R->z, R->z, E->buf[9], n);
    mpres_add(R->z, R->z, R->z, n);
    mpres_add(R->z, R->z, E->buf[7], n);
    mpres_add(R->z, R->z, E->buf[6], n);
    mpres_sub(R->z, E->buf[0], R->z, n);
    mpres_add(R->z, R->z, E->buf[8], n);
    if(mpz_sgn(R->x) == 0 && mpz_sgn(R->y) == 0 && mpz_sgn(R->z) == 0){
	// iff (X2:Y2:Z2)=(Z1:gamma^2*X1:gamma*Y1), gamma^3 = a
	fprintf(stderr, "GASP: X3, Y3 and Z3 are 0\n");
	exit(-1);
#if 0
	    // TODO: rewrite with above quantities!
	    X3:=(X1^2*Y2*Z2-X2^2*Y1*Z1) mod N;
	    // A*X1*Y2-C*X2*Z1 = A*U-C*V
	    Y3:=(Z1^2*X2*Y2-Z2^2*X1*Y1) mod N;
	    // E*Z1*X2-A*Z2*Y1 = E*V-A*W
	    Z3:=(Y1^2*X2*Z2-Y2^2*X1*Z1) mod N;
	    // C*Y1*Z2-E*Y2*X1 = C*W-E*U

	    // X3 =     Y1*(a*X1^3-Z1^3)
	    // Y3 = g^2*X1*(Z1^3-Y1^3)
	    // Z3 =   g*Z1*(Y1^3-Z1^3)

	    // can be made faster with a = aa^3, since then g = aa and we
	    // can share many things

#endif
    }
    return 1;
}

int
twisted_hessian_add(ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
{
    if(twisted_hessian_is_zero(P, E, n)){
	ell_point_set(R, Q, E, n);
	return 1;
    }
    else if(twisted_hessian_is_zero(Q, E, n)){
	ell_point_set(R, P, E, n);
	return 1;
    }
    else
	return twisted_hessian_plus(R, P, Q, E, n);
}

#if USE_ADD_SUB_CHAINS > 0
int
twisted_hessian_sub(ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
{
    int ret;

    twisted_hessian_negate(Q, E, n);
    ret = twisted_hessian_add(R, P, Q, E, n);
    twisted_hessian_negate(Q, E, n);
    return ret;
}
#endif

/* INPUT: a*x^3+y^3+1 = d*x*y
   OUTPUT: Y^2 = X^3+A*X+B
   If a=c^3, then curve isom to Hessian (c*x)^3+y^3+1=3*(d/(3*c))*(c*x)*y
   SIDE EFFECT: (x, y, c) <- (x_on_W, y_on_W, A_of_W)
 */
int
twisted_hessian_to_weierstrass(mpz_t f, mpres_t x, mpres_t y, mpres_t c, mpres_t d, mpmod_t n)
{
    int ret = ECM_NO_FACTOR_FOUND;
    mpres_t tmp;

#if DEBUG_ADD_LAWS >= 2
    printf("x_tH="); print_mpz_from_mpres(x, n); printf("\n");
    printf("y_tH="); print_mpz_from_mpres(y, n); printf("\n");
    printf("c_tH="); print_mpz_from_mpres(c, n); printf("\n");
    printf("d_tH="); print_mpz_from_mpres(d, n); printf("\n");
#endif
    mpres_init(tmp, n);
    mpres_mul_ui(tmp, c, 3, n);
    if(mpres_invert(tmp, tmp, n) == 0){
        mpres_gcd(f, tmp, n);
        ret = ECM_FACTOR_FOUND_STEP1;
    }
    else{
	mpres_mul(x, x, c, n);
	mpres_mul(c, tmp, d, n);
	/* from x^3+y^3+1=3*c*x*y to Weierstrass stuff */
	ret = hessian_to_weierstrass(f, x, y, c, n);
#if DEBUG_ADD_LAWS >= 2
	printf("A_W="); print_mpz_from_mpres(c, n); printf("\n");
	printf("x_W="); print_mpz_from_mpres(x, n); printf("\n");
	printf("y_W="); print_mpz_from_mpres(y, n); printf("\n");
#endif
    }
    mpres_clear(tmp, n);
    return ret;
}

/******************** generic ec's ********************/

void
ell_point_init(ell_point_t P, ell_curve_t E, mpmod_t n)
{
    mpres_init(P->x, n);
    mpres_init(P->y, n);
    mpres_init(P->z, n);
    if(E->type == ECM_EC_TYPE_WEIERSTRASS){
      if(E->law == ECM_LAW_AFFINE)
	mpz_set_ui(P->z, 1); /* humf */
      else if(E->law == ECM_LAW_HOMOGENEOUS)
	mpres_set_ui(P->z, 1, n);
    }
    else if(E->type == ECM_EC_TYPE_HESSIAN 
	    || E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
	mpres_set_ui(P->z, 1, n);
}

/* TODO: change this according to E->type */
void
ell_point_clear(ell_point_t P, ATTRIBUTE_UNUSED ell_curve_t E, mpmod_t n)
{
    mpres_clear(P->x, n);
    mpres_clear(P->y, n);
    mpres_clear(P->z, n);
}

#if DEBUG_ADD_LAWS >= 1
void
ell_point_print(ell_point_t P, ell_curve_t E, mpmod_t n)
{
    if(E->type == ECM_EC_TYPE_WEIERSTRASS)
	pt_w_print(P->x, P->y, P->z, E, n);
    else if(E->type == ECM_EC_TYPE_HESSIAN)
	hessian_print(P, E, n);
    else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
	twisted_hessian_print(P, E, n);
}
#endif

/* TODO: should depend on E->type... */
void
ell_point_set(ell_point_t Q, ell_point_t P,
	     ATTRIBUTE_UNUSED ell_curve_t E, ATTRIBUTE_UNUSED mpmod_t n)
{
    mpres_set(Q->x, P->x, n);
    mpres_set(Q->y, P->y, n);
    mpres_set(Q->z, P->z, n);
}

void
ell_curve_init(ell_curve_t E, int etype, int law, mpmod_t n)
{
    int i;

    E->type = etype;
    E->law = law;
    mpres_init(E->a1, n);
    mpres_init(E->a3, n);
    mpres_init(E->a2, n);
    mpres_init(E->a4, n);
    mpres_init(E->a6, n);
    mpres_set_ui(E->a1, 0, n);
    mpres_set_ui(E->a3, 0, n);
    mpres_set_ui(E->a2, 0, n);
    mpres_set_ui(E->a4, 0, n);
    mpres_set_ui(E->a6, 0, n);
    for(i = 0; i < EC_W_NBUFS; i++)
	mpres_init (E->buf[i], n);
}

void
ell_curve_init_set(ell_curve_t E, int etype, int law, mpres_t A, mpmod_t n)
{
    ell_curve_init(E, etype, law, n);
    mpres_set(E->a4, A, n);
}

void
ell_curve_set_z(ell_curve_t E, ell_curve_t zE, mpmod_t n)
{
    ell_curve_init(E, zE->type, zE->law, n);
    mpres_set_z(E->a1, zE->a1, n);
    mpres_set_z(E->a3, zE->a3, n);
    mpres_set_z(E->a2, zE->a2, n);
    mpres_set_z(E->a4, zE->a4, n);
    mpres_set_z(E->a6, zE->a6, n);
#if 0
    E->disc = zE->disc;
    if(E->disc != 0){
	mpres_init(E->sq[0], n);
	mpres_set_z(E->sq[0], zE->sq[0], n);
    }
#endif
}

void
ell_curve_clear(ell_curve_t E, mpmod_t n)
{
    int i;

    mpres_clear(E->a1, n);
    mpres_clear(E->a3, n);
    mpres_clear(E->a2, n);
    mpres_clear(E->a4, n);
    mpres_clear(E->a6, n);
    for(i = 0; i < EC_W_NBUFS; i++)
	mpres_clear (E->buf[i], n);
    /* TODO: case of sq */
}

#if DEBUG_ADD_LAWS >= 1
void
ell_curve_print(ell_curve_t E, mpmod_t n)
{
    if(E->type == ECM_EC_TYPE_WEIERSTRASS){
	printf("["); print_mpz_from_mpres(E->a1, n);
	printf(", "); print_mpz_from_mpres(E->a3, n);
	printf(", "); print_mpz_from_mpres(E->a2, n);
	printf(", "); print_mpz_from_mpres(E->a4, n);
	printf(", "); print_mpz_from_mpres(E->a6, n); printf("];\n");
    }
    else if(E->type == ECM_EC_TYPE_HESSIAN){
	printf("D:="); print_mpz_from_mpres(E->a4, n); printf(";\n");
	printf("E:=[D];\n");
    }
    else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN){
	printf("a:="); print_mpz_from_mpres(E->a4, n); printf(";\n");
	printf("d:="); print_mpz_from_mpres(E->a6, n); printf(";\n");
	printf("E:=[a, d];\n");
    }
}
#endif

/* OUTPUT: 1 if P = O_E, 0 otherwise. */
int
ell_point_is_zero(ell_point_t P, ell_curve_t E, mpmod_t n)
{
    if(E->type == ECM_EC_TYPE_WEIERSTRASS)
	return pt_w_is_zero(P->z, n);
    else if(E->type == ECM_EC_TYPE_HESSIAN)
	return hessian_is_zero(P, E, n);
    else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
	return twisted_hessian_is_zero(P, E, n);
    return ECM_ERROR;
}

void
ell_point_set_to_zero(ell_point_t P, ell_curve_t E, mpmod_t n)
{
    if(E->type == ECM_EC_TYPE_WEIERSTRASS)
	pt_w_set_to_zero(P, n);
    else if(E->type == ECM_EC_TYPE_HESSIAN)
	hessian_set_to_zero(P, E, n);
    else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
	twisted_hessian_set_to_zero(P, E, n);
}

int
ell_point_is_on_curve(ell_point_t P, ell_curve_t E, mpmod_t n)
{
    int ok = 1;

    if(ell_point_is_zero(P, E, n))
	return 1;
    if(E->type == ECM_EC_TYPE_WEIERSTRASS){
	mpres_t tmp1, tmp2;

	mpres_init(tmp1, n);
	mpres_init(tmp2, n);
	if(E->law == ECM_LAW_AFFINE){
	    /* y^2+a1*x*y+a3*y = x^3+a2*x^2+a4*x+a6? */
	    mpres_mul(tmp1, E->a1, P->x, n);
	    mpres_add(tmp1, tmp1, P->y, n);
	    mpres_add(tmp1, tmp1, E->a3, n);
	    mpres_mul(tmp1, tmp1, P->y, n);
	    
	    mpres_add(tmp2, E->a2, P->x, n);
	    mpres_mul(tmp2, tmp2, P->x, n);
	    mpres_add(tmp2, tmp2, E->a4, n);
	    mpres_mul(tmp2, tmp2, P->x, n);
	    mpres_add(tmp2, tmp2, E->a6, n);
	}
#if 0 // useless for the time being
	else{
	    /* y^2*z+a1*x*y*z+a3*y*z^2 = x^3+a2*x^2*z+a4*x*z^2+a6*z^3? */
	    /* y*z*(y+a1*x+a3*z) = ((x+a2*z)*x+a4*z^2)*x+a6*z^3? */
	    mpres_t tmp3;

	    mpres_mul(tmp1, E->a1, P->x, n);  /* a1*x */
	    mpres_add(tmp1, tmp1, P->y, n);   /* a1*x+y */
	    mpres_mul(tmp2, E->a3, P->z, n);  /* a3*z */
	    mpres_add(tmp1, tmp1, tmp2, n);   /* y+a1*x+a3*z */
	    mpres_mul(tmp1, tmp1, P->y, n);   /* y*(...) */
	    mpres_mul(tmp1, tmp1, P->z, n);   /* lhs */

	    mpres_init(tmp3, n);
	    mpres_mul(tmp2, E->a2, P->z, n);  /* a2*z */
	    mpres_add(tmp2, tmp2, P->x, n);   /* x+a2*z */
	    mpres_mul(tmp2, tmp2, P->x, n);   /* (x+a2*z)*x */
	    mpres_mul(tmp3, E->a4, P->z, n);  /* a4*z */
	    mpres_mul(tmp3, tmp3, P->z, n);   /* a4*z^2 */
	    mpres_add(tmp2, tmp2, tmp3, n);   /* (x+a2*z)*x+a4*z^2 */
	    mpres_mul(tmp2, tmp2, P->x, n);   /* (...)*x */
	    mpres_mul(tmp3, P->z, P->z, n);   /* z^2 */
	    mpres_mul(tmp3, tmp3, P->z, n);   /* z^3 */
	    mpres_mul(tmp3, tmp3, E->a6, n);  /* a6*z^3 */
	    mpres_add(tmp2, tmp2, tmp3, n);   /* rhs */
	    mpres_clear(tmp3, n);
	}
#endif
	ok = mpres_equal(tmp1, tmp2, n);

	mpres_clear(tmp1, n);
	mpres_clear(tmp2, n);
    }
    else if(E->type == ECM_EC_TYPE_HESSIAN){
	/* TODO */
    }
    else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN){
	/* TODO */
    }
    return ok;
}

#if DEBUG_ADD_LAWS >= 1
static void
ell_point_check(ell_point_t P, ell_curve_t E, mpmod_t n)
{
    if(ell_point_is_on_curve(P, E, n) == 0){
	printf("Point not on curve\n");
	printf("E:=");
	ell_curve_print(E, n);
	printf("P:=");
	pt_print(E, P, n);
	printf("\n");
	exit(-1);
    }
}
#endif

#if DEBUG_ADD_LAWS >= 1
int
ell_point_equal(ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
{
    int ret = 1;

    if(E->type == ECM_EC_TYPE_WEIERSTRASS){
	if(E->law == ECM_LAW_AFFINE)
	    return mpres_equal(P->x, Q->x, n) 
		   && mpres_equal(P->y, Q->y, n)
		   && mpres_equal(P->z, Q->z, n);
	else if(E->law == ECM_LAW_HOMOGENEOUS){
	    mpres_t tmp1, tmp2;

	    mpres_init(tmp1, n);
	    mpres_init(tmp2, n);
	    mpres_mul(tmp1, P->x, Q->z, n);
	    mpres_mul(tmp2, P->z, Q->x, n);
	    if(mpres_equal(tmp1, tmp2, n) == 0){
		printf("Px/Pz != Qx/Qz\n");
		ret = 0;
		exit(-1);
	    }
	    else{
		mpres_mul(tmp1, P->y, Q->z, n);
		mpres_mul(tmp2, P->z, Q->y, n);
		if(mpres_equal(tmp1, tmp2, n) == 0){
		    printf("Py/Pz != Qy/Qz\n");
		    ret = 0;
		    exit(-1);
		}
	    }
	    mpres_clear(tmp1, n);
	    mpres_clear(tmp2, n);
	}
    }
    return ret;
}
#endif

/* OUTPUT: 1 if everything ok, 0 otherwise */
int
ell_point_add(mpz_t f, ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
{
    if(E->type == ECM_EC_TYPE_WEIERSTRASS)
	return pt_w_add(f, R->x, R->y, R->z, P->x, P->y, P->z, 
			Q->x, Q->y, Q->z, n, E);
    else if(E->type == ECM_EC_TYPE_HESSIAN)
	return hessian_add(R, P, Q, E, n);
    else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
	return twisted_hessian_add(R, P, Q, E, n);
    else
	return ECM_ERROR;
}

#if USE_ADD_SUB_CHAINS > 0
/* R <- P-Q */
int
ell_point_sub(mpz_t f, ell_point_t R, ell_point_t P, ell_point_t Q, ell_curve_t E, mpmod_t n)
{
    if(E->type == ECM_EC_TYPE_WEIERSTRASS)
	return pt_w_sub(f, R->x, R->y, R->z, P->x, P->y, P->z,
			Q->x, Q->y, Q->z, n, E);
    else if(E->type == ECM_EC_TYPE_HESSIAN)
	return hessian_sub(R, P, Q, E, n);
    else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
	return twisted_hessian_sub(R, P, Q, E, n);
    else
	return ECM_ERROR;
}
#endif

int
ell_point_duplicate(mpz_t f, ell_point_t R, ell_point_t P, ell_curve_t E, mpmod_t n)
{
#if DEBUG_ADD_LAWS >= 2
    printf("E:=");
    ell_curve_print(E, n);
#endif
    if(E->type == ECM_EC_TYPE_WEIERSTRASS)
	return pt_w_duplicate(f, R->x, R->y, R->z, P->x, P->y, P->z, n, E);
    else if(E->type == ECM_EC_TYPE_HESSIAN)
	return hessian_duplicate(R, P, E, n);
    else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
	return twisted_hessian_duplicate(R, P, E, n);
    else
	return ECM_ERROR;
}

void
ell_point_negate(ell_point_t P, ell_curve_t E, mpmod_t n)
{
#if DEBUG_ADD_LAWS >= 2
    printf("P:="); ell_point_print(P, E, n); printf(";\n");
#endif
    if(ell_point_is_zero(P, E, n) == 0){
	if(E->type == ECM_EC_TYPE_WEIERSTRASS){
	    if(E->law == ECM_LAW_HOMOGENEOUS){
		/* FIXME: does not work for complete equation! */
		mpres_neg(P->y, P->y, n);
	    }
	    else if(E->law == ECM_LAW_AFFINE){
		/* (-P).y = -P.y-a1*P.x-a3 */
		if(mpz_sgn(E->a1) != 0
		   || mpz_sgn(E->a3) != 0
		   || mpz_sgn(E->a2) != 0){ /* FIXME */
		    printf("GROUMF\n");
		    exit(-1);
		}
		mpres_neg(P->y, P->y, n);
	    }
	}
#if USE_ADD_SUB_CHAINS > 0
	else if(E->type == ECM_EC_TYPE_HESSIAN)
	    hessian_negate(P, E, n);
	else if(E->type == ECM_EC_TYPE_TWISTED_HESSIAN)
	    twisted_hessian_negate(P, E, n);
#endif
    }
#if DEBUG_ADD_LAWS >= 2
    printf("neg(P):="); ell_point_print(P, E, n); printf(";\n");
#endif
}

/* Q <- [e]*P
   Return value: 0 if a factor is found, and the factor is in Q->x,
                 1 otherwise.
*/
int
ell_point_mul_plain (mpz_t f, ell_point_t Q, mpz_t e, ell_point_t P, ell_curve_t E, mpmod_t n)
{
  size_t l;
  int negated = 0, status = 1;
  ell_point_t P0;

  if(ell_point_is_zero(P, E, n) != 0){
      ell_point_set(Q, P, E, n);
      return 1;
  }

  if (mpz_sgn (e) == 0)
    {
      ell_point_set_to_zero(Q, E, n);
      return 1;
    }

  if (mpz_sgn (e) < 0)
    {
      negated = 1;
      mpz_neg (e, e);
      ell_point_negate(P, E, n); /* since the point is non-zero */
    }

  if (mpz_cmp_ui (e, 1) == 0){
      ell_point_set(Q, P, E, n);
      goto ell_point_mul_plain_end;
  }

  l = mpz_sizeinbase (e, 2) - 1; /* l >= 1 */

  ell_point_init(P0, E, n);
  ell_point_set(P0, P, E, n);

#if DEBUG_ADD_LAWS >= 2
  printf("P:="); ell_point_print(P, E, n); printf(";\n");
#endif
  while (l-- > 0)
    {
#if DEBUG_ADD_LAWS >= 2
	printf("P0:="); ell_point_print(P0, E, n); printf(";\n");
#endif
	if(ell_point_duplicate (f, P0, P0, E, n) == 0)
	  {
	    status = 0;
	    break;
	  }
#if DEBUG_ADD_LAWS >= 2
	printf("Rdup:="); ell_point_print(P0, E, n); printf(";\n");
	printf("dup:=ProjEcmDouble(P0, E, N); ProjEcmEqual(dup, Rdup, N);\n");
#endif
	if (mpz_tstbit (e, l))
	  {
	      if(ell_point_add (f, P0, P0, P, E, n) == 0)
	      {
		status = 0;
		break;
	      }
#if DEBUG_ADD_LAWS >= 2
	      printf("Radd:="); ell_point_print(P0, E, n); printf(";\n");
	      printf("Padd:=ProjEcmAdd(P, Rdup, E, N); ProjEcmEqual(Padd, Radd, N);\n");
#endif
	  }
    }

  ell_point_set(Q, P0, E, n);
  ell_point_clear(P0, E, n);
ell_point_mul_plain_end:

  /* Undo negation to avoid changing the caller's e value */
  if (negated){
    mpz_neg (e, e);
    ell_point_negate(P, E, n);
  }
  return status;
}

int
ell_point_mul(mpz_t f, ell_point_t Q, mpz_t e, ell_point_t P, ell_curve_t E, mpmod_t n)
{
#if 1 /* keeping it simple */
    return ell_point_mul_plain(f, Q, e, P, E, n);
#else
    return ell_point_mul_add_sub(f, Q, e, P, E, n);
#endif
}