#pragma once

Complex ixp(Complex y)
{
	return exp(I*y);
}

double wb(int n)
{
	return 2.*pi*T*n;
}
//Complex f1(Complex b, Complex c, int m)
//{
//	Complex den;
//	den = pow((b - I*wb(m)), 2) - c*c;
//	return -tanh(b/2./T)/(4.*b*den);
//}

//Complex f2(Complex b, Complex c, Complex d, int m)
//{
//	Complex den;
//	den = pow((b - I*wb(m)), 2) - c*c;
//	return -(b-d)*(b-d-I*wb(m))*tanh(b / 2. / T) / (4.*b*den);
//}
//
//Complex fD1(Complex b, Complex c, int m)
//{
//	Complex res;
//	res = f1(b,c,m)+f1(-b,c,m)+f1(c,b,m)+f1(-c,b,m);
//	return res;
//}
//
//Complex fD2(Complex b, Complex c, Complex d, int m)
//{
//	Complex res;
//	res = f2(b,c,d,m) + f2(-b,c,d,m) + f2(c,b,-d,m) + f2(-c,b,-d,m);
//	return res;
//}

double swP(Complex z, int kx, int ky, double mur)
{
	double eps[N];
	for (int px = 0; px < Nx; px++)
	{
		for (int py = 0; py < Nx; py++)
		{
			eps[Nx*px + py] = -2.*(cos(dk*px) + cos(dk*py)
				- 4.*tp*(cos(dk*px) + cos(dk*py))) - mur;
		}
	}
	//cout << "mu=" << mur << ";   U=" << U << ";   T=" << T << endl;
	int k = Nx*kx + ky;
	Complex Sig = 0.;
	for (int l = 0; l < N; l++)
	{
		for (int m = 0; m < N; m++)
		{
			double e1, e2, e3;
			e1 = eps[l]; e2 = eps[m]; e3 = eps[vsum(k, ref(l), ref(m))];
			if (abs(e2 - e3) < .000000001)
			{
				Sig += -.25*U*U / N / N * pow(cosh(e2 / (2.*T)), -2) / (z + e1);
			}
			else
			{
				Sig += -.25*U*U / N / N*(tanh(e2 / (2.*T)) - tanh(e3 / (2.*T)))
					*(tanh(e1 / (2.*T)) - 1. / tanh((e3 - e2) / (2.*T)))
					/ (z + e1 + e2 - e3);
			}
		}
	}
	Complex InG, G;
	InG = z + eps[k] + Sig;
	G = -1. / InG;
	cout << "GPARA={" << G.real() << "," << G.imag()<<"}" << endl;
	return 1./pi*G.imag();
}

Complex ss1a(Complex a, Complex b, Complex c, Complex z, double mur)
{
	Complex res,den1,den2,den3,den4;
	den1 = a*a+2.*a*b+b*b-c*c-2.*a*z-2.*b*z+z*z+2.*a*mur+2.*b*mur-2.*z*mur+mur*mur;
	den2 = a*a+b*b-2.*a*b-c*c-2.*a*mur+2.*b*mur+mur*mur+2.*a*z-2.*b*z-2.*mur*z+z*z;
	den3 = -a*a + b*b - 2.* b*c + c*c + 2.*b*mur - 2.*c*mur + mur*mur-2.*b*z + 2.*c*z -2.*mur*z + z*z;
	den4 = -a*a+b*b+2.*b*c + c*c+2.*b*mur + 2.*c*mur + mur*mur-2.*b*z-2.*c*z-2.*mur*z + z*z;
	res = -tanh((a+mur)/(2.*T))/den1/a+tanh((mur-a)/(2.*T))/den2/a-1./tanh((-b+c)/(2.*T))/den3/c
		+ 1./tanh((-b-c)/(2.*T))/den4/c;
	return tanh((b+mur)/(2.*T))/(16.*b)*res;
}

