__device__ static inline void cexp(double &re, double &im)
{
	const double e = exp(re);
	re = e * cos(im);
	im = e * sin(im);
}

__constant__ const double taylorarr24[24 * 12] = {
	// real part imaginary part , low order coefficients last
	// nsing = 0
	0.00000000000000000e-00,  3.00901111225470020e-01,
	5.00000000000000000e-01,  0.00000000000000000e-00,
	0.00000000000000000e-00, -7.52252778063675049e-01,
	-1.00000000000000000e-00,  0.00000000000000000e-00,
	0.00000000000000000e-00,  1.12837916709551257e+00,
	1.00000000000000000e-00,  0.00000000000000000e-00,
	// nsing = 1
	-2.22423508493755319e-01,  1.87966717746229718e-01,
	3.41805419240637628e-01,  3.42752593807919263e-01,
	4.66574321730757753e-01, -5.59649213591058097e-01,
	-8.05759710273191021e-01, -5.38989366115424093e-01,
	-4.88914083733395200e-01,  9.80580906465856792e-01,
	9.33757118080975970e-01,  2.82273885115127769e-01,
	// nsing = 2
	-2.60522586513312894e-01, -4.26259455096092786e-02,
	1.36549702008863349e-03,  4.39243227763478846e-01,
	6.50591493715480700e-01, -1.23422352472779046e-01,
	-3.43379903564271318e-01, -8.13862662890748911e-01,
	-7.96093943501906645e-01,  6.11271022503935772e-01,
	7.60213717643090957e-01,  4.93801903948967945e-01,
	// nsing = 3
	-1.18249853727020186e-01, -1.90471659765411376e-01,
	-2.59044664869706839e-01,  2.69333898502392004e-01,
	4.99077838344125714e-01,  2.64644800189075006e-01,
	1.26114512111568737e-01, -7.46519337025968199e-01,
	-8.47666863706379907e-01,  1.89347715957263646e-01,
	5.39641485816297176e-01,  5.97805988669631615e-01,
	// nsing = 4
	4.94825297066481491e-02, -1.71428212158876197e-01,
	-2.97766677111471585e-01,  1.60773286596649656e-02,
	1.88114210832460682e-01,  4.11734391195006462e-01,
	3.98540613293909842e-01, -4.63321903522162715e-01,
	-6.99522070542463639e-01, -1.32412024008354582e-01,
	3.33997185986131785e-01,  6.01983450812696742e-01,
	// nsing = 5
	1.18367078448232332e-01, -6.09533063579086850e-02,
	-1.74762998833038991e-01, -1.39098099222000187e-01,
	-6.71534655984154549e-02,  3.34462251996496680e-01,
	4.37429678577360024e-01, -1.59613865629038012e-01,
	-4.71863911886034656e-01, -2.92759316465055762e-01,
	1.80238737704018306e-01,  5.42834914744283253e-01,
	// nsing = 6
	8.87698096005701290e-02,  2.84339354980994902e-02,
	-3.18943083830766399e-02, -1.53946887977045862e-01,
	-1.71825061547624858e-01,  1.70734367410600348e-01,
	3.33690792296469441e-01,  3.97048587678703930e-02,
	-2.66422678503135697e-01, -3.18469797424381480e-01,
	8.48049724711137773e-02,  4.60546329221462864e-01,
	// nsing = 7
	2.99767046276705077e-02,  5.34659695701718247e-02,
	4.53131030251822568e-02, -9.37915401977138648e-02,
	-1.57982359988083777e-01,  3.82170507060760740e-02,
	1.98891589845251706e-01,  1.17546677047049354e-01,
	-1.27514335237079297e-01, -2.72741112680307074e-01,
	3.47906344595283767e-02,  3.82277517244493224e-01,
	// nsing = 8
	-7.35922494437203395e-03,  3.72011290318534610e-02,
	5.66783220847204687e-02, -3.21015398169199501e-02,
	-1.00308737825172555e-01, -2.57695148077963515e-02,
	9.67294850588435368e-02,  1.18174625238337507e-01,
	-5.21266530264988508e-02, -2.08850084114630861e-01,
	1.24443217440050976e-02,  3.19239968065752286e-01,
	// nsing = 9
	-1.66126772808035320e-02,  1.46180329587665321e-02,
	3.85927576915247303e-02,  1.18910471133003227e-03,
	-4.94003498320899806e-02, -3.93468443660139110e-02,
	3.92113167048952835e-02,  9.03306084789976219e-02,
	-1.82889636251263500e-02, -1.53816215444915245e-01,
	3.88103861995563741e-03,  2.72090310854550347e-01,
	// nsing = 10
	-1.21245068916826880e-02,  1.59080224420074489e-03,
	1.91116222508366035e-02,  1.05879549199053302e-02,
	-1.97228428219695318e-02, -3.16962067712639397e-02,
	1.34110372628315158e-02,  6.18045654429108837e-02,
	-5.52574921865441838e-03, -1.14259663804569455e-01,
	1.05534036292203489e-03,  2.37326534898818288e-01,
	// nsing = 11
	-5.96835002183177493e-03,-2.42594931567031205e-03,
	7.44753817476594184e-03,9.33450807578394386e-03,
	-6.52649522783026481e-03,-2.08165802069352019e-02,
	3.89988065678848650e-03,4.12784313451549132e-02,
	-1.44110721106127920e-03,-8.76484782997757425e-02,
	2.50210184908121337e-04,2.11131066219336647e-01,
	// nsing = 12
	-2.24505212235034193e-03,-2.38114524227619446e-03,
	2.36375918970809340e-03,5.97324040603806266e-03,
	-1.81333819936645381e-03,-1.28126250720444051e-02,
	9.69251586187208358e-04,2.83055679874589732e-02,
	-3.24986363596307374e-04,-6.97056268370209313e-02,
	5.17231862038123061e-05,1.90681117197597520e-01,
	// nsing = 13
	-6.76887607549779069e-04,-1.48589685249767064e-03,
	6.22548369472046953e-04,3.43871156746448680e-03,
	-4.26557147166379929e-04,-7.98854145009655400e-03,
	2.06644460919535524e-04,2.03107152586353217e-02,
	-6.34563929410856987e-05,-5.71425144910115832e-02,
	9.32252179140502456e-06,1.74167663785025829e-01,
	// nsing = 14
	-1.67596437777156162e-04,-8.05384193869903178e-04,
	1.37627277777023791e-04,1.97652692602724093e-03,
	-8.54392244879459717e-05,-5.23088906415977167e-03,
	3.78965577556493513e-05,1.52191559129376333e-02,
	-1.07393019498185646e-05,-4.79347862153366295e-02,
	1.46503970628861795e-06,1.60471011683477685e-01,
	// nsing = 15
	-3.45715760630978778e-05,-4.31089554210205493e-04,
	2.57350138106549737e-05,1.19449262097417514e-03,
	-1.46322227517372253e-05,-3.61303766799909378e-03,
	5.99057675687392260e-06,1.17993805017130890e-02,
	-1.57660578509526722e-06,-4.09165023743669707e-02,
	2.00739683204152177e-07,1.48879348585662670e-01,
	// nsing = 16
	-5.99735188857573424e-06,-2.42949218855805052e-04,
	4.09249090936269722e-06,7.67400152727128171e-04,
	-2.14920611287648034e-06,-2.60710519575546230e-03,
	8.17591694958640978e-07,9.38581640137393053e-03,
	-2.00910914042737743e-07,-3.54045580123653803e-02,
	2.39819738182594508e-08,1.38916449405613711e-01,
	// nsing = 17
	-8.80708505155966658e-07,-1.46479474515521504e-04,
	5.55693207391871904e-07,5.19165587844615415e-04,
	-2.71391142598826750e-07,-1.94439427580099576e-03,
	9.64641799864928425e-08,7.61536975207357980e-03,
	-2.22357616069432967e-08,-3.09762939485679078e-02,
	2.49806920458212581e-09,1.30247401712293206e-01,
	// nsing = 18
	-1.10007111030476390e-07,-9.35886150886691786e-05,
	6.46244096997824390e-08,3.65267193418479043e-04,
	-2.95175785569292542e-08,-1.48730955943961081e-03,
	9.84949251974795537e-09,6.27824679148707177e-03,
	-2.13827217704781576e-09,-2.73545766571797965e-02,
	2.26877724435352177e-10,1.22627158810895267e-01,
	// nsing = 19
	-1.17302439957657553e-08,-6.24890956722053332e-05,
	6.45231881609786173e-09,2.64799907072561543e-04,
	-2.76943921343331654e-09,-1.16094187847598385e-03,
	8.71074689656480749e-10,5.24514377390761210e-03,
	-1.78730768958639407e-10,-2.43489203319091538e-02,
	1.79658223341365988e-11,1.15870972518909888e-01,
	// nsing = 20
	-1.07084502471985403e-09,-4.31515421260633319e-05,
	5.54152563270547927e-10,1.96606443937168357e-04,
	-2.24423474431542338e-10,-9.21550077887211094e-04,
	6.67734377376211580e-11,4.43201203646827019e-03,
	-1.29896907717633162e-11,-2.18236356404862774e-02,
	1.24042409733678516e-12,1.09836276968151848e-01,
	// nsing = 21
	-8.38816525569060600e-11,-3.06091807093959821e-05,
	4.10033961556230842e-11,1.48895624771753491e-04,
	-1.57238128435253905e-11,-7.42073499862065649e-04,
	4.43938379112418832e-12,3.78197089773957382e-03,
	-8.21067867869285873e-13,-1.96793607299577220e-02,
	7.46725770201828754e-14,1.04410965521273064e-01,
	//nsing = 22
	-5.64848922712870507e-12,-2.22021942382507691e-05,
	2.61729537775838587e-12,1.14683068921649992e-04,
	-9.53316139085394895e-13,-6.05021573565916914e-04,
	2.56116039498542220e-13,3.25530796858307225e-03,
	-4.51482829896525004e-14,-1.78416955716514289e-02,
	3.91940313268087086e-15,9.95054815464739996e-02,
	//nsing = 23
	-3.27482357793897640e-13,-1.64138890390689871e-05,
	1.44278798346454523e-13,8.96362542918265398e-05,
	-5.00524303437266481e-14,-4.98699756861136127e-04,
	1.28274026095767213e-14,2.82359118537843949e-03,
	-2.16009593993917109e-15,-1.62538825704327487e-02,
	1.79368667683853708e-16,9.50473084594884184e-02
};

