/*
    PIP - Platform Independent Primitives
    Math
    Copyright (C) 2013  Ivan Pelipenko peri4ko@gmail.com

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "pimath.h"


const char Solver::methods_desc[] = "b{Methods:}\
\n  -1 - Global settings\
\n  01 - Eyler 1\
\n  02 - Eyler 2\
\n  14 - Runge-Kutta 4\
\n  23 - Adams-Bashfort-Moulton 3\
\n  24 - Adams-Bashfort-Moulton 4\
\n  32 - Polynomial Approximation 2\
\n  33 - Polynomial Approximation 3\
\n  34 - Polynomial Approximation 4\
\n  35 - Polynomial Approximation 5";

Solver::Method Solver::method_global = Solver::Eyler_2;


void Solver::solve(double u, double h) {
	switch (method) {
	case Global:
		switch (method_global) {
		case Eyler_1: solveEyler1(u, h); break;
		case Eyler_2: solveEyler2(u, h); break;
		case RungeKutta_4: solveRK4(u, h); break;
		case AdamsBashfortMoulton_2: solveABM2(u, h); break;
		case AdamsBashfortMoulton_3: solveABM3(u, h); break;
		case AdamsBashfortMoulton_4: default: solveABM4(u, h); break;
		case PolynomialApproximation_2: solvePA2(u, h); break;
		case PolynomialApproximation_3: solvePA3(u, h); break;
		case PolynomialApproximation_4: solvePA4(u, h); break;
		case PolynomialApproximation_5: solvePA5(u, h); break;
		}
		break;
	case Eyler_1: solveEyler1(u, h); break;
	case Eyler_2: solveEyler2(u, h); break;
	case RungeKutta_4: solveRK4(u, h); break;
	case AdamsBashfortMoulton_2: solveABM2(u, h); break;
	case AdamsBashfortMoulton_3: solveABM3(u, h); break;
	case AdamsBashfortMoulton_4: default: solveABM4(u, h); break;
	case PolynomialApproximation_2: solvePA2(u, h); break;
	case PolynomialApproximation_3: solvePA3(u, h); break;
	case PolynomialApproximation_4: solvePA4(u, h); break;
	case PolynomialApproximation_5: solvePA5(u, h); break;
	}
	step++;
}


void Solver::fromTF(const TransferFunction & TF) {
	if (TF.vector_An.size() >= TF.vector_Bm.size())
		size = TF.vector_An.size()-1;
	else {
		cout << "Solver error: {A} should be greater than {B}" << endl;
		return;
	}
	if (size == 0) return;

	step = 0;
	times.fill(0.);
	A.resize(size, size);
	d.resize(size + 1); d.fill(0.);
	a1.resize(size + 1); a1.fill(0.);
	b1.resize(size + 1); b1.fill(0.);
	X.resize(size); X.fill(0.);
	F.resize(5);
	for (uint i = 0; i < 5; ++i)
		F[i].resize(size), F[i].fill(0.);
	k1.resize(size); k1.fill(0.);
	k2.resize(size); k2.fill(0.);
	k3.resize(size); k3.fill(0.);
	k4.resize(size); k4.fill(0.);
	xx.resize(size); xx.fill(0.);
	XX.resize(size); XX.fill(0.);
	for (uint i = 0; i < size + 1; ++i)
		a1[size - i] = TF.vector_An[i];
	for (uint i = 0; i < TF.vector_Bm.size(); ++i)
		b1[size - i] = TF.vector_Bm[i];
	double a0 = a1[0];
	a1 /= a0;
	b1 /= a0;

	d[0] = b1[0]; // Рассчитываем вектор d
	for (uint i = 1; i < size + 1; ++i) {
		sum = 0.;
		for (uint m = 0; m < i; ++m)
			sum += a1[i - m] * d[m];
		d[i] = b1[i] - sum;
	}

	for (uint i = 0; i < size - 1; ++i) // Заполняем матрицу А
		for (uint j = 0; j < size; ++j)
				A[j][i] = (j == i + 1);
	for (uint i = 0; i < size; ++i)
		A[i][size - 1] = -a1[size - i];
	for (uint i = 0; i < size; ++i)
		d[i] = d[i + 1];
}


void Solver::solveEyler1(double u, double h) {
	/*for (uint i = 0; i < size; ++i) {
	 *		sum = 0.;
	 *		for (uint j = 0; j < size; ++j)
	 *			sum += A[j][i] * X[j];
	 *		xx[i] = sum + d[i] * u;
	}*/
	F[0] = A * X + d * u;
	X += F[0] * h;
	moveF();
}


void Solver::solveEyler2(double u, double h) {
	F[0] = A * X + d * u;
	X += (F[0] + F[1]) * h / 2.;
	moveF();
}


void Solver::solveRK4(double u, double h) {
	td = X[0] - F[0][0];
	k1 = A * X + d * u; xx = k1 * h / 2.; XX = X + xx;
	k2 = A * (XX + k1 * h / 2.) + d * (u + td/3.); xx = k2 * h / 2.; XX += xx;
	k3 = A * (XX + k2 * h / 2.) + d * (u + td*2./3.); xx = k3 * h; XX += xx;
	k4 = A * (XX + k3 * h) + d * (u + td);
	//cout << k1 << k2 << k3 << k4 << endl;
	X += (k1 + k2 * 2. + k3 * 2. + k4) * h / 6.;
	F[0] = X;
}


void Solver::solveABM2(double u, double h) {
	F[0] = A * X + d * u;
	XX = X + (F[0] * 3. - F[1]) * (h / 2.);
	F[1] = A * XX + d * u;
	X += (F[1] + F[0]) * (h / 2.);
	moveF();
}


void Solver::solveABM3(double u, double h) {
	F[0] = A * X + d * u;
	XX = X + (F[0] * 23. - F[1] * 16. + F[2] * 5.) * (h / 12.);
	F[2] = A * XX + d * u;
	X += (F[2] * 5. + F[0] * 8. - F[1]) * (h / 12.);
	moveF();
}