Complex ss1b(Complex a, Complex b, Complex c, Complex z, double mur)
{
	Complex res, den1, den2, den3, den4;
	den1 = a*a-2.*a*c-b*b+c*c - 2.*a*z+2.*c*z + z*z + 2.*a*mur-2.*c*mur - 2.*z*mur + mur*mur;
	den2 = a*a - b*b + 2.*a*c + c*c - 2.*a*mur - 2.*c*mur + mur*mur + 2.*a*z + 2.*c*z - 2.*mur*z + z*z;
	den3 = -a*a + b*b + 2.* b*c + c*c - 2.*b*mur - 2.*c*mur + mur*mur + 2.*b*z + 2.*c*z - 2.*mur*z + z*z;
	den4 = -a*a + b*b - 2.*b*c + c*c + 2.*b*mur - 2.*c*mur + mur*mur - 2.*b*z + 2.*c*z - 2.*mur*z + z*z;
	res = -tanh((a + mur)/(2.*T))/ den1 / a + tanh((mur - a)/(2.*T))/den2 / a - 1./ tanh((b + c)/(2.*T))/den3 /b
		  + 1./ tanh((-b+c) /(2.*T))/den4 /b;
	return tanh((c + mur) / (2.*T)) / (16.*c)*res;
}

Complex ss2a(Complex a, Complex b, Complex c, Complex d, Complex z, double mur)
{
	Complex res,zmu, num1, num2, num3, num4, den1 , den2,den3, den4;
	zmu = z - mur;
	
	num1 = (a-d)*(a+b+d-zmu)*tanh((a+mur)/(2.*T));
	num2 = (a+d)*(a-b-d+zmu)*tanh((-a+mur)/(2.*T));
	num3 = (c+d)*(b-c+d-zmu)/tanh((-b+c)/(2.*T));
	num4 = (c-d)*(b+c+d-zmu)/tanh((-b-c)/(2.*T));

	den1 = a*a + 2.*a*b + b*b - c*c - 2.*a*z - 2.*b*z + z*z + 2.*a*mur + 2.*b*mur - 2.*z*mur + mur*mur;
	den2 = a*a + b*b - 2.*a*b - c*c - 2.*a*mur + 2.*b*mur + mur*mur + 2.*a*z - 2.*b*z - 2.*mur*z + z*z;
	den3 = -a*a + b*b - 2.* b*c + c*c + 2.*b*mur - 2.*c*mur + mur*mur - 2.*b*z + 2.*c*z - 2.*mur*z + z*z;
	den4 = -a*a + b*b + 2.*b*c + c*c + 2.*b*mur + 2.*c*mur + mur*mur - 2.*b*z - 2.*c*z - 2.*mur*z + z*z;
	res = -num1/a/den1 + num2/a/den2+num3/c/den3 + num4/c/den4;
	return (b + d)*tanh((b+mur)/(2.*T))/(16.*b)*res;
}

Complex ss2b(Complex a, Complex b, Complex c, Complex d, Complex z, double mur)
{
	Complex res,zmu,num1,num2,num3,num4, den1, den2, den3, den4;
	zmu = z - mur;
	num1 = (a-d)*(a-c-d-zmu)*tanh((a + mur)/(2.*T));
	num2 = (a+d)*(a+c+d+zmu)*tanh((-a + mur)/(2.*T));
	num3 = (b-d)*(b+c-d+zmu)/tanh((b+c)/(2.*T));
	num4 = (b+d)*(b-c+d-zmu)/tanh((c-b)/(2.*T));
	den1 = a*a - 2.*a*c - b*b + c*c - 2.*a*z + 2.*c*z + z*z + 2.*a*mur - 2.*c*mur - 2.*z*mur + mur*mur;
	den2 = a*a - b*b + 2.*a*c + c*c - 2.*a*mur - 2.*c*mur + mur*mur + 2.*a*z + 2.*c*z - 2.*mur*z + z*z;
	den3 = -a*a + b*b + 2.* b*c + c*c - 2.*b*mur - 2.*c*mur + mur*mur + 2.*b*z + 2.*c*z - 2.*mur*z + z*z;
	den4 = -a*a + b*b - 2.*b*c + c*c + 2.*b*mur - 2.*c*mur + mur*mur - 2.*b*z + 2.*c*z - 2.*mur*z + z*z;
	res = num1/a/den1-num2/a/den2+num3/b/den3-num4/b/den4;
	return (c+d)*tanh((c+mur)/(2.*T))/(16.*c)*res;
}

