Weil-Pairing/Miller.cpp

#include"elliptic.h"


FPOINT * newfpoint(int x, int y)
{
	FPOINT * result = (FPOINT *)malloc(sizeof(fpoint));
	result->x = x;
	result->y = y;
	return result;
}

POINT * newpoint(int a, int b, int c, int d)
{
	POINT * result = (POINT *)malloc(sizeof(point));
	result->x = newfpoint(a,b);
	result->y = newfpoint(c,d);
	return result;
}

CURVE * newcurve(int A, int B)
{
	CURVE * result = (CURVE *)malloc(sizeof(curve));
	result->A = newfpoint(0,A);
	result->B = newfpoint(0,B);
	
	return result;
}

void freepoint(POINT * a)
{
	free(a->x);
	free(a->y);
}

//absolute value modulo p
int ABS(int a, int p)
{
	return (a>=0)?a%p:(p-(-a)%p)%p;
}

//extended euclid algorithm
int gcdEx(int a, int b, int *x, int *y) 
{
    if(b==0){
        *x = 1,*y = 0;
        return a;
    }
    else{
        int r = gcdEx(b, a%b, x, y);
        int t = *x;
        *x = *y;
        *y = t - a/b * *y;
        return r;
    }
}

int powermod(int a, int n, int p)
{
	int result, DB;
	
	result = 1;
	DB = a;
	
	while(n>0){
		if(n&1)result = (result*DB)%p;
		DB = (DB*DB)%p;
		n >>= 1;
	}
	
	return ABS(result,p);
}

//modulo inverse
int inver(int a, int p)
{
	int s,t;
	if(gcdEx(a,p,&s,&t)!=1)return 0;
	return ABS(s,p);
}

int randomnonq(int p)
{
	int i, k = (p-1)>>1;
	for(i = 1;i < p;i++){
		if(powermod(i,k,p) != 1)break;
	}
	return i;
}

int modsquareroot(int a, int p)
{
	if(powermod(a,(p-1)>>1,p) != 1)return -1;
	
	int r = (p-1)>>1;
	int b = randomnonq(p);
	int x = r, y = 0;

	while(!(x&1)){
		x >>= 1; y >>= 1;
		if(ABS(powermod(a,x,p)*powermod(b,y,p),p) != 1)y += r;
	}
	
	return ABS(powermod(a,(x+1)>>1,p)*powermod(b,y>>1,p),p);
}

bool millerrabin(int n, int r)
{
	int s = 0, t = n - 1;
	
	while(!(t&1)){
		s++; t >>= 1;
	}
	while(r--){
		int b = rand()%(n-1) + 2;
		int r0 = powermod(b,t,n);
		int s0 = s - 1;
		
		if(r0 == 1 || r0 == n - 1)continue;
		if(s0 < 1)return false;
		while(s0--){
			r0 = powermod(r0,2,n);
			if(r0 == n-1)break;
			if(s0 == 0)return false;
		}
	}
	
	return true;
} 

int largerandom(int n)
{
	int a = rand()%10000, b = rand()%10000;
	return (a+b*10000)%n;
}

int randomgoodprime(int n)
{
	int p = 4;
	while(!millerrabin(p,10))p = 12*largerandom(n)+11;
	return p;
}

//field element equlity judgement
bool equl(FPOINT * a, FPOINT * b)
{
	if(a->x == b->x && a->y == b->y)return true;
	return false;
}

bool equln(FPOINT * a, FPOINT * b, int p)
{
	if(a->x == ABS(-b->x,p)&&a->y == ABS(-b->y,p))return true;
	return false;
}

bool pequl(POINT * a, POINT * b)
{
	return equl(a->x,b->x)&&equl(a->y,b->y);
}

void showelement(FPOINT * p)
{
	printf("(%d,%d)\n",p->x,p->y);
}

FPOINT * fneg(FPOINT * a, int p, FPOINT * result)
{
	result->x = ABS(-a->x,p);
	result->y = ABS(-a->y,p);
	
	return result;
}

//field elements addition
FPOINT * fadd(FPOINT * a, FPOINT * b, int p, FPOINT * result)
{
	int x = ABS(a->x + b->x,p);
	int y = ABS(a->y + b->y,p);
	
	result->x = x;
	result->y = y;
	
	return result;
}

//field elements minus
FPOINT * fminus(FPOINT * a, FPOINT * b, int p, FPOINT * result)
{
	int x = ABS(a->x - b->x,p);
	int y = ABS(a->y - b->y,p);
	
	result->x = x;
	result->y = y;
	
	return result;
}

//field elements multiplication
FPOINT * fmulti(FPOINT * a, FPOINT * b, int p, FPOINT * result)
{
	int x = ABS(a->y*b->x + a->x*b->y,p);
 	int y = ABS(a->y*b->y - a->x*b->x,p);
 	
 	result->x = x;
 	result->y = y;
 	
	return result;
}

