// Test program for ed336 scalar point multiplication
// Uses Bernstein et al.s point multiplication method from Curve41417 paper
// Cache safety thanks to ed25519
// Fully Tested and debugged
// g++ -O3 ed336.cpp -o ed336
// M.Scott 23/09/2015

// 32-bit version

#include <iostream>
#include <ctime>
#include <inttypes.h>

#define WINDOW 4
#define PANES  84

using namespace std;

typedef int64_t type64;
typedef int32_t type32;

static const type32 bot28bits = 0xfffffff;

#include <stdint.h>

// w=x+y
void gadd(type32 *x,type32 *y,type32 *w)
{
	w[0]=x[0]+y[0];
	w[1]=x[1]+y[1];
	w[2]=x[2]+y[2];
	w[3]=x[3]+y[3];
	w[4]=x[4]+y[4];
	w[5]=x[5]+y[5];

	w[6]=x[6]+y[6];
	w[7]=x[7]+y[7];
	w[8]=x[8]+y[8];
	w[9]=x[9]+y[9];
	w[10]=x[10]+y[10];
	w[11]=x[11]+y[11];
}

// w=x-y
void gsub(type32 *x,type32 *y,type32 *w)
{
	w[0]=x[0]-y[0];
	w[1]=x[1]-y[1];
	w[2]=x[2]-y[2];
	w[3]=x[3]-y[3];
	w[4]=x[4]-y[4];
	w[5]=x[5]-y[5];
	w[6]=x[6]-y[6];
	w[7]=x[7]-y[7];
	w[8]=x[8]-y[8];
	w[9]=x[9]-y[9];
	w[10]=x[10]-y[10];
	w[11]=x[11]-y[11];
}

// w-=x
void gdec(type32 *x,type32 *w)
{
	w[0]-=x[0];
	w[1]-=x[1];
	w[2]-=x[2];
	w[3]-=x[3];
	w[4]-=x[4];
	w[5]-=x[5];
	w[6]-=x[6];
	w[7]-=x[7];
	w[8]-=x[8];
	w[9]-=x[9];
	w[10]-=x[10];
	w[11]-=x[11];
}

// w=x
void gcopy(type32 *x,type32 *w)
{
	w[0]=x[0];
	w[1]=x[1];
	w[2]=x[2];
	w[3]=x[3];
	w[4]=x[4];
	w[5]=x[5];
	w[6]=x[6];
	w[7]=x[7];
	w[8]=x[8];
	w[9]=x[9];
	w[10]=x[10];
	w[11]=x[11];
}

// w*=2
void gmul2(type32 *w)
{
	w[0]*=2;
	w[1]*=2;
	w[2]*=2;
	w[3]*=2;
	w[4]*=2;
	w[5]*=2;
	w[6]*=2;
	w[7]*=2;
	w[8]*=2;
	w[9]*=2;
	w[10]*=2;
	w[11]*=2;
}

// w-=2*x
void gsb2(type32 *x,type32 *w)
{
	w[0]-=2*x[0];
	w[1]-=2*x[1];
	w[2]-=2*x[2];
	w[3]-=2*x[3];
	w[4]-=2*x[4];
	w[5]-=2*x[5];
	w[6]-=2*x[6];
	w[7]-=2*x[7];
	w[8]-=2*x[8];
	w[9]-=2*x[9];
	w[10]-=2*x[10];
	w[11]-=2*x[11];
}

// reduce w - Short Coefficient Reduction
void scr(type32 *w)
{
	type32 t0,t1,t2;
	t0=w[0]&bot28bits;

	t1=w[1]+(w[0]>>28);
	w[1]=t1&bot28bits;

	t2=w[2]+(t1>>28);
	w[2]=t2&bot28bits;

	t1=w[3]+(t2>>28);
	w[3]=t1&bot28bits;

	t2=w[4]+(t1>>28);
	w[4]=t2&bot28bits;

	t1=w[5]+(t2>>28);
	w[5]=t1&bot28bits;

	t2=w[6]+(t1>>28);
	w[6]=t2&bot28bits;

	t1=w[7]+(t2>>28);
	w[7]=t1&bot28bits;

	t2=w[8]+(t1>>28);
	w[8]=t2&bot28bits;

	t1=w[9]+(t2>>28);
	w[9]=t1&bot28bits;
	
	t2=w[10]+(t1>>28);
	w[10]=t2&bot28bits;

	t1=w[11]+(t2>>28);
	w[11]=t1&bot28bits;

	w[0]=t0+3*(t1>>28);

}