__device__ void faddeeva(double zre, double zim,
	double* out_real, double* out_imag)
{
	constexpr double npi24[24] = { // precomputed values n * pi
		0.00000000000000000e+00, 3.14159265358979324e+00, 6.28318530717958648e+00,
		9.42477796076937972e+00, 1.25663706143591730e+01, 1.57079632679489662e+01,
		1.88495559215387594e+01, 2.19911485751285527e+01, 2.51327412287183459e+01,
		2.82743338823081391e+01, 3.14159265358979324e+01, 3.45575191894877256e+01,
		3.76991118430775189e+01, 4.08407044966673121e+01, 4.39822971502571053e+01,
		4.71238898038468986e+01, 5.02654824574366918e+01, 5.34070751110264851e+01,
		5.65486677646162783e+01, 5.96902604182060715e+01, 6.28318530717958648e+01,
		6.59734457253856580e+01, 6.91150383789754512e+01, 7.22566310325652445e+01,
	};
	constexpr double a24[24] = { // precomputed Fourier coefficient prefactors
		2.95408975150919338e-01, 2.75840233292177084e-01, 2.24573955224615866e-01,
		1.59414938273911723e-01, 9.86657664154541891e-02, 5.32441407876394120e-02,
		2.50521500053936484e-02, 1.02774656705395362e-02, 3.67616433284484706e-03,
		1.14649364124223317e-03, 3.11757015046197600e-04, 7.39143342960301488e-05,
		1.52794934280083635e-05, 2.75395660822107093e-06, 4.32785878190124505e-07,
		5.93003040874588103e-08, 7.08449030774820423e-09, 7.37952063581678038e-10,
		6.70217160600200763e-11, 5.30726516347079017e-12, 3.66432411346763916e-13,
		2.20589494494103134e-14, 1.15782686262855879e-15, 5.29871142946730482e-17,
	};
	const double tm = 12;

	// catch singularities in the Fourier representation At
	// z = n pi / tm , and provide a Taylor series expansion in those
	// points , and only use it when we ’re close enough to the real axis
	// that there is a chance we need it
	const double zim2 = zim * zim;
	const double maxnorm = 9.0 / 1000000.0;
	if (zim2 < maxnorm) {
		// we ’re close enough to the real axis that we need to worry about
		// singularities
		const double dnsing = tm * zre / npi24[1];
		const double dnsingmax2 = (24.0 - 1.0 / 2.0) * (24.0 - 1.0 / 2.0);
		if (dnsing * dnsing < dnsingmax2) {
			// we ’re in the interesting range of the real axis as well ...
			// deal with Re(z) < 0 so we only need N different Taylor
			// expansions ; use w(-x+iy) = conj (w(x+ iy))
			const bool negrez = zre < 0.0;
			// figure out closest singularity
			const int nsing = int(abs(dnsing) + 1.0 / 2.0);
			// and calculate just how far we are from it
			const double zmnpire = abs(zre) - (nsing * M_PI);
			const double zmnpinorm = zmnpire * zmnpire + zim2;
			// close enough to one of the singularities ?
			if (zmnpinorm < maxnorm) {
				const double* coeffs = &taylorarr24[nsing * 12];
				// calculate value of taylor expansion ...
				// ( note : there ’s no chance to vectorize this one , since
				// the value of the next iteration depend on the ones from
				// the previous iteration )
				double sumre = coeffs[0], sumim = coeffs[1];
				for (unsigned i = 1; i < 6; ++i) {
					const double re = sumre * zmnpire - sumim * zim;
					const double im = sumim * zmnpire + sumre * zim;
					sumre = re + coeffs[2 * i + 0];
					sumim = im + coeffs[2 * i + 1];
				}
				// undo the flip in real part of z if needed
				if (negrez) sumim *= -1;

				*out_real = sumre;
				*out_imag = sumim;
				return;
			}
		}
	}
	// negative Im(z) is treated by calculating for -z , and using the
	// symmetry properties of erfc (z)
	const bool negimz = zim < 0.0;
	if (negimz) {
		zre = -zre;
		zim = -zim;
	}
	const double twosqrtpi = 3.54490770181103205e+00;
	const double tmzre = tm * zre, tmzim = tm * zim;
	// calculate exp (i tm z)
	double eitmzre = -tmzim, eitmzim = tmzre;
	cexp(eitmzre, eitmzim);
	// form 1 +/ - exp (i tm z)
	const double numerarr[4] = {
		1.0 - eitmzre , -eitmzim , 1.0 + eitmzre , +eitmzim
	};
	// form tm z * (1 +/ - exp (i tm z))
	const double numertmz[4] = {
		tmzre * numerarr[0] - tmzim * numerarr[1] ,
		tmzre * numerarr[1] + tmzim * numerarr[0] ,
		tmzre * numerarr[2] - tmzim * numerarr[3] ,
		tmzre * numerarr[3] + tmzim * numerarr[2]
	};
	// common subexpressions for use inside the loop
	const double reimtmzm2 = -2.0 * tmzre * tmzim;
	const double imtmz2 = tmzim * tmzim;
	const double reimtmzm22 = reimtmzm2 * reimtmzm2;
	const double znorm = zre * zre + zim2;
	double sumre = (-a24[0] / znorm) * (numerarr[0] * zre + numerarr[1] * zim);
	double sumim = (-a24[0] / znorm) * (numerarr[1] * zre - numerarr[0] * zim);
#pragma unroll
	for (unsigned i = 0; i < 24; ++i) {
		const unsigned j = (i << 1) & 2;
		// denominator
		const double wk = imtmz2 + (npi24[i] + tmzre) * (npi24[i] - tmzre);
		// norm of denominator
		const double norm = wk * wk + reimtmzm22;
		const double f = 2.0 * tm * a24[i] / norm;
		// sum += a[i] * numer / wk
		sumre -= f * (numertmz[j] * wk + numertmz[j + 1] * reimtmzm2);
		sumim -= f * (numertmz[j + 1] * wk - numertmz[j] * reimtmzm2);
	}
	if (negimz) {
		// use erfc (-z) = 2 - erfc (z) to get good accuracy for
		// Im (z) < 0: 2 / exp (z ^2) - w(z)
		const double z2im = 2.0 * zre * zim;
		const double z2re = (zre + zim) * (zre - zim);
		double ez2re = z2re, ez2im = z2im;
		cexp(ez2re, ez2im);
		const double twoez2norm = 2.0 / (ez2re * ez2re + ez2im * ez2im);
		*out_real = twoez2norm * ez2re + sumim / twosqrtpi;
		*out_imag = -twoez2norm * ez2im - sumre / twosqrtpi;
	}
	else {
		*out_real = -sumim / twosqrtpi;
		*out_imag = sumre / twosqrtpi;
	}

}


