extb92_finite.cpp

#include "key_rate_utilsb92.hpp"
#include <fstream>
#include <iostream>


//#define USE_FLOAT128
#ifdef  USE_FLOAT128
	//#include <quadmath.h>
    //typedef __float128  long_double_t;
     typedef long double  long_double_t;
#else
    typedef long double  long_double_t;
#endif
/*
	ostream& operator<<(ostream& os, const  long_double_t& t){
		return os << (double) t; 
	}
	long_double_t fabs(const long_double_t& x){
		#ifdef USE_FLOAT128
		return fabsq(x);
		#else
		return fabs(x);
		#endif
		}

	long_double_t log(const long_double_t& x){
		#ifdef USE_FLOAT128
		return logq(x);
		#else
		return log(x);
		#endif
	}
	long_double_t sqrt(const long_double_t& x){
		#ifdef USE_FLOAT128
		return sqrtq(x);
		#else
		return sqrt(x);
		#endif
	}
	long_double_t isnan(const long_double_t& x){
		#ifdef USE_FLOAT128
		return isnanq(x);
		#else
		return isnan(x);
		#endif
	}
	long_double_t min(long double y, const long_double_t& x){
		#ifdef USE_FLOAT128
		return min((long_double_t) y, x);
		#else
		return min(y,x);
		#endif
	}*/
