/********************************************************
*														*
*		The purpose of this program is to search out possible		*
*	eigen values to the frobenius equation when operating over a		*
*	field extension.  Taken from Menezes' "Elliptic Curve Public Key	*
*	Cryptosystems" chapter 7, section 3.						*
*														*
*					Author = Mike Rosing					*
*					 date  =  Aug. 12 1999					*
*														*
********************************************************/

#include <stdio.h>
#include <stdlib.h>
#include "field2n.h"
#include "eliptic.h"
#include "multipoly.h"

extern RAMDATA ram_block[];

#define PRIMEdepth		37
#define TESTN	10

INDEX	PRIMES[32] = {3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 
					41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97};

/*  check to see if two FIELD2N values are equal.
	Returns 1 if they are, 0 if they aren't.  */

INDEX isequal( FIELD2N *a, FIELD2N *b)
{
	INDEX	i;
	ELEMENT	r;
	
	i = 0;
	r = 0;
	while ( ( !r) && ( i<MAXLONG))	
	{
		r = a->e[i] ^ b->e[i];
		i++;
	}
	if (r) return (0);
	return(1);
}

void print_FIELD( text, value)
char *text;
FIELD2N *value;
{
	INDEX i;
	
	printf("%s", text);
	SUMLOOP(i) printf ("%8x ",value->e[i]);
	printf("\n");
}

/*  Routine to multiply any MULTIPOLY by a power of x.
	This amounts to a shift in position by the power.
	Enter with z and n, returns with x^n * z.
*/
int multi_xmul( MULTIPOLY z, int n,  MULTIPOLY *xnz)
{
	FIELD2N	*src, *dst;
	int		i;
	
	xnz->degree = z.degree + n;
	if( !get_space( xnz)) return 0;
	src = Address( z);
	dst = AddressOf( xnz);
	for( i=0; i<n; i++)
	{
		null( dst);
		dst++;
	}
	multi_copy( z.degree + 1, src, dst);
}

/*  This routine makes any MULTIPOLY monic.
	Enter with pointer to object to reduce,
	returns with new ram_block in same object
	and frees up previous space.  Reduces degree
	if possible.
*/

int multi_monic( MULTIPOLY *fixit)
{
	MULTIPOLY	dup;
	FIELD2N		t1_inv, *dptr, temp;
	INDEX		i;
	
/*  copy original data and reduce degree if possible  */

	if( !multi_dup( *fixit, &dup)) return 0;
	i = dup.degree;
	dptr = Address(dup);
	while (i)
		if( !zero_check( &dptr[i])) i--;
		else break;
	dup.degree = i;

/*  use inverse of highest poewr coefficient to force a monic polynomial  */

	ring_inv( &dptr[i], &t1_inv);
	one( &dptr[i]);
	while (i)
	{
		i--;
		if( zero_check( &dptr[i]))
		{
			ring_mul( &t1_inv, &dptr[i], &temp);
			copy( &temp, &dptr[i]);
		}
	}

/*  take care of memory management  */

	free_space( fixit);
	fixit->degree = dup.degree;
	fixit->memdex = dup.memdex;
	return 1;
}

/*  core routine to square (a(x) + b(x) y).
	Used to determine sign of an eigenvalue.
	Enter with pointers to results a'(x), b'(x)
	and routine computes
	a'(x) = a^2 + b^2( x^3 + a6)
	b'(x) = x b^2
	Done modulo the curve over y, modulo g1 over x
*/

int squarey( MULTIPOLY a, MULTIPOLY b, 
			MULTIPOLY *newa, MULTIPOLY *newb,
			CURVE crv, MULTIPOLY modulus,
			MULTIPOLY *x2j_mod_g)
{
	MULTIPOLY	a2, b2, curve, temp;
	FIELD2N		*xb, *bptr;

/*  a^2 and b^2 first  */

	if( !square_modf( a, modulus, x2j_mod_g, &a2)) return 0;
	if( !square_modf( b, modulus, x2j_mod_g, &b2)) 
	{
		free_space(&a2);
		return 0;
	}

/*  create curve polynomial  */

	curve.degree = 3;
	if( !get_space( &curve)) 
	{
		free_space(&a2);
		free_space(&b2);
		return 0;
	}
	xb = Address( curve);
	copy( &crv.a6, xb);
	null( &xb[1]);
	null( &xb[2]);
	one( &xb[3]);
	if( !multi_mul( curve, b2, &temp))
	{
		free_space( &a2);
		free_space( &b2);
		free_space( &curve);
		 return 0;
	}
	if( !multi_add( a2, temp, &temp))
	{
		free_space( &a2);
		free_space( &b2);
		free_space( &curve);
		free_space( &temp);
		 return 0;
	}
	if( multi_div( temp, modulus, &temp, newa) <= 0)
	{
		free_space( &a2);
		free_space( &b2);
		free_space( &curve);
		free_space( &temp);
		return 0;
	}
	free_space( &temp);
	free_space( &curve);

/*  do y coefficient  */

	if( !multi_xmul( b2, 1, newb))
	{
		free_space( &a2);
		free_space( &b2);
		return 0;
	}
	free_space( &a2);
	free_space( &b2);
	if( newb->degree >= modulus.degree)
		if( multi_div( *newb, modulus, &temp, newb) <= 0)
			return 0;
	return 1;
}