// multiply w by a constant, w*=i

void gmuli(type32 *w,int i)
{
	type64 t;

	t=(type64)w[0]*i;
	w[0]=((type32)t)&bot28bits;

	t=(type64)w[1]*i+(t>>28);
	w[1]=((type32)t)&bot28bits;

	t=(type64)w[2]*i+(t>>28);
	w[2]=((type32)t)&bot28bits;

	t=(type64)w[3]*i+(t>>28);
	w[3]=((type32)t)&bot28bits;

	t=(type64)w[4]*i+(t>>28);
	w[4]=((type32)t)&bot28bits;

	t=(type64)w[5]*i+(t>>28);
	w[5]=((type32)t)&bot28bits;

	t=(type64)w[6]*i+(t>>28);
	w[6]=((type32)t)&bot28bits;

	t=(type64)w[7]*i+(t>>28);
	w[7]=((type32)t)&bot28bits;

	t=(type64)w[8]*i+(t>>28);
	w[8]=((type32)t)&bot28bits;

	t=(type64)w[9]*i+(t>>28);
	w[9]=((type32)t)&bot28bits;

	t=(type64)w[10]*i+(t>>28);
	w[10]=((type32)t)&bot28bits;

	t=(type64)w[11]*i+(t>>28);
	w[11]=((type32)t)&bot28bits;


	w[0]+=3*(t>>28);

}

// z=x^2

void gsqr(type32 *x,type32*z)
{
	type64 t0,t1;

	t0=2*((type64)x[0]*x[11]+(type64)x[1]*x[10]+(type64)x[2]*x[9]+(type64)x[3]*x[8]+(type64)x[4]*x[7]+(type64)x[5]*x[6]);
	z[11]=((type32) t0)&bot28bits;

	t1=(type64)x[0]*x[0]+3*(2*((type64)x[1]*x[11]+(type64)x[2]*x[10]+(type64)x[3]*x[9]+(type64)x[4]*x[8]+(type64)x[5]*x[7])+(type64)x[6]*x[6]+(t0>>28));     
	z[0]=((type32) t1)&bot28bits;

	t0=2*((type64)x[0]*x[1]+3*((type64)x[2]*x[11]+(type64)x[3]*x[10]+(type64)x[4]*x[9]+(type64)x[5]*x[8]+(type64)x[6]*x[7]))+(t1>>28);
	z[1]=((type32) t0)&bot28bits;

	t1=(type64)x[0]*(2*x[2])+(type64)x[1]*x[1]+3*(2*((type64)x[3]*x[11]+(type64)x[4]*x[10]+(type64)x[5]*x[9]+(type64)x[6]*x[8])+(type64)x[7]*x[7])+(t0>>28);
	z[2]=((type32) t1)&bot28bits;

	t0=2*((type64)x[0]*x[3]+(type64)x[1]*x[2]+3*((type64)x[4]*x[11]+(type64)x[5]*x[10]+(type64)x[6]*x[9]+(type64)x[7]*x[8]))+(t1>>28);
	z[3]=((type32) t0)&bot28bits;

	t1=2*((type64)x[0]*x[4]+(type64)x[1]*x[3])+(type64)x[2]*x[2]+3*(2*((type64)x[5]*x[11]+(type64)x[6]*x[10]+(type64)x[7]*x[9])+(type64)x[8]*x[8])+(t0>>28);
	z[4]=((type32) t1)&bot28bits;

	t0=2*((type64)x[0]*x[5]+(type64)x[1]*x[4]+(type64)x[2]*x[3]+3*((type64)x[6]*x[11]+(type64)x[7]*x[10]+(type64)x[8]*x[9]))+(t1>>28);
	z[5]=((type32) t0)&bot28bits;

	t1=2*((type64)x[0]*x[6]+(type64)x[1]*x[5]+(type64)x[2]*x[4])+(type64)x[3]*x[3]+3*(2*((type64)x[7]*x[11]+(type64)x[8]*x[10])+(type64)x[9]*x[9])+(t0>>28);
	z[6]=((type32) t1)&bot28bits;

	t0=2*((type64)x[0]*x[7]+(type64)x[1]*x[6]+(type64)x[2]*x[5]+(type64)x[3]*x[4]+3*((type64)x[8]*x[11]+(type64)x[9]*x[10]))+(t1>>28);
	z[7]=((type32) t0)&bot28bits;

	t1=2*((type64)x[0]*x[8]+(type64)x[1]*x[7]+(type64)x[2]*x[6]+(type64)x[3]*x[5])+(type64)x[4]*x[4]+3*((type64)x[9]*(2*x[11])+(type64)x[10]*x[10])+(t0>>28);
	z[8]=((type32) t1)&bot28bits;

	t0=2*((type64)x[0]*x[9]+(type64)x[1]*x[8]+(type64)x[2]*x[7]+(type64)x[3]*x[6]+(type64)x[4]*x[5]+(type64)x[10]*(3*x[11]))+(t1>>28);
	z[9]=((type32) t0)&bot28bits;

	t1=2*((type64)x[0]*x[10]+(type64)x[1]*x[9]+(type64)x[2]*x[8]+(type64)x[3]*x[7]+(type64)x[4]*x[6])+(type64)x[5]*x[5]+(type64)x[11]*(3*x[11])+(t0>>28);
	z[10]=((type32) t1)&bot28bits;

	t0=z[11]+(t1>>28);
	z[11]=t0&bot28bits;
	z[0]+=3*(t0>>28);

}