//field assignment
FPOINT * assign(FPOINT * a, FPOINT * b)
{
	a->x = b->x;
	a->y = b->y;
	return a;
}


//field element expontinal
FPOINT * fpower(FPOINT * a, int n, int p, FPOINT * result)
{
	n = n%(p*p - 1);
	
	FPOINT * DB = newfpoint(0,0);
	
	assign(DB,a);
	assign(result,ONE);

	while(n > 0){
		if(n&1)fmulti(DB,result,p,result);
		fmulti(DB,DB,p,DB);
		n >>= 1;	
	}
	
	free(DB);
	
	return result;
}

//field element inverse
FPOINT * inverse(FPOINT * a, int p, FPOINT * result)
{	
	if(a->x == 0){
		result->x = 0;
		result->y = inver(a->y,p);
		return result;
	}
	
	int x = inver(ABS(-a->y*a->y*inver(a->x,p)-a->x,p),p);
	int y = ABS(-a->y*inver(a->x,p)*x,p);
	
	result->x = x;
	result->y = y;
	
	return result;
}

//multiply by number
FPOINT * fnmulti(FPOINT * a, int b, int p, FPOINT * result)
{
	result->x = ABS(a->x*b,p);
	result->y = ABS(a->y*b,p);
	return result;
}

POINT * passign(POINT * a, POINT * b)
{
	assign(a->x,b->x);
	assign(a->y,b->y);
	return a;
}

bool testpoint(POINT * p, CURVE * c, int p1)
{
	if(pequl(p,O))return true;
	
	FPOINT * x = newfpoint(0,0);
	FPOINT * y = newfpoint(0,0);
	
	assign(x,p->x);
	
	fadd(fadd(fpower(x,3,p1,x),fmulti(c->A,p->x,p1,y),p1,x),c->B,p1,x);
	
	fpower(p->y,2,p1,y);
	
	return equl(x,y);
}

void showpoint(POINT * p)
{
	printf("\n[(%d,%d),(%d,%d)]\n",p->x->x,p->x->y,p->y->x,p->y->y);
}

POINT * pneg(POINT * a, int p, POINT * result)
{
	if(pequl(a,O)){
		passign(result,O);
		return result;
	}
	assign(result->x,a->x);
	fneg(a->y,p,result->y);
	
	return result;
}

//curve point additon
POINT * add(POINT * p1, POINT * p2, CURVE * c, int p, POINT * result)
{	
	if(pequl(p1,O))return passign(result,p2);
	
	if(pequl(p2,O))return passign(result,p1);
	
	if(equl(p1->x,p2->x)&&equln(p1->y,p2->y,p)){
		return passign(result,O);	
	}
	
	FPOINT * x, * y, *lambda;
	
	x = newfpoint(0,0);
	y = newfpoint(0,0);
	lambda = (FPOINT *)malloc(sizeof(fpoint));
	
	if(equl(p1->x,p2->x)){
		fadd(fnmulti(fpower(p1->x,2,p,lambda),3,p,lambda),c->A,p,lambda);
		fmulti(lambda,inverse(fnmulti(p1->y,2,p,x),p,x),p,lambda);
	}else{
		fminus(p2->y,p1->y,p,lambda);
		fmulti(lambda,inverse(fminus(p2->x,p1->x,p,x),p,x),p,lambda);
	}
	
	fminus(fminus(fpower(lambda,2,p,x),p1->x,p,x),p2->x,p,x);
	fminus(fmulti(fminus(p1->x,x,p,y),lambda,p,y),p1->y,p,y);
	
	assign(result->x,x);
	assign(result->y,y); 
	
	free(lambda);
	free(x);
	free(y);
	
	return result;
}

POINT * minus(POINT * p1, POINT * p2, CURVE * c, int p, POINT * result)
{
	POINT * temp = newpoint(0,0,0,0);
	add(pneg(p2,p,temp),p1,c,p,result);
	
	freepoint(temp);
	
	return result;
}

//power of points addition
POINT * ppower(POINT * a, int n, CURVE * c, int p, POINT * result)
{
	POINT * DB;

	DB = newpoint(0,0,0,0);
	
	passign(DB,a);
	passign(result,O);
	
	while(n > 0){
		if(n&1)add(DB,result,c,p,result);
		add(DB,DB,c,p,DB);
		n >>= 1;
	}
	
	freepoint(DB);
	
	return result;
}