void Solver::solveABM4(double u, double h) {
	F[0] = A * X + d * u;
	XX = X + (F[0] * 55. - F[1] * 59. + F[2] * 37. - F[3] * 9.) * (h / 24.);
	F[3] = A * XX + d * u;
	X += (F[3] * 9. + F[0] * 19. - F[1] * 5. + F[2]) * (h / 24.);
	moveF();
}


void Solver::solvePA(double u, double h, uint deg) {
	F[0] = A * X + d * u;
	M.resize(deg, deg);
	Y.resize(deg);
	pY.resize(deg);

	for (uint k = 0; k < size; ++k) {
		for (uint i = 0; i < deg; ++i) {
			td = 1.;
			ct = times[i];
			for (uint j = 0; j < deg; ++j) {
				M[j][i] = td;
				td *= ct;
			}
		}
		for (uint i = 0; i < deg; ++i)
			Y[i] = F[i][k];
		/// find polynom
		//if (step == 1) cout << M << endl << Y << endl;
		M.invert(&ok, &Y);
		//if (step == 1) cout << Y << endl;
		if (!ok) {
			solveEyler2(u, h);
			break;
		}
		/// calc last piece
		x0 = 0.;
		td = 1.;
		t = times[0];
		for (uint i = 0; i < deg; ++i) {
			x0 += Y[i] * td / (i + 1.);
			td *= t;
		}
		x0 *= t;

		x1 = 0.;
		td = 1.;
		t = times[1];
		for (uint i = 0; i < deg; ++i) {
			x1 += Y[i] * td / (i + 1.);
			td *= t;
		}
		x1 *= t;
		lp = x0 - x1;

		if (deg > 2) {
			/// calc prev piece
			x0 = 0.;
			td = 1.;
			dh = times[1] - times[2];
			if (dh != 0. && step > 1) {
				t = times[2];
				for (uint i = 0; i < deg; ++i) {
					x0 += Y[i] * td / (i + 1.);
					td *= t;
				}
				x0 *= t;
				ct = x1 - x0;
				/// calc correction
				ct -= pY[k];
			}
			/// calc final
			X[k] += lp - ct;
			pY[k] = lp;
		} else {
			X[k] += lp;
			pY[k] = lp;
		}
	}
	moveF();
}



PIFFT::PIFFT() {
	prepared = false;
}


PIVector<complexd> * PIFFT::calcFFT(const PIVector<complexd> & val) {
//    for (uint i=0; i<result.size(); i+=2)
//    {
//        result[i] = val.at(indexes[i]) + val.at(indexes[i+1]);
//        result[i+1] = val.at(indexes[i]) - val.at(indexes[i+1]);
//    }
//    return &result;
    result.clear();
	if (val.size_s() < 4) return &result;
	fftc1d(val, val.size());
    return &result;
}


PIVector<complexd> *PIFFT::calcFFTinverse(const PIVector<complexd> &val)
{
    result.clear();
    if (val.size_s() < 4) return &result;
    fftc1dinv(val, val.size());
    return &result;
}


PIVector<complexd> *PIFFT::calcHilbert(const PIVector<double> &val)
{
    result.clear();
    if (val.size_s() < 4) return &result;
    fftc1r(val, val.size());
    for (uint i=0; i<result.size()/2; i++)           result[i] = result[i]*2.;
    for (uint i=result.size()/2; i<result.size(); i++) result[i] = 0;
    fftc1dinv(result, result.size());
    return &result;
}


PIVector< complexd >* PIFFT::calcFFT(const PIVector<double> & val) {
    result.clear();
    if (val.size_s() < 4) return &result;
	fftc1r(val, val.size());
	return &result;
}


PIVector<double> PIFFT::getAmplitude() {
	PIVector<double> a;
	double tmp;
	for (uint i=0; i<result.size(); i++) {
		tmp = sqrt(result.at(i).real()*result.at(i).real()+result.at(i).imag()*result.at(i).imag());
		a.push_back(tmp);
	}
	return a;
}


void PIFFT::fftc1d(const PIVector<complexd> &a, uint n) {
	createPlan(n);
	uint i;
	PIVector<double> buf;
	buf.resize(2*n);
	for(i=0; i<n; i++) {
		buf[2*i+0] = a.at(i).real();// a->ptr.p_complex[i].x;
		buf[2*i+1] = a.at(i).imag();//a->ptr.p_complex[i].y;
	}
	ftbaseexecuteplan(&buf, 0, n, &curplan);
	result.resize(n);
	for(i=0; i<n; i++)
		result[i]=complexd(buf[2*i+0],buf[2*i+1]);
}


void PIFFT::fftc1r(const PIVector<double> & a, uint n) {
	uint i;
	if( n%2==0 ) {
		PIVector<double> buf;
		uint n2 = n/2;
		//buf.resize(n);
		buf = a;
		createPlan(n2);
		//cout << "fftr " << n2 << endl;
		ftbaseexecuteplan(&buf, 0, n2, &curplan);
		result.resize(n);
		uint idx;
		complexd hn, hmnc, v;
		for(i=0; i<=n2; i++) {
			idx = 2*(i%n2);
			hn = complexd(buf[idx+0], buf[idx+1]);
			idx = 2*((n2-i)%n2);
			hmnc = complexd(buf[idx+0], -buf[idx+1]);
			v = complexd(sin(M_PI*i/n2), cos(M_PI*i/n2));
			result[i] = ((hn + hmnc) - (v * (hn - hmnc)));
			result[i] *= 0.5;
		}
		for(i=n2+1; i<n; i++)
			result[i] = conj(result[n-i]);
	} else {
		PIVector<complexd> cbuf;
		cbuf.resize(n);
		for(i=0; i<n; i++)
			cbuf[i] = complexd(a[i], 0.);
		fftc1d(cbuf, n);
    }
}


void PIFFT::fftc1dinv(const PIVector<complexd> &a, uint n)
{
    PIVector<complexd> cbuf;
    cbuf.resize(n);
    uint i;
    for(i=0; i<n; i++)
    {
        cbuf[i] = conj(a[i]);
    }
    fftc1d(cbuf, n);
//    result.resize(n);
    for(i=0; i<n; i++)
    {
        result[i] = conj(result[i] / (double)n);
    }
}