// z=x*y - Granger's method

void gmul(type32 *x,type32 *y,type32 *z)
{
	type64 t0,t1,t2;
	type64 a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11;
	type64 b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10;
	a0=(type64)x[0]*y[0];
	a1=(type64)x[1]*y[1];
	a2=(type64)x[2]*y[2];
	a3=(type64)x[3]*y[3];
	a4=(type64)x[4]*y[4];
	a5=(type64)x[5]*y[5];
	a6=(type64)x[6]*y[6];
	a7=(type64)x[7]*y[7];
	a8=(type64)x[8]*y[8];
	a9=(type64)x[9]*y[9];
	a10=(type64)x[10]*y[10];
	a11=(type64)x[11]*y[11];

	b10=a11+a10;
	b9=b10+a9;
	b8=b9+a8;
	b7=b8+a7;
	b6=b7+a6;
	b5=b6+a5;
	b4=b5+a4;
	b3=b4+a3;
	b2=b3+a2;
	b1=b2+a1;
	b0=b1+a0;

	t0=b0-(type64)(x[0]-x[11])*(y[0]-y[11])-(type64)(x[1]-x[10])*(y[1]-y[10])-(type64)(x[2]-x[9])*(y[2]-y[9])-(type64)(x[3]-x[8])*(y[3]-y[8])-(type64)(x[4]-x[7])*(y[4]-y[7])-(type64)(x[5]-x[6])*(y[5]-y[6]);
	z[11]=((type32)t0)&bot28bits;

	t1=a0+3*(b1-(type64)(x[1]-x[11])*(y[1]-y[11])-(type64)(x[2]-x[10])*(y[2]-y[10])-(type64)(x[3]-x[9])*(y[3]-y[9])-(type64)(x[4]-x[8])*(y[4]-y[8])-(type64)(x[5]-x[7])*(y[5]-y[7])+(t0>>28));
	z[0]=((type32)t1)&bot28bits;

	a0+=a1;
	t0=a0-(type64)(x[0]-x[1])*(y[0]-y[1])+3*(b2-(type64)(x[2]-x[11])*(y[2]-y[11])-(type64)(x[3]-x[10])*(y[3]-y[10])-(type64)(x[4]-x[9])*(y[4]-y[9])-(type64)(x[5]-x[8])*(y[5]-y[8])-(type64)(x[6]-x[7])*(y[6]-y[7]))+(t1>>28);
	z[1]=((type32)t0)&bot28bits;
	a0+=a2;

	t1=a0-(type64)(x[0]-x[2])*(y[0]-y[2])+3*(b3-(type64)(x[3]-x[11])*(y[3]-y[11])-(type64)(x[4]-x[10])*(y[4]-y[10])-(type64)(x[5]-x[9])*(y[5]-y[9])-(type64)(x[6]-x[8])*(y[6]-y[8]))+(t0>>28);
	z[2]=((type32)t1)&bot28bits;
	a0+=a3;

	t0=a0-(type64)(x[0]-x[3])*(y[0]-y[3])-(type64)(x[1]-x[2])*(y[1]-y[2])+3*(b4-(type64)(x[4]-x[11])*(y[4]-y[11])-(type64)(x[5]-x[10])*(y[5]-y[10])-(type64)(x[6]-x[9])*(y[6]-y[9])-(type64)(x[7]-x[8])*(y[7]-y[8]))+(t1>>28);
	z[3]=((type32)t0)&bot28bits;
	a0+=a4;

	t1=a0-(type64)(x[0]-x[4])*(y[0]-y[4])-(type64)(x[1]-x[3])*(y[1]-y[3])+3*(b5-(type64)(x[5]-x[11])*(y[5]-y[11])-(type64)(x[6]-x[10])*(y[6]-y[10])-(type64)(x[7]-x[9])*(y[7]-y[9]))+(t0>>28);
	z[4]=((type32)t1)&bot28bits;
	a0+=a5;

	t0=a0-(type64)(x[0]-x[5])*(y[0]-y[5])-(type64)(x[1]-x[4])*(y[1]-y[4])-(type64)(x[2]-x[3])*(y[2]-y[3])+3*(b6-(type64)(x[6]-x[11])*(y[6]-y[11])-(type64)(x[7]-x[10])*(y[7]-y[10])-(type64)(x[8]-x[9])*(y[8]-y[9]))+(t1>>28);
	z[5]=((type32)t0)&bot28bits;
	a0+=a6;

	t1=a0-(type64)(x[0]-x[6])*(y[0]-y[6])-(type64)(x[1]-x[5])*(y[1]-y[5])-(type64)(x[2]-x[4])*(y[2]-y[4])+3*(b7-(type64)(x[7]-x[11])*(y[7]-y[11])-(type64)(x[8]-x[10])*(y[8]-y[10]))+(t0>>28);
	z[6]=((type32)t1)&bot28bits;
	a0+=a7;

	t0=a0-(type64)(x[0]-x[7])*(y[0]-y[7])-(type64)(x[1]-x[6])*(y[1]-y[6])-(type64)(x[2]-x[5])*(y[2]-y[5])-(type64)(x[3]-x[4])*(y[3]-y[4])+3*(b8-(type64)(x[8]-x[11])*(y[8]-y[11])-(type64)(x[9]-x[10])*(y[9]-y[10]))+(t1>>28);
	z[7]=((type32)t0)&bot28bits;
	a0+=a8;

	t1=a0-(type64)(x[0]-x[8])*(y[0]-y[8])-(type64)(x[1]-x[7])*(y[1]-y[7])-(type64)(x[2]-x[6])*(y[2]-y[6])-(type64)(x[3]-x[5])*(y[3]-y[5])+3*(b9-(type64)(x[9]-x[11])*(y[9]-y[11]))+(t0>>28);
	z[8]=((type32)t1)&bot28bits;
	a0+=a9;

	t0=a0-(type64)(x[0]-x[9])*(y[0]-y[9])-(type64)(x[1]-x[8])*(y[1]-y[8])-(type64)(x[2]-x[7])*(y[2]-y[7])-(type64)(x[3]-x[6])*(y[3]-y[6])-(type64)(x[4]-x[5])*(y[4]-y[5])+3*(b10-(type64)(x[10]-x[11])*(y[10]-y[11]))+(t1>>28);
	z[9]=((type32)t0)&bot28bits;
	a0+=a10;

	t1=a0-(type64)(x[0]-x[10])*(y[0]-y[10])-(type64)(x[1]-x[9])*(y[1]-y[9])-(type64)(x[2]-x[8])*(y[2]-y[8])-(type64)(x[3]-x[7])*(y[3]-y[7])-(type64)(x[4]-x[6])*(y[4]-y[6])+3*a11+(t0>>28);
	z[10]=((type32)t1)&bot28bits;

	t0=z[11]+(t1>>28);  
	z[11]=t0&bot28bits;  

	z[0]+=3*(t0>>28);
}

