version 1.22.0

source tree changed
detached PIConsole and PIScreen* in "pip_console" library
This commit is contained in:
2020-06-28 00:18:24 +03:00
parent 5de62b1c83
commit 42925122cb
231 changed files with 22981 additions and 22948 deletions

View File

@@ -0,0 +1,32 @@
/*! \file picompress.h
* \brief Compress class using zlib
*/
/*
PIP - Platform Independent Primitives
Compress class using zlib
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 Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PICOMPRESS_H
#define PICOMPRESS_H
#include "pibytearray.h"
PIByteArray piCompress(const PIByteArray & ba, int level = 6);
PIByteArray piDecompress(const PIByteArray & zba);
#endif // PICOMPRESS_H

243
lib/main/math/picrc.h Normal file
View File

@@ -0,0 +1,243 @@
/*! \file picrc.h
* \brief CRC checksum calculator
*/
/*
PIP - Platform Independent Primitives
CRC checksum calculator
Ivan Pelipenko peri4ko@yandex.ru
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#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];
};
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) {PIByteArray s(PIString(str).toByteArray()); return calculate(s.data(), s.size_s());}
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_16 standardCRC_16_Modbus() {return CRC_16(0x8005, 0xFFFF, 0xFFFF, false);}
inline CRC_8 standardCRC_8() {return CRC_8(0xD5, true, 0x0, 0x0);}
#endif // CRC_H

View File

@@ -0,0 +1,1209 @@
/*
PIP - Platform Independent Primitives
Evaluator designed for stream computing
Ivan Pelipenko peri4ko@yandex.ru
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "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) - arctangent
* * arcctg(x) - arccotangent
* * 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, a) - regular random with shift s and amp a
* * randomn(s, a) - normalize random with shift s and amp a
* * 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)
* * round(x) - round
*
* 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
*/
using namespace PIEvaluatorTypes;
PIEvaluatorContent::PIEvaluatorContent() {
addFunction("arcsin", 1);
addFunction("arccos", 1);
addFunction("arctg", 1);
addFunction("arcctg", 1);
addFunction("random", 2);
addFunction("randomn", 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);
addFunction("round", 1);
clearCustomVariables();
}
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() {
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<Variable>(variables[i], variables[j]);
}
}
}
BaseFunctions PIEvaluatorContent::getBaseFunction(const PIString & name) {
if (name == "sin") return bfSin;
if (name == "cos") return bfCos;
if (name == "tg") return bfTg;
if (name == "ctg") return bfCtg;
if (name == "arcsin") return bfArcsin;
if (name == "arccos") return bfArccos;
if (name == "arctg") return bfArctg;
if (name == "arcctg") return bfArcctg;
if (name == "exp") return bfExp;
if (name == "random") return bfRandom;
if (name == "randomn") return bfRandomn;
if (name == "sh") return bfSh;
if (name == "ch") return bfCh;
if (name == "th") return bfTh;
if (name == "cth") return bfCth;
if (name == "sqrt") return bfSqrt;
if (name == "sqr") return bfSqr;
if (name == "pow") return bfPow;
if (name == "abs") return bfAbs;
if (name == "ln") return bfLn;
if (name == "lg") return bfLg;
if (name == "log") return bfLog;
if (name == "im") return bfIm;
if (name == "re") return bfRe;
if (name == "arg") return bfArg;
if (name == "len") return bfLen;
if (name == "conj") return bfConj;
if (name == "sign") return bfSign;
if (name == "rad") return bfRad;
if (name == "deg") return bfDeg;
if (name == "j0") return bfJ0;
if (name == "j1") return bfJ1;
if (name == "jn") return bfJN;
if (name == "y0") return bfY0;
if (name == "y1") return bfY1;
if (name == "yn") return bfYN;
if (name == "min") return bfMin;
if (name == "max") return bfMax;
if (name == "clamp") return bfClamp;
if (name == "step") return bfStep;
if (name == "mix") return bfMix;
if (name == "defined") return bfDefined;
if (name == "round") return bfRound;
return 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("||", "|");
currentString.replaceAll(PIString::fromUTF8(""), ":");
currentString.replaceAll(PIString::fromUTF8(""), "}");
currentString.replaceAll(PIString::fromUTF8(""), "{");
currentString.replaceAll(PIString::fromUTF8(""), "&");
currentString.replaceAll(PIString::fromUTF8(""), "|");
}
void PIEvaluator::makeOutput(PIString & string) {
string.replaceAll(":", PIString::fromUTF8(""));
string.replaceAll("}", PIString::fromUTF8(""));
string.replaceAll("{", PIString::fromUTF8(""));
string.replaceAll("&", PIString::fromUTF8(""));
string.replaceAll("|", PIString::fromUTF8(""));
}
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 = etVariable;
elements[i].var_num = -666;
}
currentVariables.clear();
for (int i = 0; i < content.functionsCount(); i++) {
curfind = content.function(i).identifier;
cfunc = i;
flen = curfind.length();
fstart = 0;
while (fstart >= 0) {
fstart = tmps.find(curfind, fstart);
if (fstart < 0) break;
if (tmps[fstart + flen] != '(') {
fstart++;
continue;
}
for (int j = fstart; j < fstart + flen; j++) {
elements[j].set(etFunction, cnum, cfunc);
tmps.replace(j, 1, fc);
}
cnum++;
}
}
cnum = 0;
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(etVariable, cnum, i);
tmps.replace(j, 1, fc);
}
cnum++;
}
}
curfind = "";
cnum = 1;
for (int i = 0; i < tmps.length(); i++) {
cc = tmps[i];
switch (cpart) {
case 0:
if ((cc >= '0' && cc <= '9')) {
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) {
currentVariables.push_back(Variable("tmp" + PIString::fromNumber(cnum), curfind.toDouble()));
for (int j = i - curfind.length(); j < i; j++) {
elements[j].set(etNumber, cnum, -cnum);
tmps.replace(j, 1, fc);
}
curfind = "";
cnum++;
cpart = 0;
numFound = false;
}
}
if (cpart > 0) {
currentVariables.push_back(Variable("tmp" + PIString::fromNumber(cnum), curfind.toDouble()));
for (int j = tmps.length() - curfind.length(); j < tmps.length(); j++) {
elements[j].set(etNumber, cnum, -cnum);
tmps.replace(j, 1, fc);
}
}
cc = nc = fc;
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(etOperator, -1);
continue;
}
if (cc == '-' || cc == '+') {
elements[i].set(etOperator, -1);
if (i < tmps.length() - 1) if (elements[i + 1].type == etVariable ||
elements[i + 1].type == etFunction) continue;
if ((pc == '(' || isSign(pc) || i == 0) && i < tmps.length() - 1) {
if (elements[i + 1].type != etOperator) {
cnum = elements[i + 1].num;
elements[i].set(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;
continue;
}
}
}
if (isSign(cc)) {
elements[i].set(etOperator, -1);
continue;
}
}
/*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;
}
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 (fc == ',' || sc == ',') continue;
if (elements[i].type == etOperator && elements[ni].type == etOperator) continue;
if (fc == ')' && (elements[ni].type == etNumber || elements[ni].type == etVariable || elements[ni].type == etFunction)) needInsert = 1;
if (sc == '(' && (elements[i].type == etNumber || elements[i].type == etVariable)) needInsert = 1;
if (elements[i].type == etNumber && elements[ni].type == etNumber && elements[i].num != elements[ni].num) needInsert = 1;
if (elements[i].type == etVariable && elements[ni].type == etVariable && elements[i].num != elements[ni].num) needInsert = 1;
if ((elements[i].type == etNumber && elements[ni].type == etVariable) || (elements[i].type == etVariable && elements[ni].type == etNumber)) needInsert = 1;
if ((elements[i].type == etNumber || elements[i].type == etVariable) && elements[ni].type == etFunction) needInsert = 1;
if (elements[i].type == etFunction && elements[ni].type == etFunction && elements[i].num != elements[ni].num) needInsert = 2;
if (elements[i].type == etFunction && elements[ni].type != etFunction && sc != '(') needInsert = 2;
if (elements[pi].type == etOperator && (elements[ni].type == etFunction || elements[ni].type == etVariable) && fc == '-') needInsert = 3;
switch (needInsert) {
case 1:
currentString.insert(ni + inserted, "*");
elements.insert(ni + inserted, Element(etOperator, -1));
return true;
case 3:
currentString.insert(ni + inserted, "1*");
elements.insert(ni + inserted, Element(etOperator, -1));
return true;
}
}
return false;
}
void PIEvaluator::convert() {
int j;
Element ce, pe;
for (int i = 0; i < currentString.length(); i++) {
pe = elements[i];
if (pe.type != 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);
}
for (int i = 0; i < currentString.length(); i++) {
pe = elements[i];
if (pe.type != 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);
}
for (int i = 0; i < currentString.length(); i++) {
pe = elements[i];
if (pe.type != 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);
}
/*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; */
}
PIString PIEvaluator::preprocess(const PIString & string) {
PIString ret;
int lind;
ret = prepare(string);
convert();
instructions.clear();
variables = currentVariables;
lind = parse(currentString);
if (instructions.size() == 0) {
variables.push_back(Variable());
instructions.push_back(Instruction(oNone, PIVector<short>(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;
}
int PIEvaluator::parse(const PIString & string, int offset) {
int slen = string.length(), /*facnt,*/ farg, bcnt, k;
PIChar cc;
Element ce;
Function cfunc;
PIString sbrackets, carg;
PIVector<short> args, atmp;
PIVector<Operation> opers;
//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 etNumber:
args.push_back(ce.var_num);
continue;
case etVariable:
args.push_back(ce.var_num);
continue;
case etFunction:
i++;
cfunc = content.function(ce.var_num);
atmp.clear();
bcnt = farg = 1;
carg = "";
k = i + 1;
while (bcnt > 0) {
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(Variable());
farg = -variables.size_s();
}
instructions.push_back(Instruction(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;
continue;
case etOperator:
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(oAdd); continue;}
if (cc == '-') {opers.push_back(oSubtract); continue;}
if (cc == '*') {opers.push_back(oMultiply); continue;}
if (cc == '/') {opers.push_back(oDivide); continue;}
if (cc == '%') {opers.push_back(oResidue); continue;}
if (cc == '^') {opers.push_back(oPower); continue;}
if (cc == '=') {opers.push_back(oEqual); continue;}
if (cc == ':') {opers.push_back(oNotEqual); continue;}
if (cc == '}') {opers.push_back(oGreaterEqual); continue;}
if (cc == '{') {opers.push_back(oSmallerEqual); continue;}
if (cc == '>') {opers.push_back(oGreater); continue;}
if (cc == '<') {opers.push_back(oSmaller); continue;}
if (cc == '&') {opers.push_back(oAnd); continue;}
if (cc == '|') {opers.push_back(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;
}
int oprior = -1;
PIVector<Operation> opv;
while(1) {
operationsByPriority(++oprior, opv);
if (opv.isEmpty()) break;
for (int j = 0; j < opers.size_s(); j++) {
if (!opv.contains(opers[j])) 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(Variable());
farg = -variables.size_s();
}
}
instructions.push_back(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() {
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];
Function cf;
int fac, gac;
switch (ci.operation) {
case oNone: break;
case 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) {
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 Operation & operation) {
switch (operation) {
case oAdd: return "+";
case oSubtract: return "-";
case oMultiply: return "*";
case oDivide: return "/";
case oPower: return "^";
case oResidue: return "%";
case oEqual: return "=";
case oNotEqual: return PIString::fromUTF8("");
case oGreaterEqual: return PIString::fromUTF8("");
case oSmallerEqual: return PIString::fromUTF8("");
case oGreater: return ">";
case oSmaller: return "<";
case oAnd: return PIString::fromUTF8("");
case oOr: return PIString::fromUTF8("");
default: return "???";
}
}
void PIEvaluator::operationsByPriority(int p, PIVector<Operation> & ret) {
ret.clear();
switch (p) {
case 0: ret << oPower; break;
case 1: ret << oMultiply << oDivide << oResidue; break;
case 2: ret << oAdd << oSubtract; break;
case 3: ret << oEqual << oNotEqual << oGreaterEqual << oSmallerEqual << oGreater << oSmaller; break;
case 4: ret << oAnd; break;
case 5: ret << oOr; break;
default: break;
}
}
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 Instruction & ci) {
const Function & cfunc(content.function(ci.function));
int oi = -ci.out - 1;
complexd tmp, stmp, ttmp;
switch (cfunc.type) {
case bfSin:
tmpvars[oi].value = sin(value(ci.operators[0]));
break;
case bfCos:
tmpvars[oi].value = cos(value(ci.operators[0]));
break;
case bfTg:
tmpvars[oi].value = tan(value(ci.operators[0]));
break;
case bfCtg:
tmp = tan(value(ci.operators[0]));
if (tmp == complexd_0) tmpvars[oi].value = 0.;
else tmpvars[oi].value = complexd_1 / tmp;
break;
case bfArcsin:
tmpvars[oi].value = asinc(value(ci.operators[0]));
break;
case bfArccos:
tmpvars[oi].value = acosc(value(ci.operators[0]));
break;
case bfArctg:
tmpvars[oi].value = atanc(value(ci.operators[0]));
break;
case bfArcctg:
tmpvars[oi].value = atanc(-value(ci.operators[0])) + M_PI_2;
break;
case bfSh:
tmpvars[oi].value = sinh(value(ci.operators[0]));
break;
case bfCh:
tmpvars[oi].value = cosh(value(ci.operators[0]));
break;
case bfTh:
tmpvars[oi].value = tanh(value(ci.operators[0]));
break;
case bfCth:
tmp = tanh(value(ci.operators[0]));
if (tmp == complexd_0) tmpvars[oi].value = 0.;
else tmpvars[oi].value = complexd_1 / tmp;
break;
case bfAbs:
tmpvars[oi].value = abs(value(ci.operators[0]));
break;
case bfSqrt:
tmpvars[oi].value = sqrt(value(ci.operators[0]));
break;
case bfSqr:
tmpvars[oi].value = value(ci.operators[0]) * value(ci.operators[0]);
break;
case bfExp:
tmpvars[oi].value = exp(value(ci.operators[0]));
break;
case bfPow:
tmpvars[oi].value = pow(value(ci.operators[0]), value(ci.operators[1]));
break;
case bfLn:
tmpvars[oi].value = log(value(ci.operators[0]));
break;
case bfLg:
tmpvars[oi].value = log10(value(ci.operators[0]));
break;
case 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 bfRe:
tmpvars[oi].value = value(ci.operators[0]).real();
break;
case bfIm:
tmpvars[oi].value = value(ci.operators[0]).imag();
break;
case bfArg:
tmpvars[oi].value = arg(value(ci.operators[0]));
break;
case bfLen:
tmpvars[oi].value = abs(value(ci.operators[0]));
break;
case bfConj:
tmpvars[oi].value = conj(value(ci.operators[0]));
break;
case bfSign:
tmpvars[oi].value = value(ci.operators[0]).real() >= 0. ? complexd_1 : -complexd_1;
break;
case bfRad:
tmpvars[oi].value = value(ci.operators[0]) * complexd(deg2rad, 0.);
break;
case bfDeg:
tmpvars[oi].value = value(ci.operators[0]) * complexd(rad2deg, 0.);
break;
case bfJ0:
tmpvars[oi].value = piJ0(value(ci.operators[0]).real());
break;
case bfJ1:
tmpvars[oi].value = piJ1(value(ci.operators[0]).real());
break;
case bfJN:
tmpvars[oi].value = piJn(piRoundd(value(ci.operators[1]).real()), value(ci.operators[0]).real());
break;
case bfY0:
tmpvars[oi].value = piY0(value(ci.operators[0]).real());
break;
case bfY1:
tmpvars[oi].value = piY1(value(ci.operators[0]).real());
break;
case bfYN:
tmpvars[oi].value = piYn(piRoundd(value(ci.operators[1]).real()), value(ci.operators[0]).real());
break;
case 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 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 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 bfStep:
tmpvars[oi].value = complexd(value(ci.operators[0]).real() >= value(ci.operators[1]).real() ? complexld_1 : complexld_0);
break;
case 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 bfDefined:
tmpvars[oi].value = value(ci.operators[0]).real() > 0. ? complexd_1 : complexd_0;
break;
case bfRandom:
tmpvars[oi].value = value(ci.operators[0]) + randomd() * value(ci.operators[1]);
break;
case bfRandomn:
tmpvars[oi].value = randomn(value(ci.operators[0]).real(), value(ci.operators[1]).real());
break;
case bfRound:
tmpvars[oi].value = piRoundd(value(ci.operators[0]).real());
break;
default: break;
}
}
inline bool PIEvaluator::execInstructions() {
kvars = &(content.variables);
int oi;
complexd tmp;
tmpvars = variables;
//cout << "var count " << tmpvars.size_s() << endl;
for (int i = 0; i < instructions.size_s(); i++) {
const Instruction & 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 oAdd:
tmpvars[oi].value = value(ci.operators[0]) + value(ci.operators[1]);
break;
case oSubtract:
tmpvars[oi].value = value(ci.operators[0]) - value(ci.operators[1]);
break;
case oMultiply:
tmpvars[oi].value = value(ci.operators[0]) * value(ci.operators[1]);
break;
case 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 oResidue:
tmpvars[oi].value = residue(value(ci.operators[0]), value(ci.operators[1]));
break;
case oPower:
tmpvars[oi].value = pow(value(ci.operators[0]), value(ci.operators[1]));
break;
case oEqual:
tmpvars[oi].value = value(ci.operators[0]) == value(ci.operators[1]);
break;
case oNotEqual:
tmpvars[oi].value = value(ci.operators[0]) != value(ci.operators[1]);
break;
case oGreaterEqual:
tmpvars[oi].value = value(ci.operators[0]).real() >= value(ci.operators[1]).real();
break;
case oSmallerEqual:
tmpvars[oi].value = value(ci.operators[0]).real() <= value(ci.operators[1]).real();
break;
case oGreater:
tmpvars[oi].value = value(ci.operators[0]).real() > value(ci.operators[1]).real();
break;
case oSmaller:
tmpvars[oi].value = value(ci.operators[0]).real() < value(ci.operators[1]).real();
break;
case oAnd:
tmpvars[oi].value = value(ci.operators[0]).real() > 0. && value(ci.operators[1]).real() > 0.;
break;
case oOr:
tmpvars[oi].value = value(ci.operators[0]).real() > 0. || value(ci.operators[1]).real() > 0.;
break;
case oFunction:
execFunction(ci);
break;
case oNone:
tmpvars[oi].value = value(ci.operators[0]);
break;
default: 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;
}
PIByteArray PIEvaluator::save() const {
PIByteArray ret;
ret << content << currentVariables << unknownVars << instructions << currentString << lastError << out << correct;
return ret;
}
void PIEvaluator::load(PIByteArray ba) {
if (ba.size() <= 4) return;
ba >> content >> currentVariables >> unknownVars >> instructions >> currentString >> lastError >> out >> correct;
variables = currentVariables;
}

266
lib/main/math/pievaluator.h Normal file
View File

@@ -0,0 +1,266 @@
/*! \file pievaluator.h
* \brief Mathematic expressions calculator
*/
/*
PIP - Platform Independent Primitives
Evaluator designed for stream calculations
Ivan Pelipenko peri4ko@yandex.ru
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PIEVALUATOR_H
#define PIEVALUATOR_H
#include "pistringlist.h"
#include "pimathcomplex.h"
typedef complexd (*FuncFunc)(void * , int, complexd * );
namespace PIEvaluatorTypes {
enum PIP_EXPORT eType {etNumber, etOperator, etVariable, etFunction};
enum PIP_EXPORT Operation {oNone, oAdd, oSubtract, oMultiply, oDivide, oResidue, oPower,
oEqual, oNotEqual, oGreater, oSmaller, oGreaterEqual, oSmallerEqual,
oAnd, oOr, oFunction
};
enum PIP_EXPORT BaseFunctions {bfUnknown, bfSin, bfCos, bfTg, bfCtg,
bfArcsin, bfArccos, bfArctg, bfArcctg,
bfExp, bfRandom, bfRandomn,
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,
bfRound,
bfCustom = 0xFFFF
};
struct PIP_EXPORT Instruction {
Instruction() {out = -1; function = -1; operation = oNone;}
Instruction(Operation oper, PIVector<short> opers, short out_ind, short func = -1) {
operation = oper; operators = opers; out = out_ind; function = func;}
Operation operation;
short out;
short function;
PIVector<short> operators;
};
struct PIP_EXPORT Element {
Element() {num = 0; var_num = -1; type = etNumber;}
Element(eType new_type, short new_num, short new_var_num = -1) {set(new_type, new_num, new_var_num);}
void set(eType new_type, short new_num, short new_var_num = -1) {type = new_type; num = new_num; var_num = new_var_num;}
eType type;
short num;
short var_num;
};
struct PIP_EXPORT Function {
Function() {arguments = 0; type = bfUnknown; handler = 0;}
Function(const PIString & name, short args, BaseFunctions ftype) {identifier = name; arguments = args; type = ftype; handler = 0;}
Function(const PIString & name, short args, FuncFunc h) {identifier = name; arguments = args; type = bfCustom; handler = h;}
PIString identifier;
BaseFunctions type;
FuncFunc handler;
short arguments;
};
struct PIP_EXPORT 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;
friend inline PIByteArray & operator <<(PIByteArray & s, const PIEvaluatorContent & v);
friend inline PIByteArray & operator >>(PIByteArray & s, PIEvaluatorContent & v);
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;}
//! 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;}
//! Save to %PIByteArray evaluator state (expression, variables, errors, compiled instructions)
PIByteArray save() const;
//! Restore from %PIByteArray evaluator state (expression, variables, errors, compiled instructions)
void load(PIByteArray ba);
PIEvaluatorContent content;
private:
const PIString & prepare(const PIString & string);
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);
void operationsByPriority(int p, PIVector<PIEvaluatorTypes::Operation> & ret);
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);
PIDeque<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);}
inline PIByteArray & operator <<(PIByteArray & s, const PIEvaluatorTypes::Instruction & v) {
s << PIByteArray::RawData(&v, sizeof(v) - sizeof(v.operators)) << v.operators;
return s;
}
inline PIByteArray & operator >>(PIByteArray & s, PIEvaluatorTypes::Instruction & v) {
s >> PIByteArray::RawData(&v, sizeof(v) - sizeof(v.operators)) >> v.operators;
return s;
}
inline PIByteArray & operator <<(PIByteArray & s, const PIEvaluatorTypes::Element & v) {
s << PIByteArray::RawData(&v, sizeof(v));
return s;
}
inline PIByteArray & operator >>(PIByteArray & s, PIEvaluatorTypes::Element & v) {
s >> PIByteArray::RawData(&v, sizeof(v));
return s;
}
inline PIByteArray & operator <<(PIByteArray & s, const PIEvaluatorTypes::Variable & v) {
s << v.name << v.value;
return s;
}
inline PIByteArray & operator >>(PIByteArray & s, PIEvaluatorTypes::Variable & v) {
s >> v.name >> v.value;
return s;
}
inline PIByteArray & operator <<(PIByteArray & s, const PIEvaluatorContent & v) {
s << v.variables << v.cv_count;
return s;
}
inline PIByteArray & operator >>(PIByteArray & s, PIEvaluatorContent & v) {
s >> v.variables >> v.cv_count;
return s;
}
#endif // PIEVALUATOR_H

1893
lib/main/math/pifft.cpp Normal file
View File

@@ -0,0 +1,1893 @@
/*
PIP - Platform Independent Primitives
Class for FFT, IFFT and Hilbert transformations
Andrey Bychkov work.a.b@yandex.ru, Ivan Pelipenko peri4ko@yandex.ru
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "pifft.h"
PIFFT_double::PIFFT_double() {
}
PIVector<complexd> * PIFFT_double::calcFFT(const PIVector<complexd> & val) {
result.clear();
if (val.size_s() < 4) return &result;
fftc1d(val, val.size());
return &result;
}
PIVector<complexd> * PIFFT_double::calcFFTinverse(const PIVector<complexd> & val) {
result.clear();
if (val.size_s() < 4) return &result;
fftc1dinv(val, val.size());
return &result;
}
PIVector<complexd> * PIFFT_double::calcHilbert(const PIVector<double> & val) {
result.clear();
if (val.size_s() < 4) return &result;
fftc1r(val, val.size());
for (uint i = 0; i < result.size() / 2; i++) result[i] = result[i] * 2.;
for (uint i = result.size() / 2; i < result.size(); i++) result[i] = 0;
fftc1dinv(result, result.size());
return &result;
}
PIVector<complexd> * PIFFT_double::calcFFT(const PIVector<double> & val) {
result.clear();
if (val.size_s() < 4) return &result;
fftc1r(val, val.size());
return &result;
}
PIVector<double> PIFFT_double::getAmplitude() const {
PIVector<double> ret;
ret.resize(result.size());
double tmp;
for (uint i = 0; i < result.size(); i++) {
tmp = sqrt(result[i].real() * result[i].real() + result[i].imag() * result[i].imag());
ret[i] = tmp;
}
return ret;
}
PIVector<double> PIFFT_double::getReal() const {
PIVector<double> ret;
ret.resize(result.size());
for (uint i = 0; i < result.size(); i++)
ret[i] = result[i].real();
return ret;
}
PIVector<double> PIFFT_double::getImag() const {
PIVector<double> ret;
ret.resize(result.size());
for (uint i = 0; i < result.size(); i++)
ret[i] = result[i].imag();
return ret;
}
void PIFFT_double::fftc1d(const PIVector<complexd> & a, uint n) {
createPlan(n);
uint i;
PIVector<double> buf;
buf.resize(2 * n);
for (i = 0; i < n; i++) {
buf[2 * i + 0] = a[i].real();
buf[2 * i + 1] = a[i].imag();
}
ftbaseexecuteplan(&buf, 0, n, &curplan);
result.resize(n);
for (i = 0; i < n; i++)
result[i] = complexd(buf[2 * i + 0], buf[2 * i + 1]);
}
void PIFFT_double::fftc1r(const PIVector<double> & a, uint n) {
uint i;
if (n % 2 == 0) {
PIVector<double> buf;
uint n2 = n / 2;
buf = a;
createPlan(n2);
ftbaseexecuteplan(&buf, 0, n2, &curplan);
result.resize(n);
uint idx;
complexd hn, hmnc, v;
for (i = 0; i <= n2; i++) {
idx = 2 * (i % n2);
hn = complexd(buf[idx + 0], buf[idx + 1]);
idx = 2 * ((n2 - i) % n2);
hmnc = complexd(buf[idx + 0], -buf[idx + 1]);
v = complexd(sin(M_PI * i / n2), cos(M_PI * i / n2));
result[i] = ((hn + hmnc) - (v * (hn - hmnc)));
result[i] *= 0.5;
}
for (i = n2 + 1; i < n; i++)
result[i] = conj(result[n - i]);
} else {
PIVector<complexd> cbuf;
cbuf.resize(n);
for (i = 0; i < n; i++)
cbuf[i] = complexd(a[i], 0.);
fftc1d(cbuf, n);
}
}
void PIFFT_double::fftc1dinv(const PIVector<complexd> & a, uint n) {
PIVector<complexd> cbuf;
cbuf.resize(n);
uint i;
for (i = 0; i < n; i++) {
cbuf[i] = conj(a[i]);
}
fftc1d(cbuf, n);
for (i = 0; i < n; i++) {
result[i] = conj(result[i] / (double)n);
}
}
void PIFFT_double::createPlan(uint n) {
curplan.plan.clear();
curplan.precomputed.clear();
curplan.stackbuf.clear();
curplan.tmpbuf.clear();
if (n < 2) return;
ftbasegeneratecomplexfftplan(n, &curplan);
}
void PIFFT_double::ftbasegeneratecomplexfftplan(uint n, ftplan * plan) {
int planarraysize;
int plansize;
int precomputedsize;
int tmpmemsize;
int stackmemsize;
ae_int_t stackptr;
planarraysize = 1;
plansize = 0;
precomputedsize = 0;
stackmemsize = 0;
stackptr = 0;
tmpmemsize = 2 * n;
curplan.plan.resize(planarraysize);
int ftbase_ftbasecffttask = 0;
ftbase_ftbasegenerateplanrec(n, ftbase_ftbasecffttask, plan, &plansize, &precomputedsize, &planarraysize, &tmpmemsize, &stackmemsize, stackptr);
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_double::ftbase_ftbasegenerateplanrec(
int n,
int tasktype,
ftplan * plan,
int * plansize,
int * precomputedsize,
int * planarraysize,
int * tmpmemsize,
int * stackmemsize,
ae_int_t stackptr, int debugi) {
int k, m, n1, n2, esize, entryoffset;
int ftbase_ftbaseplanentrysize = 8;
int ftbase_ftbasecffttask = 0;
int ftbase_fftcooleytukeyplan = 0;
int ftbase_fftbluesteinplan = 1;
int ftbase_fftcodeletplan = 2;
int ftbase_fftrealcooleytukeyplan = 5;
int ftbase_fftemptyplan = 6;
if (*plansize + ftbase_ftbaseplanentrysize > (*planarraysize)) {
curplan.plan.resize(8 * (*planarraysize));
*planarraysize = 8 * (*planarraysize);
}
entryoffset = *plansize;
esize = ftbase_ftbaseplanentrysize;
*plansize = *plansize + esize;
if (n == 1) {
curplan.plan[entryoffset + 0] = esize;
curplan.plan[entryoffset + 1] = -1;
curplan.plan[entryoffset + 2] = -1;
curplan.plan[entryoffset + 3] = ftbase_fftemptyplan;
curplan.plan[entryoffset + 4] = -1;
curplan.plan[entryoffset + 5] = -1;
curplan.plan[entryoffset + 6] = -1;
curplan.plan[entryoffset + 7] = -1;
return;
}
ftbasefactorize(n, &n1, &n2);
if (n1 != 1) {
*tmpmemsize = piMax<int>(*tmpmemsize, 2 * n1 * n2);
curplan.plan[entryoffset + 0] = esize;
curplan.plan[entryoffset + 1] = n1;
curplan.plan[entryoffset + 2] = n2;
if (tasktype == ftbase_ftbasecffttask)
curplan.plan[entryoffset + 3] = ftbase_fftcooleytukeyplan;
else
curplan.plan[entryoffset + 3] = ftbase_fftrealcooleytukeyplan;
curplan.plan[entryoffset + 4] = 0;
curplan.plan[entryoffset + 5] = *plansize;
debugi++;
ftbase_ftbasegenerateplanrec(n1, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr, debugi);
curplan.plan[entryoffset + 6] = *plansize;
ftbase_ftbasegenerateplanrec(n2, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr, debugi);
curplan.plan[entryoffset + 7] = -1;
return;
} else {
if (n >= 2 && n <= 5) {
curplan.plan[entryoffset + 0] = esize;
curplan.plan[entryoffset + 1] = n1;
curplan.plan[entryoffset + 2] = n2;
curplan.plan[entryoffset + 3] = ftbase_fftcodeletplan;
curplan.plan[entryoffset + 4] = 0;
curplan.plan[entryoffset + 5] = -1;
curplan.plan[entryoffset + 6] = -1;
curplan.plan[entryoffset + 7] = *precomputedsize;
if (n == 3)
*precomputedsize = *precomputedsize + 2;
if (n == 5)
*precomputedsize = *precomputedsize + 5;
return;
} else {
k = 2 * n2 - 1;
m = ftbasefindsmooth(k);
*tmpmemsize = piMax<int>(*tmpmemsize, 2 * m);
curplan.plan[entryoffset + 0] = esize;
curplan.plan[entryoffset + 1] = n2;
curplan.plan[entryoffset + 2] = -1;
curplan.plan[entryoffset + 3] = ftbase_fftbluesteinplan;
curplan.plan[entryoffset + 4] = m;
curplan.plan[entryoffset + 5] = *plansize;
stackptr = stackptr + 2 * 2 * m;
*stackmemsize = piMax<int>(*stackmemsize, stackptr);
ftbase_ftbasegenerateplanrec(m, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr);
stackptr = stackptr - 2 * 2 * m;
curplan.plan[entryoffset + 6] = -1;
curplan.plan[entryoffset + 7] = *precomputedsize;
*precomputedsize = *precomputedsize + 2 * m + 2 * n;
return;
}
}
}
/*************************************************************************
Recurrent subroutine for precomputing FFT plans
-- ALGLIB --
Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT_double::ftbase_ftbaseprecomputeplanrec(ftplan * plan,
int entryoffset,
ae_int_t stackptr) {
int n1, n2, n, m, offs;
double v, bx, by;
int ftbase_fftcooleytukeyplan = 0;
int ftbase_fftbluesteinplan = 1;
int ftbase_fftcodeletplan = 2;
int ftbase_fhtcooleytukeyplan = 3;
int ftbase_fhtcodeletplan = 4;
int ftbase_fftrealcooleytukeyplan = 5;
if ((curplan.plan[entryoffset + 3] == ftbase_fftcooleytukeyplan || curplan.plan[entryoffset + 3] == ftbase_fftrealcooleytukeyplan) || curplan.plan[entryoffset + 3] == ftbase_fhtcooleytukeyplan) {
ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset + 5], stackptr);
ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset + 6], stackptr);
return;
}
if (curplan.plan[entryoffset + 3] == ftbase_fftcodeletplan || curplan.plan[entryoffset + 3] == ftbase_fhtcodeletplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
n = n1 * n2;
if (n == 3) {
offs = curplan.plan[entryoffset + 7];
curplan.precomputed[offs + 0] = cos(2 * M_PI / 3) - 1;
curplan.precomputed[offs + 1] = sin(2 * M_PI / 3);
return;
}
if (n == 5) {
offs = curplan.plan[entryoffset + 7];
v = 2 * M_PI / 5;
curplan.precomputed[offs + 0] = (cos(v) + cos(2 * v)) / 2 - 1;
curplan.precomputed[offs + 1] = (cos(v) - cos(2 * v)) / 2;
curplan.precomputed[offs + 2] = -sin(v);
curplan.precomputed[offs + 3] = -(sin(v) + sin(2 * v));
curplan.precomputed[offs + 4] = sin(v) - sin(2 * v);
return;
}
}
if (curplan.plan[entryoffset + 3] == ftbase_fftbluesteinplan) {
ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset + 5], stackptr);
n = curplan.plan[entryoffset + 1];
m = curplan.plan[entryoffset + 4];
offs = curplan.plan[entryoffset + 7];
for (int i = 0; i <= 2 * m - 1; i++)
curplan.precomputed[offs + i] = 0;
for (int i = 0; i < n; i++) {
bx = cos(M_PI * sqr(i) / n);
by = sin(M_PI * sqr(i) / n);
curplan.precomputed[offs + 2 * i + 0] = bx;
curplan.precomputed[offs + 2 * i + 1] = by;
curplan.precomputed[offs + 2 * m + 2 * i + 0] = bx;
curplan.precomputed[offs + 2 * m + 2 * i + 1] = by;
if (i > 0) {
curplan.precomputed[offs + 2 * (m - i) + 0] = bx;
curplan.precomputed[offs + 2 * (m - i) + 1] = by;
}
}
ftbaseexecuteplanrec(&curplan.precomputed, offs, plan, curplan.plan[entryoffset + 5], stackptr);
return;
}
}
void PIFFT_double::ftbasefactorize(int n, int * n1, int * n2) {
*n1 = *n2 = 0;
int ftbase_ftbasecodeletrecommended = 5;
if ((*n1) * (*n2) != n) {
for (int j = ftbase_ftbasecodeletrecommended; j >= 2; j--) {
if (n % j == 0) {
*n1 = j;
*n2 = n / j;
break;
}
}
}
if ((*n1) * (*n2) != n) {
for (int j = ftbase_ftbasecodeletrecommended + 1; j <= n - 1; j++) {
if (n % j == 0) {
*n1 = j;
*n2 = n / j;
break;
}
}
}
if ((*n1) * (*n2) != n) {
*n1 = 1;
*n2 = n;
}
if ((*n2) == 1 && (*n1) != 1) {
*n2 = *n1;
*n1 = 1;
}
}
/*************************************************************************
Is number smooth?
-- ALGLIB --
Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT_double::ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int * best) {
if (seed >= n) {
*best = piMini(*best, seed);
return;
}
if (leastfactor <= 2)
ftbase_ftbasefindsmoothrec(n, seed * 2, 2, best);
if (leastfactor <= 3)
ftbase_ftbasefindsmoothrec(n, seed * 3, 3, best);
if (leastfactor <= 5)
ftbase_ftbasefindsmoothrec(n, seed * 5, 5, best);
}
int PIFFT_double::ftbasefindsmooth(int n) {
int best, result;
best = 2;
while (best < n)
best = 2 * best;
ftbase_ftbasefindsmoothrec(n, 1, 2, &best);
result = best;
return result;
}
void PIFFT_double::ftbase_internalreallintranspose(PIVector<double> * a, int m, int n, int astart, PIVector<double> * buf) {
ftbase_fftirltrec(a, astart, n, buf, 0, m, m, n);
for (int i = 0; i < 2 * m * n; i++) (*a)[astart + i] = (*buf)[i];
}
void PIFFT_double::ftbase_fftirltrec(PIVector<double> * a, int astart, int astride, PIVector<double> * b, int bstart, int bstride, int m, int n) {
int idx1, idx2;
int m1, n1;
if (m == 0 || n == 0)
return;
if (piMaxi(m, n) <= 8) {
for (int i = 0; i <= m - 1; i++) {
idx1 = bstart + i;
idx2 = astart + i * astride;
for (int j = 0; j <= n - 1; j++) {
(*b)[idx1] = a->at(idx2);
idx1 = idx1 + bstride;
idx2 = idx2 + 1;
}
}
return;
}
if (n > m) {
n1 = n / 2;
if (n - n1 >= 8 && n1 % 8 != 0)
n1 = n1 + (8 - n1 % 8);
ftbase_fftirltrec(a, astart, astride, b, bstart, bstride, m, n1);
ftbase_fftirltrec(a, astart + n1, astride, b, bstart + n1 * bstride, bstride, m, n - n1);
} else {
m1 = m / 2;
if (m - m1 >= 8 && m1 % 8 != 0)
m1 = m1 + (8 - m1 % 8);
ftbase_fftirltrec(a, astart, astride, b, bstart, bstride, m1, n);
ftbase_fftirltrec(a, astart + m1 * astride, astride, b, bstart + m1, bstride, m - m1, n);
}
}
void PIFFT_double::ftbase_internalcomplexlintranspose(PIVector<double> * a, int m, int n, int astart, PIVector<double> * buf) {
ftbase_ffticltrec(a, astart, n, buf, 0, m, m, n);
for (int i = 0; i < 2 * m * n; i++)
(*a)[astart + i] = (*buf)[i];
}
void PIFFT_double::ftbase_ffticltrec(PIVector<double> * a, int astart, int astride, PIVector<double> * b, int bstart, int bstride, int m, int n) {
int idx1, idx2, m2, m1, n1;
if (m == 0 || n == 0)
return;
if (piMax<int>(m, n) <= 8) {
m2 = 2 * bstride;
for (int i = 0; i <= m - 1; i++) {
idx1 = bstart + 2 * i;
idx2 = astart + 2 * i * astride;
for (int j = 0; j <= n - 1; j++) {
(*b)[idx1 + 0] = a->at(idx2 + 0);
(*b)[idx1 + 1] = a->at(idx2 + 1);
idx1 = idx1 + m2;
idx2 = idx2 + 2;
}
}
return;
}
if (n > m) {
n1 = n / 2;
if (n - n1 >= 8 && n1 % 8 != 0)
n1 = n1 + (8 - n1 % 8);
ftbase_ffticltrec(a, astart, astride, b, bstart, bstride, m, n1);
ftbase_ffticltrec(a, astart + 2 * n1, astride, b, bstart + 2 * n1 * bstride, bstride, m, n - n1);
} else {
m1 = m / 2;
if (m - m1 >= 8 && m1 % 8 != 0)
m1 = m1 + (8 - m1 % 8);
ftbase_ffticltrec(a, astart, astride, b, bstart, bstride, m1, n);
ftbase_ffticltrec(a, astart + 2 * m1 * astride, astride, b, bstart + 2 * m1, bstride, m - m1, n);
}
}
void PIFFT_double::ftbaseexecuteplan(PIVector<double> * a, int aoffset, int n, ftplan * plan) {
ae_int_t stackptr;
stackptr = 0;
ftbaseexecuteplanrec(a, aoffset, plan, 0, stackptr);
}
/*************************************************************************
Recurrent subroutine for the FTBaseExecutePlan
Parameters:
A FFT'ed array
AOffset offset of the FFT'ed part (distance is measured in doubles)
-- ALGLIB --
Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT_double::ftbaseexecuteplanrec(PIVector<double> * a, int aoffset, ftplan * plan, int entryoffset, ae_int_t stackptr) {
int n1, n2, n, m, offs, offs1, offs2, offsa, offsb, offsp;
double hk, hnk, x, y, bx, by, v0, v1, v2, v3;
double a0x, a0y, a1x, a1y, a2x, a2y, a3x, a3y;
double t1x, t1y, t2x, t2y, t3x, t3y, t4x, t4y, t5x, t5y;
double m1x, m1y, m2x, m2y, m3x, m3y, m4x, m4y, m5x, m5y;
double s1x, s1y, s2x, s2y, s3x, s3y, s4x, s4y, s5x, s5y;
double c1, c2, c3, c4, c5;
int ftbase_fftcooleytukeyplan = 0;
int ftbase_fftbluesteinplan = 1;
int ftbase_fftcodeletplan = 2;
int ftbase_fhtcooleytukeyplan = 3;
int ftbase_fhtcodeletplan = 4;
int ftbase_fftrealcooleytukeyplan = 5;
int ftbase_fftemptyplan = 6;
PIVector<double> & tmpb(curplan.tmpbuf);
if (curplan.plan[entryoffset + 3] == ftbase_fftemptyplan)
return;
if (curplan.plan[entryoffset + 3] == ftbase_fftcooleytukeyplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n2 - 1; i++)
ftbaseexecuteplanrec(a, aoffset + i * n1 * 2, plan, curplan.plan[entryoffset + 5], stackptr);
ftbase_ffttwcalc(a, aoffset, n1, n2);
ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n1 - 1; i++)
ftbaseexecuteplanrec(a, aoffset + i * n2 * 2, plan, curplan.plan[entryoffset + 6], stackptr);
ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
return;
}
if (curplan.plan[entryoffset + 3] == ftbase_fftrealcooleytukeyplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n1 / 2 - 1; i++) {
offs = aoffset + 2 * i * n2 * 2;
for (int k = 0; k <= n2 - 1; k++)
(*a)[offs + 2 * k + 1] = (*a)[offs + 2 * n2 + 2 * k + 0];
ftbaseexecuteplanrec(a, offs, plan, curplan.plan[entryoffset + 6], stackptr);
tmpb[0] = (*a)[offs + 0];
tmpb[1] = 0;
tmpb[2 * n2 + 0] = (*a)[offs + 1];
tmpb[2 * n2 + 1] = 0;
for (int k = 1; k <= n2 - 1; k++) {
offs1 = 2 * k;
offs2 = 2 * n2 + 2 * k;
hk = (*a)[offs + 2 * k + 0];
hnk = (*a)[offs + 2 * (n2 - k) + 0];
tmpb[offs1 + 0] = 0.5 * (hk + hnk);
tmpb[offs2 + 1] = -0.5 * (hk - hnk);
hk = (*a)[offs + 2 * k + 1];
hnk = (*a)[offs + 2 * (n2 - k) + 1];
tmpb[offs2 + 0] = 0.5 * (hk + hnk);
tmpb[offs1 + 1] = 0.5 * (hk - hnk);
}
for (int k = 0; k < 2 * n2 * 2; k++) (*a)[offs + k] = tmpb[k];
}
if (n1 % 2 != 0)
ftbaseexecuteplanrec(a, aoffset + (n1 - 1)*n2 * 2, plan, curplan.plan[entryoffset + 6], stackptr);
ftbase_ffttwcalc(a, aoffset, n2, n1);
ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n2 - 1; i++)
ftbaseexecuteplanrec(a, aoffset + i * n1 * 2, plan, curplan.plan[entryoffset + 5], stackptr);
ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
return;
}
if (curplan.plan[entryoffset + 3] == ftbase_fhtcooleytukeyplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
n = n1 * n2;
ftbase_internalreallintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n2 - 1; i++)
ftbaseexecuteplanrec(a, aoffset + i * n1, plan, curplan.plan[entryoffset + 5], stackptr);
for (int i = 0; i <= n2 - 1; i++) {
for (int j = 0; j <= n1 - 1; j++) {
offsa = aoffset + i * n1;
hk = (*a)[offsa + j];
hnk = (*a)[offsa + (n1 - j) % n1];
offs = 2 * (i * n1 + j);
tmpb[offs + 0] = -0.5 * (hnk - hk);
tmpb[offs + 1] = 0.5 * (hk + hnk);
}
}
ftbase_ffttwcalc(&(curplan.tmpbuf), 0, n1, n2);
for (int j = 0; j <= n1 - 1; j++)
(*a)[aoffset + j] = tmpb[2 * j + 0] + tmpb[2 * j + 1];
if (n2 % 2 == 0) {
offs = 2 * (n2 / 2) * n1;
offsa = aoffset + n2 / 2 * n1;
for (int j = 0; j <= n1 - 1; j++)
(*a)[offsa + j] = tmpb[offs + 2 * j + 0] + tmpb[offs + 2 * j + 1];
}
for (int i = 1; i <= (n2 + 1) / 2 - 1; i++) {
offs = 2 * i * n1;
offs2 = 2 * (n2 - i) * n1;
offsa = aoffset + i * n1;
for (int j = 0; j <= n1 - 1; j++)
(*a)[offsa + j] = tmpb[offs + 2 * j + 1] + tmpb[offs2 + 2 * j + 0];
offsa = aoffset + (n2 - i) * n1;
for (int j = 0; j <= n1 - 1; j++)
(*a)[offsa + j] = tmpb[offs + 2 * j + 0] + tmpb[offs2 + 2 * j + 1];
}
ftbase_internalreallintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n1 - 1; i++)
ftbaseexecuteplanrec(a, aoffset + i * n2, plan, curplan.plan[entryoffset + 6], stackptr);
ftbase_internalreallintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
return;
}
if (curplan.plan[entryoffset + 3] == ftbase_fftcodeletplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
n = n1 * n2;
if (n == 2) {
a0x = (*a)[aoffset + 0];
a0y = (*a)[aoffset + 1];
a1x = (*a)[aoffset + 2];
a1y = (*a)[aoffset + 3];
v0 = a0x + a1x;
v1 = a0y + a1y;
v2 = a0x - a1x;
v3 = a0y - a1y;
(*a)[aoffset + 0] = v0;
(*a)[aoffset + 1] = v1;
(*a)[aoffset + 2] = v2;
(*a)[aoffset + 3] = v3;
return;
}
if (n == 3) {
offs = curplan.plan[entryoffset + 7];
c1 = curplan.precomputed[offs + 0];
c2 = curplan.precomputed[offs + 1];
a0x = (*a)[aoffset + 0];
a0y = (*a)[aoffset + 1];
a1x = (*a)[aoffset + 2];
a1y = (*a)[aoffset + 3];
a2x = (*a)[aoffset + 4];
a2y = (*a)[aoffset + 5];
t1x = a1x + a2x;
t1y = a1y + a2y;
a0x = a0x + t1x;
a0y = a0y + t1y;
m1x = c1 * t1x;
m1y = c1 * t1y;
m2x = c2 * (a1y - a2y);
m2y = c2 * (a2x - a1x);
s1x = a0x + m1x;
s1y = a0y + m1y;
a1x = s1x + m2x;
a1y = s1y + m2y;
a2x = s1x - m2x;
a2y = s1y - m2y;
(*a)[aoffset + 0] = a0x;
(*a)[aoffset + 1] = a0y;
(*a)[aoffset + 2] = a1x;
(*a)[aoffset + 3] = a1y;
(*a)[aoffset + 4] = a2x;
(*a)[aoffset + 5] = a2y;
return;
}
if (n == 4) {
a0x = (*a)[aoffset + 0];
a0y = (*a)[aoffset + 1];
a1x = (*a)[aoffset + 2];
a1y = (*a)[aoffset + 3];
a2x = (*a)[aoffset + 4];
a2y = (*a)[aoffset + 5];
a3x = (*a)[aoffset + 6];
a3y = (*a)[aoffset + 7];
t1x = a0x + a2x;
t1y = a0y + a2y;
t2x = a1x + a3x;
t2y = a1y + a3y;
m2x = a0x - a2x;
m2y = a0y - a2y;
m3x = a1y - a3y;
m3y = a3x - a1x;
(*a)[aoffset + 0] = t1x + t2x;
(*a)[aoffset + 1] = t1y + t2y;
(*a)[aoffset + 4] = t1x - t2x;
(*a)[aoffset + 5] = t1y - t2y;
(*a)[aoffset + 2] = m2x + m3x;
(*a)[aoffset + 3] = m2y + m3y;
(*a)[aoffset + 6] = m2x - m3x;
(*a)[aoffset + 7] = m2y - m3y;
return;
}
if (n == 5) {
offs = curplan.plan[entryoffset + 7];
c1 = curplan.precomputed[offs + 0];
c2 = curplan.precomputed[offs + 1];
c3 = curplan.precomputed[offs + 2];
c4 = curplan.precomputed[offs + 3];
c5 = curplan.precomputed[offs + 4];
t1x = (*a)[aoffset + 2] + (*a)[aoffset + 8];
t1y = (*a)[aoffset + 3] + (*a)[aoffset + 9];
t2x = (*a)[aoffset + 4] + (*a)[aoffset + 6];
t2y = (*a)[aoffset + 5] + (*a)[aoffset + 7];
t3x = (*a)[aoffset + 2] - (*a)[aoffset + 8];
t3y = (*a)[aoffset + 3] - (*a)[aoffset + 9];
t4x = (*a)[aoffset + 6] - (*a)[aoffset + 4];
t4y = (*a)[aoffset + 7] - (*a)[aoffset + 5];
t5x = t1x + t2x;
t5y = t1y + t2y;
(*a)[aoffset + 0] = (*a)[aoffset + 0] + t5x;
(*a)[aoffset + 1] = (*a)[aoffset + 1] + t5y;
m1x = c1 * t5x;
m1y = c1 * t5y;
m2x = c2 * (t1x - t2x);
m2y = c2 * (t1y - t2y);
m3x = -c3 * (t3y + t4y);
m3y = c3 * (t3x + t4x);
m4x = -c4 * t4y;
m4y = c4 * t4x;
m5x = -c5 * t3y;
m5y = c5 * t3x;
s3x = m3x - m4x;
s3y = m3y - m4y;
s5x = m3x + m5x;
s5y = m3y + m5y;
s1x = (*a)[aoffset + 0] + m1x;
s1y = (*a)[aoffset + 1] + m1y;
s2x = s1x + m2x;
s2y = s1y + m2y;
s4x = s1x - m2x;
s4y = s1y - m2y;
(*a)[aoffset + 2] = s2x + s3x;
(*a)[aoffset + 3] = s2y + s3y;
(*a)[aoffset + 4] = s4x + s5x;
(*a)[aoffset + 5] = s4y + s5y;
(*a)[aoffset + 6] = s4x - s5x;
(*a)[aoffset + 7] = s4y - s5y;
(*a)[aoffset + 8] = s2x - s3x;
(*a)[aoffset + 9] = s2y - s3y;
return;
}
}
if (curplan.plan[entryoffset + 3] == ftbase_fhtcodeletplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
n = n1 * n2;
if (n == 2) {
a0x = (*a)[aoffset + 0];
a1x = (*a)[aoffset + 1];
(*a)[aoffset + 0] = a0x + a1x;
(*a)[aoffset + 1] = a0x - a1x;
return;
}
if (n == 3) {
offs = curplan.plan[entryoffset + 7];
c1 = curplan.precomputed[offs + 0];
c2 = curplan.precomputed[offs + 1];
a0x = (*a)[aoffset + 0];
a1x = (*a)[aoffset + 1];
a2x = (*a)[aoffset + 2];
t1x = a1x + a2x;
a0x = a0x + t1x;
m1x = c1 * t1x;
m2y = c2 * (a2x - a1x);
s1x = a0x + m1x;
(*a)[aoffset + 0] = a0x;
(*a)[aoffset + 1] = s1x - m2y;
(*a)[aoffset + 2] = s1x + m2y;
return;
}
if (n == 4) {
a0x = (*a)[aoffset + 0];
a1x = (*a)[aoffset + 1];
a2x = (*a)[aoffset + 2];
a3x = (*a)[aoffset + 3];
t1x = a0x + a2x;
t2x = a1x + a3x;
m2x = a0x - a2x;
m3y = a3x - a1x;
(*a)[aoffset + 0] = t1x + t2x;
(*a)[aoffset + 1] = m2x - m3y;
(*a)[aoffset + 2] = t1x - t2x;
(*a)[aoffset + 3] = m2x + m3y;
return;
}
if (n == 5) {
offs = curplan.plan[entryoffset + 7];
c1 = curplan.precomputed[offs + 0];
c2 = curplan.precomputed[offs + 1];
c3 = curplan.precomputed[offs + 2];
c4 = curplan.precomputed[offs + 3];
c5 = curplan.precomputed[offs + 4];
t1x = (*a)[aoffset + 1] + (*a)[aoffset + 4];
t2x = (*a)[aoffset + 2] + (*a)[aoffset + 3];
t3x = (*a)[aoffset + 1] - (*a)[aoffset + 4];
t4x = (*a)[aoffset + 3] - (*a)[aoffset + 2];
t5x = t1x + t2x;
v0 = (*a)[aoffset + 0] + t5x;
(*a)[aoffset + 0] = v0;
m2x = c2 * (t1x - t2x);
m3y = c3 * (t3x + t4x);
s3y = m3y - c4 * t4x;
s5y = m3y + c5 * t3x;
s1x = v0 + c1 * t5x;
s2x = s1x + m2x;
s4x = s1x - m2x;
(*a)[aoffset + 1] = s2x - s3y;
(*a)[aoffset + 2] = s4x - s5y;
(*a)[aoffset + 3] = s4x + s5y;
(*a)[aoffset + 4] = s2x + s3y;
return;
}
}
if (curplan.plan[entryoffset + 3] == ftbase_fftbluesteinplan) {
n = curplan.plan[entryoffset + 1];
m = curplan.plan[entryoffset + 4];
offs = curplan.plan[entryoffset + 7];
for (int i = stackptr + 2 * n; i <= stackptr + 2 * m - 1; i++)
curplan.stackbuf[i] = 0;
offsp = offs + 2 * m;
offsa = aoffset;
offsb = stackptr;
for (int i = 0; i < n; i++) {
bx = curplan.precomputed[offsp + 0];
by = curplan.precomputed[offsp + 1];
x = (*a)[offsa + 0];
y = (*a)[offsa + 1];
curplan.stackbuf[offsb + 0] = x * bx - y * (-by);
curplan.stackbuf[offsb + 1] = x * (-by) + y * bx;
offsp = offsp + 2;
offsa = offsa + 2;
offsb = offsb + 2;
}
ftbaseexecuteplanrec(&curplan.stackbuf, stackptr, plan, curplan.plan[entryoffset + 5], stackptr + 2 * 2 * m);
offsb = stackptr;
offsp = offs;
for (int i = 0; i <= m - 1; i++) {
x = curplan.stackbuf[offsb + 0];
y = curplan.stackbuf[offsb + 1];
bx = curplan.precomputed[offsp + 0];
by = curplan.precomputed[offsp + 1];
curplan.stackbuf[offsb + 0] = x * bx - y * by;
curplan.stackbuf[offsb + 1] = -(x * by + y * bx);
offsb = offsb + 2;
offsp = offsp + 2;
}
ftbaseexecuteplanrec(&curplan.stackbuf, stackptr, plan, curplan.plan[entryoffset + 5], stackptr + 2 * 2 * m);
offsb = stackptr;
offsp = offs + 2 * m;
offsa = aoffset;
for (int i = 0; i < n; i++) {
x = curplan.stackbuf[offsb + 0] / m;
y = -curplan.stackbuf[offsb + 1] / m;
bx = curplan.precomputed[offsp + 0];
by = curplan.precomputed[offsp + 1];
(*a)[offsa + 0] = x * bx - y * (-by);
(*a)[offsa + 1] = x * (-by) + y * bx;
offsp = offsp + 2;
offsa = offsa + 2;
offsb = offsb + 2;
}
return;
}
}
/*************************************************************************
Twiddle factors calculation
-- ALGLIB --
Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT_double::ftbase_ffttwcalc(PIVector<double> * a, int aoffset, int n1, int n2) {
int n, idx, offs;
double x, y, twxm1, twy, twbasexm1, twbasey, twrowxm1, twrowy, tmpx, tmpy, v;
int ftbase_ftbaseupdatetw = 4;
n = n1 * n2;
v = -2 * M_PI / n;
twbasexm1 = -2 * sqr(sin(0.5 * v));
twbasey = sin(v);
twrowxm1 = 0;
twrowy = 0;
for (int i = 0, j = 0; i <= n2 - 1; i++) {
twxm1 = 0;
twy = 0;
for (j = 0; j <= n1 - 1; j++) {
idx = i * n1 + j;
offs = aoffset + 2 * idx;
x = (*a)[offs + 0];
y = (*a)[offs + 1];
tmpx = x * twxm1 - y * twy;
tmpy = x * twy + y * twxm1;
(*a)[offs + 0] = x + tmpx;
(*a)[offs + 1] = y + tmpy;
if (j < n1 - 1) {
if (j % ftbase_ftbaseupdatetw == 0) {
v = -2 * M_PI * i * (j + 1) / n;
twxm1 = -2 * sqr(sin(0.5 * v));
twy = sin(v);
} else {
tmpx = twrowxm1 + twxm1 * twrowxm1 - twy * twrowy;
tmpy = twrowy + twxm1 * twrowy + twy * twrowxm1;
twxm1 = twxm1 + tmpx;
twy = twy + tmpy;
}
}
}
if (i < n2 - 1) {
if (j % ftbase_ftbaseupdatetw == 0) {
v = -2 * M_PI * (i + 1) / n;
twrowxm1 = -2 * sqr(sin(0.5 * v));
twrowy = sin(v);
} else {
tmpx = twbasexm1 + twrowxm1 * twbasexm1 - twrowy * twbasey;
tmpy = twbasey + twrowxm1 * twbasey + twrowy * twbasexm1;
twrowxm1 = twrowxm1 + tmpx;
twrowy = twrowy + tmpy;
}
}
}
}
//=================================================================================================
//*************************************************************************************************
//---------------------FLOAT-----------------------FLOAT-------------------------------------------
//*************************************************************************************************
//=================================================================================================
PIFFT_float::PIFFT_float() {
}
PIVector<complexf> * PIFFT_float::calcFFT(const PIVector<complexf> & val) {
result.clear();
if (val.size_s() < 4) return &result;
fftc1d(val, val.size());
return &result;
}
PIVector<complexf> * PIFFT_float::calcFFTinverse(const PIVector<complexf> & val) {
result.clear();
if (val.size_s() < 4) return &result;
fftc1dinv(val, val.size());
return &result;
}
PIVector<complexf> * PIFFT_float::calcHilbert(const PIVector<float> & val) {
result.clear();
if (val.size_s() < 4) return &result;
fftc1r(val, val.size());
for (uint i = 0; i < result.size() / 2; i++) result[i] = result[i] * 2.f;
for (uint i = result.size() / 2; i < result.size(); i++) result[i] = 0;
fftc1dinv(result, result.size());
return &result;
}
PIVector<complexf> * PIFFT_float::calcFFT(const PIVector<float> & val) {
result.clear();
if (val.size_s() < 4) return &result;
fftc1r(val, val.size());
return &result;
}
PIVector<float> PIFFT_float::getAmplitude() const {
PIVector<float> ret;
ret.resize(result.size());
float tmp;
for (uint i = 0; i < result.size(); i++) {
tmp = sqrt(result[i].real() * result[i].real() + result[i].imag() * result[i].imag());
ret[i] = tmp;
}
return ret;
}
PIVector<float> PIFFT_float::getReal() const {
PIVector<float> ret;
ret.resize(result.size());
for (uint i = 0; i < result.size(); i++)
ret[i] = result[i].real();
return ret;
}
PIVector<float> PIFFT_float::getImag() const {
PIVector<float> ret;
ret.resize(result.size());
for (uint i = 0; i < result.size(); i++)
ret[i] = result[i].imag();
return ret;
}
void PIFFT_float::fftc1d(const PIVector<complexf> & a, uint n) {
createPlan(n);
uint i;
PIVector<float> buf;
buf.resize(2 * n);
for (i = 0; i < n; i++) {
buf[2 * i + 0] = a[i].real();
buf[2 * i + 1] = a[i].imag();
}
ftbaseexecuteplan(&buf, 0, n, &curplan);
result.resize(n);
for (i = 0; i < n; i++)
result[i] = complexf(buf[2 * i + 0], buf[2 * i + 1]);
}
void PIFFT_float::fftc1r(const PIVector<float> & a, uint n) {
uint i;
if (n % 2 == 0) {
PIVector<float> buf;
uint n2 = n / 2;
buf = a;
createPlan(n2);
ftbaseexecuteplan(&buf, 0, n2, &curplan);
result.resize(n);
uint idx;
complexf hn, hmnc, v;
for (i = 0; i <= n2; i++) {
idx = 2 * (i % n2);
hn = complexf(buf[idx + 0], buf[idx + 1]);
idx = 2 * ((n2 - i) % n2);
hmnc = complexf(buf[idx + 0], -buf[idx + 1]);
v = complexf(sin(M_PI * i / n2), cos(M_PI * i / n2));
result[i] = ((hn + hmnc) - (v * (hn - hmnc)));
result[i] *= 0.5;
}
for (i = n2 + 1; i < n; i++)
result[i] = conj(result[n - i]);
} else {
PIVector<complexf> cbuf;
cbuf.resize(n);
for (i = 0; i < n; i++)
cbuf[i] = complexf(a[i], 0.);
fftc1d(cbuf, n);
}
}
void PIFFT_float::fftc1dinv(const PIVector<complexf> & a, uint n) {
PIVector<complexf> cbuf;
cbuf.resize(n);
uint i;
for (i = 0; i < n; i++) {
cbuf[i] = conj(a[i]);
}
fftc1d(cbuf, n);
for (i = 0; i < n; i++) {
result[i] = conj(result[i] / (float)n);
}
}
void PIFFT_float::createPlan(uint n) {
curplan.plan.clear();
curplan.precomputed.clear();
curplan.stackbuf.clear();
curplan.tmpbuf.clear();
if (n < 2) return;
ftbasegeneratecomplexfftplan(n, &curplan);
}
void PIFFT_float::ftbasegeneratecomplexfftplan(uint n, ftplan * plan) {
int planarraysize;
int plansize;
int precomputedsize;
int tmpmemsize;
int stackmemsize;
ae_int_t stackptr;
planarraysize = 1;
plansize = 0;
precomputedsize = 0;
stackmemsize = 0;
stackptr = 0;
tmpmemsize = 2 * n;
curplan.plan.resize(planarraysize);
int ftbase_ftbasecffttask = 0;
ftbase_ftbasegenerateplanrec(n, ftbase_ftbasecffttask, plan, &plansize, &precomputedsize, &planarraysize, &tmpmemsize, &stackmemsize, stackptr);
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_float::ftbase_ftbasegenerateplanrec(
int n,
int tasktype,
ftplan * plan,
int * plansize,
int * precomputedsize,
int * planarraysize,
int * tmpmemsize,
int * stackmemsize,
ae_int_t stackptr, int debugi) {
int k, m, n1, n2, esize, entryoffset;
int ftbase_ftbaseplanentrysize = 4;
int ftbase_ftbasecffttask = 0;
int ftbase_fftcooleytukeyplan = 0;
int ftbase_fftbluesteinplan = 1;
int ftbase_fftcodeletplan = 2;
int ftbase_fftrealcooleytukeyplan = 5;
int ftbase_fftemptyplan = 6;
if (*plansize + ftbase_ftbaseplanentrysize > (*planarraysize)) {
curplan.plan.resize(4 * (*planarraysize));
*planarraysize = 4 * (*planarraysize);
}
entryoffset = *plansize;
esize = ftbase_ftbaseplanentrysize;
*plansize = *plansize + esize;
if (n == 1) {
curplan.plan[entryoffset + 0] = esize;
curplan.plan[entryoffset + 1] = -1;
curplan.plan[entryoffset + 2] = -1;
curplan.plan[entryoffset + 3] = ftbase_fftemptyplan;
curplan.plan[entryoffset + 4] = -1;
curplan.plan[entryoffset + 5] = -1;
curplan.plan[entryoffset + 6] = -1;
curplan.plan[entryoffset + 7] = -1;
return;
}
ftbasefactorize(n, &n1, &n2);
if (n1 != 1) {
*tmpmemsize = piMax<int>(*tmpmemsize, 2 * n1 * n2);
curplan.plan[entryoffset + 0] = esize;
curplan.plan[entryoffset + 1] = n1;
curplan.plan[entryoffset + 2] = n2;
if (tasktype == ftbase_ftbasecffttask)
curplan.plan[entryoffset + 3] = ftbase_fftcooleytukeyplan;
else
curplan.plan[entryoffset + 3] = ftbase_fftrealcooleytukeyplan;
curplan.plan[entryoffset + 4] = 0;
curplan.plan[entryoffset + 5] = *plansize;
debugi++;
ftbase_ftbasegenerateplanrec(n1, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr, debugi);
curplan.plan[entryoffset + 6] = *plansize;
ftbase_ftbasegenerateplanrec(n2, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr, debugi);
curplan.plan[entryoffset + 7] = -1;
return;
} else {
if (n >= 2 && n <= 5) {
curplan.plan[entryoffset + 0] = esize;
curplan.plan[entryoffset + 1] = n1;
curplan.plan[entryoffset + 2] = n2;
curplan.plan[entryoffset + 3] = ftbase_fftcodeletplan;
curplan.plan[entryoffset + 4] = 0;
curplan.plan[entryoffset + 5] = -1;
curplan.plan[entryoffset + 6] = -1;
curplan.plan[entryoffset + 7] = *precomputedsize;
if (n == 3)
*precomputedsize = *precomputedsize + 2;
if (n == 5)
*precomputedsize = *precomputedsize + 5;
return;
} else {
k = 2 * n2 - 1;
m = ftbasefindsmooth(k);
*tmpmemsize = piMax<int>(*tmpmemsize, 2 * m);
curplan.plan[entryoffset + 0] = esize;
curplan.plan[entryoffset + 1] = n2;
curplan.plan[entryoffset + 2] = -1;
curplan.plan[entryoffset + 3] = ftbase_fftbluesteinplan;
curplan.plan[entryoffset + 4] = m;
curplan.plan[entryoffset + 5] = *plansize;
stackptr = stackptr + 2 * 2 * m;
*stackmemsize = piMax<int>(*stackmemsize, stackptr);
ftbase_ftbasegenerateplanrec(m, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr);
stackptr = stackptr - 2 * 2 * m;
curplan.plan[entryoffset + 6] = -1;
curplan.plan[entryoffset + 7] = *precomputedsize;
*precomputedsize = *precomputedsize + 2 * m + 2 * n;
return;
}
}
}
/*************************************************************************
Recurrent subroutine for precomputing FFT plans
-- ALGLIB --
Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT_float::ftbase_ftbaseprecomputeplanrec(ftplan * plan,
int entryoffset,
ae_int_t stackptr) {
int n1, n2, n, m, offs;
float v, bx, by;
int ftbase_fftcooleytukeyplan = 0;
int ftbase_fftbluesteinplan = 1;
int ftbase_fftcodeletplan = 2;
int ftbase_fhtcooleytukeyplan = 3;
int ftbase_fhtcodeletplan = 4;
int ftbase_fftrealcooleytukeyplan = 5;
if ((curplan.plan[entryoffset + 3] == ftbase_fftcooleytukeyplan || curplan.plan[entryoffset + 3] == ftbase_fftrealcooleytukeyplan) || curplan.plan[entryoffset + 3] == ftbase_fhtcooleytukeyplan) {
ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset + 5], stackptr);
ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset + 6], stackptr);
return;
}
if (curplan.plan[entryoffset + 3] == ftbase_fftcodeletplan || curplan.plan[entryoffset + 3] == ftbase_fhtcodeletplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
n = n1 * n2;
if (n == 3) {
offs = curplan.plan[entryoffset + 7];
curplan.precomputed[offs + 0] = cos(2 * M_PI / 3) - 1;
curplan.precomputed[offs + 1] = sin(2 * M_PI / 3);
return;
}
if (n == 5) {
offs = curplan.plan[entryoffset + 7];
v = 2 * M_PI / 5;
curplan.precomputed[offs + 0] = (cos(v) + cos(2 * v)) / 2 - 1;
curplan.precomputed[offs + 1] = (cos(v) - cos(2 * v)) / 2;
curplan.precomputed[offs + 2] = -sin(v);
curplan.precomputed[offs + 3] = -(sin(v) + sin(2 * v));
curplan.precomputed[offs + 4] = sin(v) - sin(2 * v);
return;
}
}
if (curplan.plan[entryoffset + 3] == ftbase_fftbluesteinplan) {
ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset + 5], stackptr);
n = curplan.plan[entryoffset + 1];
m = curplan.plan[entryoffset + 4];
offs = curplan.plan[entryoffset + 7];
for (int i = 0; i <= 2 * m - 1; i++)
curplan.precomputed[offs + i] = 0;
for (int i = 0; i < n; i++) {
bx = cos(M_PI * sqr(i) / n);
by = sin(M_PI * sqr(i) / n);
curplan.precomputed[offs + 2 * i + 0] = bx;
curplan.precomputed[offs + 2 * i + 1] = by;
curplan.precomputed[offs + 2 * m + 2 * i + 0] = bx;
curplan.precomputed[offs + 2 * m + 2 * i + 1] = by;
if (i > 0) {
curplan.precomputed[offs + 2 * (m - i) + 0] = bx;
curplan.precomputed[offs + 2 * (m - i) + 1] = by;
}
}
ftbaseexecuteplanrec(&curplan.precomputed, offs, plan, curplan.plan[entryoffset + 5], stackptr);
return;
}
}
void PIFFT_float::ftbasefactorize(int n, int * n1, int * n2) {
*n1 = *n2 = 0;
int ftbase_ftbasecodeletrecommended = 5;
if ((*n1) * (*n2) != n) {
for (int j = ftbase_ftbasecodeletrecommended; j >= 2; j--) {
if (n % j == 0) {
*n1 = j;
*n2 = n / j;
break;
}
}
}
if ((*n1) * (*n2) != n) {
for (int j = ftbase_ftbasecodeletrecommended + 1; j <= n - 1; j++) {
if (n % j == 0) {
*n1 = j;
*n2 = n / j;
break;
}
}
}
if ((*n1) * (*n2) != n) {
*n1 = 1;
*n2 = n;
}
if ((*n2) == 1 && (*n1) != 1) {
*n2 = *n1;
*n1 = 1;
}
}
/*************************************************************************
Is number smooth?
-- ALGLIB --
Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT_float::ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int * best) {
if (seed >= n) {
*best = piMini(*best, seed);
return;
}
if (leastfactor <= 2)
ftbase_ftbasefindsmoothrec(n, seed * 2, 2, best);
if (leastfactor <= 3)
ftbase_ftbasefindsmoothrec(n, seed * 3, 3, best);
if (leastfactor <= 5)
ftbase_ftbasefindsmoothrec(n, seed * 5, 5, best);
}
int PIFFT_float::ftbasefindsmooth(int n) {
int best, result;
best = 2;
while (best < n)
best = 2 * best;
ftbase_ftbasefindsmoothrec(n, 1, 2, &best);
result = best;
return result;
}
void PIFFT_float::ftbase_internalreallintranspose(PIVector<float> * a, int m, int n, int astart, PIVector<float> * buf) {
ftbase_fftirltrec(a, astart, n, buf, 0, m, m, n);
for (int i = 0; i < 2 * m * n; i++) (*a)[astart + i] = (*buf)[i];
}
void PIFFT_float::ftbase_fftirltrec(PIVector<float> * a, int astart, int astride, PIVector<float> * b, int bstart, int bstride, int m, int n) {
int idx1, idx2;
int m1, n1;
if (m == 0 || n == 0)
return;
if (piMaxi(m, n) <= 8) {
for (int i = 0; i <= m - 1; i++) {
idx1 = bstart + i;
idx2 = astart + i * astride;
for (int j = 0; j <= n - 1; j++) {
(*b)[idx1] = a->at(idx2);
idx1 = idx1 + bstride;
idx2 = idx2 + 1;
}
}
return;
}
if (n > m) {
n1 = n / 2;
if (n - n1 >= 8 && n1 % 8 != 0)
n1 = n1 + (8 - n1 % 8);
ftbase_fftirltrec(a, astart, astride, b, bstart, bstride, m, n1);
ftbase_fftirltrec(a, astart + n1, astride, b, bstart + n1 * bstride, bstride, m, n - n1);
} else {
m1 = m / 2;
if (m - m1 >= 8 && m1 % 8 != 0)
m1 = m1 + (8 - m1 % 8);
ftbase_fftirltrec(a, astart, astride, b, bstart, bstride, m1, n);
ftbase_fftirltrec(a, astart + m1 * astride, astride, b, bstart + m1, bstride, m - m1, n);
}
}
void PIFFT_float::ftbase_internalcomplexlintranspose(PIVector<float> * a, int m, int n, int astart, PIVector<float> * buf) {
ftbase_ffticltrec(a, astart, n, buf, 0, m, m, n);
for (int i = 0; i < 2 * m * n; i++)
(*a)[astart + i] = (*buf)[i];
}
void PIFFT_float::ftbase_ffticltrec(PIVector<float> * a, int astart, int astride, PIVector<float> * b, int bstart, int bstride, int m, int n) {
int idx1, idx2, m2, m1, n1;
if (m == 0 || n == 0)
return;
if (piMax<int>(m, n) <= 8) {
m2 = 2 * bstride;
for (int i = 0; i <= m - 1; i++) {
idx1 = bstart + 2 * i;
idx2 = astart + 2 * i * astride;
for (int j = 0; j <= n - 1; j++) {
(*b)[idx1 + 0] = a->at(idx2 + 0);
(*b)[idx1 + 1] = a->at(idx2 + 1);
idx1 = idx1 + m2;
idx2 = idx2 + 2;
}
}
return;
}
if (n > m) {
n1 = n / 2;
if (n - n1 >= 8 && n1 % 8 != 0)
n1 = n1 + (8 - n1 % 8);
ftbase_ffticltrec(a, astart, astride, b, bstart, bstride, m, n1);
ftbase_ffticltrec(a, astart + 2 * n1, astride, b, bstart + 2 * n1 * bstride, bstride, m, n - n1);
} else {
m1 = m / 2;
if (m - m1 >= 8 && m1 % 8 != 0)
m1 = m1 + (8 - m1 % 8);
ftbase_ffticltrec(a, astart, astride, b, bstart, bstride, m1, n);
ftbase_ffticltrec(a, astart + 2 * m1 * astride, astride, b, bstart + 2 * m1, bstride, m - m1, n);
}
}
void PIFFT_float::ftbaseexecuteplan(PIVector<float> * a, int aoffset, int n, ftplan * plan) {
ae_int_t stackptr;
stackptr = 0;
ftbaseexecuteplanrec(a, aoffset, plan, 0, stackptr);
}
/*************************************************************************
Recurrent subroutine for the FTBaseExecutePlan
Parameters:
A FFT'ed array
AOffset offset of the FFT'ed part (distance is measured in floats)
-- ALGLIB --
Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT_float::ftbaseexecuteplanrec(PIVector<float> * a, int aoffset, ftplan * plan, int entryoffset, ae_int_t stackptr) {
int n1, n2, n, m, offs, offs1, offs2, offsa, offsb, offsp;
float hk, hnk, x, y, bx, by, v0, v1, v2, v3;
float a0x, a0y, a1x, a1y, a2x, a2y, a3x, a3y;
float t1x, t1y, t2x, t2y, t3x, t3y, t4x, t4y, t5x, t5y;
float m1x, m1y, m2x, m2y, m3x, m3y, m4x, m4y, m5x, m5y;
float s1x, s1y, s2x, s2y, s3x, s3y, s4x, s4y, s5x, s5y;
float c1, c2, c3, c4, c5;
int ftbase_fftcooleytukeyplan = 0;
int ftbase_fftbluesteinplan = 1;
int ftbase_fftcodeletplan = 2;
int ftbase_fhtcooleytukeyplan = 3;
int ftbase_fhtcodeletplan = 4;
int ftbase_fftrealcooleytukeyplan = 5;
int ftbase_fftemptyplan = 6;
PIVector<float> & tmpb(curplan.tmpbuf);
if (curplan.plan[entryoffset + 3] == ftbase_fftemptyplan)
return;
if (curplan.plan[entryoffset + 3] == ftbase_fftcooleytukeyplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n2 - 1; i++)
ftbaseexecuteplanrec(a, aoffset + i * n1 * 2, plan, curplan.plan[entryoffset + 5], stackptr);
ftbase_ffttwcalc(a, aoffset, n1, n2);
ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n1 - 1; i++)
ftbaseexecuteplanrec(a, aoffset + i * n2 * 2, plan, curplan.plan[entryoffset + 6], stackptr);
ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
return;
}
if (curplan.plan[entryoffset + 3] == ftbase_fftrealcooleytukeyplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n1 / 2 - 1; i++) {
offs = aoffset + 2 * i * n2 * 2;
for (int k = 0; k <= n2 - 1; k++)
(*a)[offs + 2 * k + 1] = (*a)[offs + 2 * n2 + 2 * k + 0];
ftbaseexecuteplanrec(a, offs, plan, curplan.plan[entryoffset + 6], stackptr);
tmpb[0] = (*a)[offs + 0];
tmpb[1] = 0;
tmpb[2 * n2 + 0] = (*a)[offs + 1];
tmpb[2 * n2 + 1] = 0;
for (int k = 1; k <= n2 - 1; k++) {
offs1 = 2 * k;
offs2 = 2 * n2 + 2 * k;
hk = (*a)[offs + 2 * k + 0];
hnk = (*a)[offs + 2 * (n2 - k) + 0];
tmpb[offs1 + 0] = 0.5 * (hk + hnk);
tmpb[offs2 + 1] = -0.5 * (hk - hnk);
hk = (*a)[offs + 2 * k + 1];
hnk = (*a)[offs + 2 * (n2 - k) + 1];
tmpb[offs2 + 0] = 0.5 * (hk + hnk);
tmpb[offs1 + 1] = 0.5 * (hk - hnk);
}
for (int k = 0; k < 2 * n2 * 2; k++) (*a)[offs + k] = tmpb[k];
}
if (n1 % 2 != 0)
ftbaseexecuteplanrec(a, aoffset + (n1 - 1)*n2 * 2, plan, curplan.plan[entryoffset + 6], stackptr);
ftbase_ffttwcalc(a, aoffset, n2, n1);
ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n2 - 1; i++)
ftbaseexecuteplanrec(a, aoffset + i * n1 * 2, plan, curplan.plan[entryoffset + 5], stackptr);
ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
return;
}
if (curplan.plan[entryoffset + 3] == ftbase_fhtcooleytukeyplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
n = n1 * n2;
ftbase_internalreallintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n2 - 1; i++)
ftbaseexecuteplanrec(a, aoffset + i * n1, plan, curplan.plan[entryoffset + 5], stackptr);
for (int i = 0; i <= n2 - 1; i++) {
for (int j = 0; j <= n1 - 1; j++) {
offsa = aoffset + i * n1;
hk = (*a)[offsa + j];
hnk = (*a)[offsa + (n1 - j) % n1];
offs = 2 * (i * n1 + j);
tmpb[offs + 0] = -0.5 * (hnk - hk);
tmpb[offs + 1] = 0.5 * (hk + hnk);
}
}
ftbase_ffttwcalc(&(curplan.tmpbuf), 0, n1, n2);
for (int j = 0; j <= n1 - 1; j++)
(*a)[aoffset + j] = tmpb[2 * j + 0] + tmpb[2 * j + 1];
if (n2 % 2 == 0) {
offs = 2 * (n2 / 2) * n1;
offsa = aoffset + n2 / 2 * n1;
for (int j = 0; j <= n1 - 1; j++)
(*a)[offsa + j] = tmpb[offs + 2 * j + 0] + tmpb[offs + 2 * j + 1];
}
for (int i = 1; i <= (n2 + 1) / 2 - 1; i++) {
offs = 2 * i * n1;
offs2 = 2 * (n2 - i) * n1;
offsa = aoffset + i * n1;
for (int j = 0; j <= n1 - 1; j++)
(*a)[offsa + j] = tmpb[offs + 2 * j + 1] + tmpb[offs2 + 2 * j + 0];
offsa = aoffset + (n2 - i) * n1;
for (int j = 0; j <= n1 - 1; j++)
(*a)[offsa + j] = tmpb[offs + 2 * j + 0] + tmpb[offs2 + 2 * j + 1];
}
ftbase_internalreallintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
for (int i = 0; i <= n1 - 1; i++)
ftbaseexecuteplanrec(a, aoffset + i * n2, plan, curplan.plan[entryoffset + 6], stackptr);
ftbase_internalreallintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
return;
}
if (curplan.plan[entryoffset + 3] == ftbase_fftcodeletplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
n = n1 * n2;
if (n == 2) {
a0x = (*a)[aoffset + 0];
a0y = (*a)[aoffset + 1];
a1x = (*a)[aoffset + 2];
a1y = (*a)[aoffset + 3];
v0 = a0x + a1x;
v1 = a0y + a1y;
v2 = a0x - a1x;
v3 = a0y - a1y;
(*a)[aoffset + 0] = v0;
(*a)[aoffset + 1] = v1;
(*a)[aoffset + 2] = v2;
(*a)[aoffset + 3] = v3;
return;
}
if (n == 3) {
offs = curplan.plan[entryoffset + 7];
c1 = curplan.precomputed[offs + 0];
c2 = curplan.precomputed[offs + 1];
a0x = (*a)[aoffset + 0];
a0y = (*a)[aoffset + 1];
a1x = (*a)[aoffset + 2];
a1y = (*a)[aoffset + 3];
a2x = (*a)[aoffset + 4];
a2y = (*a)[aoffset + 5];
t1x = a1x + a2x;
t1y = a1y + a2y;
a0x = a0x + t1x;
a0y = a0y + t1y;
m1x = c1 * t1x;
m1y = c1 * t1y;
m2x = c2 * (a1y - a2y);
m2y = c2 * (a2x - a1x);
s1x = a0x + m1x;
s1y = a0y + m1y;
a1x = s1x + m2x;
a1y = s1y + m2y;
a2x = s1x - m2x;
a2y = s1y - m2y;
(*a)[aoffset + 0] = a0x;
(*a)[aoffset + 1] = a0y;
(*a)[aoffset + 2] = a1x;
(*a)[aoffset + 3] = a1y;
(*a)[aoffset + 4] = a2x;
(*a)[aoffset + 5] = a2y;
return;
}
if (n == 4) {
a0x = (*a)[aoffset + 0];
a0y = (*a)[aoffset + 1];
a1x = (*a)[aoffset + 2];
a1y = (*a)[aoffset + 3];
a2x = (*a)[aoffset + 4];
a2y = (*a)[aoffset + 5];
a3x = (*a)[aoffset + 6];
a3y = (*a)[aoffset + 7];
t1x = a0x + a2x;
t1y = a0y + a2y;
t2x = a1x + a3x;
t2y = a1y + a3y;
m2x = a0x - a2x;
m2y = a0y - a2y;
m3x = a1y - a3y;
m3y = a3x - a1x;
(*a)[aoffset + 0] = t1x + t2x;
(*a)[aoffset + 1] = t1y + t2y;
(*a)[aoffset + 4] = t1x - t2x;
(*a)[aoffset + 5] = t1y - t2y;
(*a)[aoffset + 2] = m2x + m3x;
(*a)[aoffset + 3] = m2y + m3y;
(*a)[aoffset + 6] = m2x - m3x;
(*a)[aoffset + 7] = m2y - m3y;
return;
}
if (n == 5) {
offs = curplan.plan[entryoffset + 7];
c1 = curplan.precomputed[offs + 0];
c2 = curplan.precomputed[offs + 1];
c3 = curplan.precomputed[offs + 2];
c4 = curplan.precomputed[offs + 3];
c5 = curplan.precomputed[offs + 4];
t1x = (*a)[aoffset + 2] + (*a)[aoffset + 8];
t1y = (*a)[aoffset + 3] + (*a)[aoffset + 9];
t2x = (*a)[aoffset + 4] + (*a)[aoffset + 6];
t2y = (*a)[aoffset + 5] + (*a)[aoffset + 7];
t3x = (*a)[aoffset + 2] - (*a)[aoffset + 8];
t3y = (*a)[aoffset + 3] - (*a)[aoffset + 9];
t4x = (*a)[aoffset + 6] - (*a)[aoffset + 4];
t4y = (*a)[aoffset + 7] - (*a)[aoffset + 5];
t5x = t1x + t2x;
t5y = t1y + t2y;
(*a)[aoffset + 0] = (*a)[aoffset + 0] + t5x;
(*a)[aoffset + 1] = (*a)[aoffset + 1] + t5y;
m1x = c1 * t5x;
m1y = c1 * t5y;
m2x = c2 * (t1x - t2x);
m2y = c2 * (t1y - t2y);
m3x = -c3 * (t3y + t4y);
m3y = c3 * (t3x + t4x);
m4x = -c4 * t4y;
m4y = c4 * t4x;
m5x = -c5 * t3y;
m5y = c5 * t3x;
s3x = m3x - m4x;
s3y = m3y - m4y;
s5x = m3x + m5x;
s5y = m3y + m5y;
s1x = (*a)[aoffset + 0] + m1x;
s1y = (*a)[aoffset + 1] + m1y;
s2x = s1x + m2x;
s2y = s1y + m2y;
s4x = s1x - m2x;
s4y = s1y - m2y;
(*a)[aoffset + 2] = s2x + s3x;
(*a)[aoffset + 3] = s2y + s3y;
(*a)[aoffset + 4] = s4x + s5x;
(*a)[aoffset + 5] = s4y + s5y;
(*a)[aoffset + 6] = s4x - s5x;
(*a)[aoffset + 7] = s4y - s5y;
(*a)[aoffset + 8] = s2x - s3x;
(*a)[aoffset + 9] = s2y - s3y;
return;
}
}
if (curplan.plan[entryoffset + 3] == ftbase_fhtcodeletplan) {
n1 = curplan.plan[entryoffset + 1];
n2 = curplan.plan[entryoffset + 2];
n = n1 * n2;
if (n == 2) {
a0x = (*a)[aoffset + 0];
a1x = (*a)[aoffset + 1];
(*a)[aoffset + 0] = a0x + a1x;
(*a)[aoffset + 1] = a0x - a1x;
return;
}
if (n == 3) {
offs = curplan.plan[entryoffset + 7];
c1 = curplan.precomputed[offs + 0];
c2 = curplan.precomputed[offs + 1];
a0x = (*a)[aoffset + 0];
a1x = (*a)[aoffset + 1];
a2x = (*a)[aoffset + 2];
t1x = a1x + a2x;
a0x = a0x + t1x;
m1x = c1 * t1x;
m2y = c2 * (a2x - a1x);
s1x = a0x + m1x;
(*a)[aoffset + 0] = a0x;
(*a)[aoffset + 1] = s1x - m2y;
(*a)[aoffset + 2] = s1x + m2y;
return;
}
if (n == 4) {
a0x = (*a)[aoffset + 0];
a1x = (*a)[aoffset + 1];
a2x = (*a)[aoffset + 2];
a3x = (*a)[aoffset + 3];
t1x = a0x + a2x;
t2x = a1x + a3x;
m2x = a0x - a2x;
m3y = a3x - a1x;
(*a)[aoffset + 0] = t1x + t2x;
(*a)[aoffset + 1] = m2x - m3y;
(*a)[aoffset + 2] = t1x - t2x;
(*a)[aoffset + 3] = m2x + m3y;
return;
}
if (n == 5) {
offs = curplan.plan[entryoffset + 7];
c1 = curplan.precomputed[offs + 0];
c2 = curplan.precomputed[offs + 1];
c3 = curplan.precomputed[offs + 2];
c4 = curplan.precomputed[offs + 3];
c5 = curplan.precomputed[offs + 4];
t1x = (*a)[aoffset + 1] + (*a)[aoffset + 4];
t2x = (*a)[aoffset + 2] + (*a)[aoffset + 3];
t3x = (*a)[aoffset + 1] - (*a)[aoffset + 4];
t4x = (*a)[aoffset + 3] - (*a)[aoffset + 2];
t5x = t1x + t2x;
v0 = (*a)[aoffset + 0] + t5x;
(*a)[aoffset + 0] = v0;
m2x = c2 * (t1x - t2x);
m3y = c3 * (t3x + t4x);
s3y = m3y - c4 * t4x;
s5y = m3y + c5 * t3x;
s1x = v0 + c1 * t5x;
s2x = s1x + m2x;
s4x = s1x - m2x;
(*a)[aoffset + 1] = s2x - s3y;
(*a)[aoffset + 2] = s4x - s5y;
(*a)[aoffset + 3] = s4x + s5y;
(*a)[aoffset + 4] = s2x + s3y;
return;
}
}
if (curplan.plan[entryoffset + 3] == ftbase_fftbluesteinplan) {
n = curplan.plan[entryoffset + 1];
m = curplan.plan[entryoffset + 4];
offs = curplan.plan[entryoffset + 7];
for (int i = stackptr + 2 * n; i <= stackptr + 2 * m - 1; i++)
curplan.stackbuf[i] = 0;
offsp = offs + 2 * m;
offsa = aoffset;
offsb = stackptr;
for (int i = 0; i < n; i++) {
bx = curplan.precomputed[offsp + 0];
by = curplan.precomputed[offsp + 1];
x = (*a)[offsa + 0];
y = (*a)[offsa + 1];
curplan.stackbuf[offsb + 0] = x * bx - y * (-by);
curplan.stackbuf[offsb + 1] = x * (-by) + y * bx;
offsp = offsp + 2;
offsa = offsa + 2;
offsb = offsb + 2;
}
ftbaseexecuteplanrec(&curplan.stackbuf, stackptr, plan, curplan.plan[entryoffset + 5], stackptr + 2 * 2 * m);
offsb = stackptr;
offsp = offs;
for (int i = 0; i <= m - 1; i++) {
x = curplan.stackbuf[offsb + 0];
y = curplan.stackbuf[offsb + 1];
bx = curplan.precomputed[offsp + 0];
by = curplan.precomputed[offsp + 1];
curplan.stackbuf[offsb + 0] = x * bx - y * by;
curplan.stackbuf[offsb + 1] = -(x * by + y * bx);
offsb = offsb + 2;
offsp = offsp + 2;
}
ftbaseexecuteplanrec(&curplan.stackbuf, stackptr, plan, curplan.plan[entryoffset + 5], stackptr + 2 * 2 * m);
offsb = stackptr;
offsp = offs + 2 * m;
offsa = aoffset;
for (int i = 0; i < n; i++) {
x = curplan.stackbuf[offsb + 0] / m;
y = -curplan.stackbuf[offsb + 1] / m;
bx = curplan.precomputed[offsp + 0];
by = curplan.precomputed[offsp + 1];
(*a)[offsa + 0] = x * bx - y * (-by);
(*a)[offsa + 1] = x * (-by) + y * bx;
offsp = offsp + 2;
offsa = offsa + 2;
offsb = offsb + 2;
}
return;
}
}
/*************************************************************************
Twiddle factors calculation
-- ALGLIB --
Copyright 01.05.2009 by Bochkanov Sergey
*************************************************************************/
void PIFFT_float::ftbase_ffttwcalc(PIVector<float> * a, int aoffset, int n1, int n2) {
int n, idx, offs;
float x, y, twxm1, twy, twbasexm1, twbasey, twrowxm1, twrowy, tmpx, tmpy, v;
int ftbase_ftbaseupdatetw = 4;
n = n1 * n2;
v = -2 * M_PI / n;
twbasexm1 = -2 * sqr(sin(0.5 * v));
twbasey = sin(v);
twrowxm1 = 0;
twrowy = 0;
for (int i = 0, j = 0; i <= n2 - 1; i++) {
twxm1 = 0;
twy = 0;
for (j = 0; j <= n1 - 1; j++) {
idx = i * n1 + j;
offs = aoffset + 2 * idx;
x = (*a)[offs + 0];
y = (*a)[offs + 1];
tmpx = x * twxm1 - y * twy;
tmpy = x * twy + y * twxm1;
(*a)[offs + 0] = x + tmpx;
(*a)[offs + 1] = y + tmpy;
if (j < n1 - 1) {
if (j % ftbase_ftbaseupdatetw == 0) {
v = -2 * M_PI * i * (j + 1) / n;
twxm1 = -2 * sqr(sin(0.5 * v));
twy = sin(v);
} else {
tmpx = twrowxm1 + twxm1 * twrowxm1 - twy * twrowy;
tmpy = twrowy + twxm1 * twrowy + twy * twrowxm1;
twxm1 = twxm1 + tmpx;
twy = twy + tmpy;
}
}
}
if (i < n2 - 1) {
if (j % ftbase_ftbaseupdatetw == 0) {
v = -2 * M_PI * (i + 1) / n;
twrowxm1 = -2 * sqr(sin(0.5 * v));
twrowy = sin(v);
} else {
tmpx = twbasexm1 + twrowxm1 * twbasexm1 - twrowy * twbasey;
tmpy = twbasey + twrowxm1 * twbasey + twrowy * twbasexm1;
twrowxm1 = twrowxm1 + tmpx;
twrowy = twrowy + tmpy;
}
}
}
}

196
lib/main/math/pifft.h Normal file
View File

@@ -0,0 +1,196 @@
/*! \file pifft.h
* \brief Class for FFT, IFFT and Hilbert transformations
*/
/*
PIP - Platform Independent Primitives
Class for FFT, IFFT and Hilbert transformations
Andrey Bychkov work.a.b@yandex.ru, Ivan Pelipenko peri4ko@yandex.ru
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PIFFT_H
#define PIFFT_H
#include "pimathcomplex.h"
class PIP_EXPORT PIFFT_double
{
public:
PIFFT_double();
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() const;
PIVector<double> getReal() const;
PIVector<double> getImag() const;
private:
PIVector<complexd> result;
typedef ptrdiff_t ae_int_t;
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);
};
class PIP_EXPORT PIFFT_float
{
public:
PIFFT_float();
PIVector<complexf> * calcFFT(const PIVector<complexf> &val);
PIVector<complexf> * calcFFT(const PIVector<float> &val);
PIVector<complexf> * calcFFTinverse(const PIVector<complexf> &val);
PIVector<complexf> * calcHilbert(const PIVector<float> &val);
PIVector<float> getAmplitude() const;
PIVector<float> getReal() const;
PIVector<float> getImag() const;
private:
PIVector<complexf> result;
typedef ptrdiff_t ae_int_t;
struct ftplan {
PIVector<int> plan;
PIVector<float> precomputed;
PIVector<float> tmpbuf;
PIVector<float> stackbuf;
};
ftplan curplan;
void fftc1d(const PIVector<complexf> &a, uint n);
void fftc1r(const PIVector<float> &a, uint n);
void fftc1dinv(const PIVector<complexf> &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<float> *a, int aoffset, int n, ftplan *plan);
void ftbaseexecuteplanrec(PIVector<float> *a, int aoffset, ftplan *plan, int entryoffset, ae_int_t stackptr);
void ftbase_internalcomplexlintranspose(PIVector<float> *a, int m, int n, int astart, PIVector<float> *buf);
void ftbase_ffticltrec(PIVector<float> *a, int astart, int astride, PIVector<float> *b, int bstart, int bstride, int m, int n);
void ftbase_internalreallintranspose(PIVector<float> *a, int m, int n, int astart, PIVector<float> *buf);
void ftbase_fftirltrec(PIVector<float> *a, int astart, int astride, PIVector<float> *b, int bstart, int bstride, int m, int n);
void ftbase_ffttwcalc(PIVector<float> *a, int aoffset, int n1, int n2);
};
typedef PIFFT_double PIFFT;
typedef PIFFT_double PIFFTd;
typedef PIFFT_float PIFFTf;
#ifndef CC_VC
#define _PIFFTW_H(type) class _PIFFTW_P_##type##_ { \
public: \
_PIFFTW_P_##type##_(); \
~_PIFFTW_P_##type##_(); \
const PIVector<complex<type> > & calcFFT(const PIVector<complex<type> > & in); \
const PIVector<complex<type> > & calcFFTR(const PIVector<type> & in); \
const PIVector<complex<type> > & calcFFTI(const PIVector<complex<type> > & in); \
void preparePlan(int size, int op); \
void * impl; \
};
_PIFFTW_H(float)
_PIFFTW_H(double)
_PIFFTW_H(ldouble)
template <typename T>
class PIP_EXPORT PIFFTW
{
public:
explicit PIFFTW() {p = 0; newP(p);}
~PIFFTW() {deleteP(p);}
inline const PIVector<complex<T> > & calcFFT(const PIVector<complex<T> > & in) {return PIVector<complex<T> >().resize(in.size());}
inline const PIVector<complex<T> > & calcFFT(const PIVector<T> & in) {return PIVector<complex<T> >().resize(in.size());}
inline const PIVector<complex<T> > & calcFFTinverse(const PIVector<complex<T> > & in) {return PIVector<complex<T> >().resize(in.size());}
enum FFT_Operation {foReal, foComplex, foInverse};
inline void preparePlan(int size, FFT_Operation op) {}
private:
void operator =(const PIFFTW & );
PIFFTW(const PIFFTW &);
inline void newP(void *& _p) {}
inline void deleteP(void *& _p) {}
void * p;
};
template<> inline const PIVector<complex<float> > & PIFFTW<float>::calcFFT(const PIVector<complex<float> > & in) {return ((_PIFFTW_P_float_*)p)->calcFFT(in);}
template<> inline const PIVector<complex<float> > & PIFFTW<float>::calcFFT(const PIVector<float> & in) {return ((_PIFFTW_P_float_*)p)->calcFFTR(in);}
template<> inline const PIVector<complex<float> > & PIFFTW<float>::calcFFTinverse(const PIVector<complex<float> > & in) {return ((_PIFFTW_P_float_*)p)->calcFFTI(in);}
template<> inline void PIFFTW<float>::preparePlan(int size, FFT_Operation op) {((_PIFFTW_P_float_*)p)->preparePlan(size, op);}
template<> inline void PIFFTW<float>::newP(void *& _p) {_p = new _PIFFTW_P_float_();}
template<> inline void PIFFTW<float>::deleteP(void *& _p) {if (_p) delete (_PIFFTW_P_float_*)_p; _p = 0;}
typedef PIFFTW<float> PIFFTWf;
template<> inline const PIVector<complex<double> > & PIFFTW<double>::calcFFT(const PIVector<complex<double> > & in) {return ((_PIFFTW_P_double_*)p)->calcFFT(in);}
template<> inline const PIVector<complex<double> > & PIFFTW<double>::calcFFT(const PIVector<double> & in) {return ((_PIFFTW_P_double_*)p)->calcFFTR(in);}
template<> inline const PIVector<complex<double> > & PIFFTW<double>::calcFFTinverse(const PIVector<complex<double> > & in) {return ((_PIFFTW_P_double_*)p)->calcFFTI(in);}
template<> inline void PIFFTW<double>::preparePlan(int size, FFT_Operation op) {((_PIFFTW_P_double_*)p)->preparePlan(size, op);}
template<> inline void PIFFTW<double>::newP(void *& _p) {_p = new _PIFFTW_P_double_();}
template<> inline void PIFFTW<double>::deleteP(void *& _p) {if (_p) delete (_PIFFTW_P_double_*)_p; _p = 0;}
typedef PIFFTW<double> PIFFTWd;
template<> inline const PIVector<complex<ldouble> > & PIFFTW<ldouble>::calcFFT(const PIVector<complex<ldouble> > & in) {return ((_PIFFTW_P_ldouble_*)p)->calcFFT(in);}
template<> inline const PIVector<complex<ldouble> > & PIFFTW<ldouble>::calcFFT(const PIVector<ldouble> & in) {return ((_PIFFTW_P_ldouble_*)p)->calcFFTR(in);}
template<> inline const PIVector<complex<ldouble> > & PIFFTW<ldouble>::calcFFTinverse(const PIVector<complex<ldouble> > & in) {return ((_PIFFTW_P_ldouble_*)p)->calcFFTI(in);}
template<> inline void PIFFTW<ldouble>::preparePlan(int size, FFT_Operation op) {((_PIFFTW_P_ldouble_*)p)->preparePlan(size, op);}
template<> inline void PIFFTW<ldouble>::newP(void *& _p) {_p = new _PIFFTW_P_ldouble_();}
template<> inline void PIFFTW<ldouble>::deleteP(void *& _p) {if (_p) delete (_PIFFTW_P_ldouble_*)_p; _p = 0;}
typedef PIFFTW<ldouble> PIFFTWld;
#endif
#endif // PIFFT_H

158
lib/main/math/pigeometry.h Normal file
View File

@@ -0,0 +1,158 @@
/*! \file pigeometry.h
* \brief Geometry base class
*/
/*
PIP - Platform Independent Primitives
Geometry
Ivan Pelipenko peri4ko@yandex.ru
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PIGEOMETRY_H
#define PIGEOMETRY_H
#include "pimathbase.h"
template<typename Type>
class PIP_EXPORT PIPoint {
public:
Type x;
Type y;
PIPoint() {x = y = Type();}
PIPoint(Type x_, Type y_) {set(x_, y_);}
PIPoint<Type> & set(Type x_, Type y_) {x = x_; y = y_; return *this;}
PIPoint<Type> & move(Type x_, Type y_) {x += x_; y += y_; return *this;}
PIPoint<Type> & move(const PIPoint<Type> & p) {x += p.x; y += p.y; return *this;}
double angleRad() const {return atan2(y, x);}
int angleDeg() const {return round(atan2(y, x) * rad2deg);}
PIPoint<Type> toPolar(bool isDeg = false) const {return PIPoint<Type>(sqrt(x*x + y*y), isDeg ? angleDeg() : angleRad());}
static PIPoint<Type> fromPolar(const PIPoint<Type> & p) {return PIPoint<Type>(p.y * cos(p.x), p.y * sin(p.x));}
PIPoint<Type> operator +(const PIPoint<Type> & p) {return PIPoint<Type>(x + p.x, y + p.y);}
PIPoint<Type> operator +(const Type & p) {return PIPoint<Type>(x + p, y + p);}
PIPoint<Type> operator -(const PIPoint<Type> & p) {return PIPoint<Type>(x - p.x, y - p.y);}
PIPoint<Type> operator -(const Type & p) {return PIPoint<Type>(x - p, y - p);}
PIPoint<Type> operator -() {return PIPoint<Type>(-x, -y);}
PIPoint<Type> operator *(const Type & d) {return PIPoint<Type>(x * d, y * d);}
PIPoint<Type> operator /(const Type & d) {return PIPoint<Type>(x / d, y / d);}
bool operator ==(const PIPoint<Type> & p) const {return (x == p.x && y == p.y);}
bool operator !=(const PIPoint<Type> & p) const {return (x != p.x || y != p.y);}
};
template<typename Type>
PICout operator <<(PICout & s, const PIPoint<Type> & v) {s.setControl(0, true); s << "Point{" << v.x << ", " << v.y << "}"; s.restoreControl(); return s;}
template<typename Type>
inline PIByteArray & operator <<(PIByteArray & s, const PIPoint<Type> & v) {s << v.x << v.y; return s;}
template<typename Type>
inline PIByteArray & operator >>(PIByteArray & s, PIPoint<Type> & v) {s >> v.x >> v.y; return s;}
typedef PIPoint<int> PIPointi;
typedef PIPoint<uint> PIPointu;
typedef PIPoint<float> PIPointf;
typedef PIPoint<double> PIPointd;
template<typename Type>
class PIP_EXPORT PIRect {
public:
Type x0;
Type y0;
Type x1;
Type y1;
PIRect() {x0 = y0 = x1 = y1 = Type();}
PIRect(Type x, Type y, Type w, Type h) {set(x, y, w, h);}
PIRect(const PIPoint<Type> & tl, const PIPoint<Type> & br) {set(tl.x, tl.y, br.x, br.y);}
PIRect(const PIPoint<Type> & p0, const PIPoint<Type> & p1, const PIPoint<Type> & p2) {set(piMin<Type>(p0.x, p1.x, p2.x), piMin<Type>(p0.y, p1.y, p2.y),
piMax<Type>(p0.x, p1.x, p2.x), piMax<Type>(p0.y, p1.y, p2.y));}
PIRect<Type> & set(Type x, Type y, Type w, Type h) {x0 = x; y0 = y; x1 = x + w; y1 = y + h; return *this;}
bool pointIn(Type x, Type y) const {return (x <= x1 && x >= x0 && y <= y1 && y >= y0);}
bool pointIn(const PIPoint<Type> & p) const {return pointIn(p.x, p.y);}
bool isEmpty() const {return (x1 - x0 == 0 && y1 - y0 == 0);}
PIRect<Type> & translate(Type x, Type y) {x0 += x; x1 += x; y0 += y; y1 += y; return *this;}
PIRect<Type> & translate(const PIPoint<Type> & p) {x0 += p.x; x1 += p.x; y0 += p.y; y1 += p.y; return *this;}
PIRect<Type> translated(Type x, Type y) {PIRect<Type> r(*this); r.translate(x, y); return r;}
PIRect<Type> translated(const PIPoint<Type> & p) {PIRect<Type> r(*this); r.translate(p); return r;}
PIRect<Type> & scale(Type x, Type y) {setWidth(width() * x); setHeight(height() * y); return *this;}
PIRect<Type> & scale(const PIPoint<Type> & p) {setWidth(width() * p.x); setHeight(height() * p.y); return *this;}
PIRect<Type> scaled(Type x, Type y) {PIRect<Type> r(*this); r.scale(x, y); return r;}
PIRect<Type> scaled(const PIPoint<Type> & p) {PIRect<Type> r(*this); r.scale(p); return r;}
PIRect<Type> & normalize() {if (x0 > x1) piSwap<Type>(x0, x1); if (y0 > y1) piSwap<Type>(y0, y1); return *this;}
PIRect<Type> normalized() {PIRect<Type> r(*this); r.normalize(); return r;}
PIRect<Type> & unite(const PIRect<Type> & r) {x0 = piMin<Type>(x0, r.x0); y0 = piMin<Type>(y0, r.y0); x1 = piMax<Type>(x1, r.x1); y1 = piMax<Type>(y1, r.y1); return *this;}
PIRect<Type> united(const PIRect<Type> & rect) {PIRect<Type> r(*this); r.unite(rect); return r;}
PIRect<Type> & intersect(const PIRect<Type> & r) {x0 = piMax<Type>(x0, r.x0); y0 = piMax<Type>(y0, r.y0); x1 = piMin<Type>(x1, r.x1); y1 = piMin<Type>(y1, r.y1); if (x0 > x1 || y0 > y1) x0 = x1 = y0 = y1 = Type(0); return *this;}
PIRect<Type> intersected(const PIRect<Type> & rect) {PIRect<Type> r(*this); r.intersect(rect); return r;}
Type top() const {return y0;}
Type left() const {return x0;}
Type right() const {return x1;}
Type bottom() const {return y1;}
Type width() const {return x1 - x0;}
Type height() const {return y1 - y0;}
PIPoint<Type> topLeft() {return PIPoint<Type>(x0, y0);}
PIPoint<Type> topRigth() {return PIPoint<Type>(x1, y0);}
PIPoint<Type> bottomLeft() {return PIPoint<Type>(x0, y1);}
PIPoint<Type> bottomRight() {return PIPoint<Type>(x1, y1);}
void setTop(Type v) {y0 = v;}
void setLeft(Type v) {x0 = v;}
void setRigth(Type v) {x1 = v;}
void setBottom(Type v) {y1 = v;}
void setWidth(Type v) {x1 = x0 + v;}
void setHeight(Type v) {y1 = y0 + v;}
PIRect<Type> operator -() {return PIRect<Type>(-x0, -y0, -width(), -height());}
void operator +=(Type x) {translate(x, x);}
void operator +=(const PIPoint<Type> & p) {translate(p);}
void operator -=(Type x) {translate(-x, -x);}
void operator -=(const PIPoint<Type> & p) {translate(-p);}
void operator *=(Type p) {x0 *= p; x1 *= p; y0 *= p; y1 *= p;}
void operator /=(Type p) {x0 /= p; x1 /= p; y0 /= p; y1 /= p;}
void operator |=(const PIRect<Type> & r) {unite(r);}
void operator &=(const PIRect<Type> & r) {intersect(r);}
PIRect<Type> operator +(const PIPoint<Type> & p) {return PIRect<Type>(*this).translated(p);}
PIRect<Type> operator -(const PIPoint<Type> & p) {return PIRect<Type>(*this).translated(-p);}
PIRect<Type> operator |(const PIRect<Type> & r) {return PIRect<Type>(*this).united(r);}
PIRect<Type> operator &(const PIRect<Type> & r) {return PIRect<Type>(*this).intersected(r);}
bool operator ==(const PIRect<Type> & r) const {return (x0 == r.x0 && y0 == r.y0 && x1 == r.x1 && y1 == r.y10);}
bool operator !=(const PIRect<Type> & r) const {return (x0 != r.x0 || y0 != r.y0 || x1 != r.x1 || y1 != r.y10);}
};
template<typename Type>
PICout operator <<(PICout & s, const PIRect<Type> & v) {s.setControl(0, true); s << "Rect{" << v.x0 << ", " << v.y0 << "; " << v.x1 - v.x0 << ", " << v.y1 - v.y0 << "}"; s.restoreControl(); return s;}
template<typename Type>
inline PIByteArray & operator <<(PIByteArray & s, const PIRect<Type> & v) {s << v.x0 << v.x1 << v.y0 << v.y1; return s;}
template<typename Type>
inline PIByteArray & operator >>(PIByteArray & s, PIRect<Type> & v) {s >> v.x0 >> v.x1 >> v.y0 >> v.y1; return s;}
typedef PIRect<int> PIRecti;
typedef PIRect<uint> PIRectu;
typedef PIRect<float> PIRectf;
typedef PIRect<double> PIRectd;
#endif // PIGEOMETRY_H

View File

@@ -0,0 +1,475 @@
/*
PIP - Platform Independent Primitives
Basic mathematical functions and defines
Ivan Pelipenko peri4ko@yandex.ru
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "pimathbase.h"
#include "pitime.h"
#include <stdlib.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;
}
double randomd() {
return (double) randomi() / RAND_MAX * 2. - 1.;
}

181
lib/main/math/pimathbase.h Normal file
View File

@@ -0,0 +1,181 @@
/*! \file pimathbase.h
* \brief Basic mathematical functions and defines
*/
/*
PIP - Platform Independent Primitives
Basic mathematical functions and defines
Ivan Pelipenko peri4ko@yandex.ru
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PIMATHBASE_H
#define PIMATHBASE_H
#include "piinit.h"
#include "pivector.h"
#include "pipair.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
# include <math.h>
#else
# include <cmath>
#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.141592653589793238462643383280
#endif
#ifndef M_2PI
# define M_2PI 6.283185307179586476925286766559
#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_SQRT_PI
# define M_SQRT_PI 1.772453850905516027298167483341
#endif
#ifndef M_E
# define M_E 2.7182818284590452353602874713527
#endif
#ifndef M_LIGHT_SPEED
# define M_LIGHT_SPEED 2.99792458e+8
#endif
#ifndef M_RELATIVE_CONST
# define M_RELATIVE_CONST -4.442807633e-10;
#endif
#ifndef M_GRAVITY_CONST
# define M_GRAVITY_CONST 398600.4418e9;
#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 int pow2(const int p) {return 1 << p;}
inline int sqr(const int v) {return v * v;}
inline float 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;}
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;}
// [-1 ; 1]
double randomd();
// [-1 ; 1] normal
double randomn(double dv = 0., double sv = 1.);
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] = fabs(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

View File

@@ -0,0 +1,125 @@
/*! \file pimathcomplex.h
* \brief PIP math complex
*/
/*
PIP - Platform Independent Primitives
PIP math complex
Ivan Pelipenko peri4ko@yandex.ru, 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 Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PIMATHCOMPLEX_H
#define PIMATHCOMPLEX_H
#include <complex>
#include "pimathbase.h"
#include "pibytearray.h"
#include "pivector2d.h"
#define PIP_MATH_COMPLEX
using std::complex;
typedef complex<int> complexi;
typedef complex<short> complexs;
typedef complex<float> complexf;
typedef complex<ldouble> complexld;
#ifndef QPIEVALUATOR_COMPLEX
typedef complex<double> complexd;
const complexd complexd_i(0., 1.);
const complexd complexd_0(0.);
const complexd complexd_1(1.);
#endif
const complexld complexld_i(0., 1.);
const complexld complexld_0(0.);
const complexld complexld_1(1.);
__PICONTAINERS_SIMPLE_TYPE__(complexi)
__PICONTAINERS_SIMPLE_TYPE__(complexs)
__PICONTAINERS_SIMPLE_TYPE__(complexf)
__PICONTAINERS_SIMPLE_TYPE__(complexd)
__PICONTAINERS_SIMPLE_TYPE__(complexld)
__PIVECTOR2D_SIMPLE_TYPE__(complexi)
__PIVECTOR2D_SIMPLE_TYPE__(complexs)
__PIVECTOR2D_SIMPLE_TYPE__(complexf)
__PIVECTOR2D_SIMPLE_TYPE__(complexd)
__PIVECTOR2D_SIMPLE_TYPE__(complexld)
__PIBYTEARRAY_SIMPLE_TYPE__(complexi)
__PIBYTEARRAY_SIMPLE_TYPE__(complexs)
__PIBYTEARRAY_SIMPLE_TYPE__(complexf)
__PIBYTEARRAY_SIMPLE_TYPE__(complexd)
__PIBYTEARRAY_SIMPLE_TYPE__(complexld)
inline complexd sign(const complexd & x) {return complexd(sign(x.real()), sign(x.imag()));}
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()));}
#ifdef PIP_CXX11_SUPPORT
# define acosc acos
# define asinc asin
# define atanc atan
#else
inline complexd atanc(const complexd & c) {return complexd(0., 0.5) * 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));}
#endif
#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
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;}
//! \relatesalso PIByteArray \brief Store operator
inline PIByteArray & operator <<(PIByteArray & s, complexf v) {float t; t = v.real(); s << t; t = v.imag(); s << t; return s;}
//! \relatesalso PIByteArray \brief Store operator
inline PIByteArray & operator <<(PIByteArray & s, complexd v) {double t; t = v.real(); s << t; t = v.imag(); s << t; return s;}
//! \relatesalso PIByteArray \brief Store operator
inline PIByteArray & operator <<(PIByteArray & s, complexld v) {ldouble t; t = v.real(); s << t; t = v.imag(); s << t; return s;}
//! \relatesalso PIByteArray \brief Restore operator
inline PIByteArray & operator >>(PIByteArray & s, complexf & v) {float t0, t1; s >> t0; s >> t1; v = complexf(t0, t1); return s;}
//! \relatesalso PIByteArray \brief Restore operator
inline PIByteArray & operator >>(PIByteArray & s, complexd & v) {double t0, t1; s >> t0; s >> t1; v = complexd(t0, t1); return s;}
//! \relatesalso PIByteArray \brief Restore operator
inline PIByteArray & operator >>(PIByteArray & s, complexld & v) {ldouble t0, t1; s >> t0; s >> t1; v = complexld(t0, t1); return s;}
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 PIVector2D<double> abs(const PIVector2D<complexd> & v) {
PIVector2D<double> result(v.rows(), v.cols());
for (uint i = 0; i < v.rows(); i++)
for (uint j = 0; j < v.cols(); j++)
result[i][j] = abs(v.element(i,j));
return result;
}
#endif // PIMATHCOMPLEX_H

View File

@@ -0,0 +1,574 @@
/*! \file pimathmatrix.h
* \brief PIMathMatrix
*/
/*
PIP - Platform Independent Primitives
PIMathMatrix
Ivan Pelipenko peri4ko@yandex.ru, 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 Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PIMATHMATRIX_H
#define PIMATHMATRIX_H
#include "pimathvector.h"
#include "pimathcomplex.h"
template<typename T>
inline bool _PIMathMatrixNullCompare(const T v) {
return (piAbs(v) < T(1E-200));
}
template<>
inline bool _PIMathMatrixNullCompare<complexf >(const complexf v) {
return (abs(v) < float(1E-200));
}
template<>
inline bool _PIMathMatrixNullCompare<complexd >(const complexd v) {
return (abs(v) < double(1E-200));
}
/// 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(const PIVector<Type> & val) {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 filled(const Type & v) {_CMatrix tm; PIMM_FOR_WB(r, c) tm.m[r][c] = v; 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;}
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];}
_CMatrix & operator =(const _CMatrix & sm) {memcpy(m, sm.m, sizeof(Type) * Cols * Rows); return *this;}
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 -() const {_CMatrix tm; PIMM_FOR_WB(r, c) tm.m[r][c] = -m[r][c]; return tm;}
_CMatrix operator +(const _CMatrix & sm) const {_CMatrix tm = _CMatrix(*this); PIMM_FOR_WB(r, c) tm.m[r][c] += sm.m[r][c]; return tm;}
_CMatrix operator -(const _CMatrix & sm) const {_CMatrix tm = _CMatrix(*this); PIMM_FOR_WB(r, c) tm.m[r][c] -= sm.m[r][c]; return tm;}
_CMatrix operator *(const Type & v) const {_CMatrix tm = _CMatrix(*this); PIMM_FOR_WB(r, c) tm.m[r][c] *= v; return tm;}
_CMatrix operator /(const Type & v) const {_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 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-200)) {
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 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];
}
if (i < Cols - 1) {
if (fabs(smat.m[i+1][i+1]) < Type(1E-200)) {
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[c][r] = 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;}
#ifdef PIP_STD_IOSTREAM
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 << std::endl << " ";} s << "}"; return s;}
#endif
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 << PICoutManipulators::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 {Cols}, result is vector {Rows}
template<uint Cols, uint Rows, typename Type>
inline PIMathVectorT<Rows, Type> operator *(const PIMathMatrixT<Rows, Cols, Type> & fm,
const PIMathVectorT<Cols, 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;
}
/// Multiply vector {Rows} on matrix {Rows x Cols}, result is vector {Cols}
template<uint Cols, uint Rows, typename Type>
inline PIMathVectorT<Cols, Type> operator *(const PIMathVectorT<Rows, Type> & sv,
const PIMathMatrixT<Rows, Cols, Type> & fm) {
PIMathVectorT<Cols, Type> tv;
Type t;
for (uint j = 0; j < Cols; ++j) {
t = Type(0);
for (uint i = 0; i < Rows; ++i)
t += fm[i][j] * sv[i];
tv[j] = t;
}
return tv;
}
/// Multiply value(T) on matrix {Rows x Cols}, result is vector {Rows}
template<uint Cols, uint Rows, typename Type>
inline PIMathMatrixT<Rows, Cols, Type> operator *(const Type & x, const PIMathMatrixT<Rows, Cols, Type> & v) {
return v * x;
}
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 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 < _V2D::cols_; ++c) for (uint r = 0; r < _V2D::rows_; ++r)
#define PIMM_FOR_I(c, r) for (uint r = 0; r < _V2D::rows_; ++r) for (uint c = 0; c < _V2D::cols_; ++c)
#define PIMM_FOR_A(v) for (uint v = 0; v < _V2D::mat.size(); ++v)
#define PIMM_FOR_C(v) for (uint v = 0; v < _V2D::cols_; ++v)
#define PIMM_FOR_R(v) for (uint v = 0; v < _V2D::rows_; ++v)
template<typename Type>
class PIP_EXPORT PIMathMatrix : public PIVector2D<Type> {
typedef PIVector2D<Type> _V2D;
typedef PIMathMatrix<Type> _CMatrix;
typedef PIMathVector<Type> _CMCol;
public:
PIMathMatrix(const uint cols = 0, const uint rows = 0, const Type & f = Type()) {_V2D::resize(rows, cols, f);}
PIMathMatrix(const uint cols, const uint rows, const PIVector<Type> & val) {_V2D::resize(rows, cols); int i=0; PIMM_FOR_I(c, r) _V2D::element(r, c) = val[i++];}
PIMathMatrix(const PIVector<PIVector<Type> > & val) {if(!val.isEmpty()) {_V2D::resize(val.size(), val[0].size()); PIMM_FOR_I(c, r) _V2D::element(r, c) = val[r][c];}}
PIMathMatrix(const PIVector2D<Type> & val) {if(!val.isEmpty()) {_V2D::resize(val.rows(), val.cols()); PIMM_FOR_I(c, r) _V2D::element(r, c) = val.element(r, c);}}
static _CMatrix identity(const uint cols, const uint rows) {_CMatrix tm(cols, rows); for (uint r = 0; r < rows; ++r) for (uint c = 0; c < cols; ++c) tm.element(r, c) = (c == r ? Type(1) : Type(0)); return tm;}
static _CMatrix matrixRow(const PIMathVector<Type> & val) {return _CMatrix(val.size(), 1, val.toVector());}
static _CMatrix matrixCol(const PIMathVector<Type> & val) {return _CMatrix(1, val.size(), val.toVector());}
_CMatrix & setCol(uint index, const _CMCol & v) {PIMM_FOR_R(i) _V2D::element(i, index) = v[i]; return *this;}
_CMatrix & setRow(uint index, const _CMCol & v) {PIMM_FOR_C(i) _V2D::element(index, i) = v[i]; return *this;}
_CMatrix & swapCols(uint r0, uint r1) {PIMM_FOR_C(i) {piSwap(_V2D::element(i, r0), _V2D::element(i, r1));} return *this;}
_CMatrix & swapRows(uint c0, uint c1) {PIMM_FOR_R(i) {piSwap(_V2D::element(c0, i), _V2D::element(c1, i));} return *this;}
_CMatrix & fill(const Type & v) {PIMM_FOR_A(i) _V2D::mat[i] = v; return *this;}
bool isSquare() const {return _V2D::cols_ == _V2D::rows_;}
bool isIdentity() const {PIMM_FOR(c, r) if ((c == r) ? _V2D::element(r, c) != Type(1) : _V2D::element(r, c) != Type(0)) return false; return true;}
bool isNull() const {PIMM_FOR_A(i) if (_V2D::mat[i] != Type(0)) return false; return true;}
bool isValid() const {return !PIVector2D<Type>::isEmpty();}
_CMatrix & operator =(const PIVector<PIVector<Type> > & v) {*this = _CMatrix(v); return *this;}
bool operator ==(const _CMatrix & sm) const {PIMM_FOR_A(i) if (_V2D::mat[i] != sm.mat[i]) return false; return true;}
bool operator !=(const _CMatrix & sm) const {return !(*this == sm);}
void operator +=(const _CMatrix & sm) {PIMM_FOR_A(i) _V2D::mat[i] += sm.mat[i];}
void operator -=(const _CMatrix & sm) {PIMM_FOR_A(i) _V2D::mat[i] -= sm.mat[i];}
void operator *=(const Type & v) {PIMM_FOR_A(i) _V2D::mat[i] *= v;}
void operator /=(const Type & v) {PIMM_FOR_A(i) _V2D::mat[i] /= v;}
_CMatrix operator -() const {_CMatrix tm(*this); PIMM_FOR_A(i) tm.mat[i] = -_V2D::mat[i]; return tm;}
_CMatrix operator +(const _CMatrix & sm) const {_CMatrix tm(*this); PIMM_FOR_A(i) tm.mat[i] += sm.mat[i]; return tm;}
_CMatrix operator -(const _CMatrix & sm) const {_CMatrix tm(*this); PIMM_FOR_A(i) tm.mat[i] -= sm.mat[i]; return tm;}
_CMatrix operator *(const Type & v) const {_CMatrix tm(*this); PIMM_FOR_A(i) tm.mat[i] *= v; return tm;}
_CMatrix operator /(const Type & v) const {_CMatrix tm(*this); PIMM_FOR_A(i) tm.mat[i] /= 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 < _V2D::cols_; ++c)
for (uint r = 0; r < _V2D::rows_; ++r)
if (r == c)
ret *= m.element(r, c);
return ret;
}
Type trace(bool * ok = 0) const {
Type ret = Type(0);
if (!isSquare()) {
if (ok != 0) *ok = false;
return ret;
}
for (uint i = 0; i < _V2D::cols_; ++i) {
ret += _V2D::element(i, i);
}
if (ok != 0) *ok = true;
return ret;
}
_CMatrix & toUpperTriangular(bool * ok = 0) {
if (!isSquare()) {
if (ok != 0) *ok = false;
return *this;
}
_CMatrix smat(*this);
bool ndet;
uint crow;
Type mul;
for (uint i = 0; i < _V2D::cols_; ++i) {
ndet = true;
for (uint j = 0; j < _V2D::rows_; ++j) if (smat.element(i, j) != 0) ndet = false;
if (ndet) {
if (ok != 0) *ok = false;
return *this;
}
}
for (uint i = 0; i < _V2D::cols_; ++i) {
crow = i;
while (smat.element(i, i) == Type(0))
smat.swapRows(i, ++crow);
for (uint j = i + 1; j < _V2D::rows_; ++j) {
mul = smat.element(i, j) / smat.element(i, i);
for (uint k = i; k < _V2D::cols_; ++k) smat.element(k, j) -= mul * smat.element(k, i);
}
if (i < _V2D::cols_ - 1) {
if (_PIMathMatrixNullCompare(smat.element(i+1, i+1))) {
if (ok != 0) *ok = false;
return *this;
}
}
}
if (ok != 0) *ok = true;
_V2D::mat.swap(smat.mat);
return *this;
}
_CMatrix & invert(bool * ok = 0, _CMCol * sv = 0) {
if (!isSquare()) {
if (ok != 0) *ok = false;
return *this;
}
_CMatrix mtmp = _CMatrix::identity(_V2D::cols_, _V2D::rows_), smat(*this);
bool ndet;
uint crow;
Type mul, iddiv;
for (uint i = 0; i < _V2D::cols_; ++i) {
ndet = true;
for (uint j = 0; j < _V2D::rows_; ++j) if (smat.element(i, j) != Type(0)) ndet = false;
if (ndet) {
if (ok != 0) *ok = false;
return *this;
}
}
for (uint i = 0; i < _V2D::cols_; ++i) {
crow = i;
while (smat.element(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 < _V2D::rows_; ++j) {
mul = smat.element(i, j) / smat.element(i, i);
for (uint k = i; k < _V2D::cols_; ++k) smat.element(k, j) -= mul * smat.element(k, i);
for (uint k = 0; k < _V2D::cols_; ++k) mtmp.element(k, j) -= mul * mtmp.element(k, i);
if (sv != 0) (*sv)[j] -= mul * (*sv)[i];
}
if (i < _V2D::cols_ - 1) {
if (_PIMathMatrixNullCompare(smat.element(i+1, i+1))) {
if (ok != 0) *ok = false;
return *this;
}
}
iddiv = smat.element(i, i);
for (uint j = i; j < _V2D::cols_; ++j) smat.element(j, i) /= iddiv;
for (uint j = 0; j < _V2D::cols_; ++j) mtmp.element(j, i) /= iddiv;
if (sv != 0) (*sv)[i] /= iddiv;
}
for (uint i = _V2D::cols_ - 1; i > 0; --i) {
for (uint j = 0; j < i; ++j) {
mul = smat.element(i, j);
smat.element(i, j) -= mul;
for (uint k = 0; k < _V2D::cols_; ++k) mtmp.element(k, j) -= mul * mtmp.element(k, i);
if (sv != 0) (*sv)[j] -= mul * (*sv)[i];
}
}
if (ok != 0) *ok = true;
PIVector2D<Type>::swap(mtmp);
return *this;
}
_CMatrix inverted(bool * ok = 0) const {_CMatrix tm(*this); tm.invert(ok); return tm;}
_CMatrix transposed() const {_CMatrix tm(_V2D::rows_, _V2D::cols_); PIMM_FOR(c, r) tm.element(c, r) = _V2D::element(r, c); return tm;}
};
#ifdef PIP_STD_IOSTREAM
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.element(r, c); if (c < m.cols() - 1 || r < m.rows() - 1) s << ", ";} if (r < m.rows() - 1) s << std::endl << " ";} s << "}"; return s;}
#endif
template<typename Type>
inline PICout operator <<(PICout s, const PIMathMatrix<Type> & m) {s << "Matrix{"; for (uint r = 0; r < m.rows(); ++r) { for (uint c = 0; c < m.cols(); ++c) { s << m.element(r, c); if (c < m.cols() - 1 || r < m.rows() - 1) s << ", ";} if (r < m.rows() - 1) s << PICoutManipulators::NewLine << " ";} s << "}"; return s;}
template<typename Type>
inline PIByteArray & operator <<(PIByteArray & s, const PIMathMatrix<Type> & v) {s << (const PIVector2D<Type> &)v; return s;}
template<typename Type>
inline PIByteArray & operator >>(PIByteArray & s, PIMathMatrix<Type> & v) {s >> (PIVector2D<Type> &)v; 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.element(j, k) * sm.element(k, i);
tm.element(j, i) = 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 j = 0; j < r; ++j) {
t = Type(0);
for (uint i = 0; i < c; ++i)
t += fm.element(j, i) * sv[i];
tv[j] = t;
}
return tv;
}
/// Multiply vector {Rows} on matrix {Rows x Cols}, result is vector {Cols}
template<typename Type>
inline PIMathVector<Type> operator *(const PIMathVector<Type> & sv,
const PIMathMatrix<Type> & fm) {
uint c = fm.cols(), r = fm.rows();
PIMathVector<Type> tv(c);
Type t;
for (uint j = 0; j < c; ++j) {
t = Type(0);
for (uint i = 0; i < r; ++i)
t += fm.element(i, j) * sv[i];
tv[j] = t;
}
return tv;
}
/// Multiply value(T) on matrix {Rows x Cols}, result is vector {Rows}
template<typename Type>
inline PIMathMatrix<Type> operator *(const Type & x, const PIMathMatrix<Type> & v) {
return v * x;
}
typedef PIMathMatrix<int> PIMathMatrixi;
typedef PIMathMatrix<double> PIMathMatrixd;
template<typename T>
PIMathMatrix<complex<T> > hermitian(const PIMathMatrix<complex<T> > & m) {
PIMathMatrix<complex<T> > ret(m);
for (uint r = 0; r < ret.rows(); ++r) for (uint c = 0; c < ret.cols(); ++c) ret.element(r, c).imag(-(ret.element(r, c).imag()));
return ret.transposed();
}
#undef PIMM_FOR
#undef PIMM_FOR_I
#undef PIMM_FOR_A
#undef PIMM_FOR_C
#undef PIMM_FOR_R
#endif // PIMATHMATRIX_H

View File

@@ -0,0 +1,30 @@
/*
PIP - Platform Independent Primitives
Module includes
Ivan Pelipenko peri4ko@yandex.ru, 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 Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PIMATHMODULE_H
#define PIMATHMODULE_H
#include "pimathsolver.h"
#include "pistatistic.h"
#include "pievaluator.h"
#include "pifft.h"
#include "picrc.h"
#include "piquaternion.h"
#endif // PIMATHMODULE_H

View File

@@ -0,0 +1,242 @@
/*
PIP - Platform Independent Primitives
PIMathSolver
Ivan Pelipenko peri4ko@yandex.ru
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "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];
for (uint i = 1; i < size + 1; ++i) {
sum = 0.;
for (uint m = 0; m < i; ++m)
sum += a1[i - m] * d[m];
d[i] = b1[i] - sum;
}
for (uint i = 0; i < size - 1; ++i)
for (uint j = 0; j < size; ++j)
A[j][i] = (j == i + 1);
for (uint i = 0; i < size; ++i)
A[i][size - 1] = -a1[size - i];
for (uint i = 0; i < size; ++i)
d[i] = d[i + 1];
}
void PIMathSolver::solveEyler1(double u, double h) {
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();
}

View File

@@ -0,0 +1,89 @@
/*! \file pimathsolver.h
* \brief PIMathSolver
*/
/*
PIP - Platform Independent Primitives
PIMathSolver
Ivan Pelipenko peri4ko@yandex.ru
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PIMATHSOLVER_H
#define PIMATHSOLVER_H
#include "pimathmatrix.h"
/// Differential evaluations
struct TransferFunction {
PIVector<double> vector_Bm, vector_An;
};
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

View File

@@ -0,0 +1,252 @@
/*! \file pimathvector.h
* \brief PIMathVector
*/
/*
PIP - Platform Independent Primitives
PIMathVector
Ivan Pelipenko peri4ko@yandex.ru, 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 Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#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)
template<uint Size, typename Type = double>
class PIP_EXPORT PIMathVectorT {
typedef PIMathVectorT<Size, Type> _CVector;
public:
PIMathVectorT() {resize();}
PIMathVectorT(const PIVector<Type> & val) {resize(); PIMV_FOR(i, 0) c[i] = val[i];}
PIMathVectorT(const _CVector & st, const _CVector & fn) {resize(); set(st, fn);}
uint size() const {return Size;}
_CVector & fill(const Type & v) {PIMV_FOR(i, 0) c[i] = v; 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;}
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)));}
Type angleElevation(const _CVector & v) const {_CVector z = v - *this; double c = z.angleCos(*this); return 90.0 - acos(c) * rad2deg;}
_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;}
_CVector cross(const _CVector & v) {return (*this) * v;}
Type dot(const _CVector & v) const {return (*this) ^ v;}
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;}
_CVector & operator =(const Type & v) {PIMV_FOR(i, 0) c[i] = v; 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;}
_CVector operator &(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] *= v[i]; 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;
}
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();
return f;}
template<uint Size1, typename Type1> /// vector {Size, Type} to vector {Size1, Type1}
PIMathVectorT<Size1, Type1> turnTo() const {PIMathVectorT<Size1, Type1> tv; uint sz = piMin<uint>(Size, Size1); for (uint i = 0; i < sz; ++i) tv[i] = c[i]; return tv;}
static _CVector filled(const Type & v) {_CVector vv; PIMV_FOR(i, 0) vv[i] = v; return vv;}
private:
void resize(const Type & new_value = Type()) {for (uint i = 0; i < Size; ++i) c[i] = new_value;}
Type c[Size];
};
template<uint Size, typename Type>
inline PIMathVectorT<Size, Type> operator *(const Type & x, const PIMathVectorT<Size, Type> & v) {
return v * x;
}
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<typename T>
inline PIMathVectorT<2u, T> createVectorT2(T x, T y) {return PIMathVectorT<2u, T>(PIVector<T>() << x << y);}
template<typename T>
inline PIMathVectorT<3u, T> createVectorT3(T x, T y, T z) {return PIMathVectorT<3u, T>(PIVector<T>() << x << y << z);}
template<typename T>
inline PIMathVectorT<4u, T> createVectorT4(T x, T y, T z, T w) {return PIMathVectorT<4u, T>(PIVector<T>() << x << y << z << w);}
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;
#define createVectorT2i createVectorT2<int>
#define createVectorT3i createVectorT3<int>
#define createVectorT4i createVectorT4<int>
#define createVectorT2f createVectorT2<float>
#define createVectorT3f createVectorT3<float>
#define createVectorT4f createVectorT4<float>
#define createVectorT2d createVectorT2<double>
#define createVectorT3d createVectorT3<double>
#define createVectorT4d createVectorT4<double>
#undef PIMV_FOR
/// Vector
#define PIMV_FOR(v, s) for (uint v = s; v < c.size(); ++v)
template<typename Type>
class PIP_EXPORT PIMathVector {
typedef PIMathVector<Type> _CVector;
template<typename TypeOp> friend PIByteArray & operator <<(PIByteArray & s, const PIMathVector<TypeOp> & v);
template<typename TypeOp> friend PIByteArray & operator >>(PIByteArray & s, PIMathVector<TypeOp> & v);
public:
PIMathVector(const uint size = 0) {c.resize(size);}
PIMathVector(const PIVector<Type> & val) {c.resize(val.size()); PIMV_FOR(i, 0) c[i] = val[i];}
PIMathVector(const _CVector & st, const _CVector & fn) {c.resize(st.size()); PIMV_FOR(i, 0) c[i] = fn[i] - st[i];}
uint size() const {return c.size();}
_CVector & resize(uint size, const Type & new_value = Type()) {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 & 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 & 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 isValid() const {return !c.isEmpty();}
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) {c = v.c; return *this;}
_CVector & operator =(const Type & v) {PIMV_FOR(i, 0) c[i] = v; 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 {if (c.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;}
_CVector operator &(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] *= v[i]; return tv;}
Type operator ^(const _CVector & v) const {Type tv(0); PIMV_FOR(i, 0) tv += c[i] * v[i]; return tv;}
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();
return f;
}
template<typename Type1>
PIMathVector turnTo(uint size) const {PIMathVector<Type1> tv; uint sz = piMin<uint>(c.size(), size); for (uint i = 0; i < sz; ++i) tv[i] = c[i]; return tv;}
PIVector<Type> toVector() const {return c;}
inline Type * data() {return c.data();}
inline const Type * data() const {return c.data();}
private:
PIVector<Type> c;
};
#undef PIMV_FOR
#ifdef PIP_STD_IOSTREAM
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;}
#endif
template<typename Type>
inline PICout operator <<(PICout s, const PIMathVector<Type> & v) {s << "Vector{"; for (uint i = 0; i < v.size(); ++i) {s << v[i]; if (i < v.size() - 1) s << ", ";} s << "}"; return s;}
template<typename Type>
inline PIByteArray & operator <<(PIByteArray & s, const PIMathVector<Type> & v) {s << v.c; return s;}
template<typename Type>
inline PIByteArray & operator >>(PIByteArray & s, PIMathVector<Type> & v) {s >> v.c; return s;}
typedef PIMathVector<int> PIMathVectori;
typedef PIMathVector<double> PIMathVectord;
#endif // PIMATHVECTOR_H

View File

@@ -0,0 +1,221 @@
/*
PIP - Platform Independent Primitives
Class for quaternions
Ivan Pelipenko peri4ko@yandex.ru, 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 Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "piquaternion.h"
PIQuaternion::PIQuaternion(const PIMathVectorT3d & u, double a) {
q[0] = a;
q[1] = u[0];
q[2] = u[1];
q[3] = u[2];
}
PIMathVectorT3d PIQuaternion::eyler() const {
PIMathMatrixT33d rmat = rotationMatrix();
double angle_y = asin(rmat[0][2]), angle_x, angle_z;
double c = cos(angle_y);
if (fabs(c) < 1.E-5) {
angle_x = 0;
angle_z = atan2(rmat[1][0], rmat[1][1]);
} else {
angle_x = -atan2(-rmat[1][2] / c, rmat[2][2] / c);
angle_z = atan2(-rmat[0][1] / c, rmat[0][0] / c);
}
if (angle_z < 0) angle_z = 2*M_PI + angle_z;
return createVectorT3d(angle_x,angle_y,angle_z);
}
PIMathMatrixT33d PIQuaternion::rotationMatrix() const {
PIMathMatrixT33d ret;
double xx = q[1] * q[1];
double xy = q[1] * q[2];
double xz = q[1] * q[3];
double xw = q[1] * q[0];
double yy = q[2] * q[2];
double yz = q[2] * q[3];
double yw = q[2] * q[0];
double zz = q[3] * q[3];
double zw = q[3] * q[0];
ret[0][0] = 1. - 2. * (yy + zz);
ret[0][1] = 2. * (xy - zw);
ret[0][2] = 2. * (xz + yw);
ret[1][0] = 2. * (xy + zw);
ret[1][1] = 1. - 2. * (xx + zz);
ret[1][2] = 2. * (yz - xw);
ret[2][0] = 2. * (xz - yw);
ret[2][1] = 2. * (yz + xw);
ret[2][2] = 1. - 2. * (xx + yy);
return ret;
}
void PIQuaternion::axis(PIMathVectorT3d * ret) const {
if (!ret) return;
ret[0] = vector();
}
PIQuaternion PIQuaternion::rotated(const PIMathVectorT3d & u, double a) const {
PIQuaternion ret(*this);
ret.rotate(u, a);
return ret;
}
void PIQuaternion::rotate(const PIMathVectorT3d & u, double a) {
PIQuaternion v(u * sin(a / 2.), cos(a / 2.));
*this = (v * (*this) * v.conjugate());
}
void PIQuaternion::normalize() {
double d = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
if (d != 0.)
for (int i = 0; i < 4; ++i)
q[i] /= d;
}
PIMathMatrixT44d PIQuaternion::makeMatrix() const {
PIMathMatrixT44d ret;
ret[0][0] = q[0]; ret[0][1] = -q[1]; ret[0][2] = -q[2]; ret[0][3] = -q[3];
ret[1][0] = q[1]; ret[1][1] = q[0]; ret[1][2] = -q[3]; ret[1][3] = q[2];
ret[2][0] = q[2]; ret[2][1] = q[3]; ret[2][2] = q[0]; ret[2][3] = -q[1];
ret[3][0] = q[3]; ret[3][1] = -q[2]; ret[3][2] = q[1]; ret[3][3] = q[0];
return ret;
}
PIQuaternion PIQuaternion::fromEyler(double ax, double ay, double az) {
PIQuaternion q_heading;
PIQuaternion q_pinch;
PIQuaternion q_bank;
PIQuaternion q_tmp;
q_heading.q[0] = 0;
q_heading.q[1] = 0;
q_heading.q[2] = sin(ax / 2);
q_heading.q[3] = -cos(ax / 2);
q_pinch.q[0] = 0;
q_pinch.q[1] = sin(ay / 2);
q_pinch.q[2] = 0;
q_pinch.q[3] = -cos(ay / 2);
az = M_PI - az;
q_bank.q[0] = sin(az / 2);
q_bank.q[1] = 0;
q_bank.q[2] = 0;
q_bank.q[3] = cos(az / 2);
q_tmp = q_heading * q_pinch;
q_tmp = q_tmp * q_bank;
q_tmp.normalize();
return q_tmp;
}
PIQuaternion PIQuaternion::fromAngles(double ax, double ay, double az) {
double c1 = cos(ay/2);
double s1 = sin(ay/2);
double c2 = cos(az/2);
double s2 = sin(az/2);
double c3 = cos(ax/2);
double s3 = sin(ax/2);
double c1c2 = c1*c2;
double s1s2 = s1*s2;
double a = c1c2*c3 - s1s2*s3;
double x = c1c2*s3 + s1s2*c3;
double y = s1*c2*c3 + c1*s2*s3;
double z = c1*s2*c3 - s1*c2*s3;
PIQuaternion res;
res.q[0] = a;
res.q[1] = x;
res.q[2] = y;
res.q[3] = z;
return res;
}
PIQuaternion PIQuaternion::fromAngles2(double ax, double ay, double az) {
double om = createVectorT3d(ax,ay,az).length();
if (om == 0.) return PIQuaternion(PIMathVectorT3d(), 1.);
PIQuaternion res;
res.q[0] = cos(om/2);
res.q[1] = ax/om * sin(om/2);
res.q[2] = ay/om * sin(om/2);
res.q[3] = az/om * sin(om/2);
return res;
}
PIQuaternion PIQuaternion::fromRotationMatrix(const PIMathMatrixT33d & m) {
PIQuaternion ret;
float tr = m[0][0] + m[1][1] + m[2][2];
if (tr > 0.0f) {
ret.q[0] = tr + 1.0f;
ret.q[1] = m[2][1] - m[1][2];
ret.q[2] = m[0][2] - m[2][0];
ret.q[3] = m[1][0] - m[0][1];
} else {
if ((m[0][0] > m[1][1]) && (m[0][0] > m[2][2])) {
ret.q[0] = m[2][1] - m[1][2];
ret.q[1] = 1.0f + m[0][0] - m[1][1] - m[2][2];
ret.q[2] = m[0][1] + m[1][0];
ret.q[3] = m[0][2] + m[2][0];
} else {
if (m[1][1] > m[2][2]) {
ret.q[0] = m[0][2] - m[2][0];
ret.q[1] = m[0][1] + m[1][0];
ret.q[2] = 1.0f + m[1][1] - m[0][0] - m[2][2];
ret.q[3] = m[1][2] + m[2][1];
} else {
ret.q[0] = m[1][0] - m[0][1];
ret.q[1] = m[0][2] + m[2][0];
ret.q[2] = m[1][2] + m[2][1];
ret.q[3] = 1.0f + m[2][2] - m[0][0] - m[1][1];
}
}
}
ret.normalize();
return ret;
}
PIQuaternion operator*(const PIQuaternion & q0, const PIQuaternion & q1) {
PIMathVectorT3d v0(q0.vector()), v1(q1.vector());
double r0 = q0.q[0] * q1.q[0] - (v0^v1);
PIMathVectorT3d qv = v1*q0.q[0] + v0*q1.q[0] + v0*v1;
PIQuaternion ret;
ret.q[0] = r0;
ret.q[1] = qv[0];
ret.q[2] = qv[1];
ret.q[3] = qv[2];
return ret;
}
PIQuaternion operator *(const double &a, const PIQuaternion &q1) {
PIQuaternion ret(q1);
ret.q[0] *= a;
ret.q[1] *= a;
ret.q[2] *= a;
ret.q[3] *= a;
return ret;
}

View File

@@ -0,0 +1,66 @@
/*! \file piquaternion.h
* \brief Class for quaternions
*/
/*
PIP - Platform Independent Primitives
Class for quaternions
Ivan Pelipenko peri4ko@yandex.ru, 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 Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PIQUATERNION_H
#define PIQUATERNION_H
#include "pimathmatrix.h"
class PIQuaternion
{
friend PIQuaternion operator*(const PIQuaternion & q0, const PIQuaternion & q1);
friend PIQuaternion operator*(const double & a, const PIQuaternion & q);
public:
PIQuaternion(const PIMathVectorT3d & u = PIMathVectorT3d(), double a = 0.);
PIQuaternion conjugate() const {return PIQuaternion(-vector(), scalar());}
PIQuaternion rotated(const PIMathVectorT3d & u, double a) const;
void rotate(const PIMathVectorT3d & u, double a);
void normalize();
double & scalar() {return q[0];}
double scalar() const {return q[0];}
PIMathVectorT3d vector() const {return createVectorT3<double>(q[1], q[2], q[3]);}
PIMathVectorT3d eyler() const;
PIMathMatrixT33d rotationMatrix() const;
void axis(PIMathVectorT3d*ret) const;
static PIQuaternion fromEyler(double ax, double ay, double az);
static PIQuaternion fromRotationMatrix(const PIMathMatrixT33d & m);
static PIQuaternion fromAngles(double ax, double ay, double az);
static PIQuaternion fromAngles2(double ax, double ay, double az);
protected:
double q[4];
PIMathMatrixT44d makeMatrix() const;
};
PIQuaternion operator *(const double & a, const PIQuaternion & q);
PIQuaternion operator *(const PIQuaternion & q0, const PIQuaternion & q1);
inline PIQuaternion operator +(const PIQuaternion & q0, const PIQuaternion & q1) {return PIQuaternion(q0.vector() + q1.vector(), q0.scalar() + q1.scalar());}
inline PIQuaternion operator -(const PIQuaternion & q0, const PIQuaternion & q1) {return PIQuaternion(q0.vector() - q1.vector(), q0.scalar() - q1.scalar());}
inline PIQuaternion operator -(const PIQuaternion & q0) {return PIQuaternion(-q0.vector(), -q0.scalar());}
#endif // PIQUATERNION_H

View File

@@ -0,0 +1,85 @@
/*! \file pistatistic.h
* \brief Class for calculating math statistic in values array
*/
/*
PIP - Platform Independent Primitives
Class for calculacing math statistic in values array
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 Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#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