#include <stdio.h>
#include<string.h>
#include<xgpio.h>
#include "platform.h"
#include "xparameters.h"
#include "xio.h"
#include<string.h>
#include<stdlib.h>
#include<assert.h>



typedef unsigned long long field_t;
#define FIELD_M     7
#define FIELD_W     3
#define MSB			31
#define FIELD_P ((1UL << MSB)  - 1)
#define BYTESIZE      (MSB * FIELD_M / 8 + 1)
#define MULALGO0
#define L (2*FIELD_M+3)
#define LOOP 30
#define ITER 900
/*** field2n.h ***/

#define WORDSIZE	(sizeof(int)*8)
#define NUMBITS		113
#define TYPE2
/*#undef TYPE2 */

#ifdef TYPE2
#define field_prime	((NUMBITS<<1)+1)
#else
#define field_prime (NUMBITS+1)
#endif

#define	NUMWORD		(NUMBITS/WORDSIZE)
#define UPRSHIFT	(NUMBITS%WORDSIZE)
#define MAXLONG		(NUMWORD+1)

#define MAXBITS		(MAXLONG*WORDSIZE)
#define MAXSHIFT	(WORDSIZE-1)
#define MSB_BIG		(1L<<MAXSHIFT)

#define UPRBIT		(1L<<(UPRSHIFT-1))
#define UPRMASK		(~(-1L<<UPRSHIFT))
#define SUMLOOP(i)	for(i=0; i<MAXLONG; i++)

#define LONGWORD	(field_prime/WORDSIZE)
#define LONGSHIFT	((field_prime-1)%WORDSIZE)
#define LONGBIT		(1L<<(LONGSHIFT-1))
#define LONGMASK         (~(-1L<<LONGSHIFT))

typedef	short int INDEX;

typedef unsigned long ELEMENT;

typedef struct {
	ELEMENT 	e[MAXLONG];
}  FIELD2N;

typedef struct {
	ELEMENT e[LONGWORD+1];
} CUSTFIELD;
//---------------------------------------------------------------------------
/*** bigint.h ***/



/*  The following are used by the multiprecision integer package.
	This really is very crude.  See J. W. Crenshaw "Programmers
	Toolbox", Embedded Systems Programming Dec. 1996 for why you
	want to do this in assembler.
*/

#define	HALFSIZE	(WORDSIZE/2)
#define	HIMASK		(-1L<<HALFSIZE)
#define LOMASK		(~HIMASK)
#define CARRY		(1L<<HALFSIZE)
#define MSB_HW		(CARRY>>1)
#define	INTMAX		(4*MAXLONG-1)
#define MAXSTRING	(MAXLONG*WORDSIZE/3)

#define	INTLOOP(i)	for(i=INTMAX;i>=0;i--)

typedef struct {
	ELEMENT		hw[4*MAXLONG];
}  BIGINT;
//-----------------------------------------------------------------------


//---------------------------------------------------------------------
/*** bigint.c ***/

/*  This is a very crude large integer package.  The purpose is to teach, not to
be efficient.  See freelip or miracl for very efficient code or check out the
commercial packages like Macsyma, Maple or the other similar products.  The basis
of these routines come from J.W. Crenshaw and several articles in his "Programmers
Toolbox" in Embedded Systems Programming magazine from Dec. 1996 thru Sept. 1997
*/


/*  clear all bits in a large integer storage block.  */

void int_null(BIGINT *a)
{
	INDEX i;

	INTLOOP (i) a->hw[i] = 0;
}

/*  copy one BIGINT block to another  */

void int_copy(BIGINT *a,BIGINT *b)
{
	INDEX i;

	INTLOOP (i) b->hw[i] = a->hw[i];
	//xil_printf("here");
}

/*  for use in the far distant future, convert a packed field to a large
	integer.  Does a simple expansion.  The large integer is 4 times bigger
	to accomodate multiplication (once!).
*/

void field_to_int(FIELD2N *a,BIGINT * b)
{
	INDEX	i, j;

	int_null( b);
	for (i=NUMWORD; i>=0; i--)
	{
		j = INTMAX - ((NUMWORD - i)<<1);
		b->hw[j] = a->e[i] & LOMASK;
		j--;
		b->hw[j] = (a->e[i] & HIMASK) >> HALFSIZE;
	}
}

/*  Pack a BIGINT variable back into a FIELD2N size one.  */

void int_to_field(BIGINT * a,FIELD2N * b)
{
	INDEX i, j;

	SUMLOOP(i)
	{
		j = (i + MAXLONG) << 1;
		b->e[i] = a->hw[j+1] | (a->hw[j] << HALFSIZE);
	}
}

/*  Negate a BIGINT in place.  Each half word is complemented, then we add 1  */

void int_neg(BIGINT * a)
{
	INDEX i;

	INTLOOP(i) a->hw[i] = ~a->hw[i] & LOMASK;
	INTLOOP(i)
	{
		a->hw[i]++;
		if (a->hw[i] & LOMASK) break;
		a->hw[i] &= LOMASK;
	}
}

/*  add two BIGINTS to get a third.  c = a + b
	Unlike the polynomial or ONB math, c can be one of a or b
*/

void int_add(BIGINT * a,BIGINT * b,BIGINT * c)
{
	INDEX i;
	ELEMENT ec;

	ec = 0;
	INTLOOP (i)
	{
	/* add previous carry bit to each term  */
		ec = a->hw[i] + b->hw[i] + (ec >> HALFSIZE);
		c->hw[i] = ec & LOMASK;
	}
}

/*  subtract two BIGINTS, c = a - b == a + (-b).
	as in addition, c can point to a or b and still works
*/

void int_sub(BIGINT * a,BIGINT * b,BIGINT * c)
{
	BIGINT negb;

	int_copy( b, &negb);
	int_neg( &negb);
	int_add( a, &negb, c);
}