// Inverse x = 1/x = x^(p-2) mod p

void ginv(type32 *x)
{
	int i;
	type32 xc[12],x11[12],w[12],t1[12],t2[12];
	gsqr(x,w);       // =x^2
	gmul(w,x,xc);     // xc=x^3
	gsqr(w,t1);
	gsqr(t1,x11);      
	gmul(x11,x,t2);   // xc=x^9    
	gmul(t2,w,x11);   // x11=x^11
	
	gsqr(x11,w);
	gmul(w,t2,t1);     // t1=x^31

	gcopy(t1,w);
	gsqr(w,t2);
	gsqr(t2,w);
	gsqr(w,t2);
	gsqr(t2,w);
	gsqr(w,t2);
	gmul(t2,t1,x);   // =x^(2^10-1)
	gcopy(x,t2);

	gsqr(t2,w);
	gsqr(w,t2);
	gmul(t2,xc,w);  
	gcopy(w,xc);    // xc=x^(2^12-1)
	
	gsqr(t2,w);
	gsqr(w,t2);
	gsqr(t2,w);
	gsqr(w,t2);
	gsqr(t2,w);
	gsqr(w,t2);
	gsqr(t2,w);
	gsqr(w,t2);

	gmul(t2,x,t1);    // t1=x^(2^20-1)

	gcopy(t1,x);
	for (i=0;i<10;i++)
	{
		gsqr(t1,w);
		gsqr(w,t1);
	}
	gmul(t1,x,t2);   //t2=x^(2^40-1)

	gcopy(t2,x);
	for (i=0;i<20;i++)
	{
		gsqr(t2,w);
		gsqr(w,t2);
	}
	gmul(t2,x,t1);   //t1=x^(2^80-1)

	gcopy(t1,x);
	for (i=0;i<40;i++)
	{
		gsqr(t1,w);
		gsqr(w,t1);
	}
	gmul(t1,x,t2);   //t2=x^(2^160-1)

	gcopy(t2,x);
	for (i=0;i<80;i++)
	{
		gsqr(t2,w);
		gsqr(w,t2);
	}

	gmul(t2,x,t1);   //t1=x^(2^320-1)
	for (i=0;i<6;i++)
	{
		gsqr(t1,w);
		gsqr(w,t1);
	}

	gmul(t1,xc,t2);
	for (i=0;i<2;i++)
	{
		gsqr(t2,w);
		gsqr(w,t2);
	}
	
	gmul(t2,x11,x);
}