/* == function kernels == */
/*
Lorentzian->lorentz() A or S according to the amplitude/surface leading parameter

computed as:
where dx = x - x0, so g = gamma and g = HWHM, a- amplitude, s-surface

a g^2 / (dx^2 + g^2),
s/pi * g / (dx^2 + g^2)

*/

__global__ void LorentzS(const double *xi, double *yi, const double x0, const double s, const double w0, const double d0, const double xLimit, const double aLim, const double bLim)
{
	int myIndex = blockIdx.x * blockDim.x + threadIdx.x;
	if ((xLim == 0) || (abs(xi[myIndex] - (x0 + d0)) <= abs(xLim))) {

		yi[myIndex] += ((s * M_1_PI) * w0 / ((w0 * w0) + (xi[myIndex] - (x0 + d0))*(xi[myIndex] - (x0 + d0)))) - (aLim * xi[myIndex] + bLim);
	}
}

/*
Gaussian->gauss() A or S...
computed as

... wD = HWHM

*/

__global__ void GaussS(const double *xi, double *yi, const double x0, const double s, const double wD, const double d0, const double xLimit, const double aLim, const double bLim)
{
	int myIndex = blockIdx.x * blockDim.x + threadIdx.x;
	if ((xLim == 0) || (abs(xi[myIndex] - (x0 + d0)) <= abs(xLim))) {

		yi[myIndex] += (s * (M_SqrtLn2oPi / wD) * exp(-(xi[myIndex] - (x0 + d0))*(xi[myIndex] - (x0 + d0)) / (2 * M_1o2Ln2 * wD * wD))) - (aLim * xi[myIndex] + bLim);
	}
}