Complex SS1(Complex a, Complex b, Complex c, Complex z, double mur)
{
	Complex res,den1,den2,zmu,num1,num2a,num2b,num2c,num2d;
	if (abs(b - c) < .00000001)
	{
		zmu = z - mur;
		num1 = pow(tanh((b+mur)/(2.*T)),2)+pow(tanh((-b+mur)/(2.*T)),2) - 2.;
		den1 = a*a-pow(z-mur,2); 
		num2a = tanh((b + mur)/(2.*T))-tanh((-b+mur)/(2.*T));
		num2b = a*(a*a - 4.*b*b - zmu*zmu) / tanh(b / T);
		num2c = (a*a - 4.*b*b+zmu*zmu)*sinh(a/T)+2.*a*zmu*sinh(mur / T);
		num2d = cosh((a - mur) / (2.*T))*cosh((a + mur) / (2.*T));
	    den2 = a*(a-2.*b-zmu)*(a+2.*b-zmu)*(a-2.*b+zmu)*(a+2.*b+zmu);
	    res = 1./(16.*b*b)*(num1/den1 + 2.*num2a*(num2b-b*num2c/num2d)/den2);
	}
	else
	{
		res=ss1a(a, b, c, z, mur) + ss1a(a, -b, c, z, mur) + ss1b(a, b, c, z,mur) + ss1b(a, b, -c, z,mur);
	}
	return res;
}

Complex SS2(Complex a, Complex b, Complex c, Complex d, Complex z, double mur)
{
	Complex res,zmu,den2, num1a,num1b, num2;
	if (abs(b - c) < .00000001)
	{
		zmu = z-mur;
		num1a = ((2.* b + d - zmu)) / (-a*a + pow(2.* b - zmu, 2))
			   - ((2.* b - d + zmu)) / (-a*a + pow(2.* b + zmu, 2));
		num1b =-2.*b*(a+d)*tanh((-a + mur)/(2.*T))/((a-2.*b+zmu)*(a+2.*b+zmu)) 
			   +2.*b*(-a+d)*tanh((a + mur)/(2.*T))/((a-2.*b-zmu)*(a+2.*b-zmu));
		num2 = (d-zmu)*(pow((b-d)/cosh((-b+mur)/(2.*T)),2)+pow((b+d)/cosh((b+mur)/(2.*T)),2));
		den2 = -a*a + zmu*zmu;
		res = ((b + d)*(b - d)*(num1a/tanh(b/T)+num1b/a)*(tanh((-b+mur)/(2.*T))-tanh((b+mur)/(2.*T)))
			      -num2/den2)/(16.*b*b);
	}
	else
	{
		res = ss2a(a, b, c,d, z,mur) + ss2a(a, -b, c,d, z,mur) + ss2b(a,b,c,d,z,mur) + ss2b(a,b,-c,d,z,mur);
	}
	return res;
}

MKL_Complex16* StdToMkl(Complex* complexPtr)
{
	return reinterpret_cast<MKL_Complex16*>(complexPtr);
}

int mx(int w)
{
	return (w + 12 * Nx) % Nx;
}

int my(int w)
{
	return (w + 12 * Nx) % Nx;
}

int cy(int w)
{
	return w % Nx;
}

int cx(int w)
{
	return (w - cy(w)) / Nx;
}

int vsum(int a, int b)
{
	return Nx*mx(cx(a) + cx(b)) + my(cy(a) + cy(b) );
}

int vsum(int a, int b, int c)
{
	return vsum(vsum(a, b), c);
}

int mxAF(int w)
{
	return (w + 20 * Nx) % (Nx/2);
}