// Point Structure

typedef struct {
type32 x[12];
type32 y[12];
type32 z[12];
type32 t[12];
} ECp;

// P+=P

void double_1(ECp *P)
{
	type32 a[12],b[12],e[12],f[12],g[12],h[12];
	gsqr(P->x,a);
	gsqr(P->y,b);
	gcopy(P->t,e);
	gmul2(e);
	gadd(a,b,g);
	gcopy(g,f); f[0]-=2;
	gsub(a,b,h);
	gmul(e,f,P->x); 
	gmul(g,h,P->y); 
	gsqr(g,P->z);
	gsb2(g,P->z);
	gmul(e,h,P->t);  
	scr(P->z);
}

void double_2(ECp *P)
{
	type32 a[12],b[12],c[12],e[12],f[12],g[12],h[12];
	gsqr(P->x,a);
	gsqr(P->y,b); 
	gsqr(P->z,c); gmul2(c); 
	gadd(P->x,P->y,g); gsqr(g,e); gdec(a,e); gdec(b,e); 
	gadd(a,b,g); 
	gsub(g,c,f); scr(f);
	gsub(a,b,h); 
	gmul(e,f,P->x); 
	gmul(g,h,P->y);
	gmul(f,g,P->z); 
}

void double_3(ECp *P)
{
	type32 a[12],b[12],c[12],e[12],f[12],g[12],h[12];
	gsqr(P->x,a);
	gsqr(P->y,b);
	gsqr(P->z,c); gmul2(c);
	gadd(P->x,P->y,g); gsqr(g,e); gdec(a,e); gdec(b,e); 
	gadd(a,b,g);
	gsub(g,c,f); scr(f);
	gsub(a,b,h);
	gmul(e,f,P->x); 
	gmul(g,h,P->y); 
	gmul(f,g,P->z);
	gmul(e,h,P->t); 
}

//P+=Q;

