3.10.2013 - PIPeer release, PIConsole now can work as server and remote client. Remote console test program in directory "remote_console"

2013-10-03 16:04:02 +04:00
parent 9111640ca8
commit 4b90f2818e
56 changed files with 6422 additions and 673 deletions
--- a/pimath.cpp
+++ b/pimath.cpp
@@ -20,6 +20,443 @@
 #include "pimath.h"


+#ifndef PIP_MATH_J0
+double j0(const double & v) {
+	double x = v;
+	double xsq;
+	double nn;
+	double pzero;
+	double qzero;
+	double p1;
+	double q1;
+	double result;
+	if (x < 0) x = -x;
+	if (x > 8.) {
+		double xsq_;
+		double p2;
+		double q2;
+		double p3;
+		double q3;
+		xsq_ = 64. / (x * x);
+		p2 = 0.0;
+		p2 = 2485.271928957404011288128951 + xsq_ * p2;
+		p2 = 153982.6532623911470917825993 + xsq_ * p2;
+		p2 = 2016135.283049983642487182349 + xsq_ * p2;
+		p2 = 8413041.456550439208464315611 + xsq_ * p2;
+		p2 = 12332384.76817638145232406055 + xsq_ * p2;
+		p2 = 5393485.083869438325262122897 + xsq_ * p2;
+		q2 = 1.0;
+		q2 = 2615.700736920839685159081813 + xsq_ * q2;
+		q2 = 156001.7276940030940592769933 + xsq_ * q2;
+		q2 = 2025066.801570134013891035236 + xsq_ * q2;
+		q2 = 8426449.050629797331554404810 + xsq_ * q2;
+		q2 = 12338310.22786324960844856182 + xsq_ * q2;
+		q2 = 5393485.083869438325560444960 + xsq_ * q2;
+		p3 = -0.0;
+		p3 = -4.887199395841261531199129300  +xsq_ * p3;
+		p3 = -226.2630641933704113967255053  +xsq_ * p3;
+		p3 = -2365.956170779108192723612816  +xsq_ * p3;
+		p3 = -8239.066313485606568803548860  +xsq_ * p3;
+		p3 = -10381.41698748464093880530341  +xsq_ * p3;
+		p3 = -3984.617357595222463506790588  +xsq_ * p3;
+		q3 = 1.0;
+		q3 = 408.7714673983499223402830260 + xsq_ * q3;
+		q3 = 15704.89191515395519392882766 + xsq_ * q3;
+		q3 = 156021.3206679291652539287109 + xsq_ * q3;
+		q3 = 533291.3634216897168722255057 + xsq_ * q3;
+		q3 = 666745.4239319826986004038103 + xsq_ * q3;
+		q3 = 255015.5108860942382983170882 + xsq_ * q3;
+		pzero = p2 / q2;
+		qzero = 8. * p3 / q3 / x;
+		nn = x- M_PI / 4.;
+		result = sqrt(2. / M_PI / x) * (pzero * cos(nn) - qzero * sin(nn));
+		return result;
+	}
+	xsq = x * x;
+	p1 = 26857.86856980014981415848441;
+	p1 = -40504123.71833132706360663322 + xsq * p1;
+	p1 = 25071582855.36881945555156435 + xsq * p1;
+	p1 = -8085222034853.793871199468171 + xsq * p1;
+	p1 = 1434354939140344.111664316553 + xsq * p1;
+	p1 = -136762035308817138.6865416609 + xsq * p1;
+	p1 = 6382059341072356562.289432465 + xsq * p1;
+	p1 = -117915762910761053603.8440800 + xsq * p1;
+	p1 = 493378725179413356181.6813446 + xsq * p1;
+	q1 = 1.;
+	q1 = 1363.063652328970604442810507 + xsq * q1;
+	q1 = 1114636.098462985378182402543 + xsq * q1;
+	q1 = 669998767.2982239671814028660 + xsq * q1;
+	q1 = 312304311494.1213172572469442 + xsq * q1;
+	q1 = 112775673967979.8507056031594 + xsq * q1;
+	q1 = 30246356167094626.98627330784 + xsq * q1;
+	q1 = 5428918384092285160.200195092 + xsq * q1;
+	q1 = 493378725179413356211.3278438 + xsq * q1;
+	return p1 / q1;
+}
+#endif
+
+
+#ifndef PIP_MATH_J1
+double j1(const double & v) {
+	double x = v;
+	double s;
+	double xsq;
+	double nn;
+	double pzero;
+	double qzero;
+	double p1;
+	double q1;
+	double result;
+	s = sign(x);
+	if (x < 0)
+		x = -x;
+	if (x > 8.) {
+		double xsq_;
+		double p2;
+		double q2;
+		double p3;
+		double q3;
+		xsq_ = 64.0 / (x * x);
+		p2 = -1611.616644324610116477412898;
+		p2 = -109824.0554345934672737413139 + xsq_ * p2;
+		p2 = -1523529.351181137383255105722 + xsq_ * p2;
+		p2 = -6603373.248364939109255245434 + xsq_ * p2;
+		p2 = -9942246.505077641195658377899 + xsq_ * p2;
+		p2 = -4435757.816794127857114720794 + xsq_ * p2;
+		q2 = 1.0;
+		q2 = -1455.009440190496182453565068 + xsq_ * q2;
+		q2 = -107263.8599110382011903063867 + xsq_ * q2;
+		q2 = -1511809.506634160881644546358 + xsq_ * q2;
+		q2 = -6585339.479723087072826915069 + xsq_ * q2;
+		q2 = -9934124.389934585658967556309 + xsq_ * q2;
+		q2 = -4435757.816794127856828016962 + xsq_ * q2;
+		p3 = 35.26513384663603218592175580;
+		p3 = 1706.375429020768002061283546 + xsq_ * p3;
+		p3 = 18494.26287322386679652009819 + xsq_ * p3;
+		p3 = 66178.83658127083517939992166 + xsq_ * p3;
+		p3 = 85145.16067533570196555001171 + xsq_ * p3;
+		p3 = 33220.91340985722351859704442 + xsq_ * p3;
+		q3 = 1.0;
+		q3 = 863.8367769604990967475517183 + xsq_ * q3;
+		q3 = 37890.22974577220264142952256 + xsq_ * q3;
+		q3 = 400294.4358226697511708610813 + xsq_ * q3;
+		q3 = 1419460.669603720892855755253 + xsq_ * q3;
+		q3 = 1819458.042243997298924553839 + xsq_ * q3;
+		q3 = 708712.8194102874357377502472 + xsq_ * q3;
+		pzero = p2 / q2;
+		qzero = 8 * p3 / q3 / x;
+		nn = x - 3 * M_PI / 4;
+		result = sqrt(2 / M_PI / x) * (pzero * cos(nn) - qzero * sin(nn));
+		if (s < 0)
+			result = -result;
+		return result;
+	}
+	xsq = sqr(x);
+	p1 = 2701.122710892323414856790990;
+	p1 = -4695753.530642995859767162166 + xsq * p1;
+	p1 = 3413234182.301700539091292655 + xsq * p1;
+	p1 = -1322983480332.126453125473247 + xsq * p1;
+	p1 = 290879526383477.5409737601689 + xsq * p1;
+	p1 = -35888175699101060.50743641413 + xsq * p1;
+	p1 = 2316433580634002297.931815435 + xsq * p1;
+	p1 = -66721065689249162980.20941484 + xsq * p1;
+	p1 = 581199354001606143928.050809 + xsq * p1;
+	q1 = 1.0;
+	q1 = 1606.931573481487801970916749 + xsq * q1;
+	q1 = 1501793.594998585505921097578 + xsq * q1;
+	q1 = 1013863514.358673989967045588 + xsq * q1;
+	q1 = 524371026216.7649715406728642 + xsq * q1;
+	q1 = 208166122130760.7351240184229 + xsq * q1;
+	q1 = 60920613989175217.46105196863 + xsq * q1;
+	q1 = 11857707121903209998.37113348 + xsq * q1;
+	q1 = 1162398708003212287858.529400 + xsq * q1;
+	result = s * x * p1 / q1;
+	return result;
+}
+#endif
+
+
+#ifndef PIP_MATH_JN
+double jn(const int & n, const double & v) {
+	double x = v;
+	double pkm2;
+	double pkm1;
+	double pk;
+	double xk;
+	double r;
+	double ans;
+	int k;
+	int sg;
+	double result;
+	if (n < 0) {
+		n = -n;
+		if (n % 2 == 0)
+			sg = 1;
+		else
+			sg = -1;
+	} else
+		sg = 1;
+	if (x < 0) {
+		if (n % 2 != 0)
+			sg = -sg;
+		x = -x;
+	}
+	if (n == 0) {
+		result = sg * j0(x);
+		return result;
+	}
+	if (n == 1) {
+		result = sg * j1(x);
+		return result;
+	}
+	if (n == 2) {
+		if (x == 0)
+			result = 0;
+		else
+			result = sg * (2.0 * j1(x) / x - j0(x));
+		return result;
+	}
+	if (x < 1E-16) {
+		result = 0;
+		return result;
+	}
+	k = 53;
+	pk = 2 * (n + k);
+	ans = pk;
+	xk = x * x;
+	do {
+		pk = pk - 2.0;
+		ans = pk - xk / ans;
+		k = k - 1;
+	} while (k != 0);
+	ans = x / ans;
+	pk = 1.0;
+	pkm1 = 1.0 / ans;
+	k = n - 1;
+	r = 2 * k;
+	do {
+		pkm2 = (pkm1 * r - pk * x) / x;
+		pk = pkm1;
+		pkm1 = pkm2;
+		r = r - 2.0;
+		k = k - 1;
+	} while (k != 0);
+	if (fabs(pk) > fabs(pkm1))
+		ans = j1(x) / pk;
+	else
+		ans = j0(x) / pkm1;
+	result = sg * ans;
+	return result;
+}
+#endif
+
+
+#ifndef PIP_MATH_Y0
+double y0(const double & v) {
+	double x = v;
+	double nn;
+	double xsq;
+	double pzero;
+	double qzero;
+	double p4;
+	double q4;
+	double result;
+	if (x > 8.) {
+		double xsq_;
+		double p2;
+		double q2;
+		double p3;
+		double q3;
+		xsq_ = 64.0 / (x * x);
+		p2 = 0.0;
+		p2 = 2485.271928957404011288128951 + xsq_ * p2;
+		p2 = 153982.6532623911470917825993 + xsq_ * p2;
+		p2 = 2016135.283049983642487182349 + xsq_ * p2;
+		p2 = 8413041.456550439208464315611 + xsq_ * p2;
+		p2 = 12332384.76817638145232406055 + xsq_ * p2;
+		p2 = 5393485.083869438325262122897 + xsq_ * p2;
+		q2 = 1.0;
+		q2 = 2615.700736920839685159081813 + xsq_ * q2;
+		q2 = 156001.7276940030940592769933 + xsq_ * q2;
+		q2 = 2025066.801570134013891035236 + xsq_ * q2;
+		q2 = 8426449.050629797331554404810 + xsq_ * q2;
+		q2 = 12338310.22786324960844856182 + xsq_ * q2;
+		q2 = 5393485.083869438325560444960 + xsq_ * q2;
+		p3 = -0.0;
+		p3 = -4.887199395841261531199129300 + xsq_ * p3;
+		p3 = -226.2630641933704113967255053 + xsq_ * p3;
+		p3 = -2365.956170779108192723612816 + xsq_ * p3;
+		p3 = -8239.066313485606568803548860 + xsq_ * p3;
+		p3 = -10381.41698748464093880530341 + xsq_ * p3;
+		p3 = -3984.617357595222463506790588 + xsq_ * p3;
+		q3 = 1.0;
+		q3 = 408.7714673983499223402830260 + xsq_ * q3;
+		q3 = 15704.89191515395519392882766 + xsq_ * q3;
+		q3 = 156021.3206679291652539287109 + xsq_ * q3;
+		q3 = 533291.3634216897168722255057 + xsq_ * q3;
+		q3 = 666745.4239319826986004038103 + xsq_ * q3;
+		q3 = 255015.5108860942382983170882 + xsq_ * q3;
+		pzero = p2 / q2;
+		qzero = 8 * p3 / q3 / x;
+		nn = x - M_PI / 4;
+		result = sqrt(2 / M_PI / x) * (pzero * sin(nn) + qzero * cos(nn));
+		return result;
+	}
+	xsq = sqr(x);
+	p4 = -41370.35497933148554125235152;
+	p4 = 59152134.65686889654273830069 + xsq * p4;
+	p4 = -34363712229.79040378171030138 + xsq * p4;
+	p4 = 10255208596863.94284509167421 + xsq * p4;
+	p4 = -1648605817185729.473122082537 + xsq * p4;
+	p4 = 137562431639934407.8571335453 + xsq * p4;
+	p4 = -5247065581112764941.297350814 + xsq * p4;
+	p4 = 65874732757195549259.99402049 + xsq * p4;
+	p4 = -27502866786291095837.01933175 + xsq * p4;
+	q4 = 1.0;
+	q4 = 1282.452772478993804176329391 + xsq * q4;
+	q4 = 1001702.641288906265666651753 + xsq * q4;
+	q4 = 579512264.0700729537480087915 + xsq * q4;
+	q4 = 261306575504.1081249568482092 + xsq * q4;
+	q4 = 91620380340751.85262489147968 + xsq * q4;
+	q4 = 23928830434997818.57439356652 + xsq * q4;
+	q4 = 4192417043410839973.904769661 + xsq * q4;
+	q4 = 372645883898616588198.9980 + xsq * q4;
+	result = p4 / q4 + 2 / M_PI * j0(x) * log(x);
+	return result;
+}
+#endif
+
+
+#ifndef PIP_MATH_Y1
+double y1(const double & v) {
+	double x = v;
+	double nn;
+	double xsq;
+	double pzero;
+	double qzero;
+	double p4;
+	double q4;
+	double result;
+	if (x > 8.) {
+		double xsq_;
+		double p2;
+		double q2;
+		double p3;
+		double q3;
+		xsq_ = 64.0 / (x * x);
+		p2 = -1611.616644324610116477412898;
+		p2 = -109824.0554345934672737413139 + xsq_ * p2;
+		p2 = -1523529.351181137383255105722 + xsq_ * p2;
+		p2 = -6603373.248364939109255245434 + xsq_ * p2;
+		p2 = -9942246.505077641195658377899 + xsq_ * p2;
+		p2 = -4435757.816794127857114720794 + xsq_ * p2;
+		q2 = 1.0;
+		q2 = -1455.009440190496182453565068 + xsq_ * q2;
+		q2 = -107263.8599110382011903063867 + xsq_ * q2;
+		q2 = -1511809.506634160881644546358 + xsq_ * q2;
+		q2 = -6585339.479723087072826915069 + xsq_ * q2;
+		q2 = -9934124.389934585658967556309 + xsq_ * q2;
+		q2 = -4435757.816794127856828016962 + xsq_ * q2;
+		p3 = 35.26513384663603218592175580;
+		p3 = 1706.375429020768002061283546 + xsq_ * p3;
+		p3 = 18494.26287322386679652009819 + xsq_ * p3;
+		p3 = 66178.83658127083517939992166 + xsq_ * p3;
+		p3 = 85145.16067533570196555001171 + xsq_ * p3;
+		p3 = 33220.91340985722351859704442 + xsq_ * p3;
+		q3 = 1.0;
+		q3 = 863.8367769604990967475517183 + xsq_ * q3;
+		q3 = 37890.22974577220264142952256 + xsq_ * q3;
+		q3 = 400294.4358226697511708610813 + xsq_ * q3;
+		q3 = 1419460.669603720892855755253 + xsq_ * q3;
+		q3 = 1819458.042243997298924553839 + xsq_ * q3;
+		q3 = 708712.8194102874357377502472 + xsq_ * q3;
+		pzero = p2 / q2;
+		qzero = 8 * p3 / q3 / x;
+		nn = x - 3 * M_PI / 4;
+		result = sqrt(2 / M_PI / x) * (pzero * sin(nn) + qzero * cos(nn));
+		return result;
+	}
+	xsq = sqr(x);
+	p4 = -2108847.540133123652824139923;
+	p4 = 3639488548.124002058278999428 + xsq * p4;
+	p4 = -2580681702194.450950541426399 + xsq * p4;
+	p4 = 956993023992168.3481121552788 + xsq * p4;
+	p4 = -196588746272214065.8820322248 + xsq * p4;
+	p4 = 21931073399177975921.11427556 + xsq * p4;
+	p4 = -1212297555414509577913.561535 + xsq * p4;
+	p4 = 26554738314348543268942.48968 + xsq * p4;
+	p4 = -99637534243069222259967.44354 + xsq * p4;
+	q4 = 1.0;
+	q4 = 1612.361029677000859332072312 + xsq * q4;
+	q4 = 1563282.754899580604737366452 + xsq * q4;
+	q4 = 1128686837.169442121732366891 + xsq * q4;
+	q4 = 646534088126.5275571961681500 + xsq * q4;
+	q4 = 297663212564727.6729292742282 + xsq * q4;
+	q4 = 108225825940881955.2553850180 + xsq * q4;
+	q4 = 29549879358971486742.90758119 + xsq * q4;
+	q4 = 5435310377188854170800.653097 + xsq * q4;
+	q4 = 508206736694124324531442.4152 + xsq * q4;
+	result = x * p4 / q4 + 2 / M_PI * (j1(x) * log(x) - 1 / x);
+	return result;
+}
+#endif
+
+
+#ifndef PIP_MATH_YN
+double yn(const int & n, const double & v) {
+	int i;
+	double x = v;
+	double a;
+	double b;
+	double tmp;
+	double s;
+	double result;
+	s = 1;
+	if (n < 0) {
+		n = -n;
+		if (n % 2 != 0)
+			s = -1;
+	}
+	if (n == 0) {
+		result = y0(x);
+		return result;
+	}
+	if (n == 1) {
+		result = s * y1(x);
+		return result;
+	}
+	a = y0(x);
+	b = y1(x);
+	for (i = 1; i <= n - 1; i++) {
+		tmp = b;
+		b = 2 * i / x * b - a;
+		a = tmp;
+	}
+	result = s * b;
+	return result;
+}
+#endif
+
+
+double randomn(double dv, double sv) {
+	static bool agen = false;
+	double s = 2., v0 = 0., v1 = 0.;
+	if (agen) {
+		agen = false;
+		v1 = v1 * sqrt(-2 * log(s) / s);
+		return v1 * sv + dv;
+	}
+	while (s > 1. || s == 0.) {
+		v0 = randomd();
+		v1 = randomd();
+		s = v0*v0 + v1*v1;
+	}
+	v0 = v0 * sqrt(-2 * log(s) / s);
+	return v0 * sv + dv;
+}
+
+
+
 const char Solver::methods_desc[] = "b{Methods:}\
 \n  -1 - Global settings\
 \n  01 - Eyler 1\
@@ -70,7 +507,7 @@ 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;
+		piCout << "Solver error: {A} should be greater than {B}";
 		return;
 	}
 	if (size == 0) return;
