#include <time.h>
#include <cstdlib>
#include <math.h>
#include <iostream>
#include <sstream>
#include <fstream>
#include "find_curves.h"
#include <omp.h>

//OMP_SET_NUM_THREADS(4);

using namespace std;

void Find_Curves::smallestM(const long& b, int& m, const long& p){ //helper function for shanks, might need modification
	while( modExp(b, 1 << m, p) != 1){
		m++;
	}
}
Find_Curves::Find_Curves(){//standard constructor
	curve_order=0;
} 
Find_Curves::Find_Curves(const int& n, const vector<long>& _prime_array){ //constructor that takes an integer n. n should be prime
	curve_order=n;	
	prime_array=_prime_array;
}
long Find_Curves::modExp(long a, long b, const long& p){ 
	// computes a^b mod p by successive squaring and bit manipulation
	// a,b are arguments, p is the prime modulus
	long rval=1;
		while(b>0){
			if((b & 1) == 1){
				rval = (rval * a) % p;
			}
			b= b >> 1;
			a= (a * a) % p;
		}
		return rval;
}
long Find_Curves::jacobi(const long& a, const long& p){ 
	// computes the jacobi symbol (a / p)
	//arguments a,p. p should be a prime.
	if((a % p)==0){
		return 0;
	}
	long rval=1;
	long mod8=0;
	long temp=0;
	
	long p_temp=p;
	long a_temp= (a%p);
	while(a_temp!=0){
		while(a_temp%2==0){ //pulls out factors of 2 and computes jacobi(2/n)
			a_temp= a_temp>>1;
			mod8= (p_temp % 8);
			if( mod8==3 || mod8==5){
				if(rval==1){
					rval=-1;
				}else{
					rval=1;
				}
			}
		}
		temp=a_temp;
		a_temp=p_temp;
		p_temp=temp;

		if((a_temp%4)==3 && (p_temp%4) ==3){ //applies quadratic reciprocity
			if(rval ==1){
				rval=-1;
			}else{
				rval=1;
			}
		}
		a_temp= a_temp%p_temp;
	}
	return rval;
}
long Find_Curves::shanks(const long& a, const long& p){  
	// computes sqrt(a) mod p
	// p should be prime

	/* This method deserves some explanation. The name is derived from an algorithm of Tonelli and Shanks for computing the square root mod p, for p prime. The expected running time is O(ln^4 p). 
	(see "A Course in Computational Algebraic Number Theory" by Henri Cohen, page 32 for reference). This algorithm begins after the first two if() statements, which catch the more trivial cases to enhance overall
	efficiency. The first case ((p%4) == 3) applies for half of all primes, while the second case ((p%8) ==5 && modExp(a, (p-1)/4, p)==1) applies for half of all primes not remaining after the first test.

	*/
	srand(time(NULL));

	if((p%4)==3){ //easier for p = 3 mod 4 (complexity is O((lg p)^3) bit operations, see "Algorithmic Number Theory, Volume 1" page 155, published 1997, written by Bach and Shallit
		return modExp(a, (p+1)/4, p);
	}
	if((p%8) ==5 && modExp(a, (p-1)/4, p)==1){ //sometimes easier for p= 5 mod 8
		return modExp(a, (p+3)/8, p);
	}

	//Beginning of the algorithm by Tonelli and Shanks.
	long e= 0; 
	long q= p-1;

	while((q%2)== 0 && q>0){
		e++;
		q= q>>1;
	}
	long n=2;
	while(jacobi(n,p) != -1){
		n= (long) (rand()%p);
	}
	long z= modExp(n, q, p);
	long y, r, x, b, t;
	y=z;
	r=e;
	x= modExp(a, (q-1)/2, p);
	b= (a*x*x)%p;
	x=(a*x)%p;
	while(b%p != 1){
		int m=1;
		smallestM(b, m, p);
		t=modExp(y, 1<<(r-m-1), p);
		y=(t*t)%p;
		r=m;
		x=(x*t)%p;
		b=(b*y)%p;
	}
	return x;
}
bool Find_Curves::isEC(const long& a, const long& b, const long& p){  
	//checks for multiple roots in equation
	// Arguments: a,b are the coefficients of the curve in Weierstrass form (the one mentioned previously). p should, as usual, be prime.
	// This equation might be deprecated to the elliptic_curves class, which could provide its own method to check if an equation is singular
	// Alternatively, this might be altered to take an elliptic curve as an argument instead of three longs.
	if(((4*a*a*a + 27*b*b) %p) == 0){ // i.e. if the discriminant of the elliptic curve is zero, return false (its a singular curve, hence not elliptic). Else, return true.
		return false;
	}else{
		return true;
	}
}
long Find_Curves::ecOrder(const long& a, const long& b, const long& p){ 
	long order=0;
	for(long x= 0; x < p; x++){
		long temp= 0;
		temp= (((x*x) % p) *x) % p;
		temp= (temp+ (a* x)) %p;
		temp= (temp+b) % p;
		order += (1+jacobi(temp, p));
	}
	return order+1;
}
void Find_Curves::findGenerators(const long& a, const long& b, const long& p){
	for(long x= 0; x < p; x++){
		long temp= 0;
		temp= (((x*x) % p) *x) % p;
		temp= (temp+ (a* x)) %p;
		temp= (temp+b) % p;
		if(jacobi(temp, p)== 1){
			long y= shanks(temp, p);
			curve_list.push_back(Elliptic_Curve(a, b, p, x, y));
		}
	}
}
void Find_Curves::findECs(const long& p){
	#pragma omp parallel for collapse(2)
	for(long a=0; a< p; a++){
		for(long b=1; b<p; b++){
			if(isEC(a,b,p) && ecOrder(a, b, p) == curve_order){
				#pragma omp critical(modify_curve_list)
				{
				findGenerators(a, b, p);
				}
			}
		}
	}
} 
void Find_Curves::ECsOfOrderCurveOrder(){
	int i=0;
	while(i < prime_array.size() && ((prime_array[i]+1)+2*(long) sqrt(prime_array[i])) <= curve_order){
		i++;
	}
	for(int j=i; prime_array[j] <=curve_order; j++){
		findECs(prime_array[j]);
	}
}

void Find_Curves::printCurves(){
	while(!curve_list.empty()){
		cout << curve_list.back().to_serial() << endl;
		curve_list.pop_back();
	}
}


void Find_Curves::write_to_file(){

        ostringstream result;
        result << curve_order;
	string temp= result.str()+ "_curve_list.txt";

        ofstream curve_file;
        curve_file.open(temp.c_str());
	while(!curve_list.empty()){
	        curve_file << curve_list.back().to_serial() << endl;
		curve_list.pop_back();
	}
        curve_file.close();

}