/*
Voigt peak

Vs(x; wD,w0) = SQRT(ln2/PI) * Re[w(z)]/wD
z = SQRT(ln2) * (x/wD) , SQRT(ln2) * (w0/wD)
*/

__global__ void VoigtS(const double *xi, double *yi, const double x0, const double s, const double wD, const double w0, const double d0, const double xLimit, const double aLim, const double bLim)
{
	int myIndex = blockIdx.x * blockDim.x + threadIdx.x;
	if ((xLim == 0) || (abs(xi[myIndex] - (x0 + d0)) <= abs(xLim))) {

		double rx = (M_SqrtLn2 * (xi[myIndex] - (x0 + d0)) / wD);
		double ry = ((w0 / wD) * M_SqrtLn2);
		double vr, vi;
		faddeeva(rx, ry, &vr, &vi);
		yi[myIndex] += s * vr * (1 / wD) * M_SqrtLn2oPi - (aLim * xi[myIndex] + bLim);
	}
}
/*
VoigtH peak

Vs(x; wD,w0) = SQRT(ln2/PI) * Re[w(z)]/wD
z = SQRT(ln2) * (x/wD) , SQRT(ln2) * (w0/wD)
...the w is calculated with the Humlicek algorithm
*/

__global__ void VoigtHS(const double *xi, double *yi, const double x0, const double s, const double wD, const double w0, const double d0, const double xLimit, const double aLim, const double bLim)
{
	int myIndex = blockIdx.x * blockDim.x + threadIdx.x;
	if ((xLim == 0) || (abs(xi[myIndex] - (x0 + d0)) <= abs(xLim))) {

		double z_real = M_SqrtLn2 * (xi[myIndex] - (x0 + d0)) / wD;
		double z_img = (w0 / wD) * M_SqrtLn2;
		double vrz, vri;

		faddeeva(z_real, z_img, &vrz, &vri);
		yi[myIndex] += s * vrz * (1 / wD) * M_SqrtLn2oPi - (aLim * xi[myIndex] + bLim);
	}
}

