/*
	PIP - Platform Independent Primitives
	Class for FFT, IFFT and Hilbert transformations
	Andrey Bychkov work.a.b@yandex.ru, Ivan Pelipenko peri4ko@yandex.ru

	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU Lesser 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 Lesser General Public License for more details.

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

#include "pifft.h"


PIFFT_double::PIFFT_double() {
}


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


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


PIVector<complexd> * PIFFT_double::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_double::calcFFT(const PIVector<double> & val) {
	result.clear();
	if (val.size_s() < 4) return &result;
	fftc1r(val, val.size());
	return &result;
}


PIVector<double> PIFFT_double::getAmplitude() const {
	PIVector<double> ret;
	ret.resize(result.size());
	double tmp;
	for (uint i = 0; i < result.size(); i++) {
		tmp = sqrt(result[i].real() * result[i].real() + result[i].imag() * result[i].imag());
		ret[i] = tmp;
	}
	return ret;
}


PIVector<double> PIFFT_double::getReal() const {
	PIVector<double> ret;
	ret.resize(result.size());
	for (uint i = 0; i < result.size(); i++)
		ret[i] = result[i].real();
	return ret;
}


PIVector<double> PIFFT_double::getImag() const {
	PIVector<double> ret;
	ret.resize(result.size());
	for (uint i = 0; i < result.size(); i++)
		ret[i] = result[i].imag();
	return ret;
}


void PIFFT_double::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[i].real();
		buf[2 * i + 1] = a[i].imag();
	}
	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_double::fftc1r(const PIVector<double> & a, uint n) {
	uint i;
	if (n % 2 == 0) {
		PIVector<double> buf;
		uint n2 = n / 2;
		buf = a;
		createPlan(n2);
		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_double::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);
	for (i = 0; i < n; i++) {
		result[i] = conj(result[i] / (double)n);
	}
}


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


void PIFFT_double::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);
	assert(stackptr==0);
	curplan.stackbuf.resize(piMax<int>(stackmemsize, 1)); //ae_vector_set_length(&curplan.stackbuf, ae_maxint(stackmemsize, 1));
	curplan.tmpbuf.resize(piMax<int>(tmpmemsize, 1)); //ae_vector_set_length(&(curplan.tmpbuf), ae_maxint(tmpmemsize, 1));
	curplan.precomputed.resize(piMax<int>(precomputedsize, 1)); //ae_vector_set_length(&curplan.precomputed, ae_maxint(precomputedsize, 1));
	stackptr = 0;
	ftbase_ftbaseprecomputeplanrec(plan, 0, stackptr);
	assert(stackptr==0);
}


/*************************************************************************
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_double::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<int>(*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<int>(*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<int>(*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_double::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_double::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_double::ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int * best) {
	if (seed >= n) {
		*best = piMini(*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_double::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_double::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_double::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 (piMaxi(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_double::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_double::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<int>(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_double::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_double::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 k = 0; k < 2 * n2 * 2; k++) (*a)[offs + k] = tmpb[k];
		}
		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_double::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;
			}
		}
	}
}


//=================================================================================================
//*************************************************************************************************
//---------------------FLOAT-----------------------FLOAT-------------------------------------------
//*************************************************************************************************
//=================================================================================================


PIFFT_float::PIFFT_float() {
}


PIVector<complexf> * PIFFT_float::calcFFT(const PIVector<complexf> & val) {
	result.clear();
	if (val.size_s() < 4) return &result;
	fftc1d(val, val.size());
	return &result;
}


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


PIVector<complexf> * PIFFT_float::calcHilbert(const PIVector<float> & 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.f;
	for (uint i = result.size() / 2; i < result.size(); i++) result[i] = 0;
	fftc1dinv(result, result.size());
	return &result;
}


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


PIVector<float> PIFFT_float::getAmplitude() const {
	PIVector<float> ret;
	ret.resize(result.size());
	float tmp;
	for (uint i = 0; i < result.size(); i++) {
		tmp = sqrt(result[i].real() * result[i].real() + result[i].imag() * result[i].imag());
		ret[i] = tmp;
	}
	return ret;
}


PIVector<float> PIFFT_float::getReal() const {
	PIVector<float> ret;
	ret.resize(result.size());
	for (uint i = 0; i < result.size(); i++)
		ret[i] = result[i].real();
	return ret;
}


PIVector<float> PIFFT_float::getImag() const {
	PIVector<float> ret;
	ret.resize(result.size());
	for (uint i = 0; i < result.size(); i++)
		ret[i] = result[i].imag();
	return ret;
}


void PIFFT_float::fftc1d(const PIVector<complexf> & a, uint n) {
	createPlan(n);
	uint i;
	PIVector<float> buf;
	buf.resize(2 * n);
	for (i = 0; i < n; i++) {
		buf[2 * i + 0] = a[i].real();
		buf[2 * i + 1] = a[i].imag();
	}
	ftbaseexecuteplan(&buf, 0, n, &curplan);
	result.resize(n);
	for (i = 0; i < n; i++)
		result[i] = complexf(buf[2 * i + 0], buf[2 * i + 1]);
}


void PIFFT_float::fftc1r(const PIVector<float> & a, uint n) {
	uint i;
	if (n % 2 == 0) {
		PIVector<float> buf;
		uint n2 = n / 2;
		buf = a;
		createPlan(n2);
		ftbaseexecuteplan(&buf, 0, n2, &curplan);
		result.resize(n);
		uint idx;
		complexf hn, hmnc, v;
		for (i = 0; i <= n2; i++) {
			idx = 2 * (i % n2);
			hn = complexf(buf[idx + 0], buf[idx + 1]);
			idx = 2 * ((n2 - i) % n2);
			hmnc = complexf(buf[idx + 0], -buf[idx + 1]);
			v = complexf(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<complexf> cbuf;
		cbuf.resize(n);
		for (i = 0; i < n; i++)
			cbuf[i] = complexf(a[i], 0.);
		fftc1d(cbuf, n);
	}
}


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


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


void PIFFT_float::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);
	assertm(stackptr==0, "Internal error in FTBaseGenerateComplexFFTPlan");
	curplan.stackbuf.resize(piMax<int>(stackmemsize, 1)); //ae_vector_set_length(&curplan.stackbuf, ae_maxint(stackmemsize, 1));
	curplan.tmpbuf.resize(piMax<int>(tmpmemsize, 1)); //ae_vector_set_length(&(curplan.tmpbuf), ae_maxint(tmpmemsize, 1));
	curplan.precomputed.resize(piMax<int>(precomputedsize, 1)); //ae_vector_set_length(&curplan.precomputed, ae_maxint(precomputedsize, 1));
	stackptr = 0;
	ftbase_ftbaseprecomputeplanrec(plan, 0, stackptr);
	assertm(stackptr==0, "Internal error in FTBaseGenerateComplexFFTPlan");
}


/*************************************************************************
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_float::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 = 4;
	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(4 * (*planarraysize));
		*planarraysize = 4 * (*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<int>(*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<int>(*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<int>(*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_float::ftbase_ftbaseprecomputeplanrec(ftplan * plan,
int entryoffset,
ae_int_t stackptr) {
	int n1, n2, n, m, offs;
	float 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_float::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_float::ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int * best) {
	if (seed >= n) {
		*best = piMini(*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_float::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_float::ftbase_internalreallintranspose(PIVector<float> * a, int m, int n, int astart, PIVector<float> * 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_float::ftbase_fftirltrec(PIVector<float> * a, int astart, int astride, PIVector<float> * b, int bstart, int bstride, int m, int n) {
	int idx1, idx2;
	int m1, n1;
	if (m == 0 || n == 0)
		return;
	if (piMaxi(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_float::ftbase_internalcomplexlintranspose(PIVector<float> * a, int m, int n, int astart, PIVector<float> * 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_float::ftbase_ffticltrec(PIVector<float> * a, int astart, int astride, PIVector<float> * b, int bstart, int bstride, int m, int n) {
	int idx1, idx2, m2, m1, n1;
	if (m == 0 || n == 0)
		return;
	if (piMax<int>(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_float::ftbaseexecuteplan(PIVector<float> * 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 floats)

-- ALGLIB --
	Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT_float::ftbaseexecuteplanrec(PIVector<float> * a, int aoffset, ftplan * plan, int entryoffset, ae_int_t stackptr) {
	int n1, n2, n, m, offs, offs1, offs2, offsa, offsb, offsp;
	float hk, hnk, x, y, bx, by, v0, v1, v2, v3;
	float a0x, a0y, a1x, a1y, a2x, a2y, a3x, a3y;
	float t1x, t1y, t2x, t2y, t3x, t3y, t4x, t4y, t5x, t5y;
	float m1x, m1y, m2x, m2y, m3x, m3y, m4x, m4y, m5x, m5y;
	float s1x, s1y, s2x, s2y, s3x, s3y, s4x, s4y, s5x, s5y;
	float 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<float> & 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 k = 0; k < 2 * n2 * 2; k++) (*a)[offs + k] = tmpb[k];
		}
		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_float::ftbase_ffttwcalc(PIVector<float> * a, int aoffset, int n1, int n2) {
	int n, idx, offs;
	float 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;
			}
		}
	}
}