double swAF(Complex z, int Kx, int Ky, double pg, double mur)
{
	//gamma and x
	ComplexVectorAF  gam0AA, gam0AB, gam0BA, gam0BB, gam1AA, gam1AB, gam1BA, gam1BB;
	Complex Sig0AA = 0., Sig0AB = 0., Sig0BA = 0., Sig0BB = 0.,
		Sig1AA = 0., Sig1AB = 0., Sig1BA = 0., Sig1BB = 0.;
	double x[N];

	for (int px = 0; px < Nx / 2; px++)
	{
		for (int py = 0; py < Nx; py++)
		{
			int p = Nx*px + py;
			double ep, epp;
			ep = -4.* tp*cos(dk* px)*cos(dk*(px - py));
			//epp = -2.*tpp*(cos(4.*pi / Nx * px) + cos(4.*pi / Nx * (px - py)));
			gam0AA[p] = gam1BB[p] = -z -ep + mur + pg;
			gam0AB[p] = gam1AB[p] = 1. + ixp(2.*dk*px) + ixp(dk*py) + ixp(dk*(2.*px - py));
			gam0BA[p] = gam1BA[p] = 1. + ixp(-2.*dk*px) + ixp(-dk*py) + ixp(-dk*(2.*px - py));
			gam0BB[p] = gam1AA[p] = -z - ep + mur - pg;
			x[p] = sqrt(pg*pg + 4.*pow(cos(dk*px) + cos(dk*(px - py)),2));
		}
	}
	
	int kx = mxAF(Kx); int ky = my(Kx + Ky);
	int k = Nx * kx + ky;
	//self energy
	for (int px = 0; px < Nx / 2; px++)
	{
		for (int py = 0; py < Nx; py++)
		{
			int p = Nx*px + py;
			for (int qx = 0; qx < Nx / 2; qx++)
			{
				for (int qy = 0; qy < Nx; qy++)
				{
					int q = Nx*qx + qy;
					int rx = kx - px - qx; int ry = ky - py - qy;
					int r = Nx*mxAF(rx) + my(ry);
					Complex phase, cphase;
					phase = (1. + ixp(2.*dk*px) + ixp(dk*py) + ixp(dk*(2.*px - py)))
						*(1. + ixp(2.*dk*qx) + ixp(dk*qy) + ixp(dk*(2.*qx - qy)))
						*(1. + ixp(2.*dk*rx) + ixp(dk*ry) + ixp(dk*(2.*rx - ry)));
					cphase = phase.real() - phase.imag()*I;
					Sig0AA += pow(2.*U / N, 2)*SS2(x[p], x[q], x[r], -pg, z,mur);
					Sig0AB += pow(2.*U / N, 2)*phase*SS1(x[p], x[q], x[r], z,mur);
					Sig0BA += pow(2.*U / N, 2)*cphase*SS1(x[p], x[q], x[r], z,mur);
					Sig0BB += pow(2.*U / N, 2)*SS2(x[p], x[q], x[r], pg, z,mur);

					Sig1AA += pow(2.*U / N, 2)*SS2(x[p], x[q], x[r], pg, z,mur);
					Sig1AB += pow(2.*U / N, 2)*phase*SS1(x[p], x[q], x[r], z,mur);
					Sig1BA += pow(2.*U / N, 2)*cphase*SS1(x[p], x[q], x[r], z,mur);
					Sig1BB += pow(2.*U / N, 2)*SS2(x[p], x[q], x[r], -pg, z,mur);
				}
			}
		}
	}
	Complex ig0AA, ig1AA, ig0AB, ig1AB, ig0BA, ig1BA, ig0BB, ig1BB,
		g0AA, g1AA, g0AB, g1AB, g0BA, g1BA, g0BB, g1BB, det;
	ig0AA = gam0AA[k] + Sig0AA; ig1AA = gam1AA[k] + Sig1AA;
	ig0AB = gam0AB[k] + Sig0AB; ig1AB = gam1AB[k] + Sig1AB;
	ig0BA = gam0BA[k] + Sig0BA; ig1BA = gam1BA[k] + Sig1BA;
	ig0BB = gam0BB[k] + Sig0BB; ig1BB = gam1BB[k] + Sig1BB;

	det = ig0AA*ig0BB - ig0AB*ig0BA;
	g0AA = ig0BB / det; g1AA = ig1BB / det;
	g0AB = -ig0AB / det; g1AB = -ig1AB / det;
	g0BA = -ig0BA / det; g1BA = -ig1BA / det;
	g0BB = ig0AA / det; g1BB = ig1AA / det;

	//symmetrization
	Complex Gsym;
	Gsym = .25*(g0AA + g1AA + g0BB + g1BB + ixp(-dk*Kx)*(g0AB + g1AB) + ixp(dk*Kx)*(g0BA + g1BA));
	cout<<"GsymAF="<<Gsym.real()<<"+I("<< Gsym.imag()<<")"<<endl;
	return( 1./pi*Gsym.imag() );
}

void print(string_view s)
{
	auto t = std::time(nullptr);
	auto tm = *std::localtime(&t);
	std::ostringstream oss;
	oss << std::put_time(&tm, "%d-%m-%Y %H-%M-%S");
	auto str = oss.str();

	cout << str << ": " << s << endl;
}