/*  multiply two BIGINTs to get a third.
	Do NOT attempt to do 2 multiplies in a row without a division in between.
	You may get an overflow and there is no provision in this code to return
	an error condition for that.  See more advanced packages for correct way
	to do this.  c can *not* be one of a or b, it must be a separate storage
	location.
*/

void int_mul(BIGINT * a,BIGINT * b,BIGINT * c)
{
	ELEMENT		ea, eb, mul;
	INDEX		i, j, k;
	BIGINT		sum;

	int_null(c);

	for ( i = INTMAX; i > INTMAX/2; i--)
	{
		ea = a->hw[i];
		int_null( &sum);
		for ( j = INTMAX; j > INTMAX/2; j--)
		{
			eb = b->hw[j];
			k = i + j - INTMAX;
			mul = ea * eb + sum.hw[k];
			sum.hw[k] = mul & LOMASK;
			sum.hw[k-1] = mul >> HALFSIZE;
		}
		int_add( &sum, c, c);
	}
}

/*  unsigned divide.  Input full sized numerator (top),
	half sized denominator (bottom).
	Output half sized quotient and half sized remainder.
	Exceptionally crude but works ok for basics, error
	conditions return zero results.
*/

void int_div(BIGINT * top,BIGINT * bottom,BIGINT * quotient,BIGINT * remainder)
{
	BIGINT d, e;
	ELEMENT mask;
	INDEX	l, m, n, i, j;

/*  first step, initialize counters to most significant
	bit position in top and bottom.
*/
	int_copy( top, &d);
	int_copy( bottom, &e);
	l = (INTMAX + 1) * HALFSIZE;
	for( i=0; i<=INTMAX; i++)
	{
		if (!d.hw[i]) l -= HALFSIZE;
		else break;
	}

	mask = 1L << (HALFSIZE-1);
	for ( j=0; j<HALFSIZE; j++)
	{
		if ( !(d.hw[i] & mask))
		{
			l--;
			mask >>= 1;
		}
		else break;
	}

/*  same thing for bottom, compute msb position  */

	m = (INTMAX + 1) * HALFSIZE;
	for( i=0; i<=INTMAX; i++)
	{
		if (!e.hw[i]) m -= HALFSIZE;
		else break;
	}

	mask = 1L << (HALFSIZE-1);
	for ( j=0; j<HALFSIZE; j++)
	{
		if ( !(e.hw[i] & mask))
		{
			m--;
			mask >>= 1;
		}
		else break;
	}

/*  check for error inputs, does not check for zero, so is
	actually incorrect.
*/
	if (!m)			/*  x/1 = x */
	{
		int_copy( top, quotient);
		int_null( remainder);
	}

	if (!l | (l<m))			/*  1/x = 0 */
	{
		int_null( quotient);
		int_copy( bottom, remainder);
	}

/*  next step, shift bottom over to align msb with top msb  */

	n = l - m;
	i = n;
	while ( i > HALFSIZE )
	{
		for (j=0; j<INTMAX; j++) e.hw[j] = e.hw[j+1];
		i -= HALFSIZE;
		e.hw[INTMAX] = 0;
	}

	mask = 0;
	while ( i > 0 )
	{
		INTLOOP (j)
		{
			e.hw[j] = (e.hw[j] << 1) | mask;
			mask = e.hw[j] & CARRY ? 1 : 0;
			e.hw[j] &= LOMASK;
		}
		i--;
	}

/*  main division loop.  check to see if we can subtract shifted bottom
	from what's left on top.  If we can, set that bit in quotient and do
	subtract.  if we can't, just shift bottom right and repeat until only
	remainder is left.
*/

	int_null( quotient);
	while ( n>=0)
	{
		i = INTMAX - l/HALFSIZE;
		j = INTMAX - n/HALFSIZE;
		while ( (d.hw[i] == e.hw[i]) && ( i<INTMAX) ) i++;
		if ( d.hw[i] >= e.hw[i] )
		{
			int_sub( &d, &e, &d);
			mask = 1L << ( n%HALFSIZE );
			quotient->hw[j] |= mask;
		}
		INTLOOP(j)
		{
			if (j) mask = ( e.hw[j-1] & 1) ? CARRY : 0;
			else mask = 0;
			e.hw[j] = (e.hw[j] | mask) >> 1;
		}
		n--;
		l--;
	}

	//int_copy( &d, remainder);
	INTLOOP (i) remainder->hw[i] = d.hw[i];
	//xil_printf("here");
}
int pow(int a,int b)
{
	int temp,i;
	temp=1;
	for(i=1;i<=b;i++)
		temp=temp*a;
	return temp;
}
void str_to_bigint(char *instring,BIGINT *b,int base)
{
	char *pointer;
	int i,x,j,p,sum,loop,y;
	int temp[4];
	int length;
	length = strlen(instring);
	INTLOOP (i)
		b->hw[i]=0;
	loop=0;j=0;y=length-1;
	pointer=instring;
	while(y>=0)
	{
		pointer=instring+y;
		if(*pointer >= 48 && *pointer <= 57)
			temp[j] = *pointer - 48;
		else if(*pointer >=65 && *pointer <= 70)
			temp[j] = (*pointer - 65)+10;
		if((j + 1)%4 == 0)
		{
			p = 0;
			sum = 0;
			for(x=0; x<4 ;x++ )
			{
				p = pow(base,x);
				sum = sum + temp[x]*p;
			}
			//printf("sum -->%d\n",sum);
			b->hw[INTMAX-(loop++)]=sum;
			for(i=0;i<4;i++)
				temp[i]=0;
		}
		y--;
		j = (j+1)%4;
	}
	if(j>0)
	{
		p=0; sum=0;
		for(x=0; x<4 ; x++ )
		{
			p = pow(base,x);
			sum = sum + temp[x]*p;
		}
		b->hw[INTMAX-loop]=sum;
	}

}