void add_1(ECp *Q,ECp *P)
{
	type32 a[12],b[12],c[12],d[12],e[12],f[12],g[12],h[12];
	gmul(P->x,Q->x,a);
	gmul(P->y,Q->y,b);
	gmul(P->t,Q->t,c);
	gsub(P->z,c,f);  /* reverse sign if d is negative */
	gadd(P->z,c,g);
	gsub(b,a,h);
	gadd(P->x,P->y,c); gadd(Q->x,Q->y,d); gmul(c,d,e); gdec(a,e); gdec(b,e);
	gmul(e,f,P->x); 
	gmul(g,h,P->y); 
	gmul(f,g,P->z); 
	gmul(e,h,P->t); 
}

void add_2(ECp *Q,ECp *P)
{
	type32 a[12],b[12],c[12],d[12],e[12],f[12],g[12],h[12];
	gmul(P->x,Q->x,a);
	gmul(P->y,Q->y,b);
	gmul(P->t,Q->t,c);
	gmul(P->z,Q->z,d);
	gsub(d,c,f);  /* reverse sign if d is negative */
	gadd(d,c,g);
	gsub(b,a,h);
	gadd(P->x,P->y,c); gadd(Q->x,Q->y,d); gmul(c,d,e); gdec(a,e); gdec(b,e);
	gmul(e,f,P->x); 
	gmul(g,h,P->y); 
	gmul(f,g,P->z); 
}

//P=0

void inf(ECp *P)
{
	int i;
	for (i=0;i<12;i++)
		P->x[i]=P->y[i]=P->z[i]=P->t[i]=0;
	P->y[0]=P->z[0]=1;
}

// Initialise P

void init(type32 *x,type32 *y,ECp *P)
{
	int i;
	for (i=0;i<12;i++)
	{
		P->x[i]=x[i];
		P->y[i]=y[i];
		P->z[i]=0;
	}
	P->z[0]=1;
	gmul(x,y,P->t);
}

//P=Q

void copy(ECp *Q,ECp *P)
{
	int i;
	for (i=0;i<12;i++)
	{
		P->x[i]=Q->x[i];
		P->y[i]=Q->y[i];
		P->z[i]=Q->z[i];
		P->t[i]=Q->t[i];
	}
}

// P=-Q

void neg(ECp *Q,ECp *P)
{
	int i;
	for (i=0;i<12;i++)
	{
		P->x[i]=-Q->x[i]; 
		P->y[i]=Q->y[i];
		P->z[i]=Q->z[i];
		P->t[i]=-Q->t[i]; 
	}
}

/* Make Affine */

void norm(ECp *P)
{
	type32 w[12],t[12];
	gcopy(P->z,w);
	ginv(w);
	gmul(P->x,w,t); scr(t); gcopy(t,P->x);
	gmul(P->y,w,t); scr(t); gcopy(t,P->y);
	gmul(P->z,w,t); scr(t); gcopy(t,P->z);
	gmul(P->t,w,t); scr(t); gcopy(t,P->t);
}

/* Precomputation */

void precomp(ECp *P,ECp W[])
{
	inf(&W[0]);
	copy(P,&W[1]); gmuli(W[1].t,11111);
	copy(P,&W[2]); double_1(&W[2]);
	copy(&W[2],&W[3]); add_1(&W[1],&W[3]);

	copy(&W[2],&W[4]); double_3(&W[4]);
	copy(&W[4],&W[5]); add_1(&W[1],&W[5]);
	copy(&W[3],&W[6]); double_3(&W[6]);
	copy(&W[6],&W[7]); add_1(&W[1],&W[7]);
	copy(&W[4],&W[8]); double_3(&W[8]);

/* premultiply t parameter by curve constant */

	gmuli(W[2].t,11111);
	gmuli(W[3].t,11111);
	gmuli(W[4].t,11111);
	gmuli(W[5].t,11111);
	gmuli(W[6].t,11111);
	gmuli(W[7].t,11111);
	gmuli(W[8].t,11111);
}

/* Window of width 4 */

void window(ECp *Q,ECp *P)
{
	double_2(P);
	double_2(P);
	double_2(P);
	double_3(P);
	add_2(Q,P);
}