@@ -327,7 +764,7 @@ void PIFFT::fftc1d(const PIVector<complexd> &a, uint n) {

 void PIFFT::fftc1r(const PIVector<double> & a, uint n) {
 	uint i;
-	if( n%2==0 ) {
+	if( n%2==0) {
 		PIVector<double> buf;
 		uint n2 = n/2;
 		//buf.resize(n);
@@ -455,14 +892,14 @@ void PIFFT::ftbase_ftbasegenerateplanrec(
 	int ftbase_fftcodeletplan = 2;
 	int ftbase_fftrealcooleytukeyplan = 5;
 	int ftbase_fftemptyplan = 6;
-	if( *plansize+ftbase_ftbaseplanentrysize>(*planarraysize) ) {
+	if( *plansize+ftbase_ftbaseplanentrysize>(*planarraysize)) {
 		curplan.plan.resize(8*(*planarraysize));
 		*planarraysize = 8*(*planarraysize);
 	}
 	entryoffset = *plansize;
 	esize = ftbase_ftbaseplanentrysize;
 	*plansize = *plansize+esize;
-	if( n==1 ) {
+	if( n==1) {
 		curplan.plan[entryoffset+0] = esize;
 		curplan.plan[entryoffset+1] = -1;
 		curplan.plan[entryoffset+2] = -1;
@@ -474,12 +911,12 @@ void PIFFT::ftbase_ftbasegenerateplanrec(
 		return;
 	}
 	ftbasefactorize(n, &n1, &n2);
-	if( n1!=1 ) {
+	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 )
+		if( tasktype==ftbase_ftbasecffttask)
 			curplan.plan[entryoffset+3] = ftbase_fftcooleytukeyplan;
 		else
 			curplan.plan[entryoffset+3] = ftbase_fftrealcooleytukeyplan;
@@ -501,9 +938,9 @@ void PIFFT::ftbase_ftbasegenerateplanrec(
 			curplan.plan[entryoffset+5] = -1;
 			curplan.plan[entryoffset+6] = -1;
 			curplan.plan[entryoffset+7] = *precomputedsize;
-			if( n==3 )
+			if( n==3)
 				*precomputedsize = *precomputedsize+2;
-			if( n==5 )
+			if( n==5)
 				*precomputedsize = *precomputedsize+5;
 			return;
 		} else {
@@ -547,22 +984,22 @@ void PIFFT::ftbase_ftbaseprecomputeplanrec(ftplan* plan,
 	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 ) {
+	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 ) {
+	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 ) {
+		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 ) {
+		if( n==5) {
 			offs = curplan.plan[entryoffset+7];
 			v = 2*M_PI/5;
 			curplan.precomputed[offs+0] = (cos(v)+cos(2*v))/2-1;
@@ -573,7 +1010,7 @@ void PIFFT::ftbase_ftbaseprecomputeplanrec(ftplan* plan,
 			return;
 		}
 	}
-	if( curplan.plan[entryoffset+3]==ftbase_fftbluesteinplan ) {
+	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];
@@ -587,7 +1024,7 @@ void PIFFT::ftbase_ftbaseprecomputeplanrec(ftplan* plan,
 			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 ) {
+			if( i>0) {
 				curplan.precomputed[offs+2*(m-i)+0] = bx;
 				curplan.precomputed[offs+2*(m-i)+1] = by;
 			}
@@ -601,29 +1038,29 @@ void PIFFT::ftbase_ftbaseprecomputeplanrec(ftplan* plan,
 void PIFFT::ftbasefactorize(int n, int* n1, int* n2) {
 	*n1 = *n2 = 0;
 	int ftbase_ftbasecodeletrecommended = 5;
-	if( (*n1)*(*n2)!=n ) {
+	if( (*n1)*(*n2)!=n) {
 		for(int j=ftbase_ftbasecodeletrecommended; j>=2; j--) {
-			if( n%j==0 ) {
+			if( n%j==0) {
 				*n1 = j;
 				*n2 = n/j;
 				break;
 			}
 		}
 	}
-	if( (*n1)*(*n2)!=n ) {
+	if( (*n1)*(*n2)!=n) {
 		for(int j=ftbase_ftbasecodeletrecommended+1; j<=n-1; j++) {
-			if( n%j==0 ) {
+			if( n%j==0) {
 				*n1 = j;
 				*n2 = n/j;
 				break;
 			}
 		}
 	}
-	if( (*n1)*(*n2)!=n ) {
+	if( (*n1)*(*n2)!=n) {
 		*n1 = 1;
 		*n2 = n;
 	}
-	if( (*n2)==1 && (*n1)!=1 ) {
+	if( (*n2)==1 && (*n1)!=1) {
 		*n2 = *n1;
 		*n1 = 1;
 	}
@@ -637,15 +1074,15 @@ Is number smooth?
 	Copyright 01.05.2009 by Bochkanov Sergey
 *************************************************************************/
 void PIFFT::ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int* best) {
-	if( seed>=n ) {
+	if( seed>=n) {
        *best = piMin<int>(*best, seed);
 		return;
 	}
-	if( leastfactor<=2 )
+	if( leastfactor<=2)
 		ftbase_ftbasefindsmoothrec(n, seed*2, 2, best);
-	if( leastfactor<=3 )
+	if( leastfactor<=3)
 		ftbase_ftbasefindsmoothrec(n, seed*3, 3, best);
-	if( leastfactor<=5 )
+	if( leastfactor<=5)
 		ftbase_ftbasefindsmoothrec(n, seed*5, 5, best);
 }

@@ -670,9 +1107,9 @@ void PIFFT::ftbase_internalreallintranspose(PIVector<double>* a, int m, int n, 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 )
+	if( m==0||n==0)
 		return;
-    if( piMax<int>(m, n)<=8 ) {
+    if( piMax<int>(m, n)<=8) {
 		for(int i=0; i<=m-1; i++) {
 			idx1 = bstart+i;
 			idx2 = astart+i*astride;
@@ -684,15 +1121,15 @@ void PIFFT::ftbase_fftirltrec(PIVector<double>* a, int astart, int astride, PIVe
 		}
 		return;
 	}
-	if( n>m ) {
+	if( n>m) {
 		n1 = n/2;
-		if( n-n1>=8&&n1%8!=0 )
+		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 )
+		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);
@@ -709,9 +1146,9 @@ void PIFFT::ftbase_internalcomplexlintranspose(PIVector<double>* a, int m, int n

 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 )
+	if( m==0||n==0)
 		return;
-    if( piMax<int>(m, n)<=8 ) {
+    if( piMax<int>(m, n)<=8) {
 		m2 = 2*bstride;
 		for(int i=0; i<=m-1; i++) {
 			idx1 = bstart+2*i;
@@ -725,15 +1162,15 @@ void PIFFT::ftbase_ffticltrec(PIVector<double>* a, int astart, int astride, PIVe
 		}
 		return;
 	}
-	if( n>m ) {
+	if( n>m) {
 		n1 = n/2;
-		if( n-n1>=8&&n1%8!=0 )
+		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 )
+		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);
@@ -775,9 +1212,9 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 	int ftbase_fftemptyplan = 6;
 	PIVector<double> & tmpb(curplan.tmpbuf);

-	if( curplan.plan[entryoffset+3]==ftbase_fftemptyplan )
+	if( curplan.plan[entryoffset+3]==ftbase_fftemptyplan)
 		return;
-	if( curplan.plan[entryoffset+3]==ftbase_fftcooleytukeyplan ) {
+	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));
@@ -790,7 +1227,7 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 		ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
 		return;
 	}
-	if( curplan.plan[entryoffset+3]==ftbase_fftrealcooleytukeyplan ) {
+	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));
@@ -817,7 +1254,7 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 			}
 			for (int i=0; i<2*n2*2; i++) (*a)[offs+i] = tmpb[i];
 		}
-		if( n1%2!=0 )
+		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));
@@ -826,7 +1263,7 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 		ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
 		return;
 	}
-	if( curplan.plan[entryoffset+3]==ftbase_fhtcooleytukeyplan ) {
+	if( curplan.plan[entryoffset+3]==ftbase_fhtcooleytukeyplan) {
 		n1 = curplan.plan[entryoffset+1];
 		n2 = curplan.plan[entryoffset+2];
 		n = n1*n2;
@@ -846,7 +1283,7 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 		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 ) {
+		if( n2%2==0) {
 			offs = 2*(n2/2)*n1;
 			offsa = aoffset+n2/2*n1;
 			for(int j=0; j<=n1-1; j++)
@@ -868,11 +1305,11 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 		ftbase_internalreallintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
 		return;
 	}
-	if( curplan.plan[entryoffset+3]==ftbase_fftcodeletplan ) {
+	if( curplan.plan[entryoffset+3]==ftbase_fftcodeletplan) {
 		n1 = curplan.plan[entryoffset+1];
 		n2 = curplan.plan[entryoffset+2];
 		n = n1*n2;
-		if( n==2 ) {
+		if( n==2) {
 			a0x = (*a)[aoffset+0];
 			a0y = (*a)[aoffset+1];
 			a1x = (*a)[aoffset+2];
@@ -887,7 +1324,7 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 			(*a)[aoffset+3] = v3;
 			return;
 		}
-		if( n==3 ) {
+		if( n==3) {
 			offs = curplan.plan[entryoffset+7];
 			c1 = curplan.precomputed[offs+0];
 			c2 = curplan.precomputed[offs+1];
@@ -919,7 +1356,7 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 			(*a)[aoffset+5] = a2y;
 			return;
 		}
-		if( n==4 ) {
+		if( n==4) {
 			a0x = (*a)[aoffset+0];
 			a0y = (*a)[aoffset+1];
 			a1x = (*a)[aoffset+2];
@@ -946,7 +1383,7 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 			(*a)[aoffset+7] = m2y-m3y;
 			return;
 		}
-		if( n==5 ) {
+		if( n==5) {
 			offs = curplan.plan[entryoffset+7];
 			c1 = curplan.precomputed[offs+0];
 			c2 = curplan.precomputed[offs+1];
@@ -996,18 +1433,18 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 			return;
 		}
 	}
-	if( curplan.plan[entryoffset+3]==ftbase_fhtcodeletplan ) {
+	if( curplan.plan[entryoffset+3]==ftbase_fhtcodeletplan) {
 		n1 = curplan.plan[entryoffset+1];
 		n2 = curplan.plan[entryoffset+2];
 		n = n1*n2;
-		if( n==2 ) {
+		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 ) {
+		if( n==3) {
 			offs = curplan.plan[entryoffset+7];
 			c1 = curplan.precomputed[offs+0];
 			c2 = curplan.precomputed[offs+1];
@@ -1024,7 +1461,7 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 			(*a)[aoffset+2] = s1x+m2y;
 			return;
 		}
-		if( n==4 ) {
+		if( n==4) {
 			a0x = (*a)[aoffset+0];
 			a1x = (*a)[aoffset+1];
 			a2x = (*a)[aoffset+2];
@@ -1039,7 +1476,7 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 			(*a)[aoffset+3] = m2x+m3y;
 			return;
 		}
-		if( n==5 ) {
+		if( n==5) {
 			offs = curplan.plan[entryoffset+7];
 			c1 = curplan.precomputed[offs+0];
 			c2 = curplan.precomputed[offs+1];
@@ -1067,7 +1504,7 @@ void PIFFT::ftbaseexecuteplanrec(PIVector<double>* a, int aoffset, ftplan* plan,
 			return;
 		}
 	}
-	if( curplan.plan[entryoffset+3]==ftbase_fftbluesteinplan ) {
+	if( curplan.plan[entryoffset+3]==ftbase_fftbluesteinplan) {
 		n = curplan.plan[entryoffset+1];
 		m = curplan.plan[entryoffset+4];
 		offs = curplan.plan[entryoffset+7];
@@ -1148,8 +1585,8 @@ void PIFFT::ftbase_ffttwcalc(PIVector<double> * a, int aoffset, int n1, int n2)
 			tmpy = x*twy+y*twxm1;
 			(*a)[offs+0] = x+tmpx;
 			(*a)[offs+1] = y+tmpy;
-			if( j<n1-1 ) {
-				if( j%ftbase_ftbaseupdatetw==0 ) {
+			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);
@@ -1162,8 +1599,8 @@ void PIFFT::ftbase_ffttwcalc(PIVector<double> * a, int aoffset, int n1, int n2)
 			}
 		}

-		if( i<n2-1 ) {
-			if( j%ftbase_ftbaseupdatetw==0 ) {
+		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);
@@ -1176,52 +1613,3 @@ void PIFFT::ftbase_ffttwcalc(PIVector<double> * a, int aoffset, int n1, int n2)
 		}
 	}
 }
-
-
-
-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;
-}