int is_neg(BIGINT *a)
{
	int isneg;
	isneg=(a->hw[0]& 0x8000) >> (HALFSIZE-1);
	return isneg;
}

int is_zero(BIGINT *c)
{
	int i;
	int checkifzero = 1;
	INTLOOP (i)
	{
		if(c->hw[i]!=0)
		{
			checkifzero = 0;
			break;
		}
	}
	return checkifzero;
}
int int_comp(BIGINT *a,BIGINT *b)
{
	BIGINT c;
	int cmp;
	int_sub(a,b,&c);
	if(is_neg(&c)== 1)
		cmp = -1;
	else if(is_zero(&c))
		cmp=0;
	else
		cmp=1;
	return cmp;
}

ELEMENT get_int(BIGINT *a)
{
	ELEMENT result;
	result = a->hw[INTMAX] + (65536 * a->hw[INTMAX-1]);
	return result;
}


//--------------------------------------------------------------------

long pbytwo = (FIELD_P + 1) >> 1;;

static field_t con_one[FIELD_M] = {1, 0x68768cde, 0x5f6ca339, 0x5b9f5d37, 0x354d412e, 0x5f57d0ba, 0x47d860c5};
/*****constant for frobenius map procedure*/
/***               w^floor(ip/m)                 ****/
static field_t constant[FIELD_M][FIELD_M]={{1,1,1,1,1,1,1},
                                           {1,0x68768cde ,0x5f6ca339 ,0x5b9f5d37 ,0x354d412e ,0x5f57d0ba ,0x47d860c5},
                                           {1,0x5f6ca339, 0x354d412e, 0x47d860c5, 0x68768cde, 0x5b9f5d37, 0x5f57d0ba},
                                           {1,0x5b9f5d37, 0x47d860c5, 0x5f6ca339, 0x5f57d0ba, 0x68768cde, 0x354d412e},
                                           {1,0x354d412e, 0x68768cde, 0x5f57d0ba, 0x5f6ca339, 0x47d860c5, 0x5b9f5d37},
                                           {1,0x5f57d0ba, 0x5b9f5d37, 0x68768cde, 0x47d860c5, 0x354d412e, 0x5f6ca339},
                                           {1,0x47d860c5, 0x5f57d0ba, 0x354d412e, 0x5b9f5d37, 0x5f6ca339, 0x68768cde}};
/*****constant for frobenius map procedure*/
/***               ip mod m                 ****/
static int modvalue[FIELD_M][FIELD_M]={{0,0,0,0,0,0,0},
                                       {0,1,2,3,4,5,6},
                                       {0,1,2,3,4,5,6},
                                       {0,1,2,3,4,5,6},
                                       {0,1,2,3,4,5,6},
                                       {0,1,2,3,4,5,6},
                                       {0,1,2,3,4,5,6}};
/***       a= - 3,b= - 212       ***/
field_t CURVECONST_A[FIELD_M] ={2147483644,0,0,0,0,0,0};
field_t CURVECONST_B[FIELD_M] ={ 2147483435,0,0,0,0,0,0};

//-------Base Point **/
field_t CURVE_PT_X[FIELD_M] ={ 1265161714,1539753404,708434454,1832286394,1731473315,1007441496,1511235565};
field_t CURVE_PT_Y[FIELD_M] ={373113286, 1537269188,1290843001,1831071987,1155799427,2060698802,573726986};
//------(- Base Point)**/
field_t CURVE_mPT_X[FIELD_M] ={ 1265161714,1539753404,708434454,1832286394,1731473315,1007441496,1511235565};
field_t CURVE_mPT_Y[FIELD_M] ={1774370361, 610214459,856640646,316411660,991684220,86784845,1573756661};

typedef struct _point_t{
	field_t x[FIELD_M];
	field_t y[FIELD_M];
	field_t z[FIELD_M];   // used only in projective coordinate system, else ignored
}point_t;

point_t * phipoint[FIELD_M]={NULL,NULL,NULL,NULL,NULL,NULL,NULL};
point_t * phipointn[FIELD_M]={NULL,NULL,NULL,NULL,NULL,NULL,NULL};

FILE * fp1,*fp2;
field_t init_point[2][7]={
{571919611,	1353644485,	241420551,	809195974,	1381361500,	1291363926,	1811897767},
{364754508,	635192694,	697479810,	782863320,	785828434,	1563635474,	750235896}
};



unsigned int timestamp(void)
{
    unsigned int bottom;
    unsigned int top;
    asm volatile("xorl %%eax,%%eax\n cpuid \n" ::: "%eax", "%ebx", "%ecx", "%edx"); // flush pipeline
    asm volatile("rdtsc\n" : "=a" (bottom), "=d" (top) );                           // read rdtsc
    asm volatile("xorl %%eax,%%eax\n cpuid \n" ::: "%eax", "%ebx", "%ecx", "%edx"); // flush pipeline again
    return bottom;
}
//FILE * px,* py;
//--------------------------------------

////////////////////////////////////////
void copy_element(field_t *d, field_t *s)
{
	memcpy(d, s, FIELD_M*sizeof(field_t));
}

void init_element(field_t a[5])
{
	memset(a, 0, FIELD_M*sizeof(field_t));
}

static void fillzeros(point_t *dest)
{
	memset(dest, 0, sizeof(point_t));
}

void copypoint(point_t *dest, point_t *src)
{
	memcpy(dest, src, sizeof(point_t));
}

field_t baseadd(field_t a, field_t b)
{
	field_t c1 = a + b;
	field_t c2 = c1 - FIELD_P;


	return ((c1 >= FIELD_P) ? c2 : c1);
}