/*  subroutine to compute y^2^NUMBITS.
	requires curve data, modulus, and pointer to modular squaring table as
	well as place to store results.  since y^2 = (x^3 + a6) + xy, result is
	of the form a(x) + b(x) * y.
	Returns 1 on sucess, 0 on space failure.
*/

int yqmodg( MULTIPOLY modulus, MULTIPOLY *x2j_mod_g,
			CURVE crv, MULTIPOLY *yqx, MULTIPOLY *yqy)
{
	MULTIPOLY	xcoef, ycoef, newx, newy;
	FIELD2N		*ptr;
	INDEX		i;
	
/*  create initial value of y = 0 + 1 * y  */

	xcoef.degree = 0;
	if( !get_space( &xcoef) ) return 0;
	ptr = Address( xcoef);
	null( ptr);
	ycoef.degree = 0;
	if( !get_space( &ycoef))
	{
		free_space( &xcoef);
		return 0;
	}
	ptr = Address( ycoef);
	one( ptr);
//	for( i=0; i<NUMBITS; i++)
	for( i=0; i<TESTN; i++)
	{
		if( !squarey( xcoef, ycoef, &newx, &newy, crv, modulus, x2j_mod_g))
		{
			free_space( &xcoef);
			free_space( &ycoef);
			return 0;
		}
		free_space( &xcoef);
		free_space( &ycoef);
		multi_dup( newx, &xcoef);
		multi_dup( newy, &ycoef);
		free_space( &newx);
		free_space( &newy);
	}
	multi_dup( xcoef, yqx);
	multi_dup( ycoef, yqy);
	free_space( &xcoef);
	free_space( &ycoef);
	return 1;
}
		