void PIFFT::createPlan(uint n) {
	curplan.plan.clear();
	curplan.precomputed.clear();
	curplan.stackbuf.clear();
	curplan.tmpbuf.clear();
	if (n<2) return;
	ftbasegeneratecomplexfftplan(n, &curplan);
	prepared = true;
}


void PIFFT::ftbasegeneratecomplexfftplan(uint n, ftplan* plan) {
	int planarraysize;
	int plansize;
	int precomputedsize;
	int tmpmemsize;
	int stackmemsize;
	ae_int_t stackptr;
	planarraysize = 1;
	plansize = 0;
	precomputedsize = 0;
	stackmemsize = 0;
	stackptr = 0;
	tmpmemsize = 2*n;
	curplan.plan.resize(planarraysize);
	int ftbase_ftbasecffttask = 0;
	ftbase_ftbasegenerateplanrec(n, ftbase_ftbasecffttask, plan, &plansize, &precomputedsize, &planarraysize, &tmpmemsize, &stackmemsize, stackptr);
	if (stackptr!=0) { return;}//ae_assert(stackptr==0, "Internal error in FTBaseGenerateComplexFFTPlan: stack ptr!");
    curplan.stackbuf.resize(piMax(stackmemsize,1));//ae_vector_set_length(&curplan.stackbuf, ae_maxint(stackmemsize, 1));
    curplan.tmpbuf.resize(piMax(tmpmemsize,1));//ae_vector_set_length(&(curplan.tmpbuf), ae_maxint(tmpmemsize, 1));
    curplan.precomputed.resize(piMax(precomputedsize,1));//ae_vector_set_length(&curplan.precomputed, ae_maxint(precomputedsize, 1));
	stackptr = 0;
	ftbase_ftbaseprecomputeplanrec(plan, 0, stackptr);
	if (stackptr!=0) { return;}//ae_assert(stackptr==0, "Internal error in FTBaseGenerateComplexFFTPlan: stack ptr!");
}