POINT * randompoint(CURVE * c, int p)
{	
	POINT * result = newpoint(0,0,0,0);
	
	int x, y;
	
	if((p-1)%3 && equl(c->A,ZERO)){
		y = largerandom(p);
		int r = inver(3,p-1);
		
		x = powermod(y*y-1, r, p);
	}else{
		while(1){
			x = largerandom(p);
			y = ABS(powermod(x,3,p) + c->A->y*x + c->B->y,p);
			
			if((y=modsquareroot(y,p)) != -1)break;
		}
	}
	
	result->x = newfpoint(0,x);
	result->y = newfpoint(0,y);
	
	return result;
}

FPOINT * primitroot(int p)
{
	int x, y;
	
	y = inver(2,p);
	
	if((x = modsquareroot(-3,p)) != -1){
		y = ABS(x*y - y,p);
		return newfpoint(0,y);
	}
	
	x = modsquareroot(ABS(y*y - y + 1,p),p);
	
	return newfpoint(x,ABS(-y,p)); 
}

POINT * phi(POINT * a, int p, POINT * result)
{	
	fmulti(a->x,primitroot(p),p,result->x);
	assign(result->y, a->y);
	
	return result;
}

FPOINT * evalueline(POINT * a, POINT * b, POINT * in, int p, CURVE * c, FPOINT * result)
{
	FPOINT * temp1, *temp2;
	
	temp1 = newfpoint(0,0); temp2 = newfpoint(0,0);
	
	if(pequl(a,b)){
		fadd(fnmulti(fpower(a->x,2,p,temp1),3,p,temp1),c->A,p,temp1);
		fminus(in->x,a->x,p,temp2);
		fmulti(temp1,temp2,p,result);
		fmulti(fnmulti(a->y,2,p,temp1),fminus(in->y,a->y,p,temp2),p,temp1);
		fminus(result,temp1,p,result);
		
	}else{
		fmulti(fminus(in->x,b->x,p,temp1),fminus(a->y,b->y,p,temp2),p,result);
		fmulti(fminus(a->x,b->x,p,temp1),fminus(in->y,b->y,p,temp2),p,temp1);
		fminus(result,temp1,p,result);
	}
	
	free(temp1); free(temp2);
	
	return result;
}

bool miller(POINT * a, POINT * b, CURVE * c, int m, int p, FPOINT * f)
{
	FPOINT * temp = newfpoint(0,0);
	POINT * t = newpoint(0,0,0,0), * tp = newpoint(0,0,0,0);
	
	assign(f,ONE);
	passign(t,a);
	
	int i = 0, array[(int)logb((double)m)+1];
	
	while(m){
		if(m&1)array[i] = 1;
		else
			array[i] = 0;
		m >>= 1;
		i++;
	}
	
	for(int j = i - 1;j > 0; j--){
		fmulti(f,f,p,f); fmulti(f,evalueline(t,t,b,p,c,temp),p,f);
		if(equl(f,ZERO))return false;
		add(t,t,c,p,t);
		evalueline(t,pneg(t,p,tp),b,p,c,temp);
		if(equl(temp,ZERO))return false;
		fmulti(inverse(temp,p,temp),f,p,f);
		if(array[i] == 1){
			fmulti(f,evalueline(t,a,b,p,c,temp),p,f);
			if(equl(f,ZERO))return false;
			add(t,a,c,p,t);
			evalueline(t,pneg(t,p,tp),b,p,c,temp);
			if(equl(temp,ZERO))return false;
			fmulti(inverse(temp,p,temp),f,p,f);	
		}
	}
	return true;
}

void init()
{
	ONE = newfpoint(0,1);
	ZERO = newfpoint(0,0);
	O = newpoint(-1,-1,-1,-1);
	srand((int)time(0));
}

int main()
{
	init();
	int p=23;
	FPOINT * test = newfpoint(0,14);
	FPOINT * test1;
	
	CURVE * c = newcurve(0,1); 
	
	POINT * P1, * P2, * temp = newpoint(0,0,0,0);
	
	//add(P,P,c,p,P);
	
	P1 = randompoint(c,p);//newpoint(0,2,0,7);
	P2 = randompoint(c,p);//newpoint(0,16,0,125);
	
	miller(P1,P2,c,10,p,test);
	
	showelement(test);
	
	//passign(P2,P1);
	
	showpoint(P1);
	showpoint(P2);
	
	ppower(P1,9,c,p,temp);
	
	showpoint(temp);
	
	phi(temp,p,temp);
	
	showpoint(temp);
	//minus(P,P,c,p,P);

	//pneg(P,p,P);
	
	//passign(P,O);
	
	printf("%d\n",testpoint(temp,c,p));
	
	//int a = 4;
	
	
	//showelement(inverse(test,p,test1));
	
	//test1 = primitroot(p);
	
	//showelement(test1);
	
	//fpower(test1,3,p,test1);
	
	//showelement(test1);
	
	//printf("%d\n",inver(14,p));
	//printf("%d\n", modsquareroot(a,p));
	
	printf("%d\n",millerrabin(10,10));
	
	return 0;
}