git-svn-id: svn://db.shs.com.ru/pip@4 12ceb7fc-bf1f-11e4-8940-5bc7170c53b5

src/math/picrc.h

@@ -0,0 +1,244 @@
+/*! \file picrc.h
+ * \brief CRC checksum calculator
+*/
+/*
+    PIP - Platform Independent Primitives
+    CRC checksum calculator
+    Copyright (C) 2014  Ivan Pelipenko peri4ko@gmail.com
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef PICRC_H
+#define PICRC_H
+
+#include "pistring.h"
+
+template <int L>
+class PIP_EXPORT uint_cl {
+public:
+	uint_cl() {for (int i = 0; i < L / 8; ++i) data_[i] = 0;}
+	uint_cl(const uint_cl<L> & v) {for (int i = 0; i < L / 8; ++i) data_[i] = v.data_[i];}
+	uint_cl(uchar v) {for (int i = 0; i < L / 8; ++i) data_[i] = (i == 0 ? v : 0);}
+	uint_cl(char v) {for (int i = 0; i < L / 8; ++i) data_[i] = (i == 0 ? v : 0);}
+	uint_cl(ushort v) {int l = piMin<uint>(L / 8, sizeof(v)); memcpy(data_, &v, l); for (int i = l; i < L / 8; ++i) data_[i] = 0;}
+	uint_cl(short v) {int l = piMin<uint>(L / 8, sizeof(v)); memcpy(data_, &v, l); for (int i = l; i < L / 8; ++i) data_[i] = 0;}
+	uint_cl(uint v) {int l = piMin<uint>(L / 8, sizeof(v)); memcpy(data_, &v, l); for (int i = l; i < L / 8; ++i) data_[i] = 0;}
+	uint_cl(int v) {int l = piMin<uint>(L / 8, sizeof(v)); memcpy(data_, &v, l); for (int i = l; i < L / 8; ++i) data_[i] = 0;}
+	uint_cl(ulong v) {int l = piMin<uint>(L / 8, sizeof(v)); memcpy(data_, &v, l); for (int i = l; i < L / 8; ++i) data_[i] = 0;}
+	uint_cl(long v) {int l = piMin<uint>(L / 8, sizeof(v)); memcpy(data_, &v, l); for (int i = l; i < L / 8; ++i) data_[i] = 0;}
+	uint_cl(ullong v) {int l = piMin<uint>(L / 8, sizeof(v)); memcpy(data_, &v, l); for (int i = l; i < L / 8; ++i) data_[i] = 0;}
+	uint_cl(llong v) {int l = piMin<uint>(L / 8, sizeof(v)); memcpy(data_, &v, l); for (int i = l; i < L / 8; ++i) data_[i] = 0;}
+
+	operator bool() {for (int i = 0; i < L / 8; ++i) if (data_[i] > 0) return true; return false;}
+	operator char() {return (char)data_[0];}
+	operator short() {short t(0); int l = piMin<uint>(L / 8, sizeof(t)); memcpy(&t, data_, l); return t;}
+	operator int() {int t(0); int l = piMin<uint>(L / 8, sizeof(t)); memcpy(&t, data_, l); return t;}
+	operator long() {long t(0); int l = piMin<uint>(L / 8, sizeof(t)); memcpy(&t, data_, l); return t;}
+	operator llong() {llong t(0); int l = piMin<uint>(L / 8, sizeof(t)); memcpy(&t, data_, l); return t;}
+	operator uchar() {return data_[0];}
+	operator ushort() {ushort t(0); int l = piMin<uint>(L / 8, sizeof(t)); memcpy(&t, data_, l); return t;}
+	operator uint() {uint t(0); int l = piMin<uint>(L / 8, sizeof(t)); memcpy(&t, data_, l); return t;}
+	operator ulong() {ulong t(0); int l = piMin<uint>(L / 8, sizeof(t)); memcpy(&t, data_, l); return t;}
+	operator ullong() {ullong t(0); int l = piMin<uint>(L / 8, sizeof(t)); memcpy(&t, data_, l); return t;}
+
+	uint_cl<L> operator +(const uint_cl<L> & v) {
+		uint_cl<L> t;
+		uint cv;
+		bool ov = false;
+		for (int i = 0; i < L / 8; ++i) {
+			cv = v.data_[i] + data_[i];
+			if (ov) ++cv;
+			ov = cv > 255;
+			t.data_[i] = ov ? cv - 256 : cv;
+		}
+		return t;
+	}
+
+	uint_cl<L> operator &(const uint_cl<L> & v) const {uint_cl<L> t; for (int i = 0; i < L / 8; ++i) t.data_[i] = v.data_[i] & data_[i]; return t;}
+	uint_cl<L> operator &(const uchar & v) const {return *this & uint_cl<L>(v);}
+	uint_cl<L> operator &(const ushort & v) const {return *this & uint_cl<L>(v);}
+	uint_cl<L> operator &(const uint & v) const {return *this & uint_cl<L>(v);}
+	uint_cl<L> operator &(const ulong & v) const {return *this & uint_cl<L>(v);}
+	uint_cl<L> operator &(const ullong & v) const {return *this & uint_cl<L>(v);}
+	uint_cl<L> operator &(const char & v) const {return *this & uint_cl<L>(v);}
+	uint_cl<L> operator &(const short & v) const {return *this & uint_cl<L>(v);}
+	uint_cl<L> operator &(const int & v) const {return *this & uint_cl<L>(v);}
+	uint_cl<L> operator &(const long & v) const {return *this & uint_cl<L>(v);}
+	uint_cl<L> operator &(const llong & v) const {return *this & uint_cl<L>(v);}
+
+	uint_cl<L> operator |(const uint_cl<L> & v) const {uint_cl<L> t; for (int i = 0; i < L / 8; ++i) t.data_[i] = v.data_[i] | data_[i]; return t;}
+	uint_cl<L> operator |(const uchar & v) const {return *this | uint_cl<L>(v);}
+	uint_cl<L> operator |(const ushort & v) const {return *this | uint_cl<L>(v);}
+	uint_cl<L> operator |(const uint & v) const {return *this | uint_cl<L>(v);}
+	uint_cl<L> operator |(const ulong & v) const {return *this | uint_cl<L>(v);}
+	uint_cl<L> operator |(const ullong & v) const {return *this | uint_cl<L>(v);}
+	uint_cl<L> operator |(const char & v) const {return *this | uint_cl<L>(v);}
+	uint_cl<L> operator |(const short & v) const {return *this | uint_cl<L>(v);}
+	uint_cl<L> operator |(const int & v) const {return *this | uint_cl<L>(v);}
+	uint_cl<L> operator |(const long & v) const {return *this | uint_cl<L>(v);}
+	uint_cl<L> operator |(const llong & v) const {return *this | uint_cl<L>(v);}
+
+	uint_cl<L> operator ^(const uint_cl<L> & v) const {uint_cl<L> t; for (int i = 0; i < L / 8; ++i) t.data_[i] = v.data_[i] ^ data_[i]; return t;}
+	uint_cl<L> operator ^(const uchar & v) const {return *this ^ uint_cl<L>(v);}
+	uint_cl<L> operator ^(const ushort & v) const {return *this ^ uint_cl<L>(v);}
+	uint_cl<L> operator ^(const uint & v) const {return *this ^ uint_cl<L>(v);}
+	uint_cl<L> operator ^(const ulong & v) const {return *this ^ uint_cl<L>(v);}
+	uint_cl<L> operator ^(const ullong & v) const {return *this ^ uint_cl<L>(v);}
+	uint_cl<L> operator ^(const char & v) const {return *this ^ uint_cl<L>(v);}
+	uint_cl<L> operator ^(const short & v) const {return *this ^ uint_cl<L>(v);}
+	uint_cl<L> operator ^(const int & v) const {return *this ^ uint_cl<L>(v);}
+	uint_cl<L> operator ^(const long & v) const {return *this ^ uint_cl<L>(v);}
+	uint_cl<L> operator ^(const llong & v) const {return *this ^ uint_cl<L>(v);}
+
+	bool operator <(const uint_cl<L> & v) const {for (int i = 0; i < L / 8; ++i) {if (v.data_[i] > data_[i]) return true; if (v.data_[i] < data_[i]) return false;} return false;}
+	bool operator <=(const uint_cl<L> & v) const {for (int i = 0; i < L / 8; ++i) {if (v.data_[i] > data_[i]) return true; if (v.data_[i] < data_[i]) return false;} return true;}
+	bool operator >(const uint_cl<L> & v) const {for (int i = 0; i < L / 8; ++i) {if (v.data_[i] < data_[i]) return true; if (v.data_[i] > data_[i]) return false;} return false;}
+	bool operator >=(const uint_cl<L> & v) const {for (int i = 0; i < L / 8; ++i) {if (v.data_[i] < data_[i]) return true; if (v.data_[i] > data_[i]) return false;} return true;}
+	bool operator ==(const uint_cl<L> & v) const {for (int i = 0; i < L / 8; ++i) if (v.data_[i] != data_[i]) return false; return true;}
+	bool operator !=(const uint_cl<L> & v) const {for (int i = 0; i < L / 8; ++i) if (v.data_[i] != data_[i]) return true; return false;}
+	bool operator <=(const uint_cl<8> & v1) {return (*(uchar*)data()) <= (*(uchar*)v1.data());}
+
+	uint_cl<L> operator >>(const int & c) const {
+		uint_cl<L> t;
+		int l = L - c;
+		bool bit;
+		if (l <= 0) return t;
+		for (int i = 0; i < l; ++i) {
+			bit = 1 & (data_[(i + c) / 8] >> ((i + c) % 8));
+			if (bit) t.data_[i / 8] |= (1 << (i % 8));
+			else t.data_[i / 8] &= ~(1 << (i % 8));
+		}
+		return t;
+	}
+	uint_cl<L> operator >>(const uint & c) const {return (*this << (int)c);}
+	uint_cl<L> operator <<(const int & c) const {
+		uint_cl<L> t;
+		int l = L - c;
+		bool bit;
+		if (l <= 0) return t;
+		for (int i = c; i < L; ++i) {
+			bit = 1 & (data_[(i - c) / 8] >> ((i - c) % 8));
+			if (bit) t.data_[i / 8] |= (1 << (i % 8));
+			else t.data_[i / 8] &= ~(1 << (i % 8));
+		}
+		return t;
+	}
+	uint_cl<L> operator <<(const uint & c) const {return (*this >> (int)c);}
+
+	uint_cl<L> & inverse() const {for (int i = 0; i < L / 8; ++i) data_[i] = ~data_[i]; return *this;}
+	uint_cl<L> inversed() const {uint_cl<L> t(*this); for (int i = 0; i < L / 8; ++i) t.data_[i] = ~t.data_[i]; return t;}
+	uint_cl<L> reversed() const {
+		uint_cl<L> t;
+		bool bit;
+		for (int i = 0; i < L; ++i) {
+			bit = 1 & (data_[(L - i - 1) / 8] >> ((L - i - 1) % 8));
+			if (bit) t.data_[i / 8] |= (1 << (i % 8));
+			else t.data_[i / 8] &= ~(1 << (i % 8));
+		}
+		return t;
+	}
+
+	const uchar * data() const {return data_;}
+	uchar * data() {return data_;}
+	uint length() const {return L / 8;}
+
+private:
+	uchar data_[L / 8];
+
+};
+
+template <uint L>
+inline std::ostream & operator <<(std::ostream & s, const uint_cl<L> & v) {std::ios::fmtflags f = s.flags(); s << std::hex; for (uint i = 0; i < v.length(); ++i) {s << int(v.data()[i]); if (v.data()[i] < 0x10) s << '0'; s << ' ';} s.flags(f); return s;}
+
+inline uchar reverseByte(uchar b) {
+	uchar ret = 0;
+	bool bit;
+	for (int i = 0; i < 8; ++i) {
+		bit = 1 & (b >> (7 - i));
+		if (bit) ret |= (1 << i);
+	}
+	return ret;
+}
+
+template <uint L, typename N = uint_cl<L> >
+class PIP_EXPORT PICRC {
+public:
+	PICRC(const N & poly = N()) {poly_ = poly; reverse_poly = true; init_ = inversed(N(0)); out_ = inversed(N(0)); reverse_before_xor = reverse_data = false; initTable();}
+	PICRC(const N & poly, bool reverse_poly_, const N & initial, const N & out_xor) {poly_ = poly; reverse_poly = reverse_poly_; init_ = initial; out_ = out_xor; reverse_before_xor = reverse_data = false; initTable();}
+
+	void setInitial(const N & v) {init_ = v;}
+	void setOutXor(const N & v) {out_ = v;}
+	void setReversePolynome(bool yes) {reverse_poly = yes; initTable();}
+	void setReverseOutBeforeXOR(bool yes) {reverse_before_xor = yes;}
+	void setReverseDataBytes(bool yes) {reverse_data = yes;}
+
+	void initTable() {
+		N tmp, pol = reverse_poly ? reversed(poly_) : poly_;
+		//cout << std::hex << "poly " << (uint)N(poly_) << " -> " << (uint)N(pol) << endl;
+		for (int i = 0; i < 256; ++i) {
+			tmp = uchar(i);
+			for (int j = 0; j < 8; ++j)
+				tmp = ((tmp & 1) ? ((tmp >> 1) ^ pol) : (tmp >> 1));
+			table[i] = tmp;
+		}
+
+	}
+
+	N calculate(const void * data, int size) {
+		N crc = init_;
+		uchar * data_ = (uchar * )data, cb;
+		//cout << "process " << size << endl;
+		uchar nTemp;
+		for (int i = 0; i < size; ++i) {
+			cb = data_[i];
+			if (reverse_data) cb = reverseByte(cb);
+			nTemp = cb ^ uchar(crc);
+			crc = crc >> 8;
+			crc = crc ^ table[nTemp];
+		}
+		if (reverse_before_xor) crc = reversed(crc);
+		return crc ^ out_;
+
+	}
+	N calculate(const PIByteArray & d) {return calculate(d.data(), d.size());}
+	N calculate(const char * str) {string s(str); return calculate((void * )s.data(), s.size());}
+
+private:
+	inline N reversed(const N & v) {return v.reversed();}
+	inline N inversed(const N & v) {return v.inversed();}
+	
+	N table[256];
+	N poly_, init_, out_;
+	bool reverse_poly, reverse_before_xor, reverse_data;
+
+};
+
+template <> inline uchar  PICRC<8, uchar>::reversed(const uchar & v) {return reverseByte(v);}
+template <> inline ushort PICRC<16, ushort>::reversed(const ushort & v) {return uint_cl<16>(v).reversed();}
+template <> inline uint   PICRC<32, uint>::reversed(const uint & v) {return uint_cl<32>(v).reversed();}
+template <> inline uchar  PICRC<8, uchar>::inversed(const uchar & v) {return ~v;}
+template <> inline ushort PICRC<16, ushort>::inversed(const ushort & v) {return ~v;}
+template <> inline uint   PICRC<32, uint>::inversed(const uint & v) {return ~v;}
+
+typedef PICRC<32, uint>   CRC_32;
+typedef PICRC<24>         CRC_24;
+typedef PICRC<16, ushort> CRC_16;
+typedef PICRC<8, uchar>   CRC_8;
+
+inline CRC_32 standardCRC_32() {return CRC_32(0x04C11DB7, true, 0xFFFFFFFF, 0xFFFFFFFF);}
+inline CRC_16 standardCRC_16() {return CRC_16(0x8005, true, 0x0, 0x0);}
+inline CRC_8 standardCRC_8() {return CRC_8(0xD5, true, 0x0, 0x0);}
+
+#endif // CRC_H

src/math/pievaluator.cpp

@@ -0,0 +1,1254 @@
+/*
+    PIP - Platform Independent Primitives
+    Evaluator designed for stream computing
+    Copyright (C) 2014  Ivan Pelipenko peri4ko@gmail.com
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "pievaluator.h"
+
+
+/*! \class PIEvaluator
+ * \brief This class provide mathematical evaluations of custom expression
+ * 
+ * \section PIEvaluator_sec0 Synopsis
+ * %PIEvaluator developed for stream evaluations of once set expression.
+ * It`s create internal list of instructions on function \a check() and
+ * executes very fast on function \a evaluate(). Once given expression
+ * can be evaluated any times with different variable values. Evaluator
+ * supports many common mathematic functions described below. Also it`s
+ * automatic puts unnecessarily signs and bracets. Processed expression
+ * can be obtains with function \a expression(). If there is an error
+ * in expression you can get it with function \a error(). Last evaluated
+ * result you can get with function \a lastResult().
+ * \section PIEvaluator_sec1 Using
+ * First you should set your variables with function \a setVariable().
+ * Next give your expression with function \a check() and check for error
+ * with functions \a isCorrect() and \a error(). If expression is correct
+ * you can get processed expression with function \a expression() and
+ * evaluate it with function \a evaluate(). You can change variable values
+ * without rechecking expression.
+ * 
+ * \section PIEvaluator_sec2 Functions
+ * %PIEvaluator supports arithmetical operations with complex numbers, this
+ * is their list in priority order:
+ * * ^ (power)
+ * * * (multiply)
+ * * / (divide)
+ * * % (residue)
+ * * + (add)
+ * * - (subtract)
+ * 
+ * In addition there are compare and logical operations:
+ * * == (equal)
+ * * != (not equal)
+ * * > (greater)
+ * * < (smaller)
+ * * >= (greater or equal)
+ * * <= (smaller or equal)
+ * * && (and)
+ * * || (or)
+ * 
+ * Compare and logical functions works with real operators part and returns 0 or 1.
+ * 
+ * Mathematical functions:
+ * * sin(x) - sine
+ * * cos(x) - cosine
+ * * tg(x) - tangent
+ * * ctg(x) - cotangent
+ * * arcsin(x) - arcsine
+ * * arccos(x) - arccosine
+ * * arctg(x) -arccotangent
+ * * arcctg(x) - arctangent
+ * * sh(x) - hyperbolical sine
+ * * ch(x) - hyperbolical cosine
+ * * th(x) - hyperbolical tangent
+ * * cth(x) - hyperbolical cotangent
+ * * sqr(x) - square
+ * * sqrt(x) - square root
+ * * abs(x) - absolute value
+ * * sign(x) - sign of real part (-1 or 1)
+ * * exp(x) - exponent
+ * * pow(x, p) - x in power p
+ * * ln(x) - natural logarithm
+ * * lg(x) - decimal logarithm
+ * * log(x, b) - logarithm of x with base b
+ * * im(x) - imaginary part of complex number
+ * * re(x) - real part of complex number
+ * * arg(x) - argument of complex number
+ * * len(x) - length of complex number
+ * * conj(x) - length of complex number
+ * * rad(x) - convert degrees to radians
+ * * deg(x) - convert radians to degrees
+ * * j0(x) - Bessel function first kind order 0
+ * * j1(x) - Bessel function first kind order 1
+ * * jn(x, n) - Bessel function first kind order n
+ * * y0(x) - Bessel function second kind order 0
+ * * y1(x) - Bessel function second kind order 1
+ * * yn(x, n) - Bessel function second kind order n
+ * * random(s, f) - regular random number in range [s, f]
+ * * min(x0, x1, ...) - minimum of x0, x1, ...
+ * * max(x0, x1, ...) - maximum of x0, x1, ...
+ * * clamp(x, a, b) - trim x on range [a, b]
+ * * step(x, s) - 0 if x < s, else 1
+ * * mix(x, a, b) - interpolate between a and b linear for x (a * (1 - x) + b * x)
+ * 
+ * There are some built-in constans:
+ * * i (imaginary 1)
+ * * e
+ * * pi
+ * 
+ * All trigonometric functions takes angle in radians.
+ * 
+ * \section PIEvaluator_sec3 Example
+ * \snippet pievaluator.cpp main
+ */
+
+
+PIEvaluatorContent::PIEvaluatorContent() {
+	addFunction("arcsin", 1);
+	addFunction("arccos", 1);
+	addFunction("arctg", 1);
+	addFunction("arcctg", 1);
+	addFunction("random", 2);
+	addFunction("sin", 1);
+	addFunction("cos", 1);
+	addFunction("ctg", 1);
+	addFunction("tg", 1);
+	addFunction("exp", 1);
+	addFunction("cth", 1);
+	addFunction("sh", 1);
+	addFunction("ch", 1);
+	addFunction("th", 1);
+	addFunction("sqrt", 1);
+	addFunction("sqr", 1);
+	addFunction("pow", 2);
+	addFunction("abs", 1);
+	addFunction("ln", 1);
+	addFunction("lg", 1);
+	addFunction("log", 2);
+	addFunction("im", 1);
+	addFunction("re", 1);
+	addFunction("arg", 1);
+	addFunction("len", 1);
+	addFunction("conj", 1);
+	addFunction("sign", 1);
+	addFunction("rad", 1);
+	addFunction("deg", 1);
+	addFunction("j0", 1);
+	addFunction("j1", 1);
+	addFunction("jn", 2);
+	addFunction("y0", 1);
+	addFunction("y1", 1);
+	addFunction("yn", 2);
+	addFunction("min", -2); // (x0,x1,...)
+	addFunction("max", -2); // (x0,x1,...)
+	addFunction("clamp", 3); // (x,a,b) =  x < a ? a : (x > b ? b : x)
+	addFunction("step", 2); // (x,s) =     x >= s ? 1. : 0.  (1 if 'x' >= 's', else 0)
+	addFunction("mix", 3); // (x,a,b) =    a*(1.-x) + b*x    (interpolate between 'a' and 'b' linear for 'x')
+	addFunction("defined", 1);
+	clearCustomVariables();
+	//addVariable("n", 0.);
+	//addVariable("x1", 123);
+}
+
+
+bool PIEvaluatorContent::setVariableValue(int index, complexd new_value) {
+	if (index < 0 || index >= variables.size_s()) return false;
+	variables[index].value = new_value;
+	return true;
+}
+
+
+bool PIEvaluatorContent::setVariableName(int index, const PIString & new_name) {
+	if (index < 0 || index >= variables.size_s()) return false;
+	variables[index].name = new_name;
+	return true;
+}
+
+
+void PIEvaluatorContent::clearCustomVariables() {
+	variables.clear();
+	addVariable("i", complexd_i);
+	addVariable("pi", atan(1.) * 4.);
+	addVariable("e", exp(1.));
+	cv_count = variables.size();
+}
+
+
+void PIEvaluatorContent::sortVariables() {
+	//PIEvaluatorTypes::Variable tv;
+	for (uint i = 0; i < variables.size(); i++) {
+		for (uint j = variables.size() - 1; j > i; j--) {
+			if (variables[j].name.length() <= variables[i].name.length()) continue;
+			piSwap<PIEvaluatorTypes::Variable>(variables[i], variables[j]);
+		}
+	}
+	/*
+	 *	qDebug() << "---";
+	 *	for (int i = 0; i < variables.size(); i++) {
+	 *		qDebug() << variables[i].name;
+	 }
+	 */
+}
+
+
+PIEvaluatorTypes::BaseFunctions PIEvaluatorContent::getBaseFunction(const PIString & name) {
+	if (name == "sin") return PIEvaluatorTypes::bfSin;
+	if (name == "cos") return PIEvaluatorTypes::bfCos;
+	if (name == "tg") return PIEvaluatorTypes::bfTg;
+	if (name == "ctg") return PIEvaluatorTypes::bfCtg;
+	if (name == "arcsin") return PIEvaluatorTypes::bfArcsin;
+	if (name == "arccos") return PIEvaluatorTypes::bfArccos;
+	if (name == "arctg") return PIEvaluatorTypes::bfArctg;
+	if (name == "arcctg") return PIEvaluatorTypes::bfArcctg;
+	if (name == "exp") return PIEvaluatorTypes::bfExp;
+	if (name == "random") return PIEvaluatorTypes::bfRandom;
+	if (name == "sh") return PIEvaluatorTypes::bfSh;
+	if (name == "ch") return PIEvaluatorTypes::bfCh;
+	if (name == "th") return PIEvaluatorTypes::bfTh;
+	if (name == "cth") return PIEvaluatorTypes::bfCth;
+	if (name == "sqrt") return PIEvaluatorTypes::bfSqrt;
+	if (name == "sqr") return PIEvaluatorTypes::bfSqr;
+	if (name == "pow") return PIEvaluatorTypes::bfPow;
+	if (name == "abs") return PIEvaluatorTypes::bfAbs;
+	if (name == "ln") return PIEvaluatorTypes::bfLn;
+	if (name == "lg") return PIEvaluatorTypes::bfLg;
+	if (name == "log") return PIEvaluatorTypes::bfLog;
+	if (name == "im") return PIEvaluatorTypes::bfIm;
+	if (name == "re") return PIEvaluatorTypes::bfRe;
+	if (name == "arg") return PIEvaluatorTypes::bfArg;
+	if (name == "len") return PIEvaluatorTypes::bfLen;
+	if (name == "conj") return PIEvaluatorTypes::bfConj;
+	if (name == "sign") return PIEvaluatorTypes::bfSign;
+	if (name == "rad") return PIEvaluatorTypes::bfRad;
+	if (name == "deg") return PIEvaluatorTypes::bfDeg;
+	if (name == "j0") return PIEvaluatorTypes::bfJ0;
+	if (name == "j1") return PIEvaluatorTypes::bfJ1;
+	if (name == "jn") return PIEvaluatorTypes::bfJN;
+	if (name == "y0") return PIEvaluatorTypes::bfY0;
+	if (name == "y1") return PIEvaluatorTypes::bfY1;
+	if (name == "yn") return PIEvaluatorTypes::bfYN;
+	if (name == "min") return PIEvaluatorTypes::bfMin;
+	if (name == "max") return PIEvaluatorTypes::bfMax;
+	if (name == "clamp") return PIEvaluatorTypes::bfClamp;
+	if (name == "step") return PIEvaluatorTypes::bfStep;
+	if (name == "mix") return PIEvaluatorTypes::bfMix;
+	if (name == "defined") return PIEvaluatorTypes::bfDefined;
+	return PIEvaluatorTypes::bfUnknown;
+}
+
+const PIString & PIEvaluator::prepare(const PIString & string) {
+	currentString = string.trimmed();
+	if (currentString.isEmpty()) currentString = "0";
+	replaceOperators();
+	removeSpaces();
+	checkBrackets();
+	while (fillElements()) checkBrackets();
+	while (setSignes()) fillElements();
+	removeJunk();
+	findUnknownVariables();
+	return currentString;
+}
+
+
+void PIEvaluator::removeSpaces() {
+	PIString tmps = currentString;
+	for (int i = 0; i < tmps.length(); i++) {
+		if (tmps[i] == ' ' || tmps[i] == '\t') {
+			tmps.remove(i, 1);
+			i--;
+		}
+	}
+	currentString = tmps;
+}
+
+
+void PIEvaluator::removeJunk() {
+	PIChar cc;
+	bool junk = true;
+	int bcnt;
+	while (junk) {
+		if (currentString.left(1) != "(" || currentString.right(1) != ")") return;
+		bcnt = 1;
+		junk = false;
+		for (int i = 1; i < currentString.length(); i++) {
+			cc = currentString[i];
+			if (cc == '(') bcnt++;
+			if (cc == ')') bcnt--;
+			if (bcnt == 0) {
+				if (i == currentString.length() - 1) {
+					currentString = currentString.mid(1, currentString.length() - 2);
+					elements.pop_front();
+					elements.pop_back();
+					junk = true;
+					break;
+				} else break;
+			}
+		}
+	}
+}
+
+
+void PIEvaluator::replaceOperators() {
+	currentString.replaceAll("==", "=");
+	currentString.replaceAll("!=", ":");
+	currentString.replaceAll(">=", "}");
+	currentString.replaceAll("<=", "{");
+	currentString.replaceAll("&&", "&");
+	currentString.replaceAll("||", "|");
+}
+
+
+void PIEvaluator::makeOutput(PIString & string) {
+	string.replaceAll(":", "≠");
+	string.replaceAll("}", "≥");
+	string.replaceAll("{", "≤");
+	string.replaceAll("&", "⋀");
+	string.replaceAll("|", "⋁");
+}
+
+
+void PIEvaluator::findUnknownVariables() {
+	PIString cvar;
+	unknownVars.clear();
+	for (int i = 0; i < currentString.length(); i++) {
+		if (elements[i].var_num == -666) cvar += currentString[i];
+		else {
+			if (cvar.length() == 0) continue;
+			unknownVars << cvar;
+			cvar = "";
+		}
+	}
+	if (cvar.length() > 0) unknownVars << cvar;
+	unknownVars.removeDuplicates();
+}
+
+
+bool PIEvaluator::isSign(const PIChar & ch) {
+	return ch == '+' || ch == '-' ||
+		   ch == '*' || ch == '/' ||
+		   ch == '%' || ch == '^' ||
+		   ch == '=' || ch == ':' ||
+		   ch == '>' || ch == '<' ||
+		   ch == '}' || ch == '{' ||
+		   ch == '&' || ch == '|';
+}
+
+
+void PIEvaluator::checkBrackets() {
+	PIString tmps = currentString;
+	PIChar fc, sc;
+	int bcnt = 0, bpos = 0, inserted = 0;
+	currentString = tmps;
+	for (int i = 0; i < tmps.length(); i++) {
+		if (tmps[i] == '(') {
+			if (bcnt == 0) bpos = i;
+			bcnt++;
+		}
+		if (tmps[i] == ')') {
+			if (bcnt == 0) {
+				currentString.insert(bpos + inserted, "(");
+				inserted++;
+			} else bcnt--;
+		}
+	}
+	if (bcnt > 0) currentString += PIString(bcnt, ')');
+	tmps = currentString;
+	for (int i = 0; i < tmps.length() - 1; i++) {
+		fc = tmps[i].toLower();
+		sc = tmps[i + 1].toLower();
+		if ((fc == ')' && sc == '(') ||
+			(fc == ')' && sc >= '0' && sc <= '9') ||
+			(fc == ')' && sc >= 'a' && sc <= 'z') ) tmps.insert(++i, '*');
+	}
+	currentString = tmps;
+}
+
+
+bool PIEvaluator::fillElements() {
+	int fstart, flen, cnum = 0, cpart = 0, cfunc;
+	PIChar cc, nc, pc, fc = '!';
+	bool numFound = false;
+	PIString curfind, tmps = currentString;
+	elements.resize(tmps.length());
+	for (uint i = 0; i < elements.size(); i++) {
+		elements[i].type = PIEvaluatorTypes::etVariable;
+		elements[i].var_num = -666;
+	}
+	currentVariables.clear();
+	//qDebug().nospace() << "search for functions ...";
+	for (int i = 0; i < content.functionsCount(); i++) {
+		curfind = content.function(i).identifier;
+		cfunc = i; //(int)content.function(i).type;
+		flen = curfind.length();
+		fstart = 0;
+		while (fstart >= 0) {
+			fstart = tmps.find(curfind, fstart);
+			if (fstart < 0) break;
+			if (tmps[fstart + flen] != '(') {
+				//currentString.insert(fstart + flen, "(");
+				//return true;
+				fstart++;
+				continue;
+			}
+			for (int j = fstart; j < fstart + flen; j++) {
+				elements[j].set(PIEvaluatorTypes::etFunction, cnum, cfunc);
+				tmps.replace(j, 1, fc);
+			}
+			cnum++;
+		}
+	}
+	cnum = 0;
+	//qDebug().nospace() << "search for variables ...";
+	for (int i = 0; i < content.variablesCount(); i++) {
+		curfind = content.variable(i).name;
+		flen = curfind.length();
+		fstart = 0;
+		while (fstart >= 0) {
+			fstart = tmps.find(curfind, fstart);
+			if (fstart < 0) break;
+			for (int j = fstart; j < fstart + flen; j++) {
+				elements[j].set(PIEvaluatorTypes::etVariable, cnum, i);
+				tmps.replace(j, 1, fc);
+			}
+			cnum++;
+		}
+	}
+	curfind = "";
+	cnum = 1;
+	//qDebug().nospace() << "search for numbers ...";
+	for (int i = 0; i < tmps.length(); i++) {
+		cc = tmps[i];
+		/*if (cc == " " || cc == "(" || cc == ")") {
+			curfind = "";
+			cpart = 0;
+			numFound = false;
+			continue;
+		}*/
+		switch (cpart) {
+		case 0:
+			if ((cc >= '0' && cc <= '9')) {// || cc == '-' || cc == '+') {
+				curfind += cc;
+				cpart = 1;
+				continue;
+			}
+			if (cc == '.') {
+				curfind += cc;
+				cpart = 2;
+				continue;
+			}
+			if (cc == 'E') {
+				curfind += cc;
+				cpart = 3;
+				continue;
+			}
+			break;
+		case 1:
+			if (cc >= '0' && cc <= '9') {
+				curfind += cc;
+				continue;
+			}
+			if (cc == '.') {
+				curfind += cc;
+				cpart = 2;
+				continue;
+			}
+			if (cc == 'E') {
+				curfind += cc;
+				cpart = 3;
+				continue;
+			}
+			numFound = true;
+			break;
+		case 2:
+			if (cc >= '0' && cc <= '9') {
+				curfind += cc;
+				continue;
+			}
+			if (cc == 'E') {
+				curfind += cc;
+				cpart = 3;
+				continue;
+			}
+			numFound = true;
+			break;
+		case 3:
+			if ((cc >= '0' && cc <= '9') || cc == '-' || cc == '+') {
+				curfind += cc;
+				cpart = 4;
+				continue;
+			}
+			numFound = true;
+			break;
+		case 4:
+			if (cc >= '0' && cc <= '9') {
+				curfind += cc;
+				continue;
+			}
+			numFound = true;
+			break;
+		}
+		if (numFound) {
+			//qDebug().nospace() << "add " << cnum << ": " << curfind << " = " << curfind.toDouble();
+			currentVariables.push_back(PIEvaluatorTypes::Variable("tmp" + PIString::fromNumber(cnum), curfind.toDouble()));
+			for (int j = i - curfind.length(); j < i; j++) {
+				elements[j].set(PIEvaluatorTypes::etNumber, cnum, -cnum);
+				tmps.replace(j, 1, fc);
+			}
+			curfind = "";
+			cnum++;
+			cpart = 0;
+			numFound = false;
+		}
+	}
+	if (cpart > 0) {
+		//qDebug().nospace() << "add " << cnum << ": " << curfind << " = " << curfind.toDouble();
+		currentVariables.push_back(PIEvaluatorTypes::Variable("tmp" + PIString::fromNumber(cnum), curfind.toDouble()));
+		for (int j = tmps.length() - curfind.length(); j < tmps.length(); j++) {
+			elements[j].set(PIEvaluatorTypes::etNumber, cnum, -cnum);
+			tmps.replace(j, 1, fc);
+		}
+	}
+	cc = nc = fc;
+	//qDebug().nospace() << "search for signes ...";
+	for (int i = 0; i < tmps.length(); i++) {
+		cc = tmps[i];
+		if (i > 0) pc = tmps[i - 1];
+		else pc = fc;
+		if (i < tmps.length() - 1) nc = tmps[i + 1];
+		else nc = fc;
+		if (cc == '(' || cc == ')' || cc == ',') {
+			elements[i].set(PIEvaluatorTypes::etOperator, -1);
+			continue;
+		}
+		if (cc == '-' || cc == '+') {
+			elements[i].set(PIEvaluatorTypes::etOperator, -1);
+			if (i < tmps.length() - 1) if (elements[i + 1].type == PIEvaluatorTypes::etVariable ||
+										   elements[i + 1].type == PIEvaluatorTypes::etFunction) continue;
+			if ((pc == '(' || isSign(pc) || i == 0) && i < tmps.length() - 1) {
+				if (elements[i + 1].type != PIEvaluatorTypes::etOperator) {
+					cnum = elements[i + 1].num;
+					elements[i].set(PIEvaluatorTypes::etNumber, cnum);
+					tmps.replace(i, 1, fc);
+					///cout << "found sign " << cc << " :" << cnum - 1 << endl;
+					if (cc == '-' && currentVariables.size_s() >= cnum)
+						currentVariables[cnum - 1].value = -currentVariables[cnum - 1].value;
+					//i++;
+					continue;
+				}
+			}
+		}
+		if (isSign(cc)) {
+			elements[i].set(PIEvaluatorTypes::etOperator, -1);
+			continue;
+		}
+	}
+	/*
+	qDebug().nospace() << tmps;
+	cout << " ";
+	for (int i = 0; i < elements.size(); i++) {
+		switch (elements[i].type) {
+		case etFunction: cout << "f"; break;
+		case etNumber: cout << "n"; break;
+		case etOperator: cout << "o"; break;
+		case etVariable: cout << "v"; break;
+		}
+	}
+	cout << endl;
+	*/
+	return false;
+	//for (int i = 0; i < currentVariables.size(); i++) qDebug() << "var " << i << ": " << currentVariables[i].value.real();
+}
+
+
+bool PIEvaluator::setSignes() {
+	int inserted = 0, ni, pi = 0, needInsert = 0;
+	PIChar fc, sc, pc;
+	PIString tmps = currentString;
+	for (int i = 0; i < tmps.length() - 1; i++) {
+		needInsert = 0;
+		ni = i + 1;
+		if (i > 0) pi = i - 1;
+		fc = tmps[i].toLower();
+		sc = tmps[ni].toLower();
+		pc = tmps[pi].toLower();
+		//if (elements[i].type == etOperator || elements[ni].type == etVariable) continue;
+		if (fc == ',' || sc == ',') continue;
+		if (elements[i].type == PIEvaluatorTypes::etOperator && elements[ni].type == PIEvaluatorTypes::etOperator) continue;
+		if (fc == ')' && (elements[ni].type == PIEvaluatorTypes::etNumber || elements[ni].type == PIEvaluatorTypes::etVariable || elements[ni].type == PIEvaluatorTypes::etFunction)) needInsert = 1;
+		if (sc == '(' && (elements[i].type == PIEvaluatorTypes::etNumber || elements[i].type == PIEvaluatorTypes::etVariable)) needInsert = 1;
+		if (elements[i].type == PIEvaluatorTypes::etNumber && elements[ni].type == PIEvaluatorTypes::etNumber && elements[i].num != elements[ni].num) needInsert = 1;
+		if (elements[i].type == PIEvaluatorTypes::etVariable && elements[ni].type == PIEvaluatorTypes::etVariable && elements[i].num != elements[ni].num) needInsert = 1;
+		if ((elements[i].type == PIEvaluatorTypes::etNumber && elements[ni].type == PIEvaluatorTypes::etVariable) || (elements[i].type == PIEvaluatorTypes::etVariable && elements[ni].type == PIEvaluatorTypes::etNumber)) needInsert = 1;
+		if ((elements[i].type == PIEvaluatorTypes::etNumber || elements[i].type == PIEvaluatorTypes::etVariable) && elements[ni].type == PIEvaluatorTypes::etFunction) needInsert = 1;
+		if (elements[i].type == PIEvaluatorTypes::etFunction && elements[ni].type == PIEvaluatorTypes::etFunction && elements[i].num != elements[ni].num) needInsert = 2;
+		if (elements[i].type == PIEvaluatorTypes::etFunction && elements[ni].type != PIEvaluatorTypes::etFunction && sc != '(') needInsert = 2;
+		if (elements[pi].type == PIEvaluatorTypes::etOperator && (elements[ni].type == PIEvaluatorTypes::etFunction || elements[ni].type == PIEvaluatorTypes::etVariable) && fc == '-') needInsert = 3;
+		switch (needInsert) {
+		case 1:
+			currentString.insert(ni + inserted, "*");
+			elements.insert(ni + inserted, PIEvaluatorTypes::Element(PIEvaluatorTypes::etOperator, -1));
+			//inserted++;
+			//i++;
+			return true;
+		/*case 2:
+			currentString.insert(ni + inserted, ")");
+			currentString.insert(ni + inserted, "(");
+			elements.insert(ni + inserted, Element(etOperator, -1));
+			elements.insert(ni + inserted, Element(etOperator, -1));
+			inserted++;
+			i++;
+			return true;*/
+		case 3:
+			currentString.insert(ni + inserted, "1*");
+			elements.insert(ni + inserted, PIEvaluatorTypes::Element(PIEvaluatorTypes::etOperator, -1));
+			//inserted;
+			//i++;
+			return true;
+		}
+	}
+	/*if (elements[tmps.length() - 1].type == etFunction) {
+		currentString.insert(tmps.length() + inserted, ")");
+		currentString.insert(tmps.length() + inserted, "(");
+		elements.insert(tmps.length() + inserted, Element(etOperator, -1));
+		elements.insert(tmps.length() + inserted, Element(etOperator, -1));
+		return true;
+	}*/
+	return false;
+}
+
+
+void PIEvaluator::convert() {
+	int j;
+	PIEvaluatorTypes::Element ce, pe;
+	for (int i = 0; i < currentString.length(); i++) {
+		pe = elements[i];
+		if (pe.type != PIEvaluatorTypes::etFunction) continue;
+		j = i + 1;
+		while (j < currentString.length()) {
+			ce = elements[j];
+			if (ce != pe) break;
+			j++;
+		}
+		currentString.replace(i, j - i, " ");
+		for (int k = i + 1; k < j; k++) elements.remove(i);
+		//i++;
+	}
+	for (int i = 0; i < currentString.length(); i++) {
+		pe = elements[i];
+		if (pe.type != PIEvaluatorTypes::etNumber) continue;
+		j = i + 1;
+		while (j < currentString.length()) {
+			ce = elements[j];
+			if (ce != pe) break;
+			j++;
+		}
+		currentString.replace(i, j - i, " ");
+		for (int k = i + 1; k < j; k++) elements.remove(i);
+		//i++;
+	}
+	for (int i = 0; i < currentString.length(); i++) {
+		pe = elements[i];
+		if (pe.type != PIEvaluatorTypes::etVariable) continue;
+		j = i + 1;
+		while (j < currentString.length()) {
+			ce = elements[j];
+			if (ce != pe) break;
+			j++;
+		}
+		currentString.replace(i, j - i, " ");
+		for (int k = i + 1; k < j; k++) elements.remove(i);
+		//i++;
+	}
+	/*qDebug().nospace() << currentString;
+	cout << " ";
+	for (int i = 0; i < elements.size(); i++) {
+		switch (elements[i].type) {
+		case etFunction: cout << "f"; break;
+		case etNumber: cout << "n"; break;
+		case etOperator: cout << "o"; break;
+		case etVariable: cout << "v"; break;
+		}
+	}
+	cout << endl;*/
+}
+
+
+
+const PIString & PIEvaluator::preprocess(const PIString & string) {
+	static PIString ret;
+	int lind;
+	ret = prepare(string);
+	convert();
+	instructions.clear();
+	//qDebug() << preproc->currentString;
+	variables = currentVariables;
+	lind = parse(currentString);
+	if (instructions.size() == 0) {
+		variables.push_back(PIEvaluatorTypes::Variable());
+		instructions.push_back(PIEvaluatorTypes::Instruction(PIEvaluatorTypes::oNone, PIVector<int>(1, lind), -variables.size_s()));
+	}
+	kvars = &(content.variables);
+	/*
+	cout << endl << "variables:" << endl;
+	for (int i = 0; i < variables.size(); i++)
+		cout << i << " value = " << variables[i].value << endl;
+
+	cout << endl << "instructions:" << endl;
+	for (int i = 0; i < instructions.size(); i++) {
+		cout << i << endl;
+		cout << " operation " << instructions[i].operation << endl;
+		cout << " operators: ";
+		for (int j = 0; j < instructions[i].operators.size(); j++)
+			cout << instructions[i].operators[j] << "; ";
+		cout << endl << " function " << instructions[i].function << endl;
+		cout << " out " << instructions[i].out << endl;
+	}
+	*/
+	makeOutput(ret);
+	return ret;
+}
+
+
+PIEvaluatorTypes::Operation PIEvaluator::operationInOrder(const int & index) {
+	switch (index) {
+	case 0: return PIEvaluatorTypes::oPower;
+	case 1: return PIEvaluatorTypes::oMultiply;
+	case 2: return PIEvaluatorTypes::oDivide;
+	case 3: return PIEvaluatorTypes::oResidue;
+	case 4: return PIEvaluatorTypes::oAdd;
+	case 5: return PIEvaluatorTypes::oSubtract;
+	case 6: return PIEvaluatorTypes::oEqual;
+	case 7: return PIEvaluatorTypes::oNotEqual;
+	case 8: return PIEvaluatorTypes::oGreaterEqual;
+	case 9: return PIEvaluatorTypes::oSmallerEqual;
+	case 10: return PIEvaluatorTypes::oGreater;
+	case 11: return PIEvaluatorTypes::oSmaller;
+	case 12: return PIEvaluatorTypes::oAnd;
+	case 13: return PIEvaluatorTypes::oOr;
+	default: return PIEvaluatorTypes::oNone;
+	}
+}
+
+
+int PIEvaluator::parse(const PIString & string, int offset) {
+	int slen = string.length(), /*facnt,*/ farg, bcnt, k;
+	PIChar cc;
+	PIEvaluatorTypes::Element ce;
+	PIEvaluatorTypes::Function cfunc;
+	PIEvaluatorTypes::Operation coper;
+	PIString sbrackets, carg;
+	PIVector<int> args, atmp;
+	PIVector<PIEvaluatorTypes::Operation> opers;
+
+	///qDebug() << "to parse :" + string;
+	///cout << "           "; for (int i = 0; i < slen; i++) cout << preproc->elements[i + offset].type; cout << endl;
+
+	for (int i = 0; i < slen; i++) {
+		ce = elements[i + offset];
+		cc = string[i];
+		switch (ce.type) {
+		case PIEvaluatorTypes::etNumber:
+			args.push_back(ce.var_num);
+			continue;
+		case PIEvaluatorTypes::etVariable:
+			args.push_back(ce.var_num);
+			continue;
+		case PIEvaluatorTypes::etFunction:
+			i++;
+			cfunc = content.function(ce.var_num);
+			//facnt = cfunc.arguments;
+			atmp.clear();
+			bcnt = farg = 1;
+			///qDebug() << "function: " + cfunc.identifier;
+			//for (int k = 0; k < facnt; k++) {
+			bcnt = 1;
+			carg = "";
+			k = i + 1;
+			//if (string.size_s() <= k || k < 0) return -666;
+			while (bcnt > 0) {
+				//if (k < facnt - 1) fcomma = string.indexOf(',', j);
+				cc = string[k];
+				switch (cc.toAscii()) {
+				case '(': bcnt++; break;
+				case ')':
+					bcnt--;
+					if (bcnt == 0) {
+						///qDebug() << "arument: " << carg;
+						atmp.push_back(parse(carg, k + offset - carg.length()));
+						k++;
+						carg = "";
+						if (atmp.size_s() > 0) if (atmp.back() < 0 && farg > 0) farg = atmp.back();
+						continue;
+					}
+					break;
+				case ',':
+					if (bcnt == 1) {
+						///qDebug() << "arument: " << carg;
+						atmp.push_back(parse(carg, k + offset - carg.length()));
+						k++;
+						carg = "";
+						if (atmp.size_s() > 0) if (atmp.back() < 0 && farg > 0) farg = atmp.back();
+						continue;
+					}
+					break;
+				}
+				carg += cc;
+				k++;
+			}
+			i = k - 1;
+			if (farg > 0) {
+				variables.push_back(PIEvaluatorTypes::Variable());
+				farg = -variables.size_s();
+			}
+			instructions.push_back(PIEvaluatorTypes::Instruction(PIEvaluatorTypes::oFunction, atmp, farg, ce.var_num));
+			args.push_back(farg);
+			//for (int i = 0; i < args.size_s(); i++) cout << preproc->currentVariables[-args[i]].value << endl;
+			//i = j + 1;
+			continue;
+		case PIEvaluatorTypes::etOperator:
+			//qDebug() << "operator: " << cc;
+			if (cc == '(') {
+				sbrackets = inBrackets(string.right(slen - i));
+				args.push_back(parse(sbrackets, i + offset + 1));
+				i += sbrackets.length() + 1;
+				continue;
+			}
+			if (cc == '+') {opers.push_back(PIEvaluatorTypes::oAdd); continue;}
+			if (cc == '-') {opers.push_back(PIEvaluatorTypes::oSubtract); continue;}
+			if (cc == '*') {opers.push_back(PIEvaluatorTypes::oMultiply); continue;}
+			if (cc == '/') {opers.push_back(PIEvaluatorTypes::oDivide); continue;}
+			if (cc == '%') {opers.push_back(PIEvaluatorTypes::oResidue); continue;}
+			if (cc == '^') {opers.push_back(PIEvaluatorTypes::oPower); continue;}
+			if (cc == '=') {opers.push_back(PIEvaluatorTypes::oEqual); continue;}
+			if (cc == ':') {opers.push_back(PIEvaluatorTypes::oNotEqual); continue;}
+			if (cc == '}') {opers.push_back(PIEvaluatorTypes::oGreaterEqual); continue;}
+			if (cc == '{') {opers.push_back(PIEvaluatorTypes::oSmallerEqual); continue;}
+			if (cc == '>') {opers.push_back(PIEvaluatorTypes::oGreater); continue;}
+			if (cc == '<') {opers.push_back(PIEvaluatorTypes::oSmaller); continue;}
+			if (cc == '&') {opers.push_back(PIEvaluatorTypes::oAnd); continue;}
+			if (cc == '|') {opers.push_back(PIEvaluatorTypes::oOr); continue;}
+		}
+	}
+	/*
+	cout << "stack: " << endl << "args: ";
+	for (int i = 0; i < args.size_s(); i++) cout << args[i] << ", ";
+	cout << endl << "opers: ";
+	for (int i = 0; i < opers.size_s(); i++) cout << opers[i] << ", ";
+	*/
+	if (opers.size_s() == 0) {
+		if (args.size_s() > 0) return args.back();
+		else return -666;
+	}
+	for (int i = 0; i < PIEvaluatorTypes::operationCount; i++) {
+		coper = operationInOrder(i);
+		for (int j = 0; j < opers.size_s(); j++) {
+			if (coper == PIEvaluatorTypes::oDivide || coper == PIEvaluatorTypes::oMultiply) {
+				if (opers[j] != PIEvaluatorTypes::oDivide && opers[j] != PIEvaluatorTypes::oMultiply) continue;
+			} else {
+				if (opers[j] != coper) continue;
+			}
+			atmp.clear();
+			if (j < args.size_s() && j >= 0) atmp.push_back(args[j]);
+			else atmp.push_back(-666);
+			if (j + 1 < args.size_s() && j >= -1) atmp.push_back(args[j + 1]);
+			else atmp.push_back(-666);
+			farg = 1;
+			if (atmp[0] < 0) farg = atmp[0];
+			else {
+				if (atmp[1] < 0) farg = atmp[1];
+				else {
+					variables.push_back(PIEvaluatorTypes::Variable());
+					farg = -variables.size_s();
+				}
+			}
+			instructions.push_back(PIEvaluatorTypes::Instruction(opers[j], atmp, farg));
+			if (j >= 0 && j < args.size_s()) {
+				args.remove(j);
+				if (j < args.size_s()) args[j] = farg;
+			}
+			opers.remove(j);
+			j--;
+		}
+	}
+	return instructions.back().out;
+	///cout << endl;
+}
+
+
+bool PIEvaluator::check() {
+	PIEvaluatorTypes::Instruction ci;
+	bool error;
+	if (unknownVars.size_s() > 0) {
+		lastError = "Unknown variables: \"" + unknownVars.join("\", \"") + "\"";
+		return false;
+	}
+	for (int i = 0; i < instructions.size_s(); i++) {
+		error = false;
+		ci = instructions[i];
+		PIEvaluatorTypes::Function cf;
+		int fac, gac;
+		switch (ci.operation) {
+		case PIEvaluatorTypes::oNone: break;
+		case PIEvaluatorTypes::oFunction:
+			cf = content.function(ci.function);
+			fac = cf.arguments;
+			gac = ci.operators.size_s();
+			for (int j = 0; j < ci.operators.size_s(); j++) {
+				if (ci.operators[j] == -666) { //(ci.operators[j] < -variables.size_s() || ci.operators[j] >= kvars->size()) {
+					error = true;
+					gac--;
+				}
+			}
+			if (fac > 0) {
+				if (gac != fac) {
+					lastError = "Invalid arguments count for function \"" + cf.identifier +
+								"\", expected " + PIString::fromNumber(fac) + " but " +
+								PIString::fromNumber(gac) + " given";
+					return false;
+				}
+				if (error) {
+					lastError = "Invalid at least one of function \"" + cf.identifier + "\" argument";
+					return false;
+				}
+			}
+			if (fac < 0) {
+				if (gac < -fac) {
+					lastError = "Invalid arguments count for function \"" + cf.identifier +
+								"\", expected at least " + PIString::fromNumber(-fac) + " but " +
+								PIString::fromNumber(gac) + " given";
+					return false;
+				}
+				if (error) {
+					lastError = "Invalid at least one of function \"" + cf.identifier + "\" argument";
+					return false;
+				}
+			}
+			break;
+		default:
+			if (ci.operators[0] == -666 || ci.operators[1] == -666) error = true;
+			if (ci.operators.size_s() != 2 || error) {
+				lastError = "Invalid arguments count for operation \" " + operationChar(ci.operation) + " \"";
+				return false;
+			}
+			break;
+		}
+		if (ci.out < -variables.size_s()) {
+			lastError = "Invalid variable index \"" + PIString::fromNumber(ci.out) + "\"";
+			return false;
+		}
+		for (int j = 0; j < ci.operators.size_s(); j++) {
+			if (ci.operators[j] < -variables.size_s() || ci.operators[j] >= kvars->size_s()) {
+				lastError = "Invalid variable index \"" + PIString::fromNumber(ci.operators[j]) + "\"";
+				return false;
+			}
+		}
+	}
+	return true;
+}
+
+
+PIString PIEvaluator::inBrackets(const PIString & string) {
+	int slen = string.length(), bcnt = 0;
+	PIChar cc;
+	for (int i = 0; i < slen; i++) {
+		cc = string[i];
+		if (cc == '(') bcnt++;
+		if (cc == ')') {
+			bcnt--;
+			if (bcnt == 0) return string.mid(1, i - 1);
+		}
+	}
+	return PIString();
+}
+
+
+PIString PIEvaluator::operationChar(const PIEvaluatorTypes::Operation & operation) {
+	switch (operation) {
+	case PIEvaluatorTypes::oAdd: return "+";
+	case PIEvaluatorTypes::oSubtract: return "-";
+	case PIEvaluatorTypes::oMultiply: return "*";
+	case PIEvaluatorTypes::oDivide: return "/";
+	case PIEvaluatorTypes::oPower: return "^";
+	case PIEvaluatorTypes::oResidue: return "%";
+	case PIEvaluatorTypes::oEqual: return "=";
+	case PIEvaluatorTypes::oNotEqual: return ("≠");
+	case PIEvaluatorTypes::oGreaterEqual: return ("≥");
+	case PIEvaluatorTypes::oSmallerEqual: return ("≤");
+	case PIEvaluatorTypes::oGreater: return ">";
+	case PIEvaluatorTypes::oSmaller: return "<";
+	case PIEvaluatorTypes::oAnd: return ("⋀");
+	case PIEvaluatorTypes::oOr: return ("⋁");
+	default: return "???";
+	}
+}
+
+
+inline complexd PIEvaluator::residue(const complexd & f, const complexd & s) {
+	complexd ret;
+	if (s.real() != 0.) ret = complexd(f.real() - ((int)(f.real() / s.real())) * s.real(), 0.);
+	if (s.imag() != 0.) ret = complexd(ret.real(), f.imag() - ((int)(f.imag() / s.imag())) * s.imag());
+	return ret;
+}
+
+
+inline void PIEvaluator::execFunction(const PIEvaluatorTypes::Instruction & ci) {
+	PIEvaluatorTypes::Function cfunc = content.function(ci.function);
+	int oi = -ci.out - 1;
+	complexd tmp, stmp, ttmp;
+	//qDebug() << "function " << (int)cfunc.type;
+	switch (cfunc.type) {
+	case PIEvaluatorTypes::bfSin:
+		tmpvars[oi].value = sin(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfCos:
+		tmpvars[oi].value = cos(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfTg:
+		tmpvars[oi].value = tan(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfCtg:
+		tmp = tan(value(ci.operators[0]));
+		if (tmp == complexd_0) tmpvars[oi].value = 0.;
+		else tmpvars[oi].value = complexd_1 / tmp;
+		break;
+	case PIEvaluatorTypes::bfArcsin:
+		tmpvars[oi].value = asinc(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfArccos:
+		tmpvars[oi].value = acosc(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfArctg:
+		tmpvars[oi].value = atanc(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfArcctg:
+		tmp = atanc(value(ci.operators[0]));
+		if (tmp == complexd_0) tmpvars[oi].value = 0.;
+		else tmpvars[oi].value = complexd_1 / tmp;
+		break;
+	case PIEvaluatorTypes::bfSh:
+		tmpvars[oi].value = sinh(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfCh:
+		tmpvars[oi].value = cosh(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfTh:
+		tmpvars[oi].value = tanh(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfCth:
+		tmp = tanh(value(ci.operators[0]));
+		if (tmp == complexd_0) tmpvars[oi].value = 0.;
+		else tmpvars[oi].value = complexd_1 / tmp;
+		break;
+	case PIEvaluatorTypes::bfAbs:
+		tmpvars[oi].value = abs(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfSqrt:
+		tmpvars[oi].value = sqrt(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfSqr:
+		tmpvars[oi].value = value(ci.operators[0]) * value(ci.operators[0]);
+		break;
+	case PIEvaluatorTypes::bfExp:
+		tmpvars[oi].value = exp(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfPow:
+		tmpvars[oi].value = pow(value(ci.operators[0]), value(ci.operators[1]));
+		break;
+	case PIEvaluatorTypes::bfLn:
+		tmpvars[oi].value = log(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfLg:
+		tmpvars[oi].value = log10(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfLog:
+		tmp = log(value(ci.operators[1]));
+		if (tmp == complexd_0) tmpvars[oi].value = 0.;
+		else tmpvars[oi].value = log(value(ci.operators[0])) / tmp;
+		break;
+	case PIEvaluatorTypes::bfRe:
+		tmpvars[oi].value = value(ci.operators[0]).real();
+		break;
+	case PIEvaluatorTypes::bfIm:
+		tmpvars[oi].value = value(ci.operators[0]).imag();
+		break;
+	case PIEvaluatorTypes::bfArg:
+		tmpvars[oi].value = arg(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfLen:
+		tmpvars[oi].value = abs(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfConj:
+		tmpvars[oi].value = conj(value(ci.operators[0]));
+		break;
+	case PIEvaluatorTypes::bfSign:
+		tmpvars[oi].value = value(ci.operators[0]).real() >= 0. ? complexd_1 : -complexd_1;
+		break;
+	case PIEvaluatorTypes::bfRad:
+		tmpvars[oi].value = value(ci.operators[0]) * complexd(deg2rad, 0.);
+		break;
+	case PIEvaluatorTypes::bfDeg:
+		tmpvars[oi].value = value(ci.operators[0]) * complexd(rad2deg, 0.);
+		break;
+	case PIEvaluatorTypes::bfJ0:
+		tmpvars[oi].value = piJ0(value(ci.operators[0]).real());
+		break;
+	case PIEvaluatorTypes::bfJ1:
+		tmpvars[oi].value = piJ1(value(ci.operators[0]).real());
+		break;
+	case PIEvaluatorTypes::bfJN:
+		tmpvars[oi].value = piJn(piRoundd(value(ci.operators[1]).real()), value(ci.operators[0]).real());
+		break;
+	case PIEvaluatorTypes::bfY0:
+		tmpvars[oi].value = piY0(value(ci.operators[0]).real());
+		break;
+	case PIEvaluatorTypes::bfY1:
+		tmpvars[oi].value = piY1(value(ci.operators[0]).real());
+		break;
+	case PIEvaluatorTypes::bfYN:
+		tmpvars[oi].value = piYn(piRoundd(value(ci.operators[1]).real()), value(ci.operators[0]).real());
+		break;
+	case PIEvaluatorTypes::bfMin:
+		tmp = value(ci.operators[0]);
+		for (int i = 1; i < ci.operators.size_s(); ++i) {
+			stmp = value(ci.operators[i]);
+			tmp = complexd(piMind(tmp.real(), stmp.real()), piMind(tmp.imag(), stmp.imag()));
+		}
+		tmpvars[oi].value = tmp;
+		break;
+	case PIEvaluatorTypes::bfMax:
+		tmp = value(ci.operators[0]);
+		for (int i = 1; i < ci.operators.size_s(); ++i) {
+			stmp = value(ci.operators[i]);
+			tmp = complexd(piMaxd(tmp.real(), stmp.real()), piMaxd(tmp.imag(), stmp.imag()));
+		}
+		tmpvars[oi].value = tmp;
+		break;
+	case PIEvaluatorTypes::bfClamp:
+		tmp = value(ci.operators[0]);
+		stmp = value(ci.operators[1]);
+		ttmp = value(ci.operators[2]);
+		tmpvars[oi].value = complexd(piClampd(tmp.real(), stmp.real(), ttmp.real()), piClampd(tmp.imag(), stmp.imag(), ttmp.imag()));
+		break;
+	case PIEvaluatorTypes::bfStep:
+		tmpvars[oi].value = complexd(value(ci.operators[0]).real() >= value(ci.operators[1]).real() ? complexld_1 : complexld_0);
+		break;
+	case PIEvaluatorTypes::bfMix:
+		tmp = value(ci.operators[0]);
+		stmp = value(ci.operators[1]);
+		ttmp = value(ci.operators[2]);
+		tmpvars[oi].value = stmp.real() * (1. - tmp.real()) + ttmp.real() * tmp.real();
+		break;
+	case PIEvaluatorTypes::bfDefined:
+		tmpvars[oi].value = value(ci.operators[0]).real() > 0. ? complexd_1 : complexd_0;
+		break;
+	case PIEvaluatorTypes::bfRandom:
+		tmp = static_cast<ldouble>(rand()) / RAND_MAX;
+		stmp = value(ci.operators[1]) - value(ci.operators[0]);
+		tmpvars[oi].value = value(ci.operators[0]) + tmp * stmp;
+		break;
+	default: break;
+	}
+}
+
+
+inline bool PIEvaluator::execInstructions() {
+	PIEvaluatorTypes::Instruction ci;
+	int oi;
+	complexd tmp;
+	tmpvars = variables;
+	//cout << "var count " << tmpvars.size_s() << endl;
+	for (int i = 0; i < instructions.size_s(); i++) {
+		ci = instructions[i];
+		oi = -ci.out - 1;
+		//cout << value(ci.operators[0]) << operationChar(ci.operation) << value(ci.operators[1]) << ", " << oi << endl;
+		switch (ci.operation) {
+		case PIEvaluatorTypes::oAdd:
+			tmpvars[oi].value = value(ci.operators[0]) + value(ci.operators[1]);
+			break;
+		case PIEvaluatorTypes::oSubtract:
+			tmpvars[oi].value = value(ci.operators[0]) - value(ci.operators[1]);
+			break;
+		case PIEvaluatorTypes::oMultiply:
+			tmpvars[oi].value = value(ci.operators[0]) * value(ci.operators[1]);
+			break;
+		case PIEvaluatorTypes::oDivide:
+			tmp = value(ci.operators[1]);
+			if (tmp == complexd(0., 0.)) tmpvars[oi].value = 0.;
+			else tmpvars[oi].value = value(ci.operators[0]) / tmp;
+			break;
+		case PIEvaluatorTypes::oResidue:
+			tmpvars[oi].value = residue(value(ci.operators[0]), value(ci.operators[1]));
+			break;
+		case PIEvaluatorTypes::oPower:
+			tmpvars[oi].value = pow(value(ci.operators[0]), value(ci.operators[1]));
+			break;
+		case PIEvaluatorTypes::oEqual:
+			tmpvars[oi].value = value(ci.operators[0]) == value(ci.operators[1]);
+			break;
+		case PIEvaluatorTypes::oNotEqual:
+			tmpvars[oi].value = value(ci.operators[0]) != value(ci.operators[1]);
+			break;
+		case PIEvaluatorTypes::oGreaterEqual:
+			tmpvars[oi].value = value(ci.operators[0]).real() >= value(ci.operators[1]).real();
+			break;
+		case PIEvaluatorTypes::oSmallerEqual:
+			tmpvars[oi].value = value(ci.operators[0]).real() <= value(ci.operators[1]).real();
+			break;
+		case PIEvaluatorTypes::oGreater:
+			tmpvars[oi].value = value(ci.operators[0]).real() > value(ci.operators[1]).real();
+			break;
+		case PIEvaluatorTypes::oSmaller:
+			tmpvars[oi].value = value(ci.operators[0]).real() < value(ci.operators[1]).real();
+			break;
+		case PIEvaluatorTypes::oAnd:
+			tmpvars[oi].value = value(ci.operators[0]).real() > 0. && value(ci.operators[1]).real() > 0.;
+			break;
+		case PIEvaluatorTypes::oOr:
+			tmpvars[oi].value = value(ci.operators[0]).real() > 0. || value(ci.operators[1]).real() > 0.;
+			break;
+		case PIEvaluatorTypes::oFunction:
+			execFunction(ci);
+			break;
+		case PIEvaluatorTypes::oNone:
+			tmpvars[oi].value = value(ci.operators[0]);
+			break;
+		}
+	}
+	if (!instructions.isEmpty())
+		out = value(instructions.back().out);
+	return true;
+}
+
+
+bool PIEvaluator::check(const PIString & string) {
+	currentString = preprocess(string);
+	correct = check();
+	if (!correct) {
+		instructions.clear();
+		return false;
+	}
+	lastError = "Correct";
+	return true;
+}
+
+
+complexd PIEvaluator::evaluate() {
+	if (!execInstructions()) out = 0.;
+	if (fabs(out.real()) < 1E-300) out = complexd(0., out.imag());
+	if (fabs(out.imag()) < 1E-300) out = complexd(out.real(), 0.);
+	return out;
+}

src/math/pievaluator.h

@@ -0,0 +1,227 @@
+/*! \file pievaluator.h
+ * \brief Mathematic expressions calculator
+*/
+/*
+    PIP - Platform Independent Primitives
+    Evaluator designed for stream calculations
+    Copyright (C) 2014  Ivan Pelipenko peri4ko@gmail.com
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef PIEVALUATOR_H
+#define PIEVALUATOR_H
+
+#include "pistring.h"
+#include "pimath.h"
+
+typedef complexd (*FuncFunc)(void * , int, complexd * );
+
+namespace PIEvaluatorTypes {
+	static const int operationCount = 14;
+	
+	enum eType {etNumber, etOperator, etVariable, etFunction};
+	enum Operation {oNone, oAdd, oSubtract, oMultiply, oDivide, oResidue, oPower,
+					oEqual, oNotEqual, oGreater, oSmaller, oGreaterEqual, oSmallerEqual,
+					oAnd, oOr, oFunction
+	};
+	enum BaseFunctions {bfUnknown, bfSin, bfCos, bfTg, bfCtg,
+						bfArcsin, bfArccos, bfArctg, bfArcctg,
+						bfExp, bfRandom, bfSh, bfCh, bfTh, bfCth,
+						bfSqrt, bfSqr, bfPow, bfAbs,
+						bfLn, bfLg, bfLog, bfSign,
+						bfIm, bfRe, bfArg, bfLen, bfConj,
+						bfRad, bfDeg, bfJ0, bfJ1, bfJN,
+						bfY0, bfY1, bfYN, bfMin, bfMax,
+						bfClamp, bfStep, bfMix, bfDefined,
+						bfCustom = 0xFFFF
+	};
+	
+	struct Instruction {
+		Instruction() {;}
+		Instruction(Operation oper, PIVector<int> opers, int out_ind, int func = -1) {
+			operation = oper; operators = opers; out = out_ind; function = func;}
+		Operation operation;
+		PIVector<int> operators;
+		int out;
+		int function;
+	};
+	struct Element {
+		Element() {;}
+		Element(eType new_type, int new_num, int new_var_num = -1) {set(new_type, new_num, new_var_num);}
+		void set(eType new_type, int new_num, int new_var_num = -1) {type = new_type; num = new_num; var_num = new_var_num;}
+		eType type;
+		int num;
+		int var_num;
+	};
+	struct Function {
+		Function() {arguments = 0; type = bfUnknown; handler = 0;}
+		Function(const PIString & name, int args, BaseFunctions ftype) {identifier = name; arguments = args; type = ftype; handler = 0;}
+		Function(const PIString & name, int args, FuncFunc h) {identifier = name; arguments = args; type = bfCustom; handler = h;}
+		PIString identifier;
+		BaseFunctions type;
+		FuncFunc handler;
+		int arguments;
+	};
+	struct Variable {
+		Variable() {value = 0.;}
+		Variable(const PIString & var_name, complexd val) {name = var_name; value = val;}
+		PIString name;
+		complexd value;
+	};
+};
+/*
+	≠	:
+	≥	}
+	≤	{
+	⋀	&
+	⋁	|
+*/
+
+class PIP_EXPORT PIEvaluatorContent
+{
+	friend class PIEvaluator;
+public:
+	PIEvaluatorContent();
+	~PIEvaluatorContent() {;}
+	
+	void addFunction(const PIString & name, int args = 1) {functions.push_back(PIEvaluatorTypes::Function(name, args, getBaseFunction(name)));}
+	void addVariable(const PIString & name, const complexd & val = 0.) {variables.push_back(PIEvaluatorTypes::Variable(name, val)); sortVariables();}
+	void addCustomFunction(const PIString & name, int args_count, FuncFunc func) {functions << PIEvaluatorTypes::Function(name, args_count, func);}
+	int functionsCount() const {return functions.size();}
+	int variablesCount() const {return variables.size();}
+	int customVariablesCount() const {return variables.size() - cv_count;}
+	int findFunction(const PIString & name) const {for (uint i = 0; i < functions.size(); i++) if (functions[i].identifier == name) return i; return -1;}
+	int findVariable(const PIString & var_name) const {for (uint i = 0; i < variables.size(); i++) if (variables[i].name == var_name) return i; return -1;}
+	PIEvaluatorTypes::Function function(int index) {if (index < 0 || index >= functions.size_s()) return PIEvaluatorTypes::Function(); return functions[index];}
+	PIEvaluatorTypes::Variable variable(int index) {if (index < 0 || index >= variables.size_s()) return PIEvaluatorTypes::Variable(); return variables[index];}
+	PIEvaluatorTypes::Function function(const PIString & name) {return function(findFunction(name));}
+	PIEvaluatorTypes::Variable variable(const PIString & name) {return variable(findVariable(name));}
+	PIEvaluatorTypes::Variable customVariable(int index) {if (index < cv_count || index >= variables.size_s() + cv_count) return PIEvaluatorTypes::Variable(); return variables[index + cv_count];}
+	bool setVariableValue(int index, complexd new_value);
+	bool setVariableName(int index, const PIString & new_name);
+	bool setVariableValue(const PIString & var_name, const complexd & new_value) {return setVariableValue(findVariable(var_name), new_value);}
+	bool setVariableName(const PIString & var_name, const PIString & new_name) {return setVariableName(findVariable(var_name), new_name);}
+	void removeVariable(int index) {variables.remove(index);}
+	void removeVariable(const PIString & var_name) {removeVariable(findVariable(var_name));}
+	void clearCustomVariables();
+	void sortVariables();
+	PIEvaluatorTypes::BaseFunctions getBaseFunction(const PIString & name);
+	
+private:
+	PIVector<PIEvaluatorTypes::Function> functions;
+	PIVector<PIEvaluatorTypes::Variable> variables;
+	int cv_count;
+	
+};
+
+
+class PIP_EXPORT PIEvaluator
+{
+public:
+	
+	//! Constructs an empty evaluator
+	PIEvaluator() {correct = false; data_ = 0;}
+	
+	~PIEvaluator() {;}
+	
+	
+	//! Returns custom data
+	void * data() {return data_;}
+	
+	//! Set custom data to "_data"
+	void setData(void * _data) {data_ = _data;}
+
+	
+	//! Check mathematical expression and parse it to list of instructions
+	bool check(const PIString & string);
+	
+	//! Returns true if expression was checked succesfully
+	bool isCorrect() const {return correct;}
+	
+	//! Set variable value with name "name" to value "value". Add variable if it doesn`t exists
+	int setVariable(const PIString & name, complexd value = 0.) {if (content.findVariable(name) < 0) content.addVariable(name, value); else content.setVariableValue(name, value); return content.findVariable(name);}
+
+	//! Set variable value with index "index" to value "value". Don`t add variable if it doesn`t exists
+	void setVariable(int index, complexd value = 0.) {if (index >= 0 && index < content.variablesCount()) content.setVariableValue(index, value);}
+
+	void setCustomVariableValue(int index, complexd value = 0.) {content.variables[index + content.cv_count].value = value;}
+/*
+	//! Add function "name" with arguments count "args_count" and handler "func". Three arguments will be passed to handler: \a data(), "args_count" and array of input values.
+	void addFunction(const PIString & name, int args_count, FuncFunc func) {content.addCustomFunction(name, args_count, func);}
+*/
+	//! Evaluate last successfully checked with function \a check() expression and returns result
+	complexd evaluate();
+
+	//! Remove variable with name "name"
+	void removeVariable(const PIString & name) {content.removeVariable(name);}
+
+	//! Remove all manually added variables
+	void clearCustomVariables() {content.clearCustomVariables();}
+
+	//! Returns index of variable with name "name"
+	int variableIndex(const PIString & name) const {return content.findVariable(name);}
+	
+	//! Returns all unknown variables founded in last expression passed to \a check() function
+	const PIStringList & unknownVariables() const {return unknownVars;}
+	
+	//! Returns processed last expression passed to \a check() function
+	const PIString & expression() const {return currentString;}
+	
+	//! Returns last error description occured in \a check() function
+	const PIString & error() const {return lastError;}
+	
+	//! Returns last result of \a evaluate()
+	const complexd & lastResult() const {return out;}
+	
+	PIEvaluatorContent content;
+	
+private:
+	const PIString & prepare(const PIString & string);
+	const PIString & preprocess(const PIString & string);
+	int parse(const PIString & string, int offset = 0);
+	void convert();
+	void checkBrackets();
+	void removeSpaces();
+	void findUnknownVariables();
+	void removeJunk();
+	void replaceOperators();
+	void makeOutput(PIString & string);
+	bool fillElements();
+	bool setSignes();
+	bool isSign(const PIChar & ch);
+	PIString inverse(const PIString & string) {int len = string.length(); PIString s; for (int i = 0; i < len; i++) s += string[len - i - 1]; return s;}
+	bool check();
+	bool execInstructions();
+	PIString inBrackets(const PIString & string);
+	PIString operationChar(const PIEvaluatorTypes::Operation & operation);
+	PIEvaluatorTypes::Operation operationInOrder(const int & index);
+	complexd value(const int & index) {if (index < 0) return tmpvars[-index - 1].value; else return kvars->at(index).value;}
+	inline complexd residue(const complexd & f, const complexd & s);
+	inline void execFunction(const PIEvaluatorTypes::Instruction & ci);
+	
+	PIVector<PIEvaluatorTypes::Element> elements;
+	PIVector<PIEvaluatorTypes::Variable> currentVariables, variables, tmpvars, * kvars;
+	PIVector<PIEvaluatorTypes::Instruction> instructions;
+	PIStringList unknownVars;
+	PIString currentString, lastError;
+	complexd out;
+	bool correct;
+	void * data_;
+};
+
+inline bool operator ==(PIEvaluatorTypes::Element e1, PIEvaluatorTypes::Element e2) {return (e1.type == e2.type && e1.num == e2.num);}
+inline bool operator !=(PIEvaluatorTypes::Element e1, PIEvaluatorTypes::Element e2) {return (e1.type != e2.type || e1.num != e2.num);}
+
+#endif // PIEVALUATOR_H

src/math/pifft.cpp

@@ -0,0 +1,936 @@
+/*
+	PIP - Platform Independent Primitives
+	Class for FFT, IFFT and Hilbert transformations
+	Copyright (C) 2014  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;
+			}
+		}
+	}
+}

src/math/pifft.h

@@ -0,0 +1,77 @@
+/*! \file pifft.h
+ * \brief Class for FFT, IFFT and Hilbert transformations
+*/
+/*
+	PIP - Platform Independent Primitives
+	Class for FFT, IFFT and Hilbert transformations
+	Copyright (C) 2014  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/>.
+*/
+
+#ifndef PIFFT_H
+#define PIFFT_H
+
+#include "pimathbase.h"
+
+class PIP_EXPORT PIFFT
+{
+public:
+	PIFFT();
+	
+	PIVector<complexd> * calcFFT(const PIVector<complexd> &val);
+	PIVector<complexd> * calcFFT(const PIVector<double> &val);
+	PIVector<complexd> * calcFFTinverse(const PIVector<complexd> &val);
+	PIVector<complexd> * calcHilbert(const PIVector<double> &val);
+	PIVector<double> getAmplitude();
+	
+private:
+	PIVector<complexd> result;
+	bool prepared;
+	typedef ptrdiff_t ae_int_t;
+	void calc_coefs(uint cnt2);
+	void calc_indexes(uint cnt2, uint deep2);
+	complexd coef(uint n, uint k);
+
+	struct ftplan {
+		PIVector<int> plan;
+		PIVector<double> precomputed;
+		PIVector<double> tmpbuf;
+		PIVector<double> stackbuf;
+	};
+
+	ftplan curplan;
+
+	void fftc1d(const PIVector<complexd> &a, uint n);
+	void fftc1r(const PIVector<double> &a, uint n);
+	void fftc1dinv(const PIVector<complexd> &a, uint n);
+
+	void createPlan(uint n);
+	void ftbasegeneratecomplexfftplan(uint n, ftplan *plan);
+	void ftbase_ftbasegenerateplanrec(int n, int tasktype, ftplan *plan, int *plansize, int *precomputedsize, int *planarraysize, int *tmpmemsize, int *stackmemsize, ae_int_t stackptr, int debugi=0);
+	void ftbase_ftbaseprecomputeplanrec(ftplan *plan, int entryoffset, ae_int_t stackptr);
+	void ftbasefactorize(int n, int *n1, int *n2);
+	void ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int *best);
+	int ftbasefindsmooth(int n);
+	void ftbaseexecuteplan(PIVector<double> *a, int aoffset, int n, ftplan *plan);
+	void ftbaseexecuteplanrec(PIVector<double> *a, int aoffset, ftplan *plan, int entryoffset, ae_int_t stackptr);
+	void ftbase_internalcomplexlintranspose(PIVector<double> *a, int m, int n, int astart, PIVector<double> *buf);
+	void ftbase_ffticltrec(PIVector<double> *a, int astart, int astride, PIVector<double> *b, int bstart, int bstride, int m, int n);
+	void ftbase_internalreallintranspose(PIVector<double> *a, int m, int n, int astart, PIVector<double> *buf);
+	void ftbase_fftirltrec(PIVector<double> *a, int astart, int astride, PIVector<double> *b, int bstart, int bstride, int m, int n);
+	void ftbase_ffttwcalc(PIVector<double> *a, int aoffset, int n1, int n2);
+	
+};
+
+#endif // PIFFT_H

src/math/pimath.h

@@ -0,0 +1,30 @@
+/*! \file pimath.h
+ * \brief Many mathematical functions and classes
+*/
+/*
+	PIP - Platform Independent Primitives
+	Many mathematical functions and classes
+	Copyright (C) 2014  Ivan Pelipenko peri4ko@gmail.com
+
+	This program is free software: you can redistribute it and/or modify
+	it under the terms of the GNU General Public License as published by
+	the Free Software Foundation, either version 3 of the License, or
+	(at your option) any later version.
+
+	This program is distributed in the hope that it will be useful,
+	but WITHOUT ANY WARRANTY; without even the implied warranty of
+	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+	GNU General Public License for more details.
+
+	You should have received a copy of the GNU General Public License
+	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef PIMATH_H
+#define PIMATH_H
+
+#include "pimathsolver.h"
+#include "pistatistic.h"
+#include "pifft.h"
+
+#endif // PIMATH_H

src/math/pimathbase.cpp

@@ -0,0 +1,468 @@
+/*
+	PIP - Platform Independent Primitives
+	Basic mathematical functions and defines
+	Copyright (C) 2014  Ivan Pelipenko peri4ko@gmail.com
+
+	This program is free software: you can redistribute it and/or modify
+	it under the terms of the GNU General Public License as published by
+	the Free Software Foundation, either version 3 of the License, or
+	(at your option) any later version.
+
+	This program is distributed in the hope that it will be useful,
+	but WITHOUT ANY WARRANTY; without even the implied warranty of
+	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+	GNU General Public License for more details.
+
+	You should have received a copy of the GNU General Public License
+	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "pimathbase.h"
+
+
+double piJ0(const double & v) {
+#ifndef PIP_MATH_J0
+	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;
+#else
+	return j0(v);
+#endif
+}
+
+
+double piJ1(const double & v) {
+#ifndef PIP_MATH_J1
+	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;
+#else
+	return j1(v);
+#endif
+}
+
+
+double piJn(int n, const double & v) {
+#ifndef PIP_MATH_JN
+	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 * piJ0(x);
+		return result;
+	}
+	if (n == 1) {
+		result = sg * piJ1(x);
+		return result;
+	}
+	if (n == 2) {
+		if (x == 0)
+			result = 0;
+		else
+			result = sg * (2.0 * piJ1(x) / x - piJ0(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 = piJ1(x) / pk;
+	else
+		ans = piJ0(x) / pkm1;
+	result = sg * ans;
+	return result;
+#else
+	return jn(n, v);
+#endif
+}
+
+
+double piY0(const double & v) {
+#ifndef PIP_MATH_Y0
+	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 * piJ0(x) * log(x);
+	return result;
+#else
+	return y0(v);
+#endif
+}
+
+
+double piY1(const double & v) {
+#ifndef PIP_MATH_Y1
+	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 * (piJ1(x) * log(x) - 1 / x);
+	return result;
+#else
+	return y1(v);
+#endif
+}
+
+
+double piYn(int n, const double & v) {
+#ifndef PIP_MATH_YN
+	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 = piY0(x);
+		return result;
+	}
+	if (n == 1) {
+		result = s * piY1(x);
+		return result;
+	}
+	a = piY0(x);
+	b = piY1(x);
+	for (i = 1; i <= n - 1; i++) {
+		tmp = b;
+		b = 2 * i / x * b - a;
+		a = tmp;
+	}
+	result = s * b;
+	return result;
+#else
+	return yn(n, v);
+#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;
+}

src/math/pimathbase.h

@@ -0,0 +1,189 @@
+/*! \file pimathbase.h
+ * \brief Basic mathematical functions and defines
+*/
+/*
+	PIP - Platform Independent Primitives
+	Basic mathematical functions and defines
+	Copyright (C) 2014  Ivan Pelipenko peri4ko@gmail.com
+
+	This program is free software: you can redistribute it and/or modify
+	it under the terms of the GNU General Public License as published by
+	the Free Software Foundation, either version 3 of the License, or
+	(at your option) any later version.
+
+	This program is distributed in the hope that it will be useful,
+	but WITHOUT ANY WARRANTY; without even the implied warranty of
+	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+	GNU General Public License for more details.
+
+	You should have received a copy of the GNU General Public License
+	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef PIMATHBASE_H
+#define PIMATHBASE_H
+
+#include "piinit.h"
+#include "pibytearray.h"
+#ifdef QNX
+#  undef PIP_MATH_J0
+#  undef PIP_MATH_J1
+#  undef PIP_MATH_JN
+#  undef PIP_MATH_Y0
+#  undef PIP_MATH_Y1
+#  undef PIP_MATH_YN
+#endif
+
+#ifndef M_LN2
+#  define M_LN2 0.69314718055994530942
+#endif
+#ifndef M_LN10
+#  define M_LN10 2.30258509299404568402
+#endif
+#ifndef M_SQRT2
+#  define M_SQRT2 1.41421356237309514547
+#endif
+#ifndef M_SQRT3
+#  define M_SQRT3 1.73205080756887719318
+#endif
+#ifndef M_1_SQRT2
+#  define M_1_SQRT2 0.70710678118654746172
+#endif
+#ifndef M_1_SQRT3
+#  define M_1_SQRT3 0.57735026918962584208
+#endif
+#ifndef M_PI
+#  define M_PI 3.14159265358979323846
+#endif
+#ifndef M_2PI
+#  define M_2PI 6.28318530717958647692
+#endif
+#ifndef M_PI_3
+#  define M_PI_3 1.04719755119659774615
+#endif
+#ifndef M_2PI_3
+#  define M_2PI_3 2.0943951023931954923
+#endif
+#ifndef M_180_PI
+#  define M_180_PI 57.2957795130823208768
+#endif
+#ifndef M_PI_180
+#  define M_PI_180 1.74532925199432957692e-2
+#endif
+#ifndef M_E
+#  define M_E 2.7182818284590452353602874713527
+#endif
+#ifndef M_LIGHT_SPEED
+#  define M_LIGHT_SPEED 2.99792458e+8
+#endif
+
+const double deg2rad = M_PI_180;
+const double rad2deg = M_180_PI;
+
+inline int sign(const float & x) {return (x < 0.) ? -1 : (x > 0. ? 1 : 0);}
+inline int sign(const double & x) {return (x < 0.) ? -1 : (x > 0. ? 1 : 0);}
+inline complexd sign(const complexd & x) {return complexd(sign(x.real()), sign(x.imag()));}
+inline int pow2(const int p) {return 1 << p;}
+inline double sqr(const int v) {return v * v;}
+inline double sqr(const float & v) {return v * v;}
+inline double sqr(const double & v) {return v * v;}
+inline double sinc(const double & v) {if (v == 0.) return 1.; double t = M_PI * v; return sin(t) / t;}
+inline complexd round(const complexd & c) {return complexd(piRound<double>(c.real()), piRound<double>(c.imag()));}
+inline complexd floor(const complexd & c) {return complexd(floor(c.real()), floor(c.imag()));}
+inline complexd ceil(const complexd & c) {return complexd(ceil(c.real()), ceil(c.imag()));}
+inline complexd atanc(const complexd & c) {return -complexd(-0.5, 1.) * log((complexd_1 + complexd_i * c) / (complexd_1 - complexd_i * c));}
+inline complexd asinc(const complexd & c) {return -complexd_i * log(complexd_i * c + sqrt(complexd_1 - c * c));}
+inline complexd acosc(const complexd & c) {return -complexd_i * log(c + complexd_i * sqrt(complexd_1 - c * c));}
+#ifdef CC_GCC
+#  if CC_GCC_VERSION <= 0x025F
+inline complexd tan(const complexd & c) {return sin(c) / cos(c);}
+inline complexd tanh(const complexd & c) {return sinh(c) / cosh(c);}
+inline complexd log2(const complexd & c) {return log(c) / M_LN2;}
+inline complexd log10(const complexd & c) {return log(c) / M_LN10;}
+#  endif
+#endif
+double piJ0(const double & v);
+double piJ1(const double & v);
+double piJn(int n, const double & v);
+double piY0(const double & v);
+double piY1(const double & v);
+double piYn(int n, const double & v);
+inline double toDb(double val) {return 10. * log10(val);}
+inline double fromDb(double val) {return pow(10., val / 10.);}
+inline double toRad(double deg) {return deg * M_PI_180;}
+inline double toDeg(double rad) {return rad * M_180_PI;}
+
+template<typename T>
+inline PICout operator <<(PICout s, const complex<T> & v) {s.space(); s.setControl(0, true); s << "(" << v.real() << "; " << v.imag() << ")"; s.restoreControl(); return s;}
+
+// [-1 ; 1]
+inline double randomd() {return (double)random() / RAND_MAX * 2. - 1.;}
+// [-1 ; 1] normal
+double randomn(double dv = 0., double sv = 1.);
+
+inline PIVector<double> abs(const PIVector<complexd> & v) {
+	PIVector<double> result;
+	result.resize(v.size());
+	for (uint i = 0; i < v.size(); i++)
+		result[i] = abs(v[i]);
+	return result;
+}
+inline PIVector<double> abs(const PIVector<double> & v) {
+	PIVector<double> result;
+	result.resize(v.size());
+	for (uint i = 0; i < v.size(); i++)
+		result[i] = abs(v[i]);
+	return result;
+}
+
+
+template <typename T>
+bool OLS_Linear(const PIVector<PIPair<T, T> > & input, T * out_a, T * out_b) {
+	if (input.size_s() < 2)
+		return false;
+	int n = input.size_s();
+	T a_t0 = T(), a_t1 = T(), a_t2 = T(), a_t3 = T(), a_t4 = T(), a = T(), b = T();
+	for (int i = 0; i < n; ++i) {
+		const PIPair<T, T> & cv(input[i]);
+		a_t0 += cv.first * cv.second;
+		a_t1 += cv.first;
+		a_t2 += cv.second;
+		a_t3 += cv.first * cv.first;
+	}
+	a_t4 = n * a_t3 - a_t1 * a_t1;
+	if (a_t4 != T())
+		a = (n * a_t0 - a_t1 * a_t2) / a_t4;
+	b = (a_t2 - a * a_t1) / n;
+	if (out_a != 0) *out_a = a;
+	if (out_b != 0) *out_b = b;
+	return true;
+}
+
+
+template <typename T>
+bool WLS_Linear(const PIVector<PIPair<T, T> > & input, const PIVector<T> & weights, T * out_a, T * out_b) {
+	if (input.size_s() < 2)
+		return false;
+	if (input.size_s() != weights.size_s())
+		return false;
+	int n = input.size_s();
+	T a_t0 = T(), a_t1 = T(), a_t2 = T(), a_t3 = T(), a_t4 = T(), a_n = T(), a = T(), b = T();
+	for (int i = 0; i < n; ++i) {
+		T cp = weights[i];
+		const PIPair<T, T> & cv(input[i]);
+		a_t0 += cv.first * cv.second * cp;
+		a_t1 += cv.first * cp;
+		a_t2 += cv.second * cp;
+		a_t3 += cv.first * cv.first * cp;
+		a_n += cp;
+	}
+	a_t4 = a_n * a_t3 - a_t1 * a_t1;
+	if (a_t4 != T())
+		a = (a_n * a_t0 - a_t1 * a_t2) / a_t4;
+	b = (a_t2 - a * a_t1) / a_n;
+	if (out_a != 0) *out_a = a;
+	if (out_b != 0) *out_b = b;
+	return true;
+}
+
+#endif // PIMATHBASE_H

src/math/pimathmatrix.h

@@ -0,0 +1,504 @@
+/*! \file pimathmatrix.h
+ * \brief PIMathMatrix
+*/
+/*
+	PIP - Platform Independent Primitives
+	PIMathMatrix
+	Copyright (C) 2014  Ivan Pelipenko peri4ko@gmail.com
+
+	This program is free software: you can redistribute it and/or modify
+	it under the terms of the GNU General Public License as published by
+	the Free Software Foundation, either version 3 of the License, or
+	(at your option) any later version.
+
+	This program is distributed in the hope that it will be useful,
+	but WITHOUT ANY WARRANTY; without even the implied warranty of
+	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+	GNU General Public License for more details.
+
+	You should have received a copy of the GNU General Public License
+	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef PIMATHMATRIX_H
+#define PIMATHMATRIX_H
+
+#include "pimathvector.h"
+
+/// Matrix templated
+
+#define PIMM_FOR(r, c) for (uint c = 0; c < Cols; ++c) { for (uint r = 0; r < Rows; ++r) {
+#define PIMM_FOR_WB(r, c) for (uint c = 0; c < Cols; ++c) for (uint r = 0; r < Rows; ++r) // without brakes
+#define PIMM_FOR_I(r, c) for (uint r = 0; r < Rows; ++r) { for (uint c = 0; c < Cols; ++c) {
+#define PIMM_FOR_I_WB(r, c) for (uint r = 0; r < Rows; ++r) for (uint c = 0; c < Cols; ++c) // without brakes
+#define PIMM_FOR_C(v) for (uint v = 0; v < Cols; ++v)
+#define PIMM_FOR_R(v) for (uint v = 0; v < Rows; ++v)
+
+#pragma pack(push, 1)
+template<uint Rows, uint Cols = Rows, typename Type = double>
+class PIP_EXPORT PIMathMatrixT {
+	typedef PIMathMatrixT<Rows, Cols, Type> _CMatrix;
+	typedef PIMathMatrixT<Cols, Rows, Type> _CMatrixI;
+	typedef PIMathVectorT<Rows, Type> _CMCol;
+	typedef PIMathVectorT<Cols, Type> _CMRow;
+public:
+	PIMathMatrixT() {resize(Rows, Cols);}
+	PIMathMatrixT(Type fval, ...) {resize(Rows, Cols); va_list vl; va_start(vl, fval); PIMM_FOR_I_WB(r, c) m[r][c] = (r + c == 0 ? fval : va_arg(vl, Type)); va_end(vl);}
+	PIMathMatrixT(const PIVector<Type> & val) {resize(Rows, Cols); int i = 0; PIMM_FOR_I_WB(r, c) m[r][c] = val[i++];}
+	//PIMathMatrixT(const _CMatrix & o) {resize(Rows, Cols); int i = 0; PIMM_FOR_I_WB(r, c) m[r][c] = val[i++];}
+
+	static _CMatrix identity() {_CMatrix tm = _CMatrix(); PIMM_FOR_WB(r, c) tm.m[r][c] = (c == r ? Type(1) : Type(0)); return tm;}
+	static _CMatrix rotation(double angle) {return _CMatrix();}
+	static _CMatrix rotationX(double angle) {return _CMatrix();}
+	static _CMatrix rotationY(double angle) {return _CMatrix();}
+	static _CMatrix rotationZ(double angle) {return _CMatrix();}
+	static _CMatrix scaleX(double factor) {return _CMatrix();}
+	static _CMatrix scaleY(double factor) {return _CMatrix();}
+	static _CMatrix scaleZ(double factor) {return _CMatrix();}
+
+	uint cols() const {return Cols;}
+	uint rows() const {return Rows;}
+	_CMCol col(uint index) {_CMCol tv; PIMM_FOR_R(i) tv[i] = m[i][index]; return tv;}
+	_CMRow row(uint index) {_CMRow tv; PIMM_FOR_C(i) tv[i] = m[index][i]; return tv;}
+	_CMatrix & setCol(uint index, const _CMCol & v) {PIMM_FOR_R(i) m[i][index] = v[i]; return *this;}
+	_CMatrix & setRow(uint index, const _CMRow & v) {PIMM_FOR_C(i) m[index][i] = v[i]; return *this;}
+	_CMatrix & swapRows(uint r0, uint r1) {Type t; PIMM_FOR_C(i) {t = m[r0][i]; m[r0][i] = m[r1][i]; m[r1][i] = t;} return *this;}
+	_CMatrix & swapCols(uint c0, uint c1) {Type t; PIMM_FOR_R(i) {t = m[i][c0]; m[i][c0] = m[i][c1]; m[i][c1] = t;} return *this;}
+	_CMatrix & fill(const Type & v) {PIMM_FOR_WB(r, c) m[r][c] = v; return *this;}
+	//inline _CMatrix & set(Type fval, ...) {m[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) m[i] = va_arg(vl, Type); va_end(vl); return *this;}
+	//inline void normalize() {Type tv = length(); if (tv == Type(1)) return; PIMV_FOR(i, 0) m[i] /= tv;}
+	bool isSquare() const {return cols() == rows();}
+	bool isIdentity() const {PIMM_FOR_WB(r, c) if ((c == r) ? m[r][c] != Type(1) : m[r][c] != Type(0)) return false; return true;}
+	bool isNull() const {PIMM_FOR_WB(r, c) if (m[r][c] != Type(0)) return false; return true;}
+
+	Type & at(uint row, uint col) {return m[row][col];}
+	Type at(uint row, uint col) const {return m[row][col];}
+	Type * operator [](uint row) {return m[row];}
+	const Type * operator [](uint row) const {return m[row];}
+	void operator =(const _CMatrix & sm) {memcpy(m, sm.m, sizeof(Type) * Cols * Rows);}
+	bool operator ==(const _CMatrix & sm) const {PIMM_FOR_WB(r, c) if (m[r][c] != sm.m[r][c]) return false; return true;}
+	bool operator !=(const _CMatrix & sm) const {return !(*this == sm);}
+	void operator +=(const _CMatrix & sm) {PIMM_FOR_WB(r, c) m[r][c] += sm.m[r][c];}
+	void operator -=(const _CMatrix & sm) {PIMM_FOR_WB(r, c) m[r][c] -= sm.m[r][c];}
+	void operator *=(const Type & v) {PIMM_FOR_WB(r, c) m[r][c] *= v;}
+	void operator /=(const Type & v) {PIMM_FOR_WB(r, c) m[r][c] /= v;}
+	_CMatrix operator -() {_CMatrix tm; PIMM_FOR_WB(r, c) tm.m[r][c] = -m[r][c]; return tm;}
+	_CMatrix operator +(const _CMatrix & sm) {_CMatrix tm = _CMatrix(*this); PIMM_FOR_WB(r, c) tm.m[r][c] += sm.m[r][c]; return tm;}
+	_CMatrix operator -(const _CMatrix & sm) {_CMatrix tm = _CMatrix(*this); PIMM_FOR_WB(r, c) tm.m[r][c] -= sm.m[r][c]; return tm;}
+	_CMatrix operator *(const Type & v) {_CMatrix tm = _CMatrix(*this); PIMM_FOR_WB(r, c) tm.m[r][c] *= v; return tm;}
+	_CMatrix operator /(const Type & v) {_CMatrix tm = _CMatrix(*this); PIMM_FOR_WB(r, c) tm.m[r][c] /= v; return tm;}
+	
+	Type determinant(bool * ok = 0) const {
+		_CMatrix m(*this);
+		bool k;
+		Type ret = Type(0);
+		m.toUpperTriangular(&k);
+		if (ok) *ok = k;
+		if (!k) return ret;
+		ret = Type(1);
+		for (uint c = 0; c < Cols; ++c)
+			for (uint r = 0; r < Rows; ++r)
+				if (r == c)
+					ret *= m[r][c];
+		return ret;
+	}
+
+	_CMatrix & toUpperTriangular(bool * ok = 0) {
+		if (Cols != Rows) {
+			if (ok != 0) *ok = false;
+			return *this;
+		}
+		_CMatrix smat(*this);
+		bool ndet;
+		uint crow;
+		Type mul;
+		for (uint i = 0; i < Cols; ++i) {
+			ndet = true;
+			for (uint j = 0; j < Rows; ++j) if (smat.m[i][j] != 0) ndet = false;
+			if (ndet) {
+				if (ok != 0) *ok = false;
+				return *this;
+			}
+			for (uint j = 0; j < Cols; ++j) if (smat.m[j][i] != 0) ndet = false;
+			if (ndet) {
+				if (ok != 0) *ok = false;
+				return *this;
+			}
+		}
+		for (uint i = 0; i < Cols; ++i) {
+			crow = i;
+			while (smat.m[i][i] == Type(0))
+				smat.swapRows(i, ++crow);
+			for (uint j = i + 1; j < Rows; ++j) {
+				mul = smat.m[i][j] / smat.m[i][i];
+				for (uint k = i; k < Cols; ++k) smat.m[k][j] -= mul * smat.m[k][i];
+			}
+			if (i < Cols - 1) {
+				if (fabs(smat.m[i+1][i+1]) < Type(1E-100)) {
+					if (ok != 0) *ok = false;
+					return *this;
+				}
+			}
+		}
+		if (ok != 0) *ok = true;
+		memcpy(m, smat.m, sizeof(Type) * Cols * Rows);
+		return *this;
+	}
+
+	_CMatrix & invert(bool * ok = 0) {
+		if (Cols != Rows) {
+			if (ok != 0) *ok = false;
+			return *this;
+		}
+		_CMatrix mtmp = _CMatrix::identity(), smat(*this);
+		bool ndet;
+		uint crow;
+		Type mul, iddiv;
+		for (uint i = 0; i < Cols; ++i) {
+			ndet = true;
+			for (uint j = 0; j < Rows; ++j) if (smat.m[i][j] != 0) ndet = false;
+			if (ndet) {
+				if (ok != 0) *ok = false;
+				return *this;
+			}
+			for (uint j = 0; j < Cols; ++j) if (smat.m[j][i] != 0) ndet = false;
+			if (ndet) {
+				if (ok != 0) *ok = false;
+				return *this;
+			}
+		}
+		for (uint i = 0; i < Cols; ++i) {
+			crow = i;
+			while (smat.m[i][i] == Type(0)) {
+				++crow;
+				smat.swapRows(i, crow);
+				mtmp.swapRows(i, crow);
+			}
+			for (uint j = i + 1; j < Rows; ++j) {
+				mul = smat.m[i][j] / smat.m[i][i];
+				for (uint k = i; k < Cols; ++k) smat.m[k][j] -= mul * smat.m[k][i];
+				for (uint k = 0; k < Cols; ++k) mtmp.m[k][j] -= mul * mtmp.m[k][i];
+			}
+			//cout << i << endl << smat << endl;
+			if (i < Cols - 1) {
+				if (fabs(smat.m[i+1][i+1]) < Type(1E-100)) {
+					if (ok != 0) *ok = false;
+					return *this;
+				}
+			}
+			iddiv = smat.m[i][i];
+			for (uint j = i; j < Cols; ++j) smat.m[j][i] /= iddiv;
+			for (uint j = 0; j < Cols; ++j) mtmp.m[j][i] /= iddiv;
+		}
+		for (uint i = Cols - 1; i > 0; --i) {
+			for (uint j = 0; j < i; ++j) {
+				mul = smat.m[i][j];
+				smat.m[i][j] -= mul;
+				for (uint k = 0; k < Cols; ++k) mtmp.m[k][j] -= mtmp.m[k][i] * mul;
+			}
+		}
+		if (ok != 0) *ok = true;
+		memcpy(m, mtmp.m, sizeof(Type) * Cols * Rows);
+		return *this;
+	}
+	_CMatrix inverted(bool * ok = 0) const {_CMatrix tm(*this); tm.invert(ok); return tm;}
+	_CMatrixI transposed() const {_CMatrixI tm; PIMM_FOR_WB(r, c) tm[r][c] = m[r][c]; return tm;}
+
+private:
+	void resize(uint rows_, uint cols_, const Type & new_value = Type()) {r_ = rows_; c_ = cols_; PIMM_FOR_WB(r, c) m[r][c] = new_value;}
+	int c_, r_;
+	Type m[Rows][Cols];
+
+};
+#pragma pack(pop)
+
+
+template<> inline PIMathMatrixT<2u, 2u> PIMathMatrixT<2u, 2u>::rotation(double angle) {double c = cos(angle), s = sin(angle); PIMathMatrixT<2u, 2u> tm; tm[0][0] = tm[1][1] = c; tm[0][1] = -s; tm[1][0] = s; return tm;}
+template<> inline PIMathMatrixT<2u, 2u> PIMathMatrixT<2u, 2u>::scaleX(double factor) {PIMathMatrixT<2u, 2u> tm; tm[0][0] = factor; tm[1][1] = 1.; return tm;}
+template<> inline PIMathMatrixT<2u, 2u> PIMathMatrixT<2u, 2u>::scaleY(double factor) {PIMathMatrixT<2u, 2u> tm; tm[0][0] = 1.; tm[1][1] = factor; return tm;}
+
+template<> inline PIMathMatrixT<3u, 3u> PIMathMatrixT<3u, 3u>::rotationX(double angle) {double c = cos(angle), s = sin(angle); PIMathMatrixT<3u, 3u> tm; tm[0][0] = 1.; tm[1][1] = tm[2][2] = c; tm[2][1] = s; tm[1][2] = -s; return tm;}
+template<> inline PIMathMatrixT<3u, 3u> PIMathMatrixT<3u, 3u>::rotationY(double angle) {double c = cos(angle), s = sin(angle); PIMathMatrixT<3u, 3u> tm; tm[1][1] = 1.; tm[0][0] = tm[2][2] = c; tm[2][0] = -s; tm[0][2] = s; return tm;}
+template<> inline PIMathMatrixT<3u, 3u> PIMathMatrixT<3u, 3u>::rotationZ(double angle) {double c = cos(angle), s = sin(angle); PIMathMatrixT<3u, 3u> tm; tm[2][2] = 1.; tm[0][0] = tm[1][1] = c; tm[1][0] = s; tm[0][1] = -s; return tm;}
+template<> inline PIMathMatrixT<3u, 3u> PIMathMatrixT<3u, 3u>::scaleX(double factor) {PIMathMatrixT<3u, 3u> tm; tm[1][1] = tm[2][2] = 1.; tm[0][0] = factor; return tm;}
+template<> inline PIMathMatrixT<3u, 3u> PIMathMatrixT<3u, 3u>::scaleY(double factor) {PIMathMatrixT<3u, 3u> tm; tm[0][0] = tm[2][2] = 1.; tm[1][1] = factor; return tm;}
+template<> inline PIMathMatrixT<3u, 3u> PIMathMatrixT<3u, 3u>::scaleZ(double factor) {PIMathMatrixT<3u, 3u> tm; tm[0][0] = tm[1][1] = 1.; tm[2][2] = factor; return tm;}
+
+template<uint Rows, uint Cols, typename Type>
+inline std::ostream & operator <<(std::ostream & s, const PIMathMatrixT<Rows, Cols, Type> & m) {s << '{'; PIMM_FOR_I(r, c) s << m[r][c]; if (c < Cols - 1 || r < Rows - 1) s << ", ";} if (r < Rows - 1) s << endl << ' ';} s << '}'; return s;}
+template<uint Rows, uint Cols, typename Type>
+inline PICout operator <<(PICout s, const PIMathMatrixT<Rows, Cols, Type> & m) {s << '{'; PIMM_FOR_I(r, c) s << m[r][c]; if (c < Cols - 1 || r < Rows - 1) s << ", ";} if (r < Rows - 1) s << NewLine << ' ';} s << '}'; return s;}
+
+/// Multiply matrices {Rows0 x CR} on {CR x Cols1}, result is {Rows0 x Cols1}
+template<uint CR, uint Rows0, uint Cols1, typename Type>
+inline PIMathMatrixT<Rows0, Cols1, Type> operator *(const PIMathMatrixT<Rows0, CR, Type> & fm,
+													const PIMathMatrixT<CR, Cols1, Type> & sm) {
+	PIMathMatrixT<Rows0, Cols1, Type> tm;
+	Type t;
+	for (uint j = 0; j < Rows0; ++j) {
+		for (uint i = 0; i < Cols1; ++i) {
+			t = Type(0);
+			for (uint k = 0; k < CR; ++k)
+				t += fm[j][k] * sm[k][i];
+			tm[j][i] = t;
+		}
+	}
+	return tm;
+}
+
+/// Multiply matrix {Rows x Cols} on vector {Rows}, result is vector {Cols}
+template<uint Cols, uint Rows, typename Type>
+inline PIMathVectorT<Cols, Type> operator *(const PIMathMatrixT<Rows, Cols, Type> & fm,
+											const PIMathVectorT<Rows, Type> & sv) {
+	PIMathVectorT<Rows, Type> tv;
+	Type t;
+	for (uint j = 0; j < Rows; ++j) {
+		t = Type(0);
+		for (uint i = 0; i < Cols; ++i)
+			t += fm[j][i] * sv[i];
+		tv[j] = t;
+	}
+	return tv;
+}
+
+typedef PIMathMatrixT<2u, 2u, int> PIMathMatrixT22i;
+typedef PIMathMatrixT<3u, 3u, int> PIMathMatrixT33i;
+typedef PIMathMatrixT<4u, 4u, int> PIMathMatrixT44i;
+typedef PIMathMatrixT<2u, 2u, double> PIMathMatrixT22d;
+typedef PIMathMatrixT<3u, 3u, double> PIMathMatrixT33d;
+typedef PIMathMatrixT<4u, 4u, double> PIMathMatrixT44d;
+
+
+template<typename Type>
+class PIMathMatrix;
+
+#undef PIMV_FOR
+#undef PIMM_FOR
+#undef PIMM_FOR_WB
+#undef PIMM_FOR_I
+#undef PIMM_FOR_I_WB
+#undef PIMM_FOR_C
+#undef PIMM_FOR_R
+
+
+
+
+
+/// Matrix
+
+#define PIMM_FOR(c, r) for (uint c = 0; c < cols_; ++c) { for (uint r = 0; r < rows_; ++r) {
+#define PIMM_FOR_WB(c, r) for (uint c = 0; c < cols_; ++c) for (uint r = 0; r < rows_; ++r) // without brakes
+#define PIMM_FOR_I(c, r) for (uint r = 0; r < rows_; ++r) { for (uint c = 0; c < cols_; ++c) {
+#define PIMM_FOR_I_WB(c, r) for (uint r = 0; r < rows_; ++r) for (uint c = 0; c < cols_; ++c) // without brakes
+#define PIMM_FOR_C(v) for (uint v = 0; v < cols_; ++v)
+#define PIMM_FOR_R(v) for (uint v = 0; v < rows_; ++v)
+
+template<typename Type>
+class PIP_EXPORT PIMathMatrix {
+	typedef PIMathMatrix<Type> _CMatrix;
+	typedef PIMathVector<Type> _CMCol;
+	typedef PIMathVector<Type> _CMRow;
+public:
+	PIMathMatrix(const uint cols = 3, const uint rows = 3) {resize(cols, rows);}
+	PIMathMatrix(const uint cols, const uint rows, Type fval, ...) {resize(cols, rows); va_list vl; va_start(vl, fval); PIMM_FOR_I_WB(c, r) m[c][r] = (r + c == 0 ? fval : va_arg(vl, Type)); va_end(vl);}
+	PIMathMatrix(const uint cols, const uint rows, const PIVector<Type> & val) {resize(cols, rows); int i = 0; PIMM_FOR_I_WB(c, r) m[c][r] = val[i++];}
+
+	static _CMatrix identity(const uint cols_, const uint rows_) {_CMatrix tm(cols_, rows_); PIMM_FOR_WB(c, r) tm.m[c][r] = (c == r ? Type(1) : Type(0)); return tm;}
+
+	uint cols() const {return cols_;}
+	uint rows() const {return rows_;}
+	_CMCol col(uint index) {_CMCol tv; PIMM_FOR_R(i) tv[i] = m[index][i]; return tv;}
+	_CMRow row(uint index) {_CMRow tv; PIMM_FOR_C(i) tv[i] = m[i][index]; return tv;}
+	_CMatrix & resize(const uint cols, const uint rows, const Type & new_value = Type()) {cols_ = cols; rows_ = rows; m.resize(cols); PIMM_FOR_C(i) m[i].resize(rows, new_value); return *this;}
+	_CMatrix & setCol(uint index, const _CMCol & v) {PIMM_FOR_R(i) m[index][i] = v[i]; return *this;}
+	_CMatrix & setRow(uint index, const _CMRow & v) {PIMM_FOR_C(i) m[i][index] = v[i]; return *this;}
+	_CMatrix & swapRows(uint r0, uint r1) {Type t; PIMM_FOR_C(i) {t = m[i][r0]; m[i][r0] = m[i][r1]; m[i][r1] = t;} return *this;}
+	_CMatrix & swapCols(uint c0, uint c1) {Type t; PIMM_FOR_R(i) {t = m[c0][i]; m[c0][i] = m[c1][i]; m[c1][i] = t;} return *this;}
+	_CMatrix & fill(const Type & v) {PIMM_FOR_WB(c, r) m[c][r] = v; return *this;}
+	//inline _CMatrix & set(Type fval, ...) {m[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) m[i] = va_arg(vl, Type); va_end(vl); return *this;}
+	//inline void normalize() {Type tv = length(); if (tv == Type(1)) return; PIMV_FOR(i, 0) m[i] /= tv;}
+	bool isSquare() const {return cols() == rows();}
+	bool isIdentity() const {PIMM_FOR_WB(c, r) if ((c == r) ? m[c][r] != Type(1) : m[c][r] != Type(0)) return false; return true;}
+	bool isNull() const {PIMM_FOR_WB(c, r) if (m[c][r] != Type(0)) return false; return true;}
+
+	Type & at(uint col, uint row) {return m[col][row];}
+	Type at(uint col, uint row) const {return m[col][row];}
+	PIVector<Type> & operator [](uint col) {return m[col];}
+	PIVector<Type> operator [](uint col) const {return m[col];}
+	void operator =(const _CMatrix & sm) {m = sm.m;}
+	bool operator ==(const _CMatrix & sm) const {PIMM_FOR_WB(c, r) if (m[c][r] != sm.m[c][r]) return false; return true;}
+	bool operator !=(const _CMatrix & sm) const {return !(*this == sm);}
+	void operator +=(const _CMatrix & sm) {PIMM_FOR_WB(c, r) m[c][r] += sm.m[c][r];}
+	void operator -=(const _CMatrix & sm) {PIMM_FOR_WB(c, r) m[c][r] -= sm.m[c][r];}
+	void operator *=(const Type & v) {PIMM_FOR_WB(c, r) m[c][r] *= v;}
+	void operator /=(const Type & v) {PIMM_FOR_WB(c, r) m[c][r] /= v;}
+	_CMatrix operator -() {_CMatrix tm(*this); PIMM_FOR_WB(c, r) tm.m[c][r] = -m[c][r]; return tm;}
+	_CMatrix operator +(const _CMatrix & sm) {_CMatrix tm(*this); PIMM_FOR_WB(c, r) tm.m[c][r] += sm.m[c][r]; return tm;}
+	_CMatrix operator -(const _CMatrix & sm) {_CMatrix tm(*this); PIMM_FOR_WB(c, r) tm.m[c][r] -= sm.m[c][r]; return tm;}
+	_CMatrix operator *(const Type & v) {_CMatrix tm(*this); PIMM_FOR_WB(c, r) tm.m[c][r] *= v; return tm;}
+	_CMatrix operator /(const Type & v) {_CMatrix tm(*this); PIMM_FOR_WB(c, r) tm.m[c][r] /= v; return tm;}
+
+	_CMatrix & toUpperTriangular(bool * ok = 0) {
+		if (cols_ != rows_) {
+			if (ok != 0) *ok = false;
+			return *this;
+		}
+		_CMatrix smat(*this);
+		bool ndet;
+		uint crow;
+		Type mul;
+		for (uint i = 0; i < cols_; ++i) {
+			ndet = true;
+			for (uint j = 0; j < rows_; ++j) if (smat.m[i][j] != 0) ndet = false;
+			if (ndet) {
+				if (ok != 0) *ok = false;
+				return *this;
+			}
+			for (uint j = 0; j < cols_; ++j) if (smat.m[j][i] != 0) ndet = false;
+			if (ndet) {
+				if (ok != 0) *ok = false;
+				return *this;
+			}
+		}
+		for (uint i = 0; i < cols_; ++i) {
+			crow = i;
+			while (smat.m[i][i] == Type(0))
+				smat.swapRows(i, ++crow);
+			for (uint j = i + 1; j < rows_; ++j) {
+				mul = smat.m[i][j] / smat.m[i][i];
+				for (uint k = i; k < cols_; ++k) smat.m[k][j] -= mul * smat.m[k][i];
+			}
+			if (i < cols_ - 1) {
+				if (fabs(smat.m[i+1][i+1]) < Type(1E-100)) {
+					if (ok != 0) *ok = false;
+					return *this;
+				}
+			}
+		}
+		if (ok != 0) *ok = true;
+		m = smat.m;
+		return *this;
+	}
+
+	_CMatrix & invert(bool * ok = 0, _CMCol * sv = 0) {
+		if (cols_ != rows_) {
+			if (ok != 0) *ok = false;
+			return *this;
+		}
+		_CMatrix mtmp = _CMatrix::identity(cols_, rows_), smat(*this);
+		bool ndet;
+		uint crow;
+		Type mul, iddiv;
+		for (uint i = 0; i < cols_; ++i) {
+			ndet = true;
+			for (uint j = 0; j < rows_; ++j) if (smat.m[i][j] != 0) ndet = false;
+			if (ndet) {
+				if (ok != 0) *ok = false;
+				return *this;
+			}
+			for (uint j = 0; j < cols_; ++j) if (smat.m[j][i] != 0) ndet = false;
+			if (ndet) {
+				if (ok != 0) *ok = false;
+				return *this;
+			}
+		}
+		for (uint i = 0; i < cols_; ++i) {
+			crow = i;
+			while (smat.m[i][i] == Type(0)) {
+				++crow;
+				smat.swapRows(i, crow);
+				mtmp.swapRows(i, crow);
+				if (sv != 0) sv->swap(i, crow);
+			}
+			for (uint j = i + 1; j < rows_; ++j) {
+				mul = smat.m[i][j] / smat.m[i][i];
+				for (uint k = i; k < cols_; ++k) smat.m[k][j] -= mul * smat.m[k][i];
+				for (uint k = 0; k < cols_; ++k) mtmp.m[k][j] -= mul * mtmp.m[k][i];
+				if (sv != 0) (*sv)[j] -= mul * (*sv)[i];
+			}
+			//cout << i << endl << smat << endl;
+			if (i < cols_ - 1) {
+				if (fabs(smat.m[i+1][i+1]) < Type(1E-100)) {
+					if (ok != 0) *ok = false;
+					return *this;
+				}
+			}
+			iddiv = smat.m[i][i];
+			for (uint j = i; j < cols_; ++j) smat.m[j][i] /= iddiv;
+			for (uint j = 0; j < cols_; ++j) mtmp.m[j][i] /= iddiv;
+			if (sv != 0) (*sv)[i] /= iddiv;
+		}
+		for (uint i = cols_ - 1; i > 0; --i) {
+			for (uint j = 0; j < i; ++j) {
+				mul = smat.m[i][j];
+				smat.m[i][j] -= mul;
+				for (uint k = 0; k < cols_; ++k) mtmp.m[k][j] -= mtmp.m[k][i] * mul;
+				if (sv != 0) (*sv)[j] -= mul * (*sv)[i];
+			}
+		}
+		if (ok != 0) *ok = true;
+		m = mtmp.m;
+		return *this;
+	}
+	_CMatrix inverted(bool * ok = 0) {_CMatrix tm(*this); tm.invert(ok); return tm;}
+	_CMatrix transposed() {_CMatrix tm(rows_, cols_); PIMM_FOR_WB(c, r) tm[r][c] = m[c][r]; return tm;}
+
+private:
+	uint cols_, rows_;
+	PIVector<PIVector<Type> > m;
+
+};
+
+template<typename Type>
+inline std::ostream & operator <<(std::ostream & s, const PIMathMatrix<Type> & m) {s << '{'; for (uint r = 0; r < m.rows(); ++r) { for (uint c = 0; c < m.cols(); ++c) { s << m[c][r]; if (c < m.cols() - 1 || r < m.rows() - 1) s << ", ";} if (r < m.rows() - 1) s << endl << ' ';} s << '}'; return s;}
+template<typename Type>
+inline PICout operator <<(PICout s, const PIMathMatrix<Type> & m) {s << '{'; for (uint r = 0; r < m.rows(); ++r) { for (uint c = 0; c < m.cols(); ++c) { s << m[c][r]; if (c < m.cols() - 1 || r < m.rows() - 1) s << ", ";} if (r < m.rows() - 1) s << NewLine << ' ';} s << '}'; return s;}
+
+/// Multiply matrices {CR x Rows0} on {Cols1 x CR}, result is {Cols1 x Rows0}
+template<typename Type>
+inline PIMathMatrix<Type> operator *(const PIMathMatrix<Type> & fm,
+									const PIMathMatrix<Type> & sm) {
+	uint cr = fm.cols(), rows0 = fm.rows(), cols1 = sm.cols();
+	PIMathMatrix<Type> tm(cols1, rows0);
+	if (fm.cols() != sm.rows()) return tm;
+	Type t;
+	for (uint j = 0; j < rows0; ++j) {
+		for (uint i = 0; i < cols1; ++i) {
+			t = Type(0);
+			for (uint k = 0; k < cr; ++k)
+				t += fm[k][j] * sm[i][k];
+			tm[i][j] = t;
+		}
+	}
+	return tm;
+}
+
+/// Multiply matrix {Cols x Rows} on vector {Cols}, result is vector {Rows}
+template<typename Type>
+inline PIMathVector<Type> operator *(const PIMathMatrix<Type> & fm,
+									const PIMathVector<Type> & sv) {
+	uint c = fm.cols(), r = fm.rows();
+	PIMathVector<Type> tv(r);
+	if (c != sv.size()) return tv;
+	Type t;
+	for (uint i = 0; i < r; ++i) {
+		t = Type(0);
+		for (uint j = 0; j < c; ++j)
+			t += fm[j][i] * sv[j];
+		tv[i] = t;
+	}
+	return tv;
+}
+
+typedef PIMathMatrix<int> PIMathMatrixi;
+typedef PIMathMatrix<double> PIMathMatrixd;
+
+#undef PIMV_FOR
+#undef PIMM_FOR
+#undef PIMM_FOR_WB
+#undef PIMM_FOR_I
+#undef PIMM_FOR_I_WB
+#undef PIMM_FOR_C
+#undef PIMM_FOR_R
+
+#endif // PIMATHMATRIX_H

src/math/pimathsolver.cpp

@@ -0,0 +1,248 @@
+/*
+	PIP - Platform Independent Primitives
+	PIMathSolver
+	Copyright (C) 2014  Ivan Pelipenko peri4ko@gmail.com
+
+	This program is free software: you can redistribute it and/or modify
+	it under the terms of the GNU General Public License as published by
+	the Free Software Foundation, either version 3 of the License, or
+	(at your option) any later version.
+
+	This program is distributed in the hope that it will be useful,
+	but WITHOUT ANY WARRANTY; without even the implied warranty of
+	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+	GNU General Public License for more details.
+
+	You should have received a copy of the GNU General Public License
+	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "pimathsolver.h"
+
+const char PIMathSolver::methods_desc[] = "b{Methods:}\
+\n  -1 - Global settings\
+\n  01 - Eyler 1\
+\n  02 - Eyler 2\
+\n  14 - Runge-Kutta 4\
+\n  23 - Adams-Bashfort-Moulton 3\
+\n  24 - Adams-Bashfort-Moulton 4\
+\n  32 - Polynomial Approximation 2\
+\n  33 - Polynomial Approximation 3\
+\n  34 - Polynomial Approximation 4\
+\n  35 - Polynomial Approximation 5";
+
+PIMathSolver::Method PIMathSolver::method_global = PIMathSolver::Eyler_2;
+
+
+void PIMathSolver::solve(double u, double h) {
+	switch (method) {
+	case Global:
+		switch (method_global) {
+		case Eyler_1: solveEyler1(u, h); break;
+		case Eyler_2: solveEyler2(u, h); break;
+		case RungeKutta_4: solveRK4(u, h); break;
+		case AdamsBashfortMoulton_2: solveABM2(u, h); break;
+		case AdamsBashfortMoulton_3: solveABM3(u, h); break;
+		case AdamsBashfortMoulton_4: default: solveABM4(u, h); break;
+		case PolynomialApproximation_2: solvePA2(u, h); break;
+		case PolynomialApproximation_3: solvePA3(u, h); break;
+		case PolynomialApproximation_4: solvePA4(u, h); break;
+		case PolynomialApproximation_5: solvePA5(u, h); break;
+		}
+		break;
+	case Eyler_1: solveEyler1(u, h); break;
+	case Eyler_2: solveEyler2(u, h); break;
+	case RungeKutta_4: solveRK4(u, h); break;
+	case AdamsBashfortMoulton_2: solveABM2(u, h); break;
+	case AdamsBashfortMoulton_3: solveABM3(u, h); break;
+	case AdamsBashfortMoulton_4: default: solveABM4(u, h); break;
+	case PolynomialApproximation_2: solvePA2(u, h); break;
+	case PolynomialApproximation_3: solvePA3(u, h); break;
+	case PolynomialApproximation_4: solvePA4(u, h); break;
+	case PolynomialApproximation_5: solvePA5(u, h); break;
+	}
+	step++;
+}
+
+
+void PIMathSolver::fromTF(const TransferFunction & TF) {
+	if (TF.vector_An.size() >= TF.vector_Bm.size())
+		size = TF.vector_An.size()-1;
+	else {
+		piCout << "PIMathSolver error: {A} should be greater than {B}";
+		return;
+	}
+	if (size == 0) return;
+
+	step = 0;
+	times.fill(0.);
+	A.resize(size, size);
+	d.resize(size + 1); d.fill(0.);
+	a1.resize(size + 1); a1.fill(0.);
+	b1.resize(size + 1); b1.fill(0.);
+	X.resize(size); X.fill(0.);
+	F.resize(5);
+	for (uint i = 0; i < 5; ++i)
+		F[i].resize(size), F[i].fill(0.);
+	k1.resize(size); k1.fill(0.);
+	k2.resize(size); k2.fill(0.);
+	k3.resize(size); k3.fill(0.);
+	k4.resize(size); k4.fill(0.);
+	xx.resize(size); xx.fill(0.);
+	XX.resize(size); XX.fill(0.);
+	for (uint i = 0; i < size + 1; ++i)
+		a1[size - i] = TF.vector_An[i];
+	for (uint i = 0; i < TF.vector_Bm.size(); ++i)
+		b1[size - i] = TF.vector_Bm[i];
+	double a0 = a1[0];
+	a1 /= a0;
+	b1 /= a0;
+
+	d[0] = b1[0]; // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> d
+	for (uint i = 1; i < size + 1; ++i) {
+		sum = 0.;
+		for (uint m = 0; m < i; ++m)
+			sum += a1[i - m] * d[m];
+		d[i] = b1[i] - sum;
+	}
+
+	for (uint i = 0; i < size - 1; ++i) // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20>
+		for (uint j = 0; j < size; ++j)
+				A[j][i] = (j == i + 1);
+	for (uint i = 0; i < size; ++i)
+		A[i][size - 1] = -a1[size - i];
+	for (uint i = 0; i < size; ++i)
+		d[i] = d[i + 1];
+}
+
+
+void PIMathSolver::solveEyler1(double u, double h) {
+	/*for (uint i = 0; i < size; ++i) {
+	 *		sum = 0.;
+	 *		for (uint j = 0; j < size; ++j)
+	 *			sum += A[j][i] * X[j];
+	 *		xx[i] = sum + d[i] * u;
+	}*/
+	F[0] = A * X + d * u;
+	X += F[0] * h;
+	moveF();
+}
+
+
+void PIMathSolver::solveEyler2(double u, double h) {
+	F[0] = A * X + d * u;
+	X += (F[0] + F[1]) * h / 2.;
+	moveF();
+}
+
+
+void PIMathSolver::solveRK4(double u, double h) {
+	td = X[0] - F[0][0];
+	k1 = A * X + d * u; xx = k1 * h / 2.; XX = X + xx;
+	k2 = A * (XX + k1 * h / 2.) + d * (u + td/3.); xx = k2 * h / 2.; XX += xx;
+	k3 = A * (XX + k2 * h / 2.) + d * (u + td*2./3.); xx = k3 * h; XX += xx;
+	k4 = A * (XX + k3 * h) + d * (u + td);
+	//cout << k1 << k2 << k3 << k4 << endl;
+	X += (k1 + k2 * 2. + k3 * 2. + k4) * h / 6.;
+	F[0] = X;
+}
+
+
+void PIMathSolver::solveABM2(double u, double h) {
+	F[0] = A * X + d * u;
+	XX = X + (F[0] * 3. - F[1]) * (h / 2.);
+	F[1] = A * XX + d * u;
+	X += (F[1] + F[0]) * (h / 2.);
+	moveF();
+}
+
+
+void PIMathSolver::solveABM3(double u, double h) {
+	F[0] = A * X + d * u;
+	XX = X + (F[0] * 23. - F[1] * 16. + F[2] * 5.) * (h / 12.);
+	F[2] = A * XX + d * u;
+	X += (F[2] * 5. + F[0] * 8. - F[1]) * (h / 12.);
+	moveF();
+}
+
+
+void PIMathSolver::solveABM4(double u, double h) {
+	F[0] = A * X + d * u;
+	XX = X + (F[0] * 55. - F[1] * 59. + F[2] * 37. - F[3] * 9.) * (h / 24.);
+	F[3] = A * XX + d * u;
+	X += (F[3] * 9. + F[0] * 19. - F[1] * 5. + F[2]) * (h / 24.);
+	moveF();
+}
+
+
+void PIMathSolver::solvePA(double u, double h, uint deg) {
+	F[0] = A * X + d * u;
+	M.resize(deg, deg);
+	Y.resize(deg);
+	pY.resize(deg);
+
+	for (uint k = 0; k < size; ++k) {
+		for (uint i = 0; i < deg; ++i) {
+			td = 1.;
+			ct = times[i];
+			for (uint j = 0; j < deg; ++j) {
+				M[j][i] = td;
+				td *= ct;
+			}
+		}
+		for (uint i = 0; i < deg; ++i)
+			Y[i] = F[i][k];
+		/// find polynom
+		//if (step == 1) cout << M << endl << Y << endl;
+		M.invert(&ok, &Y);
+		//if (step == 1) cout << Y << endl;
+		if (!ok) {
+			solveEyler2(u, h);
+			break;
+		}
+		/// calc last piece
+		x0 = 0.;
+		td = 1.;
+		t = times[0];
+		for (uint i = 0; i < deg; ++i) {
+			x0 += Y[i] * td / (i + 1.);
+			td *= t;
+		}
+		x0 *= t;
+
+		x1 = 0.;
+		td = 1.;
+		t = times[1];
+		for (uint i = 0; i < deg; ++i) {
+			x1 += Y[i] * td / (i + 1.);
+			td *= t;
+		}
+		x1 *= t;
+		lp = x0 - x1;
+
+		if (deg > 2) {
+			/// calc prev piece
+			x0 = 0.;
+			td = 1.;
+			dh = times[1] - times[2];
+			if (dh != 0. && step > 1) {
+				t = times[2];
+				for (uint i = 0; i < deg; ++i) {
+					x0 += Y[i] * td / (i + 1.);
+					td *= t;
+				}
+				x0 *= t;
+				ct = x1 - x0;
+				/// calc correction
+				ct -= pY[k];
+			}
+			/// calc final
+			X[k] += lp - ct;
+			pY[k] = lp;
+		} else {
+			X[k] += lp;
+			pY[k] = lp;
+		}
+	}
+	moveF();
+}

src/math/pimathsolver.h

@@ -0,0 +1,94 @@
+/*! \file pimathsolver.h
+ * \brief PIMathSolver
+*/
+/*
+	PIP - Platform Independent Primitives
+	PIMathSolver
+	Copyright (C) 2014  Ivan Pelipenko peri4ko@gmail.com
+
+	This program is free software: you can redistribute it and/or modify
+	it under the terms of the GNU General Public License as published by
+	the Free Software Foundation, either version 3 of the License, or
+	(at your option) any later version.
+
+	This program is distributed in the hope that it will be useful,
+	but WITHOUT ANY WARRANTY; without even the implied warranty of
+	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+	GNU General Public License for more details.
+
+	You should have received a copy of the GNU General Public License
+	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef PIMATHSOLVER_H
+#define PIMATHSOLVER_H
+
+#include "pimathmatrix.h"
+
+/// Differential evaluations
+
+struct TransferFunction { // <20><><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
+	PIVector<double> vector_Bm, vector_An;
+};
+
+// <20><><EFBFBD><EFBFBD><EFBFBD>, <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
+// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD>. <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD>. <20><>-<2D><><EFBFBD>:
+//	<09><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
+//	<09><><EFBFBD><EFBFBD><EFBFBD>-<2D><><EFBFBD><EFBFBD><EFBFBD> 4-<2D><> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
+//	<09><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>-<2D><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>-<2D><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> 2, 3, 4 <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
+class PIP_EXPORT PIMathSolver
+{
+public:
+	enum Method {Global = -1,
+				Eyler_1 = 01,
+				Eyler_2 = 02,
+				EylerKoshi = 03,
+				RungeKutta_4 = 14,
+				AdamsBashfortMoulton_2 = 22,
+				AdamsBashfortMoulton_3 = 23,
+				AdamsBashfortMoulton_4 = 24,
+				PolynomialApproximation_2 = 32,
+				PolynomialApproximation_3 = 33,
+				PolynomialApproximation_4 = 34,
+				PolynomialApproximation_5 = 35
+				};
+
+	PIMathSolver() {times.resize(4); step = 0;}
+
+	void solve(double u, double h);
+	void fromTF(const TransferFunction & TF);
+	void setMethod(Method m) {method = m;}
+	void setTime(double time) {times.pop_back(); times.push_front(time);}
+
+	void solveEyler1(double u, double h);
+	void solveEyler2(double u, double h);
+	void solveRK4(double u, double h);
+	void solveABM2(double u, double h);
+	void solveABM3(double u, double h);
+	void solveABM4(double u, double h);
+	void solvePA(double u, double h, uint deg);
+	void solvePA2(double u, double h) {if (step > 0) solvePA(u, h, 2); else solveEyler1(u, h);}
+	void solvePA3(double u, double h) {if (step > 1) solvePA(u, h, 3); else solvePA2(u, h);}
+	void solvePA4(double u, double h) {if (step > 2) solvePA(u, h, 4); else solvePA3(u, h);}
+	void solvePA5(double u, double h) {if (step > 3) solvePA(u, h, 5); else solvePA4(u, h);}
+
+	PIMathVectord X;
+	static Method method_global;
+	static const char methods_desc[];
+
+private:
+	void moveF() {for (uint i = F.size() - 1; i > 0; --i) F[i] = F[i - 1];}
+
+	PIMathMatrixd A, M;
+	PIMathVectord d, a1, b1;
+	PIMathVectord k1, k2, k3, k4, xx;
+	PIMathVectord XX, Y, pY;
+	PIVector<PIMathVectord> F;
+	PIVector<double> times;
+	uint size, step;
+	Method method;
+	double sum, td, ct, lp, dh, t, x1, x0;
+	bool ok;
+};
+
+#endif // PIMATHSOLVER_H

src/math/pimathvector.h

@@ -0,0 +1,229 @@
+/*! \file pimathvector.h
+ * \brief PIMathVector
+*/
+/*
+	PIP - Platform Independent Primitives
+	PIMathVector
+	Copyright (C) 2014  Ivan Pelipenko peri4ko@gmail.com
+
+	This program is free software: you can redistribute it and/or modify
+	it under the terms of the GNU General Public License as published by
+	the Free Software Foundation, either version 3 of the License, or
+	(at your option) any later version.
+
+	This program is distributed in the hope that it will be useful,
+	but WITHOUT ANY WARRANTY; without even the implied warranty of
+	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+	GNU General Public License for more details.
+
+	You should have received a copy of the GNU General Public License
+	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef PIMATHVECTOR_H
+#define PIMATHVECTOR_H
+
+#include "pimathbase.h"
+
+template<uint Cols, uint Rows, typename Type>
+class PIMathMatrixT;
+
+
+/// Vector templated
+
+#define PIMV_FOR(v, s) for (uint v = s; v < Size; ++v)
+
+#pragma pack(push, 1)
+template<uint Size, typename Type = double>
+class PIP_EXPORT PIMathVectorT {
+	typedef PIMathVectorT<Size, Type> _CVector;
+public:
+	PIMathVectorT() {resize(Size);}
+	//PIMathVectorT(Type val) {resize(Size); PIMV_FOR(i, 0) c[i] = val;}
+	PIMathVectorT(Type fval, ...) {resize(Size); c[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] = va_arg(vl, Type); va_end(vl);}
+	PIMathVectorT(const PIVector<Type> & val) {resize(Size); PIMV_FOR(i, 0) c[i] = val[i];}
+	PIMathVectorT(const _CVector & st, const _CVector & fn) {resize(Size); set(st, fn);}
+
+	uint size() const {return Size;}
+	_CVector & fill(const Type & v) {PIMV_FOR(i, 0) c[i] = v; return *this;}
+	_CVector & set(Type fval, ...) {c[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] = va_arg(vl, Type); va_end(vl); return *this;}
+	_CVector & set(const _CVector & st, const _CVector & fn) {PIMV_FOR(i, 0) c[i] = fn[i] - st[i]; return *this;}
+	_CVector & move(const Type & v) {PIMV_FOR(i, 0) c[i] += v; return *this;}
+	_CVector & move(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i]; return *this;}
+	_CVector & move(Type fval, ...) {c[0] += fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] += va_arg(vl, Type); va_end(vl); return *this;}
+	Type lengthSqr() const {Type tv(0); PIMV_FOR(i, 0) tv += (c[i] * c[i]); return tv;}
+	Type length() const {return sqrt(lengthSqr());}
+	Type manhattanLength() const {Type tv(0); PIMV_FOR(i, 0) tv += fabs(c[i]); return tv;}
+	Type angleCos(const _CVector & v) const {Type tv = v.length() * length(); return (tv == Type(0) ? Type(0) : ((*this) ^ v) / tv);}
+	Type angleSin(const _CVector & v) const {Type tv = angleCos(v); return sqrt(Type(1) - tv * tv);}
+	Type angleRad(const _CVector & v) const {return acos(angleCos(v));}
+	Type angleDeg(const _CVector & v) const {return toDeg(acos(angleCos(v)));}
+	_CVector projection(const _CVector & v) {Type tv = v.length(); return (tv == Type(0) ? _CVector() : v * (((*this) ^ v) / tv));}
+	_CVector & normalize() {Type tv = length(); if (tv == Type(1)) return *this; if (piAbs<Type>(tv) <= Type(1E-100)) {fill(Type(0)); return *this;} PIMV_FOR(i, 0) c[i] /= tv; return *this;}
+	_CVector normalized() {_CVector tv(*this); tv.normalize(); return tv;}
+	bool isNull() const {PIMV_FOR(i, 0) if (c[i] != Type(0)) return false; return true;}
+	bool isOrtho(const _CVector & v) const {return ((*this) ^ v) == Type(0);}
+
+	Type & at(uint index) {return c[index];}
+	Type at(uint index) const {return c[index];}
+	Type & operator [](uint index) {return c[index];}
+	Type operator [](uint index) const {return c[index];}
+	_CVector & operator =(const _CVector & v) {memcpy(c, v.c, sizeof(Type) * Size); return *this;}
+	bool operator ==(const _CVector & v) const {PIMV_FOR(i, 0) if (c[i] != v[i]) return false; return true;}
+	bool operator !=(const _CVector & v) const {return !(*this == c);}
+	void operator +=(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i];}
+	void operator -=(const _CVector & v) {PIMV_FOR(i, 0) c[i] -= v[i];}
+	void operator *=(const Type & v) {PIMV_FOR(i, 0) c[i] *= v;}
+	void operator *=(const _CVector & v) {PIMV_FOR(i, 0) c[i] *= v[i];}
+	void operator /=(const Type & v) {PIMV_FOR(i, 0) c[i] /= v;}
+	void operator /=(const _CVector & v) {PIMV_FOR(i, 0) c[i] /= v[i];}
+	_CVector operator -() const {_CVector tv; PIMV_FOR(i, 0) tv[i] = -c[i]; return tv;}
+	_CVector operator +(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] += v[i]; return tv;}
+	_CVector operator -(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] -= v[i]; return tv;}
+	_CVector operator *(const Type & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] *= v; return tv;}
+	_CVector operator /(const Type & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] /= v; return tv;}
+	_CVector operator /(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] /= v[i]; return tv;}
+	_CVector operator *(const _CVector & v) const {if (Size > 3) return _CVector(); _CVector tv; tv.fill(Type(1)); tv[0] = c[1]*v[2] - v[1]*c[2]; tv[1] = v[0]*c[2] - c[0]*v[2]; tv[2] = c[0]*v[1] - v[0]*c[1]; return tv;}
+	Type operator ^(const _CVector & v) const {Type tv(0); PIMV_FOR(i, 0) tv += c[i] * v[i]; return tv;}
+	
+	PIMathMatrixT<1, Size, Type> transposed() const {
+		PIMathMatrixT<1, Size, Type> ret;
+		PIMV_FOR(i, 0) ret[0][i] = c[i];
+		return ret;
+	}
+
+	operator PIMathMatrixT<1, Size, Type>() {return transposed();}
+	operator PIMathMatrixT<Size, 1, Type>() {
+		PIMathMatrixT<Size, 1, Type> ret;
+		PIMV_FOR(i, 0) ret[i][0] = c[i];
+		return ret;
+	}
+	Type distToLine(const _CVector & lp0, const _CVector & lp1) {
+		_CVector a(lp0, lp1), b(lp0, *this), c(lp1, *this);
+		Type f = fabs(a[0]*b[1] - a[1]*b[0]) / a.length();//, s = b.length() + c.length() - a.length();
+		return f;}
+
+	template<uint Size1, typename Type1> /// vector {Size, Type} to vector {Size1, Type1}
+	PIMathVectorT<Size1, Type1> turnTo() {PIMathVectorT<Size1, Type1> tv; uint sz = piMin<uint>(Size, Size1); for (uint i = 0; i < sz; ++i) tv[i] = c[i]; return tv;}
+
+private:
+	void resize(uint size, const Type & new_value = Type()) {s = size; for (int i = 0; i < s; ++i) c[i] = new_value;}
+	
+	int s;
+	Type c[Size];
+
+};
+#pragma pack(pop)
+
+template<uint Size, typename Type>
+inline std::ostream & operator <<(std::ostream & s, const PIMathVectorT<Size, Type> & v) {s << '{'; PIMV_FOR(i, 0) {s << v[i]; if (i < Size - 1) s << ", ";} s << '}'; return s;}
+template<uint Size, typename Type>
+inline PICout operator <<(PICout s, const PIMathVectorT<Size, Type> & v) {s << '{'; PIMV_FOR(i, 0) {s << v[i]; if (i < Size - 1) s << ", ";} s << '}'; return s;}
+template<uint Size, typename Type>
+inline bool operator ||(const PIMathVectorT<Size, Type> & f, const PIMathVectorT<Size, Type> & s) {return (f * s).isNull();}
+template<uint Size, typename Type>
+inline PIMathVectorT<Size, Type> sqrt(const PIMathVectorT<Size, Type> & v) {PIMathVectorT<Size, Type> ret; PIMV_FOR(i, 0) {ret[i] = sqrt(v[i]);} return ret;}
+template<uint Size, typename Type>
+inline PIMathVectorT<Size, Type> sqr(const PIMathVectorT<Size, Type> & v) {PIMathVectorT<Size, Type> ret; PIMV_FOR(i, 0) {ret[i] = sqr(v[i]);} return ret;}
+
+template<uint Size, typename Type>
+inline PIByteArray & operator <<(PIByteArray & s, const PIMathVectorT<Size, Type> & v) {for (uint i = 0; i < Size; ++i) s << v[i]; return s;}
+template<uint Size, typename Type>
+inline PIByteArray & operator >>(PIByteArray & s, PIMathVectorT<Size, Type> & v) {for (uint i = 0; i < Size; ++i) s >> v[i]; return s;}
+
+//template<uint Size0, typename Type0 = double, uint Size1 = Size0, typename Type1 = Type0> /// vector {Size0, Type0} to vector {Size1, Type1}
+//inline operator PIMathVectorT<Size1, Type1>(const PIMathVectorT<Size0, Type0> & v) {PIMathVectorT<Size1, Type1> tv; uint sz = piMin<uint>(Size0, Size1); for (uint i = 0; i < sz; ++i) tv[i] = v[i]; return tv;}
+
+typedef PIMathVectorT<2u, int> PIMathVectorT2i;
+typedef PIMathVectorT<3u, int> PIMathVectorT3i;
+typedef PIMathVectorT<4u, int> PIMathVectorT4i;
+typedef PIMathVectorT<2u, double> PIMathVectorT2d;
+typedef PIMathVectorT<3u, double> PIMathVectorT3d;
+typedef PIMathVectorT<4u, double> PIMathVectorT4d;
+
+#undef PIMV_FOR
+
+/// Vector
+
+#define PIMV_FOR(v, s) for (uint v = s; v < size_; ++v)
+
+template<typename Type>
+class PIP_EXPORT PIMathVector {
+	typedef PIMathVector<Type> _CVector;
+public:
+	PIMathVector(const uint size = 3) {resize(size);}
+	PIMathVector(const uint size, Type fval, ...) {resize(size); c[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] = va_arg(vl, Type); va_end(vl);}
+	PIMathVector(const PIVector<Type> & val) {resize(val.size()); PIMV_FOR(i, 0) c[i] = val[i];}
+	PIMathVector(const _CVector & st, const _CVector & fn) {resize(st.size()); PIMV_FOR(i, 0) c[i] = fn[i] - st[i];}
+
+	uint size() const {return size_;}
+	_CVector & resize(uint size, const Type & new_value = Type()) {size_ = size; c.resize(size, new_value); return *this;}
+	_CVector resized(uint size, const Type & new_value = Type()) {_CVector tv = _CVector(*this); tv.resize(size, new_value); return tv;}
+	_CVector & fill(const Type & v) {PIMV_FOR(i, 0) c[i] = v; return *this;}
+	_CVector & set(Type fval, ...) {c[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] = va_arg(vl, Type); va_end(vl); return *this;}
+	_CVector & move(const Type & v) {PIMV_FOR(i, 0) c[i] += v; return *this;}
+	_CVector & move(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i]; return *this;}
+	_CVector & move(Type fval, ...) {c[0] += fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] += va_arg(vl, Type); va_end(vl); return *this;}
+	_CVector & swap(uint fe, uint se) {piSwap<Type>(c[fe], c[se]); return *this;}
+	Type lengthSqr() const {Type tv(0); PIMV_FOR(i, 0) tv += (c[i] * c[i]); return tv;}
+	Type length() const {return sqrt(lengthSqr());}
+	Type manhattanLength() const {Type tv(0); PIMV_FOR(i, 0) tv += fabs(c[i]); return tv;}
+	Type angleCos(const _CVector & v) const {Type tv = v.length() * length(); return (tv == Type(0) ? Type(0) : ((*this) ^ v) / tv);}
+	Type angleSin(const _CVector & v) const {Type tv = angleCos(v); return sqrt(Type(1) - tv * tv);}
+	Type angleRad(const _CVector & v) const {return acos(angleCos(v));}
+	Type angleDeg(const _CVector & v) const {return toDeg(acos(angleCos(v)));}
+	_CVector projection(const _CVector & v) {Type tv = v.length(); return (tv == Type(0) ? _CVector() : v * (((*this) ^ v) / tv));}
+	_CVector & normalize() {Type tv = length(); if (tv == Type(1)) return *this; if (piAbs<Type>(tv) <= Type(1E-100)) {fill(Type(0)); return *this;} PIMV_FOR(i, 0) c[i] /= tv; return *this;}
+	_CVector normalized() {_CVector tv(*this); tv.normalize(); return tv;}
+	bool isNull() const {PIMV_FOR(i, 0) if (c[i] != Type(0)) return false; return true;}
+	bool isOrtho(const _CVector & v) const {return ((*this) ^ v) == Type(0);}
+
+	Type & at(uint index) {return c[index];}
+	Type at(uint index) const {return c[index];}
+	Type & operator [](uint index) {return c[index];}
+	Type operator [](uint index) const {return c[index];}
+	void operator =(const _CVector & v) {c = v.c;}
+	bool operator ==(const _CVector & v) const {PIMV_FOR(i, 0) if (c[i] != v[i]) return false; return true;}
+	bool operator !=(const _CVector & v) const {return !(*this == c);}
+	void operator +=(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i];}
+	void operator -=(const _CVector & v) {PIMV_FOR(i, 0) c[i] -= v[i];}
+	void operator *=(const Type & v) {PIMV_FOR(i, 0) c[i] *= v;}
+	void operator *=(const _CVector & v) {PIMV_FOR(i, 0) c[i] *= v[i];}
+	void operator /=(const Type & v) {PIMV_FOR(i, 0) c[i] /= v;}
+	void operator /=(const _CVector & v) {PIMV_FOR(i, 0) c[i] /= v[i];}
+	_CVector operator -() {_CVector tv; PIMV_FOR(i, 0) tv[i] = -c[i]; return tv;}
+	_CVector operator +(const _CVector & v) {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] += v[i]; return tv;}
+	_CVector operator -(const _CVector & v) {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] -= v[i]; return tv;}
+	_CVector operator *(const Type & v) {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] *= v; return tv;}
+	_CVector operator /(const Type & v) {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] /= v; return tv;}
+	_CVector operator *(const _CVector & v) {if (size_ < 3) return _CVector(); _CVector tv; tv.fill(Type(1)); tv[0] = c[1]*v[2] - v[1]*c[2]; tv[1] = v[0]*c[2] - c[0]*v[2]; tv[2] = c[0]*v[1] - v[0]*c[1]; return tv;}
+	Type operator ^(const _CVector & v) const {Type tv(0); PIMV_FOR(i, 0) tv += c[i] * v[i]; return tv;}
+
+	//inline operator PIMathMatrix<1, Size, Type>() {return PIMathMatrix<1, Size, Type>(c);}
+	Type distToLine(const _CVector & lp0, const _CVector & lp1) {
+		_CVector a(lp0, lp1), b(lp0, *this), c(lp1, *this);
+		Type f = fabs(a[0]*b[1] - a[1]*b[0]) / a.length();//, s = b.length() + c.length() - a.length();
+		return f;}
+
+	template<typename Type1>
+	PIMathVector turnTo(uint size) {PIMathVector<Type1> tv; uint sz = piMin<uint>(size_, size); for (uint i = 0; i < sz; ++i) tv[i] = c[i]; return tv;}
+
+private:
+	uint size_;
+	PIVector<Type> c;
+
+};
+
+#undef PIMV_FOR
+
+template<typename Type>
+inline std::ostream & operator <<(std::ostream & s, const PIMathVector<Type> & v) {s << '{'; for (uint i = 0; i < v.size(); ++i) {s << v[i]; if (i < v.size() - 1) s << ", ";} s << '}'; return s;}
+template<typename Type>
+inline PICout operator <<(PICout s, const PIMathVector<Type> & v) {s << '{'; for (uint i = 0; i < v.size(); ++i) {s << v[i]; if (i < v.size() - 1) s << ", ";} s << '}'; return s;}
+
+typedef PIMathVector<int> PIMathVectori;
+typedef PIMathVector<double> PIMathVectord;
+
+//#include "pimathmatrix.h"
+
+#endif // PIMATHVECTOR_H

src/math/pistatistic.h

@@ -0,0 +1,89 @@
+/*! \file pistatistic.h
+ * \brief Class for calculating math statistic in values array
+*/
+/*
+	PIP - Platform Independent Primitives
+	Class for calculacing math statistic in values array
+	Copyright (C) 2014  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/>.
+*/
+
+#ifndef PISTATISTIC_H
+#define PISTATISTIC_H
+
+#include "pimathbase.h"
+
+template <typename T>
+class PIP_EXPORT PIStatistic {
+public:
+	PIStatistic() {mean = variance = skewness = kurtosis = T();}
+	
+	static T calculateMean(const PIVector<T> & val) {
+		T ret = T();
+		int n = val.size();
+		if (n < 1)
+			return ret;
+		for (int i = 0; i < n; i++)
+			ret += val[i];
+		return ret / n;
+	}
+	bool calculate(const PIVector<T> & val, const T & given_mean) {
+		T v = T(), v1 = T(), v2 = T(), stddev = T(), var = T();
+		int i, n = val.size();
+		if (n < 2)
+			return false;
+		mean = given_mean;
+		variance = skewness = kurtosis = T();
+		/*
+		* 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 / n;
+		var = (v1 / n - v2) / (n - 1);
+		if (var < T())
+			var = T();
+		stddev = sqrt(var);
+		/*
+		* Skewness and kurtosis
+		*/
+		if (stddev != T()) {
+			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;
+	}
+	bool calculate(const PIVector<T> & val) {return calculate(val, calculateMean(val));}
+	
+	T mean;
+	T variance;
+	T skewness;
+	T kurtosis;
+};
+
+typedef PIStatistic<int> PIStatistici;
+typedef PIStatistic<float> PIStatisticf;
+typedef PIStatistic<double> PIStatisticd;
+
+#endif // PISTATISTIC_H