/*************************************************************************
Recurrent subroutine for the FFTGeneratePlan:

PARAMETERS:
	N                   plan size
	IsReal              whether input is real or not.
						subroutine MUST NOT ignore this flag because real
						inputs comes with non-initialized imaginary parts,
						so ignoring this flag will result in corrupted output
	HalfOut             whether full output or only half of it from 0 to
						floor(N/2) is needed. This flag may be ignored if
						doing so will simplify calculations
	Plan                plan array
	PlanSize            size of used part (in integers)
	PrecomputedSize     size of precomputed array allocated yet
	PlanArraySize       plan array size (actual)
	TmpMemSize          temporary memory required size
	BluesteinMemSize    temporary memory required size

-- ALGLIB --
	Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT::ftbase_ftbasegenerateplanrec(
		int n,
		int tasktype,
		ftplan* plan,
		int* plansize,
		int* precomputedsize,
		int* planarraysize,
		int* tmpmemsize,
		int* stackmemsize,
		ae_int_t stackptr, int debugi)
{
	int k, m, n1, n2, esize, entryoffset;
	int ftbase_ftbaseplanentrysize = 8;
	int ftbase_ftbasecffttask = 0;
	int ftbase_fftcooleytukeyplan = 0;
	int ftbase_fftbluesteinplan = 1;
	int ftbase_fftcodeletplan = 2;
	int ftbase_fftrealcooleytukeyplan = 5;
	int ftbase_fftemptyplan = 6;
	if( *plansize+ftbase_ftbaseplanentrysize>(*planarraysize) ) {
		curplan.plan.resize(8*(*planarraysize));
		*planarraysize = 8*(*planarraysize);
	}
	entryoffset = *plansize;
	esize = ftbase_ftbaseplanentrysize;
	*plansize = *plansize+esize;
	if( n==1 ) {
		curplan.plan[entryoffset+0] = esize;
		curplan.plan[entryoffset+1] = -1;
		curplan.plan[entryoffset+2] = -1;
		curplan.plan[entryoffset+3] = ftbase_fftemptyplan;
		curplan.plan[entryoffset+4] = -1;
		curplan.plan[entryoffset+5] = -1;
		curplan.plan[entryoffset+6] = -1;
		curplan.plan[entryoffset+7] = -1;
		return;
	}
	ftbasefactorize(n, &n1, &n2);
	if( n1!=1 ) {
        *tmpmemsize = piMax(*tmpmemsize, 2*n1*n2);
		curplan.plan[entryoffset+0] = esize;
		curplan.plan[entryoffset+1] = n1;
		curplan.plan[entryoffset+2] = n2;
		if( tasktype==ftbase_ftbasecffttask )
			curplan.plan[entryoffset+3] = ftbase_fftcooleytukeyplan;
		else
			curplan.plan[entryoffset+3] = ftbase_fftrealcooleytukeyplan;
		curplan.plan[entryoffset+4] = 0;
		curplan.plan[entryoffset+5] = *plansize;
		debugi++;
		ftbase_ftbasegenerateplanrec(n1, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr,debugi);
		curplan.plan[entryoffset+6] = *plansize;
		ftbase_ftbasegenerateplanrec(n2, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr,debugi);
		curplan.plan[entryoffset+7] = -1;
		return;
	} else {
		if (n>=2 && n<=5) {
			curplan.plan[entryoffset+0] = esize;
			curplan.plan[entryoffset+1] = n1;
			curplan.plan[entryoffset+2] = n2;
			curplan.plan[entryoffset+3] = ftbase_fftcodeletplan;
			curplan.plan[entryoffset+4] = 0;
			curplan.plan[entryoffset+5] = -1;
			curplan.plan[entryoffset+6] = -1;
			curplan.plan[entryoffset+7] = *precomputedsize;
			if( n==3 )
				*precomputedsize = *precomputedsize+2;
			if( n==5 )
				*precomputedsize = *precomputedsize+5;
			return;
		} else {
			k = 2*n2-1;
			m = ftbasefindsmooth(k);
            *tmpmemsize = piMax(*tmpmemsize, 2*m);
			curplan.plan[entryoffset+0] = esize;
			curplan.plan[entryoffset+1] = n2;
			curplan.plan[entryoffset+2] = -1;
			curplan.plan[entryoffset+3] = ftbase_fftbluesteinplan;
			curplan.plan[entryoffset+4] = m;
			curplan.plan[entryoffset+5] = *plansize;
			stackptr = stackptr+2*2*m;
            *stackmemsize = piMax(*stackmemsize, stackptr);
			ftbase_ftbasegenerateplanrec(m, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr);
			stackptr = stackptr-2*2*m;
			curplan.plan[entryoffset+6] = -1;
			curplan.plan[entryoffset+7] = *precomputedsize;
			*precomputedsize = *precomputedsize+2*m+2*n;
			return;
		}
	}
}


/*************************************************************************
Recurrent subroutine for precomputing FFT plans

-- ALGLIB --
	Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT::ftbase_ftbaseprecomputeplanrec(ftplan* plan,
										int entryoffset,
										ae_int_t stackptr)
{
	int n1, n2, n, m, offs;
	double v, bx, by;
	int ftbase_fftcooleytukeyplan = 0;
	int ftbase_fftbluesteinplan = 1;
	int ftbase_fftcodeletplan = 2;
	int ftbase_fhtcooleytukeyplan = 3;
	int ftbase_fhtcodeletplan = 4;
	int ftbase_fftrealcooleytukeyplan = 5;
	if( (curplan.plan[entryoffset+3]==ftbase_fftcooleytukeyplan||curplan.plan[entryoffset+3]==ftbase_fftrealcooleytukeyplan)||curplan.plan[entryoffset+3]==ftbase_fhtcooleytukeyplan ) {
		ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset+5], stackptr);
		ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset+6], stackptr);
		return;
	}
	if( curplan.plan[entryoffset+3]==ftbase_fftcodeletplan||curplan.plan[entryoffset+3]==ftbase_fhtcodeletplan ) {
		n1 = curplan.plan[entryoffset+1];
		n2 = curplan.plan[entryoffset+2];
		n = n1*n2;
		if( n==3 ) {
			offs = curplan.plan[entryoffset+7];
			curplan.precomputed[offs+0] = cos(2*M_PI/3)-1;
			curplan.precomputed[offs+1] = sin(2*M_PI/3);
			return;
		}
		if( n==5 ) {
			offs = curplan.plan[entryoffset+7];
			v = 2*M_PI/5;
			curplan.precomputed[offs+0] = (cos(v)+cos(2*v))/2-1;
			curplan.precomputed[offs+1] = (cos(v)-cos(2*v))/2;
			curplan.precomputed[offs+2] = -sin(v);
			curplan.precomputed[offs+3] = -(sin(v)+sin(2*v));
			curplan.precomputed[offs+4] = sin(v)-sin(2*v);
			return;
		}
	}
	if( curplan.plan[entryoffset+3]==ftbase_fftbluesteinplan ) {
		ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset+5], stackptr);
		n = curplan.plan[entryoffset+1];
		m = curplan.plan[entryoffset+4];
		offs = curplan.plan[entryoffset+7];
		for(int i=0; i<=2*m-1; i++)
			curplan.precomputed[offs+i] = 0;
		for(int i=0; i<n; i++) {
			bx = cos(M_PI*sqr(i)/n);
			by = sin(M_PI*sqr(i)/n);
			curplan.precomputed[offs+2*i+0] = bx;
			curplan.precomputed[offs+2*i+1] = by;
			curplan.precomputed[offs+2*m+2*i+0] = bx;
			curplan.precomputed[offs+2*m+2*i+1] = by;
			if( i>0 ) {
				curplan.precomputed[offs+2*(m-i)+0] = bx;
				curplan.precomputed[offs+2*(m-i)+1] = by;
			}
		}
		ftbaseexecuteplanrec(&curplan.precomputed, offs, plan, curplan.plan[entryoffset+5], stackptr);
		return;
	}
}


void PIFFT::ftbasefactorize(int n, int* n1, int* n2) {
	*n1 = *n2 = 0;
	int ftbase_ftbasecodeletrecommended = 5;
	if( (*n1)*(*n2)!=n ) {
		for(int j=ftbase_ftbasecodeletrecommended; j>=2; j--) {
			if( n%j==0 ) {
				*n1 = j;
				*n2 = n/j;
				break;
			}
		}
	}
	if( (*n1)*(*n2)!=n ) {
		for(int j=ftbase_ftbasecodeletrecommended+1; j<=n-1; j++) {
			if( n%j==0 ) {
				*n1 = j;
				*n2 = n/j;
				break;
			}
		}
	}
	if( (*n1)*(*n2)!=n ) {
		*n1 = 1;
		*n2 = n;
	}
	if( (*n2)==1 && (*n1)!=1 ) {
		*n2 = *n1;
		*n1 = 1;
	}
}


/*************************************************************************
Is number smooth?

-- ALGLIB --
	Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT::ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int* best) {
	if( seed>=n ) {
        *best = piMin(*best, seed);
		return;
	}
	if( leastfactor<=2 )
		ftbase_ftbasefindsmoothrec(n, seed*2, 2, best);
	if( leastfactor<=3 )
		ftbase_ftbasefindsmoothrec(n, seed*3, 3, best);
	if( leastfactor<=5 )
		ftbase_ftbasefindsmoothrec(n, seed*5, 5, best);
}


int PIFFT::ftbasefindsmooth(int n) {
	int best, result;
	best = 2;
	while(best<n)
		best = 2*best;
	ftbase_ftbasefindsmoothrec(n, 1, 2, &best);
	result = best;
	return result;
}


void PIFFT::ftbase_internalreallintranspose(PIVector<double>* a, int m, int n, int astart, PIVector<double>* buf) {
	ftbase_fftirltrec(a, astart, n, buf, 0, m, m, n);
	for (int i=0; i<2*m*n; i++) (*a)[astart+i] = (*buf)[i];
}


void PIFFT::ftbase_fftirltrec(PIVector<double>* a, int astart, int astride, PIVector<double>* b, int bstart, int bstride, int m, int n) {
	int idx1, idx2;
	int m1, n1;
	if( m==0||n==0 )
		return;
    if( piMax(m, n)<=8 ) {
		for(int i=0; i<=m-1; i++) {
			idx1 = bstart+i;
			idx2 = astart+i*astride;
			for(int j=0; j<=n-1; j++) {
				(*b)[idx1] = a->at(idx2);
				idx1 = idx1+bstride;
				idx2 = idx2+1;
			}
		}
		return;
	}
	if( n>m ) {
		n1 = n/2;
		if( n-n1>=8&&n1%8!=0 )
			n1 = n1+(8-n1%8);
		ftbase_fftirltrec(a, astart, astride, b, bstart, bstride, m, n1);
		ftbase_fftirltrec(a, astart+n1, astride, b, bstart+n1*bstride, bstride, m, n-n1);
	} else {
		m1 = m/2;
		if( m-m1>=8&&m1%8!=0 )
			m1 = m1+(8-m1%8);
		ftbase_fftirltrec(a, astart, astride, b, bstart, bstride, m1, n);
		ftbase_fftirltrec(a, astart+m1*astride, astride, b, bstart+m1, bstride, m-m1, n);
	}
}


void PIFFT::ftbase_internalcomplexlintranspose(PIVector<double>* a, int m, int n, int astart, PIVector<double>* buf) {
	ftbase_ffticltrec(a, astart, n, buf, 0, m, m, n);
	for (int i=0; i<2*m*n; i++)
		(*a)[astart+i] = (*buf)[i];
}


void PIFFT::ftbase_ffticltrec(PIVector<double>* a, int astart, int astride, PIVector<double>* b, int bstart, int bstride, int m, int n) {
	int idx1, idx2, m2, m1, n1;
	if( m==0||n==0 )
		return;
    if( piMax(m, n)<=8 ) {
		m2 = 2*bstride;
		for(int i=0; i<=m-1; i++) {
			idx1 = bstart+2*i;
			idx2 = astart+2*i*astride;
			for(int j=0; j<=n-1; j++) {
				(*b)[idx1+0] = a->at(idx2+0);
				(*b)[idx1+1] = a->at(idx2+1);
				idx1 = idx1+m2;
				idx2 = idx2+2;
			}
		}
		return;
	}
	if( n>m ) {
		n1 = n/2;
		if( n-n1>=8&&n1%8!=0 )
			n1 = n1+(8-n1%8);
		ftbase_ffticltrec(a, astart, astride, b, bstart, bstride, m, n1);
		ftbase_ffticltrec(a, astart+2*n1, astride, b, bstart+2*n1*bstride, bstride, m, n-n1);
	} else {
		m1 = m/2;
		if( m-m1>=8&&m1%8!=0 )
			m1 = m1+(8-m1%8);
		ftbase_ffticltrec(a, astart, astride, b, bstart, bstride, m1, n);
		ftbase_ffticltrec(a, astart+2*m1*astride, astride, b, bstart+2*m1, bstride, m-m1, n);
	}
}


void PIFFT::ftbaseexecuteplan(PIVector<double>* a, int aoffset, int n, ftplan* plan) {
	ae_int_t stackptr;
	stackptr = 0;
	ftbaseexecuteplanrec(a, aoffset, plan, 0, stackptr);
}


/*************************************************************************
Recurrent subroutine for the FTBaseExecutePlan

Parameters:
	A           FFT'ed array
	AOffset     offset of the FFT'ed part (distance is measured in doubles)

-- ALGLIB --
	Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan, int entryoffset, ae_int_t stackptr) {
	int n1, n2, n, m, offs, offs1, offs2, offsa, offsb, offsp;
	double hk, hnk, x, y, bx, by, v0, v1, v2, v3;
	double a0x, a0y, a1x, a1y, a2x, a2y, a3x, a3y;
	double t1x, t1y, t2x, t2y, t3x, t3y, t4x, t4y, t5x, t5y;
	double m1x, m1y, m2x, m2y, m3x, m3y, m4x, m4y, m5x, m5y;
	double s1x, s1y, s2x, s2y, s3x, s3y, s4x, s4y, s5x, s5y;
	double c1, c2, c3, c4, c5;
	int ftbase_fftcooleytukeyplan = 0;
	int ftbase_fftbluesteinplan = 1;
	int ftbase_fftcodeletplan = 2;
	int ftbase_fhtcooleytukeyplan = 3;
	int ftbase_fhtcodeletplan = 4;
	int ftbase_fftrealcooleytukeyplan = 5;
	int ftbase_fftemptyplan = 6;
	PIVector<double> & tmpb(curplan.tmpbuf);

	if( curplan.plan[entryoffset+3]==ftbase_fftemptyplan )
		return;
	if( curplan.plan[entryoffset+3]==ftbase_fftcooleytukeyplan ) {
		n1 = curplan.plan[entryoffset+1];
		n2 = curplan.plan[entryoffset+2];
		ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
		for(int i=0; i<=n2-1; i++)
			ftbaseexecuteplanrec(a, aoffset+i*n1*2, plan, curplan.plan[entryoffset+5], stackptr);
		ftbase_ffttwcalc(a, aoffset, n1, n2);
		ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
		for(int i=0; i<=n1-1; i++)
			ftbaseexecuteplanrec(a, aoffset+i*n2*2, plan, curplan.plan[entryoffset+6], stackptr);
		ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
		return;
	}
	if( curplan.plan[entryoffset+3]==ftbase_fftrealcooleytukeyplan ) {
		n1 = curplan.plan[entryoffset+1];
		n2 = curplan.plan[entryoffset+2];
		ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
		for(int i=0; i<=n1/2-1; i++) {
			offs = aoffset+2*i*n2*2;
			for(int k=0; k<=n2-1; k++)
				(*a)[offs+2*k+1] = (*a)[offs+2*n2+2*k+0];
			ftbaseexecuteplanrec(a, offs, plan, curplan.plan[entryoffset+6], stackptr);
			tmpb[0] = (*a)[offs+0];
			tmpb[1] = 0;
			tmpb[2*n2+0] = (*a)[offs+1];
			tmpb[2*n2+1] = 0;
			for(int k=1; k<=n2-1; k++) {
				offs1 = 2*k;
				offs2 = 2*n2+2*k;
				hk = (*a)[offs+2*k+0];
				hnk = (*a)[offs+2*(n2-k)+0];
				tmpb[offs1+0] = 0.5*(hk+hnk);
				tmpb[offs2+1] = -0.5*(hk-hnk);
				hk = (*a)[offs+2*k+1];
				hnk = (*a)[offs+2*(n2-k)+1];
				tmpb[offs2+0] = 0.5*(hk+hnk);
				tmpb[offs1+1] = 0.5*(hk-hnk);
			}
			for (int i=0; i<2*n2*2; i++) (*a)[offs+i] = tmpb[i];
		}
		if( n1%2!=0 )
			ftbaseexecuteplanrec(a, aoffset+(n1-1)*n2*2, plan, curplan.plan[entryoffset+6], stackptr);
		ftbase_ffttwcalc(a, aoffset, n2, n1);
		ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
		for(int i=0; i<=n2-1; i++)
			ftbaseexecuteplanrec(a, aoffset+i*n1*2, plan, curplan.plan[entryoffset+5], stackptr);
		ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
		return;
	}
	if( curplan.plan[entryoffset+3]==ftbase_fhtcooleytukeyplan ) {
		n1 = curplan.plan[entryoffset+1];
		n2 = curplan.plan[entryoffset+2];
		n = n1*n2;
		ftbase_internalreallintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
		for(int i=0; i<=n2-1; i++)
			ftbaseexecuteplanrec(a, aoffset+i*n1, plan, curplan.plan[entryoffset+5], stackptr);
		for(int i=0; i<=n2-1; i++) {
			for(int j=0; j<=n1-1; j++) {
				offsa = aoffset+i*n1;
				hk = (*a)[offsa+j];
				hnk = (*a)[offsa+(n1-j)%n1];
				offs = 2*(i*n1+j);
				tmpb[offs+0] = -0.5*(hnk-hk);
				tmpb[offs+1] = 0.5*(hk+hnk);
			}
		}
		ftbase_ffttwcalc(&(curplan.tmpbuf), 0, n1, n2);
		for(int j=0; j<=n1-1; j++)
			(*a)[aoffset+j] = tmpb[2*j+0]+tmpb[2*j+1];
		if( n2%2==0 ) {
			offs = 2*(n2/2)*n1;
			offsa = aoffset+n2/2*n1;
			for(int j=0; j<=n1-1; j++)
				(*a)[offsa+j] = tmpb[offs+2*j+0]+tmpb[offs+2*j+1];
		}
		for(int i=1; i<=(n2+1)/2-1; i++) {
			offs = 2*i*n1;
			offs2 = 2*(n2-i)*n1;
			offsa = aoffset+i*n1;
			for(int j=0; j<=n1-1; j++)
				(*a)[offsa+j] = tmpb[offs+2*j+1]+tmpb[offs2+2*j+0];
			offsa = aoffset+(n2-i)*n1;
			for(int j=0; j<=n1-1; j++)
				(*a)[offsa+j] = tmpb[offs+2*j+0]+tmpb[offs2+2*j+1];
		}
		ftbase_internalreallintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
		for(int i=0; i<=n1-1; i++)
			ftbaseexecuteplanrec(a, aoffset+i*n2, plan, curplan.plan[entryoffset+6], stackptr);
		ftbase_internalreallintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
		return;
	}
	if( curplan.plan[entryoffset+3]==ftbase_fftcodeletplan ) {
		n1 = curplan.plan[entryoffset+1];
		n2 = curplan.plan[entryoffset+2];
		n = n1*n2;
		if( n==2 ) {
			a0x = (*a)[aoffset+0];
			a0y = (*a)[aoffset+1];
			a1x = (*a)[aoffset+2];
			a1y = (*a)[aoffset+3];
			v0 = a0x+a1x;
			v1 = a0y+a1y;
			v2 = a0x-a1x;
			v3 = a0y-a1y;
			(*a)[aoffset+0] = v0;
			(*a)[aoffset+1] = v1;
			(*a)[aoffset+2] = v2;
			(*a)[aoffset+3] = v3;
			return;
		}
		if( n==3 ) {
			offs = curplan.plan[entryoffset+7];
			c1 = curplan.precomputed[offs+0];
			c2 = curplan.precomputed[offs+1];
			a0x = (*a)[aoffset+0];
			a0y = (*a)[aoffset+1];
			a1x = (*a)[aoffset+2];
			a1y = (*a)[aoffset+3];
			a2x = (*a)[aoffset+4];
			a2y = (*a)[aoffset+5];
			t1x = a1x+a2x;
			t1y = a1y+a2y;
			a0x = a0x+t1x;
			a0y = a0y+t1y;
			m1x = c1*t1x;
			m1y = c1*t1y;
			m2x = c2*(a1y-a2y);
			m2y = c2*(a2x-a1x);
			s1x = a0x+m1x;
			s1y = a0y+m1y;
			a1x = s1x+m2x;
			a1y = s1y+m2y;
			a2x = s1x-m2x;
			a2y = s1y-m2y;
			(*a)[aoffset+0] = a0x;
			(*a)[aoffset+1] = a0y;
			(*a)[aoffset+2] = a1x;
			(*a)[aoffset+3] = a1y;
			(*a)[aoffset+4] = a2x;
			(*a)[aoffset+5] = a2y;
			return;
		}
		if( n==4 ) {
			a0x = (*a)[aoffset+0];
			a0y = (*a)[aoffset+1];
			a1x = (*a)[aoffset+2];
			a1y = (*a)[aoffset+3];
			a2x = (*a)[aoffset+4];
			a2y = (*a)[aoffset+5];
			a3x = (*a)[aoffset+6];
			a3y = (*a)[aoffset+7];
			t1x = a0x+a2x;
			t1y = a0y+a2y;
			t2x = a1x+a3x;
			t2y = a1y+a3y;
			m2x = a0x-a2x;
			m2y = a0y-a2y;
			m3x = a1y-a3y;
			m3y = a3x-a1x;
			(*a)[aoffset+0] = t1x+t2x;
			(*a)[aoffset+1] = t1y+t2y;
			(*a)[aoffset+4] = t1x-t2x;
			(*a)[aoffset+5] = t1y-t2y;
			(*a)[aoffset+2] = m2x+m3x;
			(*a)[aoffset+3] = m2y+m3y;
			(*a)[aoffset+6] = m2x-m3x;
			(*a)[aoffset+7] = m2y-m3y;
			return;
		}
		if( n==5 ) {
			offs = curplan.plan[entryoffset+7];
			c1 = curplan.precomputed[offs+0];
			c2 = curplan.precomputed[offs+1];
			c3 = curplan.precomputed[offs+2];
			c4 = curplan.precomputed[offs+3];
			c5 = curplan.precomputed[offs+4];
			t1x = (*a)[aoffset+2]+(*a)[aoffset+8];
			t1y = (*a)[aoffset+3]+(*a)[aoffset+9];
			t2x = (*a)[aoffset+4]+(*a)[aoffset+6];
			t2y = (*a)[aoffset+5]+(*a)[aoffset+7];
			t3x = (*a)[aoffset+2]-(*a)[aoffset+8];
			t3y = (*a)[aoffset+3]-(*a)[aoffset+9];
			t4x = (*a)[aoffset+6]-(*a)[aoffset+4];
			t4y = (*a)[aoffset+7]-(*a)[aoffset+5];
			t5x = t1x+t2x;
			t5y = t1y+t2y;
			(*a)[aoffset+0] = (*a)[aoffset+0]+t5x;
			(*a)[aoffset+1] = (*a)[aoffset+1]+t5y;
			m1x = c1*t5x;
			m1y = c1*t5y;
			m2x = c2*(t1x-t2x);
			m2y = c2*(t1y-t2y);
			m3x = -c3*(t3y+t4y);
			m3y = c3*(t3x+t4x);
			m4x = -c4*t4y;
			m4y = c4*t4x;
			m5x = -c5*t3y;
			m5y = c5*t3x;
			s3x = m3x-m4x;
			s3y = m3y-m4y;
			s5x = m3x+m5x;
			s5y = m3y+m5y;
			s1x = (*a)[aoffset+0]+m1x;
			s1y = (*a)[aoffset+1]+m1y;
			s2x = s1x+m2x;
			s2y = s1y+m2y;
			s4x = s1x-m2x;
			s4y = s1y-m2y;
			(*a)[aoffset+2] = s2x+s3x;
			(*a)[aoffset+3] = s2y+s3y;
			(*a)[aoffset+4] = s4x+s5x;
			(*a)[aoffset+5] = s4y+s5y;
			(*a)[aoffset+6] = s4x-s5x;
			(*a)[aoffset+7] = s4y-s5y;
			(*a)[aoffset+8] = s2x-s3x;
			(*a)[aoffset+9] = s2y-s3y;
			return;
		}
	}
	if( curplan.plan[entryoffset+3]==ftbase_fhtcodeletplan ) {
		n1 = curplan.plan[entryoffset+1];
		n2 = curplan.plan[entryoffset+2];
		n = n1*n2;
		if( n==2 ) {
			a0x = (*a)[aoffset+0];
			a1x = (*a)[aoffset+1];
			(*a)[aoffset+0] = a0x+a1x;
			(*a)[aoffset+1] = a0x-a1x;
			return;
		}
		if( n==3 ) {
			offs = curplan.plan[entryoffset+7];
			c1 = curplan.precomputed[offs+0];
			c2 = curplan.precomputed[offs+1];
			a0x = (*a)[aoffset+0];
			a1x = (*a)[aoffset+1];
			a2x = (*a)[aoffset+2];
			t1x = a1x+a2x;
			a0x = a0x+t1x;
			m1x = c1*t1x;
			m2y = c2*(a2x-a1x);
			s1x = a0x+m1x;
			(*a)[aoffset+0] = a0x;
			(*a)[aoffset+1] = s1x-m2y;
			(*a)[aoffset+2] = s1x+m2y;
			return;
		}
		if( n==4 ) {
			a0x = (*a)[aoffset+0];
			a1x = (*a)[aoffset+1];
			a2x = (*a)[aoffset+2];
			a3x = (*a)[aoffset+3];
			t1x = a0x+a2x;
			t2x = a1x+a3x;
			m2x = a0x-a2x;
			m3y = a3x-a1x;
			(*a)[aoffset+0] = t1x+t2x;
			(*a)[aoffset+1] = m2x-m3y;
			(*a)[aoffset+2] = t1x-t2x;
			(*a)[aoffset+3] = m2x+m3y;
			return;
		}
		if( n==5 ) {
			offs = curplan.plan[entryoffset+7];
			c1 = curplan.precomputed[offs+0];
			c2 = curplan.precomputed[offs+1];
			c3 = curplan.precomputed[offs+2];
			c4 = curplan.precomputed[offs+3];
			c5 = curplan.precomputed[offs+4];
			t1x = (*a)[aoffset+1]+(*a)[aoffset+4];
			t2x = (*a)[aoffset+2]+(*a)[aoffset+3];
			t3x = (*a)[aoffset+1]-(*a)[aoffset+4];
			t4x = (*a)[aoffset+3]-(*a)[aoffset+2];
			t5x = t1x+t2x;
			v0 = (*a)[aoffset+0]+t5x;
			(*a)[aoffset+0] = v0;
			m2x = c2*(t1x-t2x);
			m3y = c3*(t3x+t4x);
			s3y = m3y-c4*t4x;
			s5y = m3y+c5*t3x;
			s1x = v0+c1*t5x;
			s2x = s1x+m2x;
			s4x = s1x-m2x;
			(*a)[aoffset+1] = s2x-s3y;
			(*a)[aoffset+2] = s4x-s5y;
			(*a)[aoffset+3] = s4x+s5y;
			(*a)[aoffset+4] = s2x+s3y;
			return;
		}
	}
	if( curplan.plan[entryoffset+3]==ftbase_fftbluesteinplan ) {
		n = curplan.plan[entryoffset+1];
		m = curplan.plan[entryoffset+4];
		offs = curplan.plan[entryoffset+7];
		for(int i=stackptr+2*n; i<=stackptr+2*m-1; i++)
			curplan.stackbuf[i] = 0;
		offsp = offs+2*m;
		offsa = aoffset;
		offsb = stackptr;
		for(int i=0; i<n; i++) {
			bx = curplan.precomputed[offsp+0];
			by = curplan.precomputed[offsp+1];
			x = (*a)[offsa+0];
			y = (*a)[offsa+1];
			curplan.stackbuf[offsb+0] = x*bx-y*(-by);
			curplan.stackbuf[offsb+1] = x*(-by)+y*bx;
			offsp = offsp+2;
			offsa = offsa+2;
			offsb = offsb+2;
		}
		ftbaseexecuteplanrec(&curplan.stackbuf, stackptr, plan, curplan.plan[entryoffset+5], stackptr+2*2*m);
		offsb = stackptr;
		offsp = offs;
		for(int i=0; i<=m-1; i++) {
			x = curplan.stackbuf[offsb+0];
			y = curplan.stackbuf[offsb+1];
			bx = curplan.precomputed[offsp+0];
			by = curplan.precomputed[offsp+1];
			curplan.stackbuf[offsb+0] = x*bx-y*by;
			curplan.stackbuf[offsb+1] = -(x*by+y*bx);
			offsb = offsb+2;
			offsp = offsp+2;
		}
		ftbaseexecuteplanrec(&curplan.stackbuf, stackptr, plan, curplan.plan[entryoffset+5], stackptr+2*2*m);
		offsb = stackptr;
		offsp = offs+2*m;
		offsa = aoffset;
		for(int i=0; i<n; i++) {
			x = curplan.stackbuf[offsb+0]/m;
			y = -curplan.stackbuf[offsb+1]/m;
			bx = curplan.precomputed[offsp+0];
			by = curplan.precomputed[offsp+1];
			(*a)[offsa+0] = x*bx-y*(-by);
			(*a)[offsa+1] = x*(-by)+y*bx;
			offsp = offsp+2;
			offsa = offsa+2;
			offsb = offsb+2;
		}
		return;
	}
}


/*************************************************************************
Twiddle factors calculation

-- ALGLIB --
	Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT::ftbase_ffttwcalc(PIVector<double> * a, int aoffset, int n1, int n2) {
	int n, idx, offs;
	double x, y, twxm1, twy, twbasexm1, twbasey, twrowxm1, twrowy, tmpx, tmpy, v;
	int ftbase_ftbaseupdatetw = 4;
	n = n1*n2;
	v = -2*M_PI/n;
	twbasexm1 = -2*sqr(sin(0.5*v));
	twbasey = sin(v);
	twrowxm1 = 0;
	twrowy = 0;
	for(int i=0, j = 0; i<=n2-1; i++) {
		twxm1 = 0;
		twy = 0;
		for(j=0; j<=n1-1; j++) {
			idx = i*n1+j;
			offs = aoffset+2*idx;
			x = (*a)[offs+0];
			y = (*a)[offs+1];
			tmpx = x*twxm1-y*twy;
			tmpy = x*twy+y*twxm1;
			(*a)[offs+0] = x+tmpx;
			(*a)[offs+1] = y+tmpy;
			if( j<n1-1 ) {
				if( j%ftbase_ftbaseupdatetw==0 ) {
					v = -2*M_PI*i*(j+1)/n;
					twxm1 = -2*sqr(sin(0.5*v));
					twy = sin(v);
				} else {
					tmpx = twrowxm1+twxm1*twrowxm1-twy*twrowy;
					tmpy = twrowy+twxm1*twrowy+twy*twrowxm1;
					twxm1 = twxm1+tmpx;
					twy = twy+tmpy;
				}
			}
		}

		if( i<n2-1 ) {
			if( j%ftbase_ftbaseupdatetw==0 ) {
				v = -2*M_PI*(i+1)/n;
				twrowxm1 = -2*sqr(sin(0.5*v));
				twrowy = sin(v);
			} else {
				tmpx = twbasexm1+twrowxm1*twbasexm1-twrowy*twbasey;
				tmpy = twbasey+twrowxm1*twbasey+twrowy*twbasexm1;
				twrowxm1 = twrowxm1+tmpx;
				twrowy = twrowy+tmpy;
			}
		}
	}
}



PIStatistic::PIStatistic() {
	mean = 0.;
	variance = 0.;
	skewness = 0.;
	kurtosis = 0.;
}


bool PIStatistic::calculate(const PIVector<double> & val) {
	double v = 0., v1 = 0., v2 = 0., stddev = 0.;
	int i, n = val.size();
	if (n < 2)
		return false;
	/*
	* Mean
	*/
	for (i = 0; i < n; i++)
		mean += val[i];
	mean /= n;
	/*
	* Variance (using corrected two-pass algorithm)
	*/
	for (i = 0; i < n; i++)
		v1 += sqr(val[i] - mean);
	for (i = 0; i < n; i++)
		v2 += val[i] - mean;
	v2 = sqr(v2) / n;
	variance = (v1 - v2) / (n - 1);
	if(variance < 0)
		variance = 0.;
	stddev = sqrt(variance);
	/*
	* Skewness and kurtosis
	*/
	if (stddev != 0) {
		for (i = 0; i < n; i++) {
			v = (val[i] - mean) / stddev;
			v2 = sqr(v);
			skewness = skewness + v2 * v;
			kurtosis = kurtosis + sqr(v2);
		}
		skewness /= n;
		kurtosis = kurtosis / n - 3.;
	}
	return true;
}