static field_t basesub(field_t a, field_t b)
{
	field_t c;
	c = (baseadd(a, FIELD_P - b));

	return c;
}

/*#if defined(MULALGO0)
static void sanitycheck()
{
	if(FIELD_P >= (1U<<31)-1){
		print("Wrong Multiplier Configured\n");
		//exit(-1);
	}
}*/
/*
 * This function can be used when the field cardinality is less than 2^32
 */
static field_t basemul(field_t a, field_t b)
{
	unsigned long long pr, q1;

	pr = a * (unsigned long long)b;
	q1 = pr / FIELD_P;

//	printf("%016llx\n", pr);

	return (pr - (q1 * FIELD_P));
}

static field_t basemulx2(field_t a)
{
	field_t msb = a >> (MSB - 1);
	a = a << 1;
	return (a & FIELD_P) + msb;
}

/* Constant multiplication with 3 */
inline static field_t basemulx3(field_t a)
{
	field_t a2;
	field_t msb = a >> (MSB - 1);
	a2 = ((a << 1) + msb) & FIELD_P;
	return baseadd(a2 , a);
}

/* Constant multiplication with 4 */
inline static field_t basemulx4(field_t a)
{
	field_t msb = a >> (MSB - 2);
	a = a << 2;
	return (a & FIELD_P) + msb;
}

/* Constant multiplication with 8 */
inline static field_t basemulx8(field_t a)
{
	field_t msb = a >> (MSB - 3);
	a = a << 3;
	return (a & FIELD_P) + msb;
}

#define ISEVEN(x)     ((x & 1ULL) == 0)

/* Binary Euclidean Algorithm to find
   the Inverse of element a in GF(FIELD_P)
*/
static field_t baseinverse(field_t a)
{
	field_t u, v, e;
	field_t x1, x2;
	field_t Pby2;

	if(a == 0) return 0;
	Pby2 = (FIELD_P + 1) >> 1;
	u = a; v = FIELD_P;
	x1 = 1ULL; x2 = 0ULL;

	while(u != 1 && v != 1){
		while(ISEVEN(u)){
			u = u >> 1;
			if (ISEVEN(x1))
				x1 = x1 >> 1;
			else{
				x1 = ((x1 - 1) >> 1) + Pby2;
			}
		}
		while(ISEVEN(v)){
			v = v >> 1;
			if (ISEVEN(x2))
				x2 = x2 >> 1;
			else{
				x2 = ((x2 - 1) >> 1) + Pby2;
			}
		}
		if(u >= v){
			u = basesub(u, v);
			x1 = basesub(x1, x2);

		}else{
			v = basesub(v, u);
			x2 = basesub(x2,x1);
		}
	}
	if (u == 1){
		return x1 ;
	}
	else{
		return x2;
	}
}



static void oefadd(field_t *res, field_t *a, field_t *b)
{
	int i;

	for(i=0; i<FIELD_M; ++i){
		res[i] = baseadd(a[i], b[i]);
	}
}

static void oefsub(field_t *res, field_t *a, field_t *b)
{
	int i;

	for(i=0; i<FIELD_M; ++i){
		res[i] = basesub(a[i], b[i]);
	}
}

/*
	OEF Classical multiplier. Also does the modular reduction.
	Computes : res = ( a * b ) mod P
 */
static void oefmul(field_t *res, field_t *a, field_t *b)
{
	int i, i1, j;
	field_t upper, tmp;

	for(i=0; i<FIELD_M; ++i){
		res[i] = 0;
		for(j=0; j<i+1; ++j){
			tmp = basemul(a[j], b[i-j]);
			res[i] = baseadd(res[i], tmp);
		}

		i1 = i + FIELD_M;
		upper = 0;
		for(j=i1-FIELD_M+1; j<FIELD_M; ++j){
				tmp = basemul(a[j], b[i1-j]);
				upper = baseadd(upper, tmp);
		}
		res[i] = baseadd(res[i], basemul(upper, FIELD_W));
	}
}

/* OEF Square - is there a better algorithm ? */
static void oefsqr(field_t *res, field_t *a)
{
	oefmul(res, a, a);
}

/* OEF Multiplication with constant 2 */
static void oefmulx2(field_t *res, field_t *a)
{
	int i;
	for(i=0; i<FIELD_M; ++i){
		res[i] = basemulx2(a[i]);
	}
}

/* OEF Multiplication with constant 3 */
static void oefmulx3(field_t *res, field_t *a)
{
	int i;
	for(i=0; i<FIELD_M; ++i){
		res[i] = basemulx3(a[i]);
	}
}
//98081778997765566787677889999000766545677887887765548891
/* OEF Multiplication with constant 4 */
static void oefmulx4(field_t *res, field_t *a)
{
	int i;
	for(i=0; i<FIELD_M; ++i){
		res[i] = basemulx4(a[i]);
	}
}

/* OEF Multiplication with constant 8 */
static void oefmulx8(field_t *res, field_t *a)
{
	int i;
	for(i=0; i<FIELD_M; ++i){
		res[i] = basemulx8(a[i]);
	}
}

static void frobeniusmap(field_t *b, field_t *a)
{
	int i;

	for(i=0; i<FIELD_M; ++i){
		b[i] = basemul(a[i], con_one[i]);
	}
}