void output(ECp *P)
{
	cout << "x[0]= " << hex << P->x[0] << endl;
	cout << "x[1]= " << hex << P->x[1] << endl;
	cout << "x[2]= " << hex << P->x[2] << endl;
	cout << "x[3]= " << hex << P->x[3] << endl;
	cout << "x[4]= " << hex << P->x[4] << endl;
	cout << "x[5]= " << hex << P->x[5] << endl;
	cout << "x[6]= " << hex << P->x[6] << endl;
	cout << "x[7]= " << hex << P->x[7] << endl;
	cout << "x[8]= " << hex << P->x[8] << endl;
	cout << "x[9]= " << hex << P->x[9] << endl;
	cout << "x[10]= " << hex << P->x[10] << endl;
	cout << "x[11]= " << hex << P->x[11] << endl;
	cout << endl;

	cout << "y[0]= " << hex << P->y[0] << endl;
	cout << "y[1]= " << hex << P->y[1] << endl;
	cout << "y[2]= " << hex << P->y[2] << endl;
	cout << "y[3]= " << hex << P->y[3] << endl;
	cout << "y[4]= " << hex << P->y[4] << endl;
	cout << "y[5]= " << hex << P->y[5] << endl;
	cout << "y[6]= " << hex << P->y[6] << endl;
	cout << "y[7]= " << hex << P->y[7] << endl;
	cout << "y[8]= " << hex << P->y[8] << endl;
	cout << "y[9]= " << hex << P->y[9] << endl;
	cout << "y[10]= " << hex << P->y[10] << endl;
	cout << "y[11]= " << hex << P->y[11] << endl;
	cout << endl;
}

/*
Constant time table look-up - borrowed from ed25519 
*/

void fe_cmov(type32 f[],type32 g[],int ib)
{
  type32 b=ib;
  b=-b;
  f[0]^=(f[0]^g[0])&b;
  f[1]^=(f[1]^g[1])&b;
  f[2]^=(f[2]^g[2])&b;
  f[3]^=(f[3]^g[3])&b;
  f[4]^=(f[4]^g[4])&b;
  f[5]^=(f[5]^g[5])&b;
  f[6]^=(f[6]^g[6])&b;
  f[7]^=(f[7]^g[7])&b;
  f[8]^=(f[8]^g[8])&b;
  f[9]^=(f[9]^g[9])&b;
  f[10]^=(f[10]^g[10])&b;
  f[11]^=(f[11]^g[11])&b;
}

static void cmov(ECp *w,ECp *u,int b)
{
  fe_cmov(w->x,u->x,b);
  fe_cmov(w->y,u->y,b);
  fe_cmov(w->z,u->z,b);
  fe_cmov(w->t,u->t,b);
}

// return 1 if b==c, no branching
static int equal(int b,int c)
{
	int x=b^c;
	x-=1;  // if x=0, x now -1
	return ((x>>31)&1);
}

static void select(ECp *T,ECp W[],int b)
{
  ECp MT; 
  int m=b>>31;
  int babs=(b^m)-m;

  cmov(T,&W[0],equal(babs,0));  // conditional move
  cmov(T,&W[1],equal(babs,1));
  cmov(T,&W[2],equal(babs,2));
  cmov(T,&W[3],equal(babs,3));
  cmov(T,&W[4],equal(babs,4));
  cmov(T,&W[5],equal(babs,5));
  cmov(T,&W[6],equal(babs,6));
  cmov(T,&W[7],equal(babs,7));
  cmov(T,&W[8],equal(babs,8)); 

  neg(T,&MT);  // minus t
  cmov(T,&MT,m&1);
}

/* Point Multiplication - exponent is 333 bits */

void mul(int *w,ECp *P)
{
	ECp W[1+(1<<(WINDOW-1))],S[2],Q;
	int i,j,m;
	precomp(P,W);

	copy(&W[w[PANES-1]],P);  
	for (i=PANES-2;i>=0;i--)
	{
		select(&Q,W,w[i]);
		window(&Q,P);
	}
	norm(P); 
}

// #define TEST  /* define to multiply by group order */