__global__ void HTP(const double *xi, double *yi, const double x0, const double d0, const double s, const double nuVC, const double eta, const double Ylm,
	const double x0c, const gpuCplx C0Tilde, const gpuCplx C0, const gpuCplx C2, const gpuCplx iCplx, const gpuCplx OneCplx,
	const gpuCplx invC2Tilde, const gpuCplx Yroot, const gpuCplx Y, const gpuCplx mOneCplx,
	const double vMP, const double xLimit, const double aLim, const double bLim)
{
	int myIndex = blockIdx.x * blockDim.x + threadIdx.x;
	if ((xLimit == 0) || (abs(xi[myIndex] - (x0 + d0)) <= abs(xLimit))) {

		gpuCplx ix = iCplx * (x0 - xi[myIndex]);
		gpuCplx X = (ix + C0Tilde) * invC2Tilde;

		gpuCplx Zm, Zp;

		Zm = sqrt(X + Y) - Yroot;
		Zp = Zm + (2.0 * Yroot);

		if (abs(Y) < (double)1.0e-15*abs(X)) {
			Zm = sqrt(X);
			Zp = sqrt(X + Y);
		}
		if (abs(X) < (double)3.0e-8*abs(Y)) {
			Zm = (ix + C0Tilde) * x0c;
			Zp = sqrt(X + Y) + Yroot;
		}

		double real, img;

		gpuCplx wZp = iCplx * Zm;
		faddeeva(wZp.real(), wZp.imag(), &real, &img);
		gpuCplx wZm(real, img);

		wZp = iCplx * Zp;
		faddeeva(wZp.real(), wZp.imag(), &real, &img);
		wZp.real(real), wZp.imag(img);

		gpuCplx Aterm = M_SqrtPI * x0c * (wZm - wZp);
		gpuCplx Bterm = invC2Tilde * vMP * vMP * (mOneCplx + M_SqrtPI*(OneCplx / (2.0*Yroot))*(OneCplx - Zm*Zm)*wZm - M_SqrtPI*(OneCplx / (2.0*Yroot))*(OneCplx - Zp*Zp)*wZp);

		if (abs(Y) < (double)1.0e-15*abs(X)) {
			Aterm = invC2Tilde*(OneCplx / X - (double)1.5 / (X*X));
			Bterm = invC2Tilde*(mOneCplx + (OneCplx - X - (double)2.0*Y)*((OneCplx / X) - (double)1.5*(1.0 / (X*X))) + (double)2.0*M_SqrtPI*wZp);

			if (abs(sqrt(X)) <= (double)4000.0) {
				Aterm = ((double)2.0*M_SqrtPI*invC2Tilde)*(M_1oSqrtPI - wZm);
				Bterm = invC2Tilde * vMP * vMP * (mOneCplx + (double)2.0*M_SqrtPI*(OneCplx - X - (double)2.0*Y)*(M_1oSqrtPI - wZm) + (double)2.0*M_SqrtPI*wZp);
			}
		}
		if (abs(X) < (double)3.0e-8*abs(Y)) {
			Aterm = M_SqrtPI * (wZm - wZp) * x0c;
			Bterm = invC2Tilde * vMP * vMP * (mOneCplx + M_SqrtPI / ((double)2.0*Yroot) * (OneCplx - Zm*Zm) * wZm - M_SqrtPI / ((double)2.0*Yroot) * (OneCplx - Zp*Zp) * wZp);
		}

		gpuCplx Hterm = Aterm * (nuVC - eta*(C0 - 1.5*C2)) + Bterm * (eta*C2 / (vMP*vMP));

		gpuCplx PeakComplex = (OneCplx + iCplx * Ylm) * Aterm * (OneCplx / (OneCplx - Hterm));
		yi[myIndex] += s * M_1_PI * PeakComplex.real() - (aLim * xi[myIndex] + bLim);
	}
}