static void oefinverse(field_t *ainv, field_t *a)
{
  	
	field_t b1[FIELD_M], b2[FIELD_M];
	field_t cinv;
	int i;

	init_element(b1);
	copy_element(b2, a);
  	
	for(i = 0; i < FIELD_M-1; ++i){
		frobeniusmap(b1, b2);
		oefmul(b2, b1, a);
	}


	cinv = baseinverse(b2[0]);
	init_element(b2);
	b2[0] = cinv;

	oefmul(ainv, b2, b1);


#ifdef __DEBUG__
	oefmul(b1, a, ainv);
	printfield(b1);
	if (b1[0] != 1) goto fatalerror;
	for(i = 1; i < FIELD_M; ++i){
		if(b1[i] != 0) goto fatalerror;
	}
	return;
fatalerror:
	print("oefinv: Error in finding inverse (at i)\n");
	printfield(a);
	printfield(ainv);
	exit(-1);
#endif

}
static void pt_double(point_t *pt)
{
	field_t a[FIELD_M], b[FIELD_M], c[FIELD_M], d[FIELD_M];
	field_t temp1[FIELD_M], temp2[FIELD_M];
	////int seal_doub=1;


	/* A <- y1^2 */
	oefsqr(a, pt->y);           // a = y1^2

	////seal_doub=!seal_doub;
 ////   XGpio_DiscreteWrite(&pin, 1, seal_doub);

	/* B <- [[4 . x1] . A] */
	oefmulx4(temp1, pt->x);     // 4.x1
	oefmul(b, temp1, a);        // temp1.a

	/* D <- (3 . X1^2) + (CURVE_A . z1^4) */
	oefsqr(c, pt->x);           // x1^2
	oefmulx3(temp1, c);         // 3.x1^2
	oefsqr(temp2, pt->z);       // z1^2
	oefsqr(c, temp2);           // z1^4
	oefmul(temp2, CURVECONST_A, c);        // CURVE_A.z1^4
	oefadd(d, temp2, temp1);    // 3.z1^2 + CURVE_A.z1^4

	/* C <- [8 . [A^2]] */
	oefsqr(temp1, a);           // a^2
	oefmulx8(c, temp1);         // 8.a^2

	/* X3 <- [[D^2] - [2B]] */
	oefsqr(temp1, d);           // d^2
	oefmulx2(temp2, b);         // 2.b
	oefsub(pt->x, temp1, temp2); // d^2 - 2b

	/* Z3 <- [2[Y1 . Z1]] */
	oefmul(temp1, pt->y, pt->z);// y1.z1
	oefmulx2(pt->z, temp1);     // 2.y1.z1

	/* Y3 <- [[D(B - X3)] - C] */
	oefsub(temp1, b, pt->x);    // b - x3
	oefmul(temp2, temp1, d);    // d(b - x3)
	oefsub(pt->y, temp2, c);    // d(b - x3) - c

}
static void pt_add(point_t *pt1, point_t *pt2)
{
	field_t a[FIELD_M], b[FIELD_M], c[FIELD_M], d[FIELD_M], e[FIELD_M];
	field_t f[FIELD_M], g[FIELD_M], h[FIELD_M], i[FIELD_M];


	oefsqr(a, pt1->z);          // a = z1^2

	/*seal=!seal;
	    XGpio_DiscreteWrite(&pin, 1, seal);*/


	oefmul(b, pt1->z, a);       // b = z1.a
	oefmul(c, pt2->x, a);       // c = x2.a
	oefmul(d, pt2->y, b);       // d = y2.b
	oefsub(e, c, pt1->x);       // e = c - x1
	oefsub(f, d, pt1->y);       // f = d - y1
	oefsqr(g, e);               // g = e^2
	oefmul(h, g, e);            // h = g.e
	oefmul(i, g, pt1->x);       // i = x1.g

	// x3 = f^2 - (h + 2i)
	oefsqr(a, f);             // f^2
	oefsub(b, a, h);          // f^2 - h
	oefmulx2(c, i);           // 2i
	oefsub(pt1->x, b, c);     // f^2 - h - 2i

	// y3 = f(i - x3) - y1.h
	oefsub(a, i, pt1->x);     // i - x3
	oefmul(b, f, a);          // f(i-x3)
	oefmul(c, h, pt1->y);     // y1.h
	oefsub(pt1->y, b, c);     // f(i-x3) - y1.h

	// z3 = z1.e
	copy_element(d, pt1->z);
	oefmul(pt1->z, e, d); // e.z1
}
/***             Computing   Pi=phi^i(P)                    ***/
//This is precomputed before scalar multiplication starts
void computephi(field_t *a,int p,field_t * res)
{
	int i,pi_i;
	field_t b[FIELD_M];

	for(i=0;i<FIELD_M;i++)
	{
		b[i]=basemul(a[i],constant[p][i]);
		res[i]=b[i];
	}
}

//scalar K is expanded as a poynomial in phi
//finding out the coefficients ui of the equation k=Summation(ui.phi^i(P))
/******************************************************
void basephiexpansion(mpz_t k,field_t * uphi)
{
	int i;
	mpz_t x,y,v,uphi_z[L],field_p,zero;
	mpz_init(zero);
	mpz_init(x);
	mpz_init(y);
	mpz_init(v);
	mpz_init(field_p);
	mpz_set(x,k);
	mpz_set_str(zero,"0",10);
	mpz_set_str(field_p,"2147483647",10);


	for(i=0;i<L;i++)
	{
		uphi[i]=0;
		mpz_init(uphi_z[i]);
	}
	i=0;

	while(1)
	{
		if(!mpz_cmp(x,zero) && !mpz_cmp(y,zero))
			break;

		mpz_tdiv_qr(v,uphi_z[i],x,field_p);
		uphi[i]=mpz_get_ui(uphi_z[i]);
		mpz_add(x,v,y);
		mpz_neg(y,v);
		i++;
	}
	/*
	for(i=0;i<FIELD_M;i++)
	{
	    printf("\nuphi-->%d,  %lld",i,uphi[i]);
	}
}*/
//the Ui values are reduced to Di values
/************************************************
void basephiopti(field_t * u,field_t * d)
{
	int i;
	for(i=0;i<FIELD_M;i++)
		d[i]=0;
	for(i=0;i<FIELD_M;i++)
	{
		if((i+FIELD_M)<L && (i+(2*FIELD_M))<L)
			d[i]=u[i]+u[i+FIELD_M]+u[i+(2*FIELD_M)];
		else if((i+FIELD_M)<L)
			d[i]=u[i]+u[i+FIELD_M];
		else
			d[i]=u[i];
	}
}
*********************************************************/