int main()
{
	uint64_t bef,aft;
	int w[PANES];
	int ii,lpz=10000;
	type32 xs[12],ys[12];
	ECp P;

/* Base point on Curve */

	xs[0]=0xc;
	xs[1]=0x0;
	xs[2]=0x0;
	xs[3]=0x0;
	xs[4]=0x0;
	xs[5]=0x0;
	xs[6]=0x0;
	xs[7]=0x0;
	xs[8]=0x0;
	xs[9]=0x0;
	xs[10]=0x0;
	xs[11]=0x0;

	ys[1]= 0xBEC6850;
	ys[0]= 0x5FE8632;
	ys[3]= 0x5650CA0;
	ys[2]= 0x365DB12;
	ys[5]= 0x1C7EF43;
	ys[4]= 0x5B6DB5D;
	ys[7]= 0x53D1B14;
	ys[6]= 0xB46C381;
	ys[9]= 0xE18E1C1;
	ys[8]= 0x61D0078;
	ys[11]=0x0C0DC61;
	ys[10]=0x6B56502;

// random multiplier arranged in signed windows of size 6
#ifndef TEST
#if WINDOW==4

w[0]= 5; w[1]= -5; w[2]= 2; w[3]= -5; w[4]= -8;
w[5]= -3; w[6]= 5; w[7]= -2; w[8]= -4; w[9]= 3; 
w[10]= -5; w[11]= -4; w[12]= -6; w[13]= -4; w[14]= -7;
w[15]= -8; w[16]= -5; w[17]= -2; w[18]= 4; w[19]= -2;
w[20]= -7; w[21]= 2; w[22]= 7; w[23]= -5; w[24]= -8;
w[25]= -4; w[26]= -1; w[27]= -8; w[28]= 7; w[29]= -6;
w[30]= -7; w[31]= -5; w[32]= -4; w[33]= -5; w[34]= -3;
w[35]= 6; w[36]= -2; w[37]= -3; w[38]= -3; w[39]= -7;
w[40]= -4; w[41]= -2; w[42]= 3; w[43]= 7; w[44]= -1;
w[45]= 0; w[46]= -1; w[47]= -5; w[48]= -1; w[49]= -4;
w[50]= -4; w[51]= 5; w[52]= -1; w[53]= -6; w[54]= 7;
w[55]= 4; w[56]= 6; w[57]= -4; w[58]= -2; w[59]= -2; 
w[60]= 2; w[61]= -4; w[62]= 5; w[63]= 1; w[64]= 3;
w[65]= -6; w[66]= -7; w[67]= -2; w[68]= 1; w[69]= 7;
w[70]= 2; w[71]= -6; w[72]= 6; w[73]= 3; w[74]= -5;
w[75]= 6; w[76]= -7; w[77]= -6; w[78]= -4; w[79]= 6;
w[80]= -6; w[81]= -8; w[82]= -3; w[83]= 1;


#endif

// Group Order - for testing
#else

#if WINDOW==4
w[0]= 5; w[1]= 0; w[2]= -8; w[3]= -5; w[4]= 0;
w[5]= -6; w[6]= -1; w[7]= 4; w[8]= 7; w[9]= -7;
w[10]= 6; w[11]= -7; w[12]= 0; w[13]= 3; w[14]= -5;
w[15]= -8; w[16]= -5; w[17]= -5; w[18]= -4; w[19]= 4;
w[20]= -7; w[21]= 0; w[22]= -8; w[23]= -7; w[24]= -4;
w[25]= 7; w[26]= -3; w[27]= -4; w[28]= 1; w[29]= -4;
w[30]= 1; w[31]= 5; w[32]= -8; w[33]= -6; w[34]= -5;
w[35]= 0; w[36]= 6; w[37]= 1; w[38]= 4; w[39]= 1;
w[40]= 7; w[41]= 0; w[42]= 0; w[43]= 0; w[44]= 0;
w[45]= 0; w[46]= 0; w[47]= 0; w[48]= 0; w[49]= 0;
w[50]= 0; w[51]= 0; w[52]= 0; w[53]= 0; w[54]= 0;
w[55]= 0; w[56]= 0; w[57]= 0; w[58]= 0; w[59]= 0;
w[60]= 0; w[61]= 0; w[62]= 0; w[63]= 0; w[64]= 0;
w[65]= 0; w[66]= 0; w[67]= 0; w[68]= 0; w[69]= 0;
w[70]= 0; w[71]= 0; w[72]= 0; w[73]= 0; w[74]= 0;
w[75]= 0; w[76]= 0; w[77]= 0; w[78]= 0; w[79]= 0;
w[80]= 0; w[81]= 0; w[82]= 0; w[83]= 2;
#endif

#endif

	for (ii=0;ii<lpz;ii++)
	{
		init(xs,ys,&P);
		mul(w,&P);
	}

	output(&P);

	return 0;
}