__global__ void HTP_simple(const double *xi, double *yi, const double x0, const double d0, const double s, const double nuVC, const double eta, const double Ylm,
	const double x0c, const gpuCplx C0Tilde, const gpuCplx C0, const gpuCplx C2, const gpuCplx iCplx, const gpuCplx OneCplx,
	const double vMP, const double xLimit, const double aLim, const double bLim)
{
	int myIndex = blockIdx.x * blockDim.x + threadIdx.x;
	if ((xLimit == 0) || (abs(xi[myIndex] - (x0 + d0)) <= abs(xLimit))) {

		gpuCplx ix = iCplx * (x0 - xi[myIndex]);
		gpuCplx Zm = (ix + C0Tilde) * x0c;
		gpuCplx Bterm;
		double real, img;

		gpuCplx wZm = Zm * iCplx;
		faddeeva(wZm.real(), wZm.imag(), &real, &img);
		wZm.real(real), wZm.imag(img);

		gpuCplx Aterm = M_SqrtPI * x0c * wZm;
		if (abs(sqrt(Zm)) < (double)4000.0)
		{
			Bterm = M_SqrtPI * x0c * vMP * vMP *  (((OneCplx - (Zm*Zm)) * wZm) + Zm * M_1oSqrtPI);
		}
		else
		{
			Bterm = M_SqrtPI * x0c * vMP * vMP *  (wZm + ((double)0.5 / Zm) - ((double)0.75 *  M_1oSqrtPI * (Zm * Zm *Zm)));
		}

		Bterm = Aterm * (nuVC - eta*(C0 - 1.5*C2)) + Bterm * (eta*C2 / (vMP*vMP));

		wZm = (OneCplx + iCplx * Ylm) * Aterm * (OneCplx / (OneCplx - Bterm));
		yi[myIndex] += s * M_1_PI * wZm.real() - (aLim * xi[myIndex] + bLim);
	}
}

int PeakAdvCUDA(double *xs, long pointCount, double xLimit, double aLim, double bLim,
	double shape, double s, double x0, double wD, double w0, double d0,
	double w2, double d2, double nuVC, double eta, double Ylm, double* &gpu_xs, double* &gpu_ys)
{
	int ret = 0;

	//Parameters defining the constants used during computations.
	const int blockSize = 256;
	const int leftover = pointCount % blockSize;

	int groupSize = 1;
	if (leftover == 0) groupSize = 0;
	groupSize += pointCount / blockSize;

	if (gpu_xs == nullptr) {
		cudaMalloc(&gpu_xs, groupSize * blockSize * sizeof(double));
		cudaMalloc(&gpu_ys, groupSize * blockSize * sizeof(double));

		cudaMemcpy(gpu_xs, xs, pointCount * sizeof(double), cudaMemcpyHostToDevice);
		cudaMemset(gpu_ys, 0, pointCount * sizeof(double));

		ret = (int)cudaGetLastError();
	}

	cudaDeviceSynchronize();

	switch ((int)shape)
	{
	case 0:
		LorentzS << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, s, w0, d0, xLimit, aLim, bLim);
		break;

	case 1:
		GaussS << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, s, wD, d0, xLimit, aLim, bLim);
		break;

	case 2:
		VoigtS << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, s, wD, w0, d0, xLimit, aLim, bLim);
		break;

	case 3:
		VoigtHS << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, s, wD, w0, d0, xLimit, aLim, bLim);
		break;

	case 7: {
		const gpuCplx C0 = gpuCplx(w0, d0);
		const gpuCplx C2 = gpuCplx(w2, d2);

		const gpuCplx C0Tilde(((double)1.0 - eta) * (C0 - (double)1.5 * C2) + nuVC);
		const gpuCplx C2Tilde(C2 * ((double)1.0 - eta));

		const double vMP = P_CoSqrtLn2 * (wD / x0);
		const double x0c = (P_C / (vMP * x0));

		const gpuCplx iCplx((double)0.0, (double)1.0);
		const gpuCplx ZEROCplx((double)0.0, (double)0.0);
		const gpuCplx OneCplx((double)1.0, (double)0.0);
		const gpuCplx mOneCplx((double)-1.0, (double)0.0);

		if (C2Tilde != ZEROCplx) {
			const gpuCplx invC2Tilde(OneCplx / C2Tilde);
			const gpuCplx Yroot(x0 * vMP * invC2Tilde * P_1o2C);
			const gpuCplx Y(Yroot * Yroot);

			HTP << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, d0, s, nuVC, eta, Ylm, x0c, C0Tilde, C0, C2, iCplx, OneCplx,
				invC2Tilde, Yroot, Y, mOneCplx, vMP, xLimit, aLim, bLim);
		}
		else {
			if ((nuVC == 0) && (eta == 0) && (w2 == 0) && (d2 == 0))
				VoigtHS << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, s, wD, w0, d0, xLimit, aLim, bLim);
			else
				HTP_simple << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, d0, s, nuVC, eta, Ylm, x0c, C0Tilde, C0, C2, iCplx, OneCplx, vMP, xLimit, aLim, bLim);
		};
		break;
	}
	case 12:
		tester_wofz << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, d0, wD, w0);
		break;

	case 123:
		tester_fadd << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, d0, wD, w0);
		break;
	case 1234:
		tester << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, d0, wD, w0);
		break;

	default:
		break;
	}

	return ret;
}