void generateRandomly(BIGINT *k,int * r1,int * r2,BIGINT * k_1,BIGINT * k_2)
{
	BIGINT s_1,s_2,k1,k2,r_1,r_2,temp1,temp2,nem1,nem2,denom,rem;
	BIGINT a,b,c,d;
	int_null(&s_2);
	int_copy(k,&s_1);
	int_null(&a);
	int_null(&b);
	int_null(&c);
	int_null(&d);
	
	do
	{
		a.hw[INTMAX] = rand()%1024;
		b.hw[INTMAX] = rand()%1024;
		c.hw[INTMAX] = rand()%1024;
		d.hw[INTMAX] = rand()%1024;
		
		int_mul(&a,&b,&temp1);
		int_mul(&c,&d,&temp2);
		int_sub(&temp1,&temp2,&denom);
	}while(is_zero(&denom));
	
	
	int_mul(&d,&s_1,&temp1);
	int_mul(&b,&s_2,&temp2);
	int_sub(&temp1,&temp2,&nem1);
	
	int_mul(&c,&s_1,&temp1);
	int_mul(&a,&s_2,&temp2);
	int_sub(&temp2,&temp1,&nem2);
	
	int_div(&nem1,&denom,&k1,&rem);
	int_div(&nem2,&denom,&k2,&rem);
	
	int_mul(&a,&k1,&temp1);
	int_mul(&b,&k2,&temp2);
	int_add(&temp1,&temp2,k_1);
	
        int_mul(&c,&k1,&temp1);
	int_mul(&d,&k2,&temp2);
	int_add(&temp1,&temp2,k_2);
	
	int_sub(&s_1,k_1,&r_1);
	int_sub(&s_2,k_2,&r_2);
	
	*r1 = r_1.hw[INTMAX];
	*r2 = r_2.hw[INTMAX];
	
}

void expansion(BIGINT * k,long * uphi)
{
	int i,checkifneg,x1,j;
	BIGINT t,x,y,v,tv,uphi_z[L],fieldp,zero;
    int_null(&zero);
	int_null(&x);
	int_null(&y);
	int_null(&v);
	int_null(&t);
	int_null(&tv);
	int_null(&zero);
	int_null(&fieldp);
	int_copy(k,&x);

	//str_to_bigint("7FFFFFFF",&fieldp,16);
	//str_to_bigint("14485",&t,16);
	fieldp.hw[INTMAX] = 0x7FFF;
	fieldp.hw[INTMAX-1] = 0xFFFF;
	
	t.hw[INTMAX] = 0x1448 ;
	t.hw[INTMAX-1] = 0x5;
	
	
	//
	//for(i=0;i<=INTMAX;i++)
	//xil_printf("%d",fieldp.hw[INTMAX]);
	//xil_printf("%d",fieldp.hw[INTMAX-1]);

	for(i=0;i<L;i++)
	{
		uphi[i]=0;
		int_null(&uphi_z[i]);
	}
	i=0;
/*******************/

	

	while(!(is_zero(&x) && is_zero(&y)))
	{

		checkifneg=is_neg(&x);
		if(checkifneg)
			int_neg(&x);

		int_div(&x,&fieldp,&v,&uphi_z[i]);

		uphi[i]=get_int(&uphi_z[i]);

		if(checkifneg)
			uphi[i]= - uphi[i];

		if(uphi[i] > pbytwo)
			uphi[i] = uphi[i] - FIELD_P;
		else if(uphi[i] <= -pbytwo)
			uphi[i] = uphi[i] + FIELD_P;


		if(checkifneg && !is_zero(&v))
		{
			int_mul(&t,&v,&tv);
			int_neg(&tv);
		}
		else
			int_mul(&t,&v,&tv);


		int_add(&tv,&y,&x);
		int_sub(&zero,&v,&y);

      		i++;


	}

}

//scalar K is expanded as a poynomial in phi
//finding out the coefficients ui of the equation k=Summation(ui.phi^i(P))

void basephiexpansion(BIGINT *k,long *uphi,int * r_1,int *r_2)
{
	int i,checkifneg,x1,j;
	BIGINT t,x,y,v,tv,fieldp,zero;
	long  uphi_1[L],uphi_2[L];
	BIGINT k_1, k_2;
	generateRandomly(k,r_1,r_2,&k_1, &k_2);
	expansion(&k_1,uphi_1);
	expansion(&k_2,uphi_2);
	long phi2=uphi_2[FIELD_M-1];
	
	for(i=FIELD_M;i>0;i--)
	{
		uphi_2[i]=uphi_2[i-1];
	}
	uphi_2[0]=0;
	for(i=0;i<FIELD_M;i++)
	{
		uphi[i]=uphi_1[i]+uphi_2[i];
	}
	

}

	


//the Ui values are reduced to Di values

void basephiopti(long * u,long * d)
{
	int i;
	for(i=0;i<FIELD_M;i++)

		d[i]=0;

	for(i=0;i<FIELD_M;i++)

	{

		if((i+FIELD_M)<L && (i+(2*FIELD_M))<L)

			d[i]=u[i]+u[i+FIELD_M]+u[i+(2*FIELD_M)];

		else if((i+FIELD_M)<L)

			d[i]=u[i]+u[i+FIELD_M];

		else

			d[i]=u[i];

	}
	
	for(i=0;i<FIELD_M;i++)
	{
		if(d[i] > pbytwo)
		{
			d[i] = d[i] - FIELD_P;
        }
		else if(d[i] <= -pbytwo)
			d[i] = d[i] + FIELD_P;
			
	    //printf("%ld\n",d[i]);
	}
	//printf("\n");


}