main()
{
	MULTIPOLY	f[PRIMEdepth], x2j_modf_table[(PRIMEdepth*PRIMEdepth-1)/2];
	MULTIPOLY	xqx, fw_term, gcd_term;
	MULTIPOLY	x2j_modg[PRIMEdepth], yqx, yqy;
	MULTIPOLY	t1, t2, t3, t4, p1, p2, p3;
	MULTIPOLY	xf3, ff2, fff, abar, bbar, hbar, x3a6;
	INDEX		prime_search, k, wlist[2];
	INDEX		wdex, escape, i, j, w;
	CURVE		crv;
	FIELD2N		normalize, temp, *gptr, *fptr, *xptr;
	FIELD2N		c[256], sum;
	int			good, numterms;

	MULTIPOLY	xq[11], prime, tmp;
	FIELD2N		*xqptr, *primeptr;    
/*	MULTIPOLY	test1, test2, test3, test4;
	FIELD2N		*testptr, t1, t2;  
	int			x;  */
	MULTIPOLY	fr1, fr2, fr3, xm, xone;
	MULTIPOLY	term1, term2;
	MULTIPOLY	yqa[11], yqb[11], g1;
	
/*  initialize lowest level math  and memory management */

	init_ram_space();
	init_ring_math();

/*  initialize curve we are testing  */

	crv.form = 0;
	null( &crv.a2);
	null( &crv.a6);
	crv.a6.e[NUMWORD] = 0xa;

/*  create multiplication list  */

/*	copy( &crv.a6, &c[0]);
	for( i=1; i<6; i++) 
	{
		ring_mul( &c[0], &c[i-1], &c[i]);
	}

/*   check via brute force calculation of x^q.  Since all tables start with x^1
	I initialize that table entry outside and change all the for each f[i] */

	xq[0].degree = 1;
	if( !get_space( &xq[0]))
	{
		printf("ridiculous!\n");
		exit (0);
	}
	xqptr = Address( xq[0]);
	null( xqptr);
	one( &xqptr[1]);
	
/*  create list of division polynomials for this curve  */

	numterms	 = gen_division_polynomial( f, PRIMEdepth+1, crv);
	if( numterms != PRIMEdepth+1)
	{
		printf( "Failed to generate all possible division polynomials.\n");
		printf("found %d out of %d requested. \n", numterms, PRIMEdepth+1);
		exit(-1);
	}

/*  check that division polynomials and math all works  */

	for( i= 4; i<PRIMEdepth/2;  i++)
	{
		multi_mul( f[13], f[13], &term1);
		multi_mul( term1, f[23], &term1);
		multi_mul( term1, f[25], &term1);
		multi_mul( f[24], f[24], &term2);
		multi_mul( term2, f[14], &term2);
		multi_mul( term2, f[12], &term2);
		multi_add( term1, term2, &fw_term);
		multi_gcd( fw_term, f[37], &gcd_term);
		if( gcd_term.degree != f[37].degree)
		{
			printf( "failed at i = %d\n", i);
			exit(0);
		}
	}
	printf("it works!\n");
/*	it does!!	
/*  create constant x for use in test routines below  */

	xm.degree = 1024;
	if(!get_space( &xm)) exit(0);
	xptr = Address( xm);
	for (i=0; i<1024; i++) null( &xptr[i]);
	one( &xptr[1024]);
	xone.degree = 1;
	if(!get_space(&xone)) exit(-1);
	xptr = Address( xone);
	null( xptr);
	one( &xptr[1]);
	
/*  For each division polynomial of a prime compute (x^2^n + x) mod f[k].
	When the product of all the primes used is greater than 2^(n/2) we can find
	the order of the curve using the Chinese Rmeainder Theorem and the
	eigenvalues found with the Frobenius mod k.
*/

	prime_search = 8;
	while( (k=PRIMES[prime_search]) <= PRIMEdepth)
	{

/*  Do brute force calculation as check  */

/*		multi_dup( f[k], &prime);
		for( i=1; i<field_prime; i++)
		{
			if( !multi_mul( xq[i-1], xq[i-1], &tmp)) printf("ERROR at i=%d\n", i);
			if( !multi_div( tmp, prime, &tmp, &xq[i])) printf("Error at i=%d\n",i);
		}

/*  brute force divide x^1024 by f[k] and see if you get same xq[i] as above  */

//		if( !multi_div(  xm, f[k], &fr1, &fr2)) printf("Error dividing x^1024\n");
		
		if( !gen_xmodf( f[k], x2j_modf_table))
		{
			printf( "Failed to generate squaring table at k = %d\n", k);
			exit(-1);
		}

/*  then compue x^2^n + x mod f[k]  */

		if( !xqmodf( f[k], x2j_modf_table, &xqx))
		{
			printf( "Failed to compute x^2^n mod f[%d}.\n", k);
			exit(-1);
		}
		xptr = Address( xqx);
		xptr[1].e[NUMWORD] ^= 1;

		free_space( &prime);
		for( i=1; i<11; i++) free_space( &xq[i]);
		printf("finished prime %d\n", k);
		prime_search++;
		continue;
		
/*  For each w in range 1 ... (k-1)/2, see if g1(x) (equation 7.4 in Menezes)
	has a common factor with division polynomial f[k].  If not, skip this w.
	If so, we have more work to do.
*/
		escape = 0;
		wdex = 0;
		wlist[0] = 0;
		wlist[1] = 0;
		for( w=1; w<= k/2; w++)
		{
/*  compare various ways of computing the same thing to see if there is
	a bug in higher level math.
*/

/*			multi_mul( f[w-1], f[w+1], &fr1);
			multi_div( fr1, f[w], &fr3, &fr1);
			multi_div( fr1, f[w], &fr3, &fr1);
			multi_add( fr1, xone, &fr1);
			multi_div( fr1, f[k], &fr3, &fr1);
			multi_gcd( xq[NUMBITS], fr1, &fr2);
			
/*  original method  */
			
			multi_mul( f[w], f[w], &fw_term);
			multi_mul( fw_term, xqx, &fw_term);
			multi_mul( f[w-1], f[w+1], &gcd_term);
			multi_add( gcd_term, fw_term, &fw_term);
			multi_gcd( fw_term, f[k], &gcd_term);
			if( gcd_term.degree < 1) continue;  /*  gcd == 1, no eigen value  */

/*  if we get here, g1 is a polynomial, not a constant.  Check for special conditions.  */

			printf("degree of g1 = %d\n", gcd_term.degree);
			if( gcd_term.degree == k-1)
			{

/*  then both +w and -w are eigen values of the Frobenius.  */

				printf( "Eigenvalues at k = %d are %d and %d.\n", k, w, -w);
				escape = 1;
			}
			if( gcd_term.degree == f[k].degree)
			{

/*  check to see if g1 == f[k] for all coefficients.  */

				gptr = Address(gcd_term);
				fptr = Address( f[k]);
				if( !isequal( &gptr[gcd_term.degree] , &fptr[gcd_term.degree]))
				{
/*  normalize coefficients so highest powers match  */

					ring_inv( &gptr[gcd_term.degree], &temp);
					ring_mul( &temp, &fptr[gcd_term.degree], &normalize);
					for( i=0; i<=gcd_term.degree; i++)
					{
						ring_mul( &normalize, &gptr[i], &temp);
						copy( &temp, &gptr[i]);
					}
				}

/*  check each term for no match (if all terms match we are equal)  */

				for( i=0; (i<gcd_term.degree) && (!escape); i++)
				{
					SUMLOOP(j)
					{
						if( fptr[i].e[j] != gptr[i].e[j])
						{
							escape = -1;
							wlist[wdex] = w;
							wdex++;
							break;
						}
					}
				}
			
/*  if g1 == f[k] then both eigen values are the same and we only get one 
	useful output.
*/

				if( !escape)
				{
					printf( "One eigenvalue at k = %d is %d\n", k, w);
					escape = 1;
				}
			}

/*  once an eigenvalue is found we need to determine its sign.
	This is done modulo the g1 Polynomial (gcd_term).  First step 
	is to build squaring table modulo g1.
*/
			if( !multi_monic( &gcd_term)) exit(-11);
			if( !gen_xmodf( gcd_term, x2j_modg))
			{
				printf("Failed to generate g1 squaring table\n");
				exit(-1);
			}

/*  next finde y^q mod g1  */

			if( !yqmodg( gcd_term, x2j_modg, crv, &yqx, &yqy))
			{
				printf("Failed to compute y^q mod g1\n");
				exit(-1);
			}
/*  plus y to get y^q + y  */

			xptr = Address( yqy);
			xptr->e[NUMWORD] ^= 1L;

/*  combine all y terms separately from x terms for equation 7.5 
	in Menezes.
	yterms = xf[w]^3( y^q + y)_y + f[w-1]f[w]f[w+1]
	xterms = xf[w]^3*y^q_x + x^2 f[w-1]f[w]f[w+1] + f[w-2]f[w+1]^2
*/
/*  first y term = t2, xf3 = x*f[w]^3  */

			if( !multi_mul( f[w], f[w], &xf3)) exit(-2);
			if( !multi_mul( f[w], xf3, &xf3)) exit(-2);
			if( !multi_mul( xone, xf3, &xf3)) exit(-2);
			if( multi_div( xf3, gcd_term, &t2, &t3) <= 0) exit(-2);
			free_space( &t2);
			if( !multi_mul( yqy, t3, &t2)) exit(-2);
			free_space(&t3);
			if( multi_div( t2, gcd_term, &t1, &t2) <= 0) exit(-2);
			free_space( &t1);

/*  second y term = t3, fff = f[w-1]f[w]f[w+1]  */

			if( !multi_mul( f[w], f[w-1], &fff)) exit(-3);
			if( !multi_mul( f[w+1], fff, &fff)) exit(-3);
			if( multi_div( fff, gcd_term, &t4, &t3) <= 0) exit(-3);
			free_space( &t4);
			if( !multi_add( t2, t3, &bbar)) exit(-3);
	
/*  compute middle term ff2 = f[w-2]f[w+1]^2  */

			if( !multi_mul( f[w-2], f[w+1], &ff2)) exit(-4);
			if( !multi_mul( f[w+1], ff2, &ff2)) exit(-4);
			if( multi_div( ff2, gcd_term, &p2, &ff2) <= 0) exit(-4);
			free_space( &p2);

/*  compute a_bar in Menezes  */

			if( !multi_mul( xf3, yqx, &p2)) exit(-5);
			
			/*  do x^2 * fff  */
			if( !multi_xmul( fff, 2, &p3)) exit(-5);
			
			if( !multi_add( p3, p2, &p2)) exit(-5);
			free_space( &p3);
			if( !multi_add( ff2, p2, &abar)) exit(-5);
			free_space( &p2);
			if( multi_div( abar, gcd_term, &p2, &abar) <= 0) exit(-5);
			free_space( &p2);

			free_space( &t1);
			free_space( &t2);
			free_space( &t3);
			free_space( &p1);
			free_space( &p2);
			x3a6.degree = 3;	/*  create x^3 + a6  */
			if( !get_space( &x3a6)) exit(-6);
			xptr = Address( x3a6);
			copy( &crv.a6, xptr);
			null( &xptr[1]);
			null( &xptr[2]);
			one( &xptr[3]);
			
/*  compute y^q by brute force.  Take modulo g1 after all terms found to
	see if modulo calculations done correctly.
*/

			yqa[0].degree = 0;
			if( !get_space( &yqa[0])) exit(-10);
			xptr = Address( yqa[0]);
			null( xptr);
			yqb[0].degree = 0;
			if( !get_space( &yqb[0])) exit(-10);
			xptr = Address( yqb[0]);
			one( xptr);
			for( i=0; i<NUMBITS; i++)
			{
				if( !multi_mul( yqa[i], yqa[i], &p1)) exit(-10);
				if( !multi_mul( yqb[i], yqb[i], &p2)) exit(-10);
				if( !multi_mul( p2, x3a6, &t1)) exit(-10);
				if( !multi_add( t1, p1, &yqa[i+1])) exit(-10);
				if( !multi_xmul( p2, 1, &yqb[i+1])) exit(-10);
				free_space( &t1);
				free_space(&p1);
				free_space( &p2);
			}
			if( multi_div( yqa[TESTN], gcd_term, &t1, &t2) <= 0 ) exit(-10);
			if( multi_div( yqb[TESTN], gcd_term, &p1, &p2) <= 0) exit(-10);
			/*  stop here and compare t2 with yqx and p2 with yqy  */
			
			/*  t2 = bbar^2  */
			if( !square_modf( bbar, gcd_term, x2j_modg, &t2)) exit(-6);
			if( !multi_mul( x3a6, t2, &t3)) exit(-6);
			free_space( &t2);
			if( !square_modf( abar, gcd_term, x2j_modg, &p2)) exit(-6);
			if( !multi_add( t3, p2, &p2)) exit(-6);
			free_space( &t3);
			/*  p2 = abar^2 + (x^3 + a6)bbar^2  */
			/*  multiply by x: t1 = x*abar*bbar */
			if( !multi_mul( abar, bbar, &p1)) exit(-6);
			if( !multi_xmul( p1, 1, &t1)) exit(-6);
			/*  add it all up  */
			if( !multi_add( p2, t1, &hbar)) exit( -6);

			free_space( &t1);
			free_space( &p1);
			free_space( &p2);
			if( !multi_gcd( hbar, gcd_term, &t1)) exit(-7);
			if( t1.degree) wlist[wdex] = -w;
			else wlist[wdex] = w;
			wdex++;
			free_space( &t1);
			free_space( &yqy);
			free_space( &yqx);
			free_space( &xf3);
			free_space( &ff2);
			free_space( &fff);
			free_space( &abar);
			free_space( &bbar);
			free_space( &x3a6);
			destroy_xmodf( gcd_term, x2j_modg);
	
			if( escape == 1) break;

/*  Either +w or -w is an eigen value, but we don't know which.  
	Add w to w list and keep searching.
*/
/*			wlist[wdex] = w;
			wdex++;
			
/*  see if both eigen values found  */

			if ( wdex > 1) break;
		} /*  end of w loop  */
		if( wdex )
		{
			if (wdex == 2)
				printf("Eigenvalues are %d and %d at k = %d\n", 
					wlist[0], wlist[1], k);
			else
				printf("Duplicate eigenvalue %d at k = %d\n",
					wlist[0], k);
			escape = 1;
		}
		prime_search++;
		if( !escape) printf("No eigenvalues for k = %d\n",k);
		free_space( &xqx);
		destroy_xmodf( f[k], x2j_modf_table);
		free_space( &fw_term);
		free_space( &gcd_term);
	}/*  end of k loop  */
	exit(0);
}/*  end of program  */