__global__ void funcHTP(const double *xi, double *yi, const double* inVector, const double *preVector, const gpuCplx *cplxVector, const gpuCplx iCplx, const gpuCplx OneCplx,
						const gpuCplx mOneCplx, const long pointCount, const long peakCount)
{
	unsigned int yIndex = blockIdx.x * blockDim.x + threadIdx.x;
	unsigned int xIndex = fmod( (float)yIndex , (float) pointCount);
	unsigned int fIndex = yIndex / pointCount;
	const unsigned int inVectorArgs = 13, preVectorArgs = 2, cplxVectorArgs = 7;

	if (yIndex < pointCount * peakCount) {
		const double s = inVector[(fIndex * inVectorArgs)], x0 = inVector[(fIndex * inVectorArgs) + 1], d0 = inVector[(fIndex * inVectorArgs) + 4], nuVC = inVector[(fIndex * inVectorArgs) + 7], vMP = preVector[(fIndex * preVectorArgs)];
		const double eta = inVector[(fIndex * inVectorArgs) + 8], Ylm = inVector[(fIndex * inVectorArgs) + 9], xLimit = inVector[(fIndex * inVectorArgs) + 10], Alim = inVector[(fIndex * inVectorArgs) + 11], Blim = inVector[(fIndex * inVectorArgs) + 12], x0c = preVector[(fIndex * preVectorArgs) + 1];
		const gpuCplx C0 = cplxVector[(fIndex * cplxVectorArgs)], C2 = cplxVector[(fIndex * cplxVectorArgs) + 1], C0Tilde = cplxVector[(fIndex * cplxVectorArgs) + 2], C2Tilde = cplxVector[(fIndex * cplxVectorArgs) + 3], invC2Tilde = cplxVector[(fIndex * cplxVectorArgs) + 4];
		const gpuCplx Yroot = cplxVector[(fIndex * cplxVectorArgs) + 5], Y = cplxVector[(fIndex * cplxVectorArgs) + 6];

		if ((xLimit == 0) || (abs(xi[xIndex] - (x0 + d0)) <= abs(xLimit))) {

			gpuCplx ix = iCplx * (x0 - xi[xIndex]);
			gpuCplx X = (ix + C0Tilde) * invC2Tilde;

			gpuCplx Zm, Zp;

			Zm = sqrt(X + Y) - Yroot;
			Zp = Zm + (2.0 * Yroot);

			if (abs(Y) < (double)1.0e-15*abs(X)) {
				Zm = sqrt(X);
				Zp = sqrt(X + Y);
			}
			if (abs(X) < (double)3.0e-8*abs(Y)) {
				Zm = (ix + C0Tilde) * x0c;
				Zp = sqrt(X + Y) + Yroot;
			}

			double real, img;

			gpuCplx wZp = iCplx * Zm;
			faddeeva(wZp.real(), wZp.imag(), &real, &img);
			gpuCplx wZm(real, img);

			wZp = iCplx * Zp;
			faddeeva(wZp.real(), wZp.imag(), &real, &img);
			wZp.real(real), wZp.imag(img);

			gpuCplx Aterm = M_SqrtPI * x0c * (wZm - wZp);
			gpuCplx Bterm = invC2Tilde * vMP * vMP * (mOneCplx + M_SqrtPI * (OneCplx / (2.0*Yroot))*(OneCplx - Zm * Zm)*wZm - M_SqrtPI * (OneCplx / (2.0*Yroot))*(OneCplx - Zp * Zp)*wZp);

			if (abs(Y) < (double)1.0e-15*abs(X)) {
				Aterm = invC2Tilde * (OneCplx / X - (double)1.5 / (X*X));
				Bterm = invC2Tilde * (mOneCplx + (OneCplx - X - (double)2.0*Y)*((OneCplx / X) - (double)1.5*(1.0 / (X*X))) + (double)2.0*M_SqrtPI*wZp);

				if (abs(sqrt(X)) <= (double)4000.0) {
					Aterm = ((double)2.0*M_SqrtPI*invC2Tilde)*(M_1oSqrtPI - wZm);
					Bterm = invC2Tilde * vMP * vMP * (mOneCplx + (double)2.0*M_SqrtPI*(OneCplx - X - (double)2.0*Y)*(M_1oSqrtPI - wZm) + (double)2.0*M_SqrtPI*wZp);
				}
			}
			if (abs(X) < (double)3.0e-8*abs(Y)) {
				Aterm = M_SqrtPI * (wZm - wZp) * x0c;
				Bterm = invC2Tilde * vMP * vMP * (mOneCplx + M_SqrtPI / ((double)2.0*Yroot) * (OneCplx - Zm * Zm) * wZm - M_SqrtPI / ((double)2.0*Yroot) * (OneCplx - Zp * Zp) * wZp);
			}

			gpuCplx Hterm = Aterm * (nuVC - eta * (C0 - 1.5*C2)) + Bterm * (eta*C2 / (vMP*vMP));

			gpuCplx PeakComplex = (OneCplx + iCplx * Ylm) * Aterm * (OneCplx / (OneCplx - Hterm));
			yi[(xIndex * peakCount) + fIndex] = s * M_1_PI * PeakComplex.real() - (Alim * xi[xIndex] + Blim);
		}
	}
}