void smul_random(int r_1,int r_2,point_t * q)
{
	point_t P[2];
	int msb[2],maxmsb,r[2];
	int i,j,count,start;
	msb[0]  = sizeof(long)*8 - 1;
	msb[1]  = sizeof(long)*8 - 1;
	///////////////////////////////////
	if(r_1<0)
	{
		copypoint(&P[0],phipointn[0]);
		r[0] = - r_1;
	}
	else
	{
		copypoint(&P[0],phipoint[0]);
		r[0] = r_1;
	}
		
		
	if(r_2<0)
	{
		copypoint(&P[1],phipointn[1]);
		r[1] = - r_2;
	}
	else
	{
		copypoint(&P[1],phipoint[1]);
		r[1] = r_2;
	}
	//////////////////////////////////
	//find msbs of each r[i]
	for(i=0;i<2;i++)
	{
		while(1)
		{
			if(r[i] >> msb[i])
				break;
			
			msb[i]--;
			if(msb[i] < 0)
			{
				count++;
				break;
			}
		}
	}
	if(count==FIELD_M)
	{
		fillzeros(q);
		printf("random k is zero!");
		exit(0);
	}
	
	if(msb[0]>msb[1])
		maxmsb = msb[0];
	else
		maxmsb = msb[1];
	

	j=maxmsb;	
	start=1;
	do{

		
		if(start > 1)
		{
			pt_double(q);
	
		}
		
		//--------------------------------------------------------
		for(i=0;i<2;i++)
		{
			if((r[i]>>j)&1)
			{
				if(start==1)
				{
					copypoint(q,&P[i]);
					start++;
				}
				else
				{
					pt_add(q,&P[i]);
				}
			
			}
		}
		j--;
	}while(j>=0);
	
}
int findmax(int *a)
{
	int i,max=-2;
	for(i=0;i<FIELD_M;i++)
	{
		if(a[i] > max)
			max=a[i];
	}
	//xil_printf("%d",max);
	return max;
}
void smul(BIGINT * k, point_t *res, field_t data_x[],field_t data_y[])
{

	int word,  i, j, found=0;
	int malloced;
	FILE * fp1;
	unsigned t1,t2,tval[1000];
	point_t * basepoint = NULL;
	int r1,r2;
	int flag[FIELD_M]={0,0,0,0,0,0,0};
	int msb[FIELD_M],maxmsb,count=0;
	long d[FIELD_M];
	int start,m1,m2;
	long u[L];
	///////////// seal=0;
	point_t q,q_1;
	
	field_t zinv[FIELD_M], zinv2[FIELD_M], zinv3[FIELD_M];

	for(i=0;i<L;i++)
		u[i]=0;
	
	malloced = 0;
	
	if(basepoint == NULL)
	{
		malloced = 1;
		basepoint = (point_t *) malloc(sizeof(point_t));
		// initilize  basepoint
		copy_element(basepoint->x, data_x);
		copy_element(basepoint->y, data_y);
	}

	init_element(basepoint->z);
	basepoint->z[0] = 1;
	
	init_element(q.x);
	init_element(q.y);
	init_element(q.z);
	
	/*****Step - 1
	/*****Base-phi expansion of K is performed*****/
	basephiexpansion(k,u,&r1,&r2);
	
	//u[] contains the base-phi expansion of k
	/*****Optimization of the Bas-phi expansion of k*****/
	basephiopti(u,d);
	/*** **d[] contains the optimized representation of k obtained from the previous base-phi representation
	*******/
	
	
	/*********Step - 2*/
	for(i=0;i<FIELD_M;i++)
	{
		if(d[i] < 0)
		{
			flag[i]=1;
			d[i] = -d[i];
		}
	}
	
	count=0;
	
	
	
	for(i=0;i<FIELD_M;i++)
	msb[i]  = sizeof(long)*8 - 1;
	//find msbs of each d_i
	for(i=0;i<FIELD_M;i++)
	{
		while(1)
		{
			if(d[i] >> msb[i])
				break;
			
			msb[i]--;
			if(msb[i] < 0)
			{
				count++;
				break;
			}
		}
		//xil_printf("msb= %d\n",msb[i]);
	}
	if(count==FIELD_M)
	{
		fillzeros(&q);
		printf("k is zero!");
		exit(0);
	}
	//find the maximum msb among the d_is
	maxmsb = findmax(msb);
	j=maxmsb;
	fp1 = fopen("countertime_1.txt","a");
	t1=timestamp();
	
	//--------------------------------------------------------
	start=1;m1=0;m2=0;
	do{

		
		if(start > 1)
		{
			pt_double(&q);
	
		}
		
		//--------------------------------------------------------
		for(i=0;i<FIELD_M;i++)
		{
			if(flag[i]==0 && (d[i]>>j)&1)
			{
				if(start==1)
				{
					copypoint(&q,phipoint[i]);
					start++;
				}
				else
				{
					pt_add(&q,phipoint[i]);
					//---------------------------------------------------
					m1++;
				}
			
			}
			else if(flag[i]==1 && (d[i]>>j)&1)
			{
				if(start==1)
				{
					copypoint(&q,phipointn[i]);
					start++;
				}
				else
				{
					pt_add(&q,phipointn[i]);
					//----------------------------------------------------
					m2++;
				}
			
			}
		}
		j--;
	}while(j>=0);
	/****Step - 3 */
	smul_random(r1,r2,&q_1);
	/****Step - 4*/
	pt_add(&q,&q_1);
	/**Final Inversion step**/
	// convert to affine coordinates
	//printf("%u",q.z[0]);
	oefinverse(zinv, q.z);          // zinv
	oefsqr(zinv2, zinv);            // zinv^2
	oefmul(zinv3, zinv2, zinv );    // zinv^3
	oefmul(res->x, q.x, zinv2);     // x.(zinv^2)
	oefmul(res->y, q.y, zinv3);     // y.(zinv^3)
	init_element(res->z);
	res->z[0] = 1;
	//printf("\n\n");
	if(malloced) free(basepoint);
	t2=timestamp();
	tval[i]=t2-t1;
	fprintf(fp1, "%u \n", tval[i]);
	fclose(fp1);


}

