/*
	PIP - Platform Independent Primitives
	Class for FFT, IFFT and Hilbert transformations
	Copyright (C) 2015  Ivan Pelipenko peri4ko@gmail.com, Andrey Bychkov work.a.b@yandex.ru

	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 "pifft.h"


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


PIVector<complexd> * PIFFT::calcFFT(const PIVector<complexd> & val) {
	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[i].real() * result[i].real() + result[i].imag() * result[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[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::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::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::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<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);
	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<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::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 = 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::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 (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::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<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::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 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::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;
			}
		}
	}
}