int FuncAdvCUDA(double *xs, double *ys, long pointCount, size_t peakCount, double shape, double *inputVector, double* &gpu_xs, double* &gpu_ys)
{
	int ret = 0;

	//Parameters defining the constants used during computations.
	const int blockSize = 256;
	const int leftover = (pointCount*peakCount) % blockSize;
	double *inVector;

	int groupSize = 1;
	if (leftover == 0) groupSize = 0;
	groupSize += (pointCount*peakCount) / blockSize;

	cudaMalloc(&gpu_xs, groupSize * blockSize * sizeof(double));
	cudaMalloc(&gpu_ys, groupSize * blockSize * sizeof(double));
	cudaMalloc(&inVector, 13 * peakCount * sizeof(double));

	cudaMemcpy(inVector, inputVector, 13 * peakCount * sizeof(double), cudaMemcpyHostToDevice);
	cudaMemcpy(gpu_xs, xs, pointCount * sizeof(double), cudaMemcpyHostToDevice);

	ret = (int)cudaGetLastError();

	switch ((int)shape)
	{
	case 0:
		//LorentzS << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, s, w0, d0, xLimit, aLim, bLim);
		break;

	case 1:
		//GaussS << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, s, wD, d0, xLimit, aLim, bLim);
		break;

	case 2:
		//VoigtS << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, s, wD, w0, d0, xLimit, aLim, bLim);
		break;

	case 3:
		//VoigtHS << <groupSize, blockSize >> > (gpu_xs, gpu_ys, x0, s, wD, w0, d0, xLimit, aLim, bLim);
		break;

	case 7: {
		const gpuCplx iCplx((double)0.0, (double)1.0);
		const gpuCplx OneCplx((double)1.0, (double)0.0);
		const gpuCplx mOneCplx((double)-1.0, (double)0.0);
		
		double *preVector, *precomputed;
		gpuCplx *cplxVector, *precomputedCplx;

		cudaMallocHost(&precomputed, 2 * peakCount * sizeof(double));
		cudaMallocHost(&precomputedCplx, 7 * peakCount * sizeof(gpuCplx));

		for (int i = 0; i < peakCount; i++) {

			const double x0 = inputVector[(i*13) + 1], w0 = inputVector[(i * 13) + 3], d0 = inputVector[(i * 13) + 4], nuVC = inputVector[(i * 13) + 7];
			const double eta = inputVector[(i * 13) + 8], wD = inputVector[(i * 13) + 2], w2 = inputVector[(i * 13) + 5], d2 = inputVector[(i * 13) + 6];

			precomputedCplx[i * 7] = gpuCplx(w0, d0);
			precomputedCplx[(i * 7) + 1] = gpuCplx(w2, d2);

			precomputedCplx[(i * 7) + 2] = gpuCplx(((double)1.0 - eta) * (precomputedCplx[i * 7] - (double)1.5 * precomputedCplx[(i * 7) + 1]) + nuVC);
			precomputedCplx[(i * 7) + 3] = gpuCplx(precomputedCplx[(i * 7) + 1] * ((double)1.0 - eta));

			precomputed[i * 2] = P_CoSqrtLn2 * (wD / x0);
			precomputed[(i * 2) + 1] = (P_C / (precomputed[i * 2] * x0));

			precomputedCplx[(i * 7) + 4] = gpuCplx(OneCplx / precomputedCplx[(i * 7) + 3]);
			precomputedCplx[(i * 7) + 5] = gpuCplx(x0 * precomputed[i * 2] * precomputedCplx[(i * 7) + 4] * P_1o2C);
			precomputedCplx[(i * 7) + 6] = gpuCplx(precomputedCplx[(i * 7) + 5] * precomputedCplx[(i * 7) + 5]);
		}

		cudaMalloc(&preVector, 2 * peakCount * sizeof(double));
		cudaMalloc(&cplxVector, 7 * peakCount * sizeof(gpuCplx));

		cudaMemcpy(preVector, precomputed, 2 * peakCount * sizeof(double), cudaMemcpyHostToDevice);
		cudaMemcpy(cplxVector, precomputedCplx, 7 * peakCount * sizeof(gpuCplx), cudaMemcpyHostToDevice);

		funcHTP << <groupSize, blockSize >> > (gpu_xs, gpu_ys, inVector, preVector, cplxVector, iCplx, OneCplx, mOneCplx, pointCount, peakCount);
		
		
		cudaFree(preVector);
		cudaFree(cplxVector);
		cudaFreeHost(precomputed);
		cudaFreeHost(precomputedCplx);
		}
		break;

	default:
		break;
	}
	
	//Final summation of calculated values.
	double *sumVector;
	cudaDeviceSynchronize();

	cudaMallocHost(&sumVector, pointCount * peakCount * sizeof(double));
	cudaMemcpy(sumVector, gpu_ys, pointCount * peakCount * sizeof(double), cudaMemcpyDeviceToHost);
	for (int i = 0; i < pointCount; i++) {
		double sum = 0;
		for (int j = 0; j < peakCount; j++) {
			sum += sumVector[(i*peakCount) + j];
		}
		
		ys[i] = sum;
	}

	cudaFreeHost(sumVector);
	cudaFreeHost(inputVector);
	cudaFree(inVector);
	cudaFree(gpu_xs);
	cudaFree(gpu_ys);
	return ret;
}

void cleanUp(double *ys, long pointCount, double *gpu_xs, double *gpu_ys)
{
	cudaDeviceSynchronize(); 
	cudaMemcpy(ys, gpu_ys, pointCount * sizeof(double), cudaMemcpyDeviceToHost);

	cudaFree(gpu_xs);
	cudaFree(gpu_ys);
}