void nop(long long int delay)
{
	int i;
	for(i=0;i<delay;i++)
		;
}
void print_pt(point_t * p)
{
	int i;
	
	for(i=0;i<FIELD_M;i++)
	   printf("%lld\t",p->x[i]);
	printf("\n");
	
	for(i=0;i<FIELD_M;i++)
	   printf("%lld\t",p->y[i]);
	printf("\n");
	
	for(i=0;i<FIELD_M;i++)
	   printf("%lld\t",p->z[i]);
	printf("\n");
}
point_2p(point_t *p,point_t * res)
{
	field_t zinv[FIELD_M], zinv2[FIELD_M], zinv3[FIELD_M];
	point_t p2;
	
	init_element(p2.x);
	init_element(p2.y);
	init_element(p2.z);
	copypoint(&p2,p);
	pt_double(&p2);
	oefinverse(zinv, p2.z);          // zinv
	oefsqr(zinv2, zinv);            // zinv^2
	oefmul(zinv3, zinv2, zinv );    // zinv^3
	oefmul(res->x, p2.x, zinv2);     // x.(zinv^2)
	oefmul(res->y, p2.y, zinv3);     // y.(zinv^3)
	init_element(res->z);
	res->z[0] = 1;
}

point_add(point_t * p1,point_t * p2,point_t * res)
{
	field_t zinv[FIELD_M], zinv2[FIELD_M], zinv3[FIELD_M];
	point_t p3;
	
	init_element(p3.x);
	init_element(p3.y);
	init_element(p3.z);
	copypoint(&p3,p1);
	pt_add(&p3,p2);
	oefinverse(zinv, p3.z);          // zinv
	oefsqr(zinv2, zinv);            // zinv^2
	oefmul(zinv3, zinv2, zinv );    // zinv^3
	oefmul(res->x, p3.x, zinv2);     // x.(zinv^2)
	oefmul(res->y, p3.y, zinv3);     // y.(zinv^3)
	init_element(res->z);
	res->z[0] = 1;

}



int main()
{
	//init_platform();

	int i,j,malloced,seal=0;
	unsigned t1,t2;
	field_t data_x[7];
	field_t data_y[7];
	unsigned int tval[1000];
	
	BIGINT big;
	int_null(&big);

	/******Space allocation for constant phi points********/
	malloced=0;
	for(i=0;i<FIELD_M;i++)
	{
		if(phipoint[i]==NULL)
		{
			phipoint[i] = (point_t *)malloc(sizeof(point_t));
			init_element(phipoint[i]->x);
			init_element(phipoint[i]->y);
			init_element(phipoint[i]->z);
			malloced=1;
		}
	}
	
	/*****Table reference multiplication based on base-phi*****/
	//Precomputation of Pi - s
	for(i=0;i<FIELD_M;i++)
	{
		computephi(CURVE_PT_X,i,phipoint[i]->x);
		computephi(CURVE_PT_Y,i,phipoint[i]->y);
		phipoint[i]->z[0]=1;
	}
	
	
	//------------------------------------------------------------------------------
	/******Space allocation for negetive constant phi points********/
	malloced=0;
	for(i=0;i<FIELD_M;i++)
	{
		if(phipointn[i]==NULL)
		{
			malloced=1;
			phipointn[i] = (point_t *)malloc(sizeof(point_t));
			init_element(phipointn[i]->x);
			init_element(phipointn[i]->y);
			init_element(phipointn[i]->z);
		}
	}
	
	//assign P0 = - P
	copy_element(phipointn[0]->x,CURVE_mPT_X);
	copy_element(phipointn[0]->y,CURVE_mPT_Y);
	init_element(phipointn[0]->z);
	phipointn[0]->z[0]=1;
	/*****Table reference multiplication based on base-phi*****/
	//Precomputation of Pi - s
	for(i=1;i<FIELD_M;i++)
	{
		computephi(CURVE_mPT_X,i,phipointn[i]->x);
		computephi(CURVE_mPT_Y,i,phipointn[i]->y);
		phipointn[i]->z[0]=1;
	}
	
	//fp1 = fopen("countertime.txt","a");
	//fp2 = fopen("ycoord.txt","a");
	
    /***********
     * P1 = Base Point
     * P2 = i * P1, i = 2,3,4,...
     * P3 = P1 + P2
     * **************/	
	point_t * p1,p2,p3;
	p1 = (point_t*)malloc(sizeof(point_t));
	copy_element(p1->x,CURVE_PT_X);
	copy_element(p1->y,CURVE_PT_Y);
	init_element(p1->z);
	p1->z[0]=1;
	
	
	
    big.hw[INTMAX]=0xFFFF;
    big.hw[INTMAX-1]=0x5555;
    big.hw[INTMAX-2]=0x8888;
    big.hw[INTMAX-3]=0xDDDD;
    
    copy_element(data_x,p1->x);
    copy_element(data_y,p1->y);
	
	for(i=1 ; i<=1000 ; i++)
	{
		//t1=timestamp();
		smul(&big,&p3,data_x,data_y);
		//t2=timestamp();
		//tval[i]=t2-t1;
		//fprintf(fp1, "%u \n", tval[i]);
	}

	
	if(malloced)
	for(i=0;i<FIELD_M;i++)
	free(phipoint[i]);
	 
	//fclose(fp1);
	//fclose(fp2);
	
	return 0;
}