//namespace Extb92Finite{

	//using namespace key_rate_utils;
	using namespace std;
	double STEP = 1e-3;
	
	double calculate_q(double alph, bool mistmatched = true){
		//no mismatched would be different -- and probably here we would have two different kinds
		auto a2 = pow(alph, 2);
		auto b2 = 1-a2;

		return ((AtoB["P0A"] - a2*AtoB["P00"] - b2*AtoB["P01"])/2) - 
			   ((AtoB["PA1"] - a2*AtoB["P01"] - b2*AtoB["P11"])/2);
		
	}

	double calculate_blam0(double E12, double alph, bool mismatched = true){
		//Calculate 
		auto a2 = pow(alph, 2);
		auto b2 = 1-a2;

		auto ret = calculate_q(alph,mismatched) - a2*AtoB["P01"] 
		+ b2 * E12;
		//cout << "L0= " << ret << endl;
		return ret;
	}

	double calc_E03(double alph, double E12, bool mismatched){
		auto a2 = pow(alph, 2);
		auto b2 = 1-a2;

		//when a is 0 we div by zero!!
		/*
		cout << "E03 lines" << endl;
		cout << a2*b2*(AtoB["P00"]+AtoB["P11"]) << endl;
		cout << b2*b2*AtoB["P10"] + a2*a2*AtoB["P01"] << endl;
		cout << b2*(AtoB["PA0"] - a2*AtoB["P00"] - b2*AtoB["P10"]) << endl;
		cout << -a2*(AtoB["P0A"] - a2*AtoB["P00"] - b2*AtoB["P01"]) << endl;
		cout << - b2*(AtoB["P1A"] - a2*AtoB["P10"] - b2*AtoB["P11"]) << endl;
		cout << a2*(AtoB["PA1"] - a2*AtoB["P10"] - b2*AtoB["P11"]) << endl;
		cout << - AtoB["PAB"] << endl;
		cout << "/" << (2*a2*b2) << endl;
*/		
		auto ret = 
				(a2*b2*(AtoB["P00"]+AtoB["P11"]) +
				 b2*b2*AtoB["P10"] + a2*a2*AtoB["P01"] +
				 b2*(AtoB["PA0"] - a2*AtoB["P00"] - b2*AtoB["P10"]) 
				- a2*(AtoB["P0A"] - a2*AtoB["P00"] - b2*AtoB["P01"]) 
				- b2*(AtoB["P1A"] - a2*AtoB["P10"] - b2*AtoB["P11"]) 
				+ a2*(AtoB["PA1"] - a2*AtoB["P10"] - b2*AtoB["P11"]) 
				- AtoB["PAB"])/(2*a2*b2) - E12;
		//cout << "E12 = " << E12 << endl;
		//cout << "E03 = " << ret<< endl;

		return ret;

	}

	double calculate_blam1(double E12, double alph, bool mismatched = true){
		//Calculate 
		auto a2 = pow(alph, 2);
		auto b2 = 1-a2;

		auto ret = calculate_q(alph, mismatched) - a2*AtoB["P01"] 
		+ b2 * calc_E03(alph, E12, mismatched);
		//cout << "L1= " << ret << endl;

		return ret;
		
	}


	
	//Need to change generate noise statistics tp add X stats for bob!
	//in finite case iteratate throuugh noise stats and change them based on uncertainty"

	double calculate_key_rate_asymptotic(double alph, double noise, bool symmetric=true, double CAD=1, bool mismatched = true){
		generateNoiseStatistics(symmetric, alph, noise);
		//add CAD
		vector<double> E0 {AtoB["P01"], 1-AtoB["P0A"]};
		vector<double> E1 {AtoB["PAB"], 1-AtoB["PA0"]};
		
		/*cout << "E0 = ";
		for (auto p: E0)
			cout << p << " , ";
		cout << endl;
		cout << "E1 = ";
		for (auto p: E1)
			cout << p << " , ";
		cout << endl;*/
		double limE12 = sqrt(AtoB["P01"]*AtoB["P10"]);
		double min_ent = 1;

		//cout << "LimE12 = " << limE12 << endl;
		double iterations = 2*limE12/STEP;
		int i = 0;
		for (double E12 = -limE12; E12 <= limE12 ; E12+=STEP){
			i++;
			vector<double> Ls{calculate_blam0(E12,alph, mismatched), calculate_blam1(E12, alph,mismatched)};
			/*cout << "Ls = [";
			for (auto i: Ls)
				cout << i << ", ";
			cout << endl;
			cout << "L1+L2= " << Ls[0] + Ls[1] << endl; */
			double s_ent = s_entropy(E0, E1, Ls, CAD);
			//cout << "[-] Completed " << i << "out of " << iterations << " iterations" << endl;
			if (s_ent > -.5 && s_ent < min_ent)
				cout << "[+] Found new best min ent at E12= " << E12 << " S(A|E) = " << s_ent << endl;
				min_ent = s_ent;
		}

		vector<double> Pij{AtoB["P0B"], AtoB["P01"], AtoB["PAB"], AtoB["PA1"]};
		vector<double> Pij_key_non_norm;
		vector<double> Pij_key;

		transform(Pij.begin(), Pij.end(), back_inserter(Pij_key_non_norm),
				  [CAD](double p) -> double {return pow(p, CAD);} );
		
		double p_acc = 0;
		for (auto p: Pij_key_non_norm)
			p_acc += p; 
	
		transform(Pij_key_non_norm.begin(), Pij_key_non_norm.end(), back_inserter(Pij_key),
				  [p_acc](double p) -> double {return p/p_acc;} );
		/*
		cout << "Pij = ";
		for (auto p: Pij)
			cout << p << " , ";
		cout << endl;
		cout << "Pij_non_norm = ";
		for (auto p: Pij_key_non_norm)
			cout << p << " , ";
		cout << endl;
		cout << p_acc << endl;
		cout << "Pij_key/^ = ";
		for (auto p: Pij_key)
			cout << p << " , ";
		cout <<endl;
		*/
		//cout << "Best min_ent found : " << min_ent << endl;
		double HAB = H_AB(Pij_key);
		cout << "H(A|B): = " << HAB << endl;
		double key_rate = min_ent - HAB;
		cout << "Key Rate for alph = " << alph << " Q = " << noise << " is " << key_rate << endl;
		
		//return key_rate;  
		double Q = noise;
		double a2 = alph*alph;
		cout << "Asymptotic key rate = " << key_rate << endl;
		cout << "Times factor = "<< (Q+(1-2*Q)*(1-a2))*.5 << endl; 
		cout << "Effective key rate = " << key_rate*(Q+(1-2*Q)*(1-a2))*.5 << endl;
		return key_rate*(Q+(1-2*Q)*(1-a2))*.5;

	}

	double calculate_key_rate_finite_inside(double alph, int CAD = 1, bool mismatched = true){
		vector<double> E0 {AtoB["P01"], 1-AtoB["P0A"]};
		vector<double> E1 {AtoB["PAB"], 1-AtoB["PA0"]};

		double limE12 = sqrt(AtoB["P01"]*AtoB["P10"]);
		double min_ent = 1;
		//cout << AtoB["P01"] << " " << AtoB["P10"] << endl;
		double iterations = 2*limE12/STEP;
		int i = 0;
		//cout << "LimE12 = " << limE12 << endl;
		for (double E12 = -limE12; E12 <= limE12 ; E12+=STEP){
			i++;
			vector<double> Ls{calculate_blam0(E12,alph, mismatched), calculate_blam1(E12, alph,mismatched)};
			
			double s_ent = s_entropy(E0, E1, Ls, CAD);
			//cout << "sent = " << s_ent << endl;
			//cout << "[-] Completed " << i << "out of " << iterations << " iterations and got " << s_ent << endl;
			if (s_ent > -.5 && s_ent < min_ent)
				//cout << "[+] Found new best min ent at E12= " << E12 << " S(A|E) = " << s_ent << endl;
				min_ent = s_ent;
		}
		//cout << "returning " << min_ent << endl;
		return min_ent;
		}

	double calculate_key_rate_finite_outside(double alph, double noise, long long N, 
		double Penc = .5, double Pzb = .5, double CAD = 1, 
		bool symmetric = true, bool mismatched = true){


		cout << "alph " << alph <<endl;
		cout << "noise "<<  noise <<endl;
		cout << "N " << N <<endl;
		cout << "Penc " << Penc <<endl;
		generateNoiseStatistics(symmetric,alph,noise);
		long long C01, C10, C0A, C1A, CA0, CAB;
		long long CA1, C0B, numkey;
		long_double_t a2 = alph*alph;
		long_double_t Q = noise;

		C01 = N*Penc*Pzb*Q/2;
		C10 = N*(1-Penc)*Pzb*Q;
		C0A = N*Penc*(1-Pzb)*(Q+(1-2*Q)*a2)/2;
		C1A = N*(1-Penc)*(1-Pzb)*(Q+(1-2*Q)*(1-a2));
		CA0 = N*Penc*Pzb*(Q+(1-2*Q)*a2)/2;
		CAB = N*Penc*(1-Pzb)*Q/2;

		CA1 = N*Penc*   Pzb * (Q+(1-2*Q)*(1-a2))/2;
		C0B = N*Penc*(1-Pzb)* (Q+(1-2*Q)*(1-a2))/2;

		numkey = CA1 + C0B;

		long_double_t B01, B10, B0A, B1A, BA0, BAB; //base probs
		long_double_t M01, M10, M0A, M1A, MA0, MAB; //offsets that lead to worst minimum
		long_double_t MM01, MM10, MM0A, MM1A, MMA0, MMAB; //offsets for the best key rate
		long_double_t best_key_ent;
		long_double_t best_eps_bar, best_eps_bar_prime;
		long_double_t best_key_rate = -float('inf');


		B01 = AtoB["P01"]; B10 = AtoB["P10"]; B0A = AtoB["P0A"];
		B1A = AtoB["P1A"]; BA0 = AtoB["PA0"]; BAB = AtoB["PAB"];

		bool skipEbar, skipEbarP, skip01, skip10, skip0A, skip1A, skipA0, skipAB;
		skipEbar = true;
		skipEbarP = true;

		skip01 = false; 
		skip10 = false;
		skip0A = false;
		skip1A = false;
		skipA0 = true;
		skipAB = true;
		if(true){
			skip01 = false; 
			skip10 = false;
			skip0A = false;
			skip1A = false;
			skipA0 = false;
			skipAB = false;
		}
		long_double_t step01, step10, step0A, step1A, stepA0, stepAB;
		int thresh = 1e7; //this tries to ensure that we arent searching massive ranges if its small N

		step01 = N > thresh ? .001 : .075; 
		step10 = N > thresh ? .001 : .075; 
		step0A = N > thresh ? .001 : .075; 
		step1A = N > thresh ? .001 : .075; 
		stepA0 = N > thresh ? .001 : .075; 
		stepAB = N > thresh ? .001 : .075;

		bool smart_step = true;

		long_double_t leakEC, eps, epsEC;
		eps = 1e-9; 	//security parameter, larger this is bigger key rate but less secure
		epsEC = 1e-10;  // probability of EC failing

		
		long_double_t eps_bar_step = 1e-10;
		long_double_t eps_bar_prime_step = 1e-10;
		int eps_count = 0;

		long_double_t Ebar_start = eps_bar_step; // where we start if we arent skipping
		long_double_t Ebar_lim = eps-epsEC;
		
		long_double_t Ebar_prime_start = eps_bar_prime_step; // where we start if we arent skipping
		long_double_t Ebar_prime_lim;
		
		if (skipEbar){
			Ebar_start = Ebar_lim - eps_bar_step; //this can be whatever we decide
			Ebar_lim = Ebar_start + eps_bar_step; //this ensures we dont iterate
		}
		

		long_double_t start01;
		long_double_t start10;
		long_double_t start0A;
		long_double_t start1A;
		long_double_t startA0;
		long_double_t startAB;

		long_double_t lim01;
		long_double_t lim10;
		long_double_t lim0A;
		long_double_t lim1A;
		long_double_t limA0;
		long_double_t limAB;
		long_double_t eps_bar;
		long_double_t eps_bar_prime;

		for (eps_bar = Ebar_start; eps_bar < Ebar_lim-eps_bar_step/2; eps_bar+=eps_bar_step){
			Ebar_prime_lim = eps_bar;
			if (skipEbarP){
				Ebar_prime_start = Ebar_prime_lim - eps_bar_prime_step; //this can be whatever we decide
				Ebar_lim = Ebar_prime_start + eps_bar_prime_step; //this ensures we dont iterate
			}

			for (eps_bar_prime = Ebar_prime_start; eps_bar_prime < Ebar_prime_lim - eps_bar_prime_step/2; eps_bar_prime+=eps_bar_prime_step){
				double min_ent = 1;
				double ent;
				eps_count+=1;
				cout << "On eps iteration "  << eps_count << endl;
				cout << "	Eps_bar = "		 << eps_bar << endl;
				cout << "	Eps_bar_prime = "   << eps_bar_prime << endl;

				//auto gamma = [eps_bar_prime] (int m) -> long_double_t 
				//	{return sqrt((2*log(1/eps_bar_prime) + 2*log(m+1))/m);};
				//auto gamma = [eps_bar_prime] (int m) -> long_double_t 
				//	{return sqrt(log(2/eps_bar_prime)/(2*m));};
				//auto gamma = [eps_bar_prime] (int m) -> long_double_t 
				//	{return sqrt(log(2/(1-pow(1-eps_bar_prime,1.0/6)))/(2*m));};
				auto gamma = [eps_bar_prime] (int m) -> long_double_t 
					{return sqrt(log(2/1.67e-10)/(2*m));};

				auto G01 = gamma(C01); if (isnan(G01)){cout << "G01 is NaN" << endl; G01=0;}
				auto G10 = gamma(C10); if (isnan(G10)){cout << "G10 is NaN" << endl; G10=0;}
				auto G0A = gamma(C0A); if (isnan(G0A)){cout << "G0A is NaN" << endl; G0A=0;}
				auto G1A = gamma(C1A); if (isnan(G1A)){cout << "G1A is NaN" << endl; G1A=0;}
				auto GA0 = gamma(CA0); if (isnan(GA0)){cout << "GA0 is NaN" << endl; GA0=0;}
				auto GAB = gamma(CAB); if (isnan(GAB)){cout << "GAB is NaN" << endl; GAB=0;}


				//these should be changed to whatever we determine
				start01 = skip01 ?  G01 : -G01; // free
				start10 = skip10 ?  G10 : -G10; // max
				start0A = skip0A ?  G0A : -G0A; // free
				start1A = skip1A ?  G1A : -G1A; // max
				startA0 = skipA0 ? -GA0 : -GA0; // min
				startAB = skipAB ?  GAB : -GAB; // max

				lim01 = skip01 ?  start01 : G01;
				lim10 = skip10 ?  start10 : G10;
				lim0A = skip0A ?  start0A : G0A;
				lim1A = skip1A ?  start1A : G1A;
				limA0 = skipA0 ?  startA0 : GA0;
				limAB = skipAB ?  startAB : GAB;


				if (smart_step){
					step01 = 2*G01/4;
					step10 = 2*G10/4;
					step0A = 2*G0A/4;
					step1A = 2*G1A/4;
					stepA0 = 2*GA0/4;
					stepAB = 2*GAB/4;

				}

				cout << "	Num__key = " << numkey << endl;
				cout << "	C01 = " << C01 << " step size= " << step01 <<" G01 = " << G01 << endl; 
				cout << "	C10 = " << C10 << " step size= " << step10 <<" G10 = " << G10 << endl; 
				cout << "	C0A = " << C0A << " step size= " << step0A <<" G0A = " << G0A << endl; 
				cout << "	C1A = " << C1A << " step size= " << step1A <<" G1A = " << G1A << endl; 
				cout << "	CA0 = " << CA0 << " step size= " << stepA0 <<" GA0 = " << GA0 << endl; 
				cout << "	CAB = " << CAB << " step size= " << stepAB <<" GAB = " << GAB << endl;

				int stat_count = 0;
				for (long_double_t P01 = start01; P01 <=lim01; P01+=step01){
					
					for (long_double_t P10 = start10; P10 <=lim10; P10+=step10){
						for (long_double_t P0A = start0A; P0A <=lim0A; P0A+=step0A){
			
							
							for (long_double_t P1A = start1A; P1A <=lim1A; P1A+=step1A){
								for (long_double_t PA0 = startA0; PA0 <=limA0; PA0+=stepA0){
									for (long_double_t PAB = startAB; PAB <=limAB; PAB+=stepAB){
										stat_count +=1;

										AtoB["P01"] = (B01 + P01) < 0 ? 0 : (B01 + P01) > 1 ? 1 : (B01 + P01);
										AtoB["P10"] = (B10 + P10) < 0 ? 0 : (B10 + P10) > 1 ? 1 : (B10 + P10);
										AtoB["P00"] = 1 - AtoB["P01"];
										AtoB["P11"] = 1 - AtoB["P10"];

										AtoB["P0A"] = (B0A + P0A) < 0 ? 0 : (B0A + P0A) > 1 ? 1 : (B0A + P0A);
										AtoB["P1A"] = (B1A + P1A) < 0 ? 0 : (B1A + P1A) > 1 ? 1 : (B1A + P1A);
										AtoB["P0B"] = 1 - AtoB["P0A"];
										AtoB["P1B"] = 1 - AtoB["P1A"];

										AtoB["PA0"] = (BA0 + PA0) < 0 ? 0 : (BA0 + PA0) > 1 ? 1 : (BA0 + PA0);
										AtoB["PA1"] = 1 - AtoB["PA0"];
										AtoB["PAB"] = (BAB + PAB) < 0 ? 0 : (BAB + PAB) > 1 ? 1 : (BAB + PAB);
										AtoB["PAA"] = 1 - AtoB["PAB"];

										

										ent = calculate_key_rate_finite_inside(alph, CAD, mismatched);
										
										if (ent < min_ent){
											min_ent = ent;
											M01 = P01;
											M10 = P10;
											M0A = P0A;
											M1A = P1A;
											MA0 = PA0;
											MAB = PAB;
										}
										if (stat_count % 10000 == 0){
											cout << "	Completed  " << stat_count << " iterations" << endl;
											
										}
										if (PAB+stepAB > limAB && PAB!=limAB){
											PAB = limAB;
											//cout << "	Forcing PAB" << endl;
										}//cout << "	PAB = " << PAB << "limit = " <<limAB <<endl;
										if (GAB==0){
											break;
										}										

									}
									if (PA0+stepA0 > limA0 && PA0!=limA0){
										PA0 = limA0;
										//cout << "	Forcing PA0" << endl;
									}//cout << "	PA0 = " << PA0 << "limit = " <<limA0 <<endl;
									if (GA0==0){
										break;
									}									

								}
								if (P1A+step1A > lim1A && P1A!=lim1A){
									P1A = lim1A;
									//cout << "	Forcing P1A" << endl;
								}//cout << "	P1A = " << P1A << "limit = " <<lim1A <<endl;
								if (G1A==0){
									break;
								}								
							}
							if (P0A+step1A > lim0A && P0A!=lim0A){
								P0A = lim0A;
								//cout << "	Forcing P0A" << endl;
							}//cout << "	P0A = " << P0A << "limit = " <<lim0A <<endl;
							if (G0A==0){
								break;
							}							

						}
						if (P10+step10 > lim10 && P10!=lim10){
							P10 = lim10;
							//cout << "	Forcing P10" << endl;
						}//cout << "	P10 = " << P10 << "limit = " <<lim10 <<endl;
						if (G10==0){
							break;
						}						

					}
					if (P01+step01 > lim01 && P01!=lim01){
						P01 = lim01;
						//cout << "	Forcing P01" << endl;
					}//cout << "	P01 = " << P01 << "limit = " <<lim01 <<endl;
					if (G01==0){
						break;
					}					

				}
				cout << "Searched stat space of size " << stat_count << endl;
				cout << "Found min ent on eps iteration "  << eps_count << endl;
				cout << "	Eps_bar = "		 << eps_bar << endl;
				cout << "	Eps_bar_prime = "   << eps_bar_prime << endl;
				cout <<"	C01 = " << C01 <<" gamma(C01) = "<< gamma(C01)<< " M01= " << M01 << endl;
				cout <<"	C10 = " << C10 <<" gamma(C10) = "<< gamma(C10)<< " M10= " << M10 << endl;
				cout <<"	C0A = " << C0A <<" gamma(C0A) = "<< gamma(C0A)<< " M0A= " << M0A << endl;
				cout <<"	C1A = " << C1A <<" gamma(C1A) = "<< gamma(C1A)<< " M1A= " << M1A << endl;
				cout <<"	CA0 = " << CA0 <<" gamma(CA0) = "<< gamma(CA0)<< " MA0= " << MA0 << endl;
				cout <<"	CAB = " << CAB <<" gamma(CAB) = "<< gamma(CAB)<< " MAB= " << MAB << endl;
				cout <<" 	Min ent = " << min_ent << endl;

				long_double_t delta = 2*safe_log2(1/((2*(eps-epsEC-eps_bar))))+7*sqrt(numkey*safe_log2(2/(eps_bar-eps_bar_prime)));
				//double key_rate = min_ent - leakEC - delta/numkey;
				long_double_t QBER;
				QBER = (AtoB["PAB"]+AtoB["P01"])/(AtoB["PAB"]+AtoB["P01"]+AtoB["PA1"]+AtoB["P0B"]);
				cout << "Q " << noise << endl;
				cout << "QBER " << QBER << endl;
				
				auto worstcaseQBER = (B01+BAB + gamma(C01)+gamma(CAB))/(B01+BAB+1-B0A+1-BA0-gamma(C01)-gamma(CAB)-gamma(C0A)-gamma(CA0));
				leakEC = 1.2*h_entropy(worstcaseQBER); //this is what they suggest in the paper -amount of info leaked during EC

				cout << "Worstcase QBER = "<< worstcaseQBER << endl;
				cout << "Normal QBER was = " << QBER << endl;;
				double key_rate = min_ent - leakEC -delta/numkey;

				cout << "for eps_bar = " << eps_bar << " and eps_bar_prime = " << eps_bar_prime << endl;
				cout << "	[-] Min Entropy is " << min_ent << endl;
				cout << "	[-] Key rate = "<<  min_ent << " - " << leakEC << 
					" - " << delta<<"/"<<numkey  << " = " << key_rate << endl;

				if (key_rate > best_key_rate){
					cout << "[+] New best Key rate = "<<  min_ent << " - " << leakEC << " - " << delta<<"/"<<numkey  << " = " << key_rate <<endl;;
					best_key_rate = key_rate;
					cout << "[+] eps_bar " << eps_bar << endl;
					cout << "[+] eps_bar_prime " << eps_bar_prime <<endl;
					best_key_ent = min_ent;
					MM01 = M01;
					MM10 = M10;
					MM0A = M0A;
					MM1A = M1A;
					MMA0 = MA0;
					MMAB = MAB;
					best_eps_bar_prime = eps_bar_prime;
					best_eps_bar = eps_bar;
				}

				if (eps_bar_prime+eps_bar_prime_step > eps_bar){
					eps_bar_prime = eps_bar-eps_bar_prime_step;
				} 
				
			}


			if (eps_bar+eps_bar_step > eps-epsEC){ //dont need this because it will never be equal} && eps_bar!=eps-epsEC){
				eps_bar = eps-epsEC - eps_bar_step;
			}
		}

		//auto gamma = [best_eps_bar_prime](int m) -> long_double_t 
		//{return sqrt((2*log(1/best_eps_bar_prime) + 2*log(m+1))/m);};
		//auto gamma = [best_eps_bar_prime] (int m) -> long_double_t 
		//			{return sqrt(log(2/best_eps_bar_prime)/(2*m));};
		//auto gamma = [eps_bar_prime] (int m) -> long_double_t 
		//			{return sqrt(log(2/(1-pow(1-eps_bar_prime,1.0/6)))/(2*m));};
		auto gamma = [eps_bar_prime] (int m) -> long_double_t 
					{return sqrt(log(2/1.67e-10)/(2*m));};

		cout << "Q = " << Q << " Alph^2 = " << a2 << " N = " << N << endl;
		cout << "Searched for best key rate with paramaters: " << endl;
		cout << "eps = " << eps << ", epsEC = " << epsEC << endl;
		cout << "Searched with step size = " << eps_bar_step << " for best eps_bar in range " <<  Ebar_start << " to " << Ebar_lim << endl;
		cout << "Searched with step size = " << eps_bar_prime_step << " for best eps_bar_prime in range " << Ebar_prime_start << " to " << Ebar_prime_lim << endl;

		cout << "Best Key rate was found at eps_bar = " << 
			best_eps_bar << " and eps_bar_prime= " << best_eps_bar_prime << endl;

		cout << "Search for P01 with steps " << step01 << " in range " << start01 << " to " << lim01 << " best was " << MM01 << " " << (fabs(MM01) == gamma(C01) ? (MM01 < 0 ? "min" : "max") : "free") << endl;
		cout << "Search for P10 with steps " << step10 << " in range " << start10 << " to " << lim10 << " best was " << MM10 << " " << (fabs(MM10) == gamma(C10) ? (MM10 < 0 ? "min" : "max") : "free") << endl;
		cout << "Search for P0A with steps " << step0A << " in range " << start0A << " to " << lim0A << " best was " << MM0A << " " << (fabs(MM0A) == gamma(C0A) ? (MM0A < 0 ? "min" : "max") : "free") << endl;
		cout << "Search for P1A with steps " << step1A << " in range " << start1A << " to " << lim1A << " best was " << MM1A << " " << (fabs(MM1A) == gamma(C1A) ? (MM1A < 0 ? "min" : "max") : "free") << endl;
		cout << "Search for PA0 with steps " << stepA0 << " in range " << startA0 << " to " << limA0 << " best was " << MMA0 << " " << (fabs(MMA0) == gamma(CA0) ? (MMA0 < 0 ? "min" : "max") : "free") << endl;
		cout << "Search for PAB with steps " << stepAB << " in range " << startAB << " to " << limAB << " best was " << MMAB << " " << (fabs(MMAB) == gamma(CAB) ? (MMAB < 0 ? "min" : "max") : "free") << endl;

		long_double_t delta = 2*safe_log2(1/((2*(eps-epsEC-best_eps_bar))))+7*sqrt(numkey*safe_log2(2/(best_eps_bar-best_eps_bar_prime)));

		cout << "Finite min_ent was " <<  best_key_ent << endl;
		cout << "Best_key Rate comp was " << best_key_ent << " - " << leakEC << 
					" - " << delta<<"/"<<numkey  << " = " << best_key_rate << endl;
		cout << "Best Key Rate is " << best_key_rate << endl;
		cout << "Numkey/N = " << (double) numkey/N ;
		cout << "Best Effective Key Rate is" << numkey * best_key_rate / N << endl;
		if (best_key_ent == 1)
			best_key_rate = -1;
		//return best_key_rate;
		return numkey * best_key_rate / N;

	}



	

		

	//};
	double noise_tolerance(double alph){
		double noise = STEP;
		while (calculate_key_rate_asymptotic(alph, noise) > 0)
			noise+=STEP/100;
		noise-=STEP/100;
		cout << "Maximum Noise tolerance for alpha = "<< alph <<" is " << noise <<endl;
		return noise;
	}
	vector<long_double_t> optimize_key_rate_finite(long_double_t start_alph, long_double_t start_Penc, long_double_t noise, long long N){
		long_double_t Penc = start_Penc;
		long_double_t alph = start_alph;
		
		long_double_t best_finite_key = -float('inf');
		long_double_t best_finite_key_Penc = -float('inf');
		long_double_t prev_finite_key_Penc = -float('inf');
		long_double_t prev_finite_key_alpha = -float('inf');

		long_double_t best_finite_alph = start_alph;
		long_double_t best_finite_alph_penc = .99L;
		long_double_t best_alph = .01;

		long_double_t bestPenc = .01;
		long_double_t key_rate;

		cout << "[*] Starting Noise = " << noise << endl;
		while(Penc < .96L){
			prev_finite_key_alpha = -float('inf');
			best_finite_key_Penc = -float('inf');
			alph = min(start_alph, best_finite_alph);
			cout << "  [*] Starting Penc = " << Penc << endl;
			while(alph < .96L){
				streambuf *old = cout.rdbuf(); // <-- save        
			    stringstream ss;

			    cout.rdbuf (ss.rdbuf());       // <-- redirect
				key_rate = calculate_key_rate_finite_outside(alph, noise, N, Penc);
				cout.rdbuf (old);              // <-- restore
				cout << "    [-] For noise = " << noise << " key rate was " << key_rate << " at alph = " << alph << " and penc = " << Penc <<endl;
				if (key_rate < prev_finite_key_alpha && best_finite_key_Penc > 0){
					cout << "    [X] Found alpha inflection point, breaking" << endl;
					break;
				}
				if(key_rate > best_finite_key_Penc){
					cout << "    [+] For noise = " << noise << " new best key rate over alpha is " << key_rate << " found at alph = " << alph << " and penc = " << Penc <<endl;
					best_finite_key_Penc = key_rate;
					best_finite_alph_penc = alph;
				}
				
				alph+=.02;
				prev_finite_key_alpha = key_rate;
			}
			cout << "  [-] For noise = " << noise << " best key rate over alpha was " << best_finite_key_Penc << " at alph = " << best_finite_alph_penc << " and penc = " << Penc <<endl;
			if (false &&best_finite_key_Penc < prev_finite_key_Penc){
				cout << "  [X] Found Penc inflection point, breaking" << endl;
				break;
			}
			if (best_finite_key_Penc > best_finite_key){
				best_finite_alph = best_finite_alph_penc;
				bestPenc = Penc;
				best_finite_key = best_finite_key_Penc;
				cout << "  [+] For noise = " << noise << " new best key rate over Penc is " << best_finite_key_Penc << " found at alph = " << best_finite_alph << " and penc = " << Penc <<endl;
				

			}
			else if (best_finite_key > 0 && best_finite_key_Penc <= 0){
				cout << "  [*] Was 0 over all alpha, previous best was " << best_finite_key << "  breaking" << endl;
				break;
			}
			prev_finite_key_Penc = best_finite_key_Penc;
			Penc +=.02;
		}/*
		if (best_finite_key < 0){
			best_finite_key = 0;
			best_finite_alph = 0;
			bestPenc = 0;
		}*/
		vector<long_double_t>  ret {best_finite_key, best_finite_alph, bestPenc};
		
		
		return ret;
	
	}
	void data_gen(long N, ofstream& data_file){
		long_double_t noise = 0.01;
		long_double_t alph = .01;
		long_double_t asym_key_rate = 0;

		long_double_t best_alph = 0;
		vector <long_double_t> best_asym_rates;
		vector <long_double_t> best_asym_effective;
		vector <long_double_t> best_asym_alph;

		data_file << "optimal effective key rate at various noise " << noise << endl;
		data_file << " N = " << N << endl;
		data_file << "Q, asym_Key_rate, Alph"<<endl;
		
		while ( noise <= .1101){
			alph = .01;
			asym_key_rate = -float('inf');
			while (alph < 1){
				streambuf *old = cout.rdbuf(); // <-- save        
			    stringstream ss;

			    cout.rdbuf (ss.rdbuf());       // <-- redirect
				long_double_t key_rate = calculate_key_rate_asymptotic(alph, noise);
			    cout.rdbuf (old);              // <-- restore

				if (key_rate > asym_key_rate){
					//cout << " 	For noise = " << noise << " new best key rate was " << key_rate << "found at alph = " << alph <<endl;

					asym_key_rate = key_rate;
					best_alph = alph;
				}
				alph+=.01;
			}

			best_asym_rates.push_back(asym_key_rate);
			best_asym_alph.push_back(best_alph);
			cout << "[+] For noise = " << noise << " Best key rate was " << asym_key_rate << " found at alph = " << best_alph <<endl;
			data_file << noise << ", " << asym_key_rate << ", " << best_alph << endl;
			noise+=.01;
		}

		vector <long_double_t> best_finite_rates;
		vector <long_double_t> best_finite_alphas;
		vector <long_double_t> best_finite_pencs;

		long_double_t best_finite_key = .01L;
		long_double_t best_finite_key_Penc = .01L;
		long_double_t prev_finite_key_Penc = .01L;
		long_double_t prev_finite_key_alpha = .01L;

		long_double_t best_finite_alph =.01L;
		long_double_t best_finite_alph_penc =.01L;

		long_double_t bestPenc =.01L;
		long_double_t Penc;
		long_double_t key_rate;

		noise = 0.01;
		data_file << endl;
		data_file << "Q, finite_key_rate, Alph, Penc" << endl;

		while (noise <= .1101){
			best_finite_key = 0;
			Penc = min(.01L, bestPenc);
			best_finite_key_Penc = 0.01L;
			prev_finite_key_Penc = 0.01L;
			cout << "[*] Starting Noise = " << noise << endl;
			
			auto ret = optimize_key_rate_finite(.05L, .05L, noise, N);
			cout << ret[0] << endl;
			
			best_finite_key = ret[0];
			best_finite_alph = ret[1];
			bestPenc = ret[2];
			
			best_finite_alphas.push_back(best_finite_alph);
			best_finite_pencs.push_back(bestPenc);
			best_finite_rates.push_back(best_finite_key);
			cout << "[+] For noise = " << noise << " best key rate was " << best_finite_key << " found at alph = " << best_finite_alph << " and penc = " << bestPenc <<endl;
			data_file << noise << ", " << best_finite_key << ", " << best_finite_alph << ", " << bestPenc << endl;
			noise+=.01;

		}

	}

	
	void data_gen2(long_double_t noise, ofstream& data_file){
		double asym_key_rate = 0;
		long_double_t alph = .01L;
		long_double_t best_alph = .99L;
		data_file << "optimal effective key rate at various N at noise " << noise << endl;
		while (alph < 1){
			streambuf *old = cout.rdbuf(); // <-- save        
		    stringstream ss;

		    cout.rdbuf (ss.rdbuf());       // <-- redirect
			long_double_t key_rate = calculate_key_rate_asymptotic(alph, noise);
		    cout.rdbuf (old);              // <-- restore

			if (key_rate > asym_key_rate){
				//cout << " 	For noise = " << noise << " new best key rate was " << key_rate << "found at alph = " << alph <<endl;

				asym_key_rate = key_rate;
				best_alph = alph;
			}
			alph+=.01;
		}

		cout << "[+] For noise = " << noise << " Best asym key rate was " << asym_key_rate << " found at alph = " << best_alph <<endl;
		data_file << "10^n, key_rate, alpha, Penc" << endl;
		data_file << "assym" << ", " << asym_key_rate << ", " << best_alph << ", 1" << endl;

		vector<long double> N_powers {1e6, 5e6, 1e7, 5e7, 1e8, 5e8, 1e9, 5e9};
		for (auto n : N_powers){
			auto N = (long) n;
			cout << "[*] Starting N= " << N << endl;
			auto ret = optimize_key_rate_finite(.05, .05L, noise, N);
			long_double_t best_finite_key = ret[0];
			long_double_t best_finite_alph = ret[1];
			long_double_t bestPenc = ret[2]; 
			cout << "[+] For N = " << N << " best key rate was " << best_finite_key << " found at alph = " << best_finite_alph << " and penc = " << bestPenc <<endl;
			data_file << N << ", " << best_finite_key << ", " << best_finite_alph << ", " << bestPenc << endl;
			

		}


	}

int main(int argc, char** argv){
	using namespace std;
	//using namespace Extb92Finite;
	//auto keyrate = calculate_key_rate_asymptotic(sqrt(atof(argv[1])), atof(argv[2]));
	//cout << " asymtptoic key rate = " << keyrate << endl;
	for (auto i: AtoB)
		cout << i.first << " " << i.second << endl;
	//noise_tolerance(atof(argv[1]));
	//calculate_key_rate_finite_outside(sqrt(atof(argv[1])), atof(argv[2]), atol(argv[3]));
	//for (auto i: AtoB)
	//	cout << i.first << " " << i.second << endl;
	//cout << "asymtptoic key rate = " << keyrate << endl;

	if(argc == 5){
		cout << "Finite Key_rate = " << calculate_key_rate_finite_outside(sqrt(atof(argv[1])), atof(argv[2]), atol(argv[3]), atof(argv[4])) << endl;

	}
	if(argc == 3){
		ofstream data_file;
		data_file.open(argv[2], ios::out);

		data_gen(atol(argv[1]), data_file);
		//data_gen2(atof(argv[1]), data_file);
	}
	if(argc == 4){
		cout << "asymptotic key_rate = " << calculate_key_rate_asymptotic(sqrt(atof(argv[1])), atof(argv[2])) << endl;
for (auto i: AtoB)
		cout << i.first << " " << i.second << endl;
	}


	return 0;
}