723 lines
32 KiB
C++
723 lines
32 KiB
C++
#ifndef PIMATH_H
|
||
#define PIMATH_H
|
||
|
||
#include "picontainers.h"
|
||
#ifndef QNX
|
||
# include <cmath>
|
||
# include <complex>
|
||
#else
|
||
# include <math.h>
|
||
# include <complex.h>
|
||
#endif
|
||
|
||
#define M_2PI 6.28318530717958647692
|
||
#define M_PI_3 1.04719755119659774615
|
||
|
||
using std::complex;
|
||
|
||
typedef complex<int> complexi;
|
||
typedef complex<float> complexf;
|
||
typedef complex<double> complexd;
|
||
typedef complex<long double> complexld;
|
||
const complexld complexld_i(0., 1.);
|
||
const complexld complexld_0(0.);
|
||
const complexld complexld_1(1.);
|
||
const complexd complexd_i(0., 1.);
|
||
const complexd complexd_0(0.);
|
||
const complexd complexd_1(1.);
|
||
|
||
const double deg2rad = atan(1.) / 45.;
|
||
const double rad2deg = 45. / atan(1.);
|
||
|
||
inline int pow2(int p) {return 1 << p;}
|
||
inline double sqr(const double & v) {return v * v;}
|
||
inline double sinc(const double & v) {double t = M_PI * v; return sin(t) / t;}
|
||
inline complexd atanc(const complexd & c) {return -complexd(-0.5, 1.) * log((complexd_1 + complexd_i * c) / (complexd_1 - complexd_i * c));}
|
||
inline complexd asinc(const complexd & c) {return -complexd_i * log(complexd_i * c + sqrt(complexd_1 - c * c));}
|
||
inline complexd acosc(const complexd & c) {return -complexd_i * log(c + complexd_i * sqrt(complexd_1 - c * c));}
|
||
#ifdef QNX
|
||
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;}
|
||
inline double j0(const double & v) {v;}
|
||
inline double j1(const double & v) {v;}
|
||
inline double jn(const int & n, const double & v) {v;}
|
||
inline double y0(const double & v) {v;}
|
||
inline double y1(const double & v) {v;}
|
||
inline double yn(const int & n, const double & v) {v;}
|
||
#endif
|
||
|
||
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 PIMathVectorT {
|
||
typedef PIMathVectorT<Size, Type> _CVector;
|
||
public:
|
||
PIMathVectorT() {resize(Size);}
|
||
//PIMathVectorT(Type val) {resize(Size); PIMV_FOR(i, 0) c[i] = val;}
|
||
PIMathVectorT(Type fval, ...) {resize(Size); c[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] = va_arg(vl, Type); va_end(vl);}
|
||
PIMathVectorT(const PIVector<Type> & val) {resize(Size); PIMV_FOR(i, 0) c[i] = val[i];}
|
||
PIMathVectorT(const _CVector & st, const _CVector & fn) {resize(Size); set(st, fn);}
|
||
|
||
inline uint size() const {return Size;}
|
||
inline _CVector & fill(const Type & v) {PIMV_FOR(i, 0) c[i] = v; return *this;}
|
||
inline _CVector & set(Type fval, ...) {c[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] = va_arg(vl, Type); va_end(vl); return *this;}
|
||
inline _CVector & set(const _CVector & st, const _CVector & fn) {PIMV_FOR(i, 0) c[i] = fn[i] - st[i]; return *this;}
|
||
inline _CVector & move(const Type & v) {PIMV_FOR(i, 0) c[i] += v; return *this;}
|
||
inline _CVector & move(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i]; return *this;}
|
||
inline _CVector & move(Type fval, ...) {c[0] += fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] += va_arg(vl, Type); va_end(vl); return *this;}
|
||
inline Type lengthSqr() const {Type tv(0); PIMV_FOR(i, 0) tv += (c[i] * c[i]); return tv;}
|
||
inline Type length() const {return sqrt(lengthSqr());}
|
||
inline Type manhattanLength() const {Type tv(0); PIMV_FOR(i, 0) tv += fabs(c[i]); return tv;}
|
||
inline Type angleCos(const _CVector & v) const {Type tv = v.length() * length(); return (tv == Type(0) ? Type(0) : ((*this) ^ v) / tv);}
|
||
inline Type angleSin(const _CVector & v) const {Type tv = angleCos(v); return sqrt(Type(1) - tv * tv);}
|
||
inline Type angleRad(const _CVector & v) const {return acos(angleCos(v));}
|
||
inline Type angleDeg(const _CVector & v) const {return acos(angleCos(v)) * rad2deg;}
|
||
inline _CVector projection(const _CVector & v) {Type tv = v.length(); return (tv == Type(0) ? _CVector() : v * (((*this) ^ v) / tv));}
|
||
inline _CVector & normalize() {Type tv = length(); if (tv == Type(1)) return *this; PIMV_FOR(i, 0) c[i] /= tv; return *this;}
|
||
inline _CVector normalized() {_CVector tv(*this); tv.normalize(); return tv;}
|
||
inline bool isNull() const {PIMV_FOR(i, 0) if (c[i] != Type(0)) return false; return true;}
|
||
inline bool isOrtho(const _CVector & v) const {return ((*this) ^ v) == Type(0);}
|
||
|
||
inline Type & at(uint index) {return c[index];}
|
||
inline Type at(uint index) const {return c[index];}
|
||
inline Type & operator [](uint index) {return c[index];}
|
||
inline Type operator [](uint index) const {return c[index];}
|
||
inline void operator =(const _CVector & v) {c = v.c;}
|
||
inline bool operator ==(const _CVector & v) const {PIMV_FOR(i, 0) if (c[i] != v[i]) return false; return true;}
|
||
inline bool operator !=(const _CVector & v) const {return !(*this == c);}
|
||
inline void operator +=(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i];}
|
||
inline void operator -=(const _CVector & v) {PIMV_FOR(i, 0) c[i] -= v[i];}
|
||
inline void operator *=(const Type & v) {PIMV_FOR(i, 0) c[i] *= v;}
|
||
inline void operator *=(const _CVector & v) {PIMV_FOR(i, 0) c[i] *= v[i];}
|
||
inline void operator /=(const Type & v) {PIMV_FOR(i, 0) c[i] /= v;}
|
||
inline void operator /=(const _CVector & v) {PIMV_FOR(i, 0) c[i] /= v[i];}
|
||
inline _CVector operator -() const {_CVector tv; PIMV_FOR(i, 0) tv[i] = -c[i]; return tv;}
|
||
inline _CVector operator +(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] += v[i]; return tv;}
|
||
inline _CVector operator -(const _CVector & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] -= v[i]; return tv;}
|
||
inline _CVector operator *(const Type & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] *= v; return tv;}
|
||
inline _CVector operator /(const Type & v) const {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] /= v; return tv;}
|
||
inline _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;}
|
||
inline Type operator ^(const _CVector & v) const {Type tv(0); PIMV_FOR(i, 0) tv += c[i] * v[i]; return tv;}
|
||
|
||
inline operator PIMathMatrixT<1, Size, Type>() {return PIMathMatrixT<1, Size, Type>(c);}
|
||
inline Type distToLine(const _CVector & lp0, const _CVector & lp1) {
|
||
_CVector a(lp0, lp1), b(lp0, *this), c(lp1, *this);
|
||
Type f = fabs(a[0]*b[1] - a[1]*b[0]) / a.length();//, s = b.length() + c.length() - a.length();
|
||
return f;}
|
||
|
||
template<uint Size1, typename Type1> /// vector {Size, Type} to vector {Size1, Type1}
|
||
inline PIMathVectorT<Size1, Type1> turnTo() {PIMathVectorT<Size1, Type1> tv; uint sz = piMin<uint>(Size, Size1); for (uint i = 0; i < sz; ++i) tv[i] = c[i]; return tv;}
|
||
|
||
private:
|
||
inline void resize(uint size, const Type & new_value = Type()) {c.resize(size, new_value);}
|
||
PIVector<Type> c;
|
||
|
||
};
|
||
|
||
template<uint Size, typename Type>
|
||
inline std::ostream & operator <<(std::ostream & s, const PIMathVectorT<Size, Type> & v) {s << '{'; PIMV_FOR(i, 0) {s << v[i]; if (i < Size - 1) s << ", ";} s << '}'; return s;}
|
||
template<uint Size, typename Type>
|
||
inline const bool operator ||(const PIMathVectorT<Size, Type> & f, const PIMathVectorT<Size, Type> & s) {return (f * s).isNull();}
|
||
|
||
//template<uint Size0, typename Type0 = double, uint Size1 = Size0, typename Type1 = Type0> /// vector {Size0, Type0} to vector {Size1, Type1}
|
||
//inline operator PIMathVectorT<Size1, Type1>(const PIMathVectorT<Size0, Type0> & v) {PIMathVectorT<Size1, Type1> tv; uint sz = piMin<uint>(Size0, Size1); for (uint i = 0; i < sz; ++i) tv[i] = v[i]; return tv;}
|
||
|
||
typedef PIMathVectorT<2u, int> PIMathVectorT2i;
|
||
typedef PIMathVectorT<3u, int> PIMathVectorT3i;
|
||
typedef PIMathVectorT<4u, int> PIMathVectorT4i;
|
||
typedef PIMathVectorT<2u, double> PIMathVectorT2d;
|
||
typedef PIMathVectorT<3u, double> PIMathVectorT3d;
|
||
typedef PIMathVectorT<4u, double> PIMathVectorT4d;
|
||
|
||
/// Matrix templated
|
||
|
||
#define PIMM_FOR(c, r) for (uint c = 0; c < Cols; ++c) { for (uint r = 0; r < Rows; ++r) {
|
||
#define PIMM_FOR_WB(c, r) for (uint c = 0; c < Cols; ++c) for (uint r = 0; r < Rows; ++r) // without brakes
|
||
#define PIMM_FOR_I(c, r) for (uint r = 0; r < Rows; ++r) { for (uint c = 0; c < Cols; ++c) {
|
||
#define PIMM_FOR_I_WB(c, r) for (uint r = 0; r < Rows; ++r) for (uint c = 0; c < Cols; ++c) // without brakes
|
||
#define PIMM_FOR_C(v) for (uint v = 0; v < Cols; ++v)
|
||
#define PIMM_FOR_R(v) for (uint v = 0; v < Rows; ++v)
|
||
|
||
template<uint Cols, uint Rows = Cols, typename Type = double>
|
||
class PIMathMatrixT {
|
||
typedef PIMathMatrixT<Cols, Rows, Type> _CMatrix;
|
||
typedef PIMathMatrixT<Rows, Cols, Type> _CMatrixI;
|
||
typedef PIMathVectorT<Rows, Type> _CMCol;
|
||
typedef PIMathVectorT<Cols, Type> _CMRow;
|
||
public:
|
||
PIMathMatrixT() {resize(Cols, Rows);}
|
||
PIMathMatrixT(Type fval, ...) {resize(Cols, Rows); va_list vl; va_start(vl, fval); PIMM_FOR_I_WB(c, r) m[c][r] = (r + c == 0 ? fval : va_arg(vl, Type)); va_end(vl);}
|
||
PIMathMatrixT(const PIVector<Type> & val) {resize(Cols, Rows); int i = 0; PIMM_FOR_I_WB(c, r) m[c][r] = val[i++];}
|
||
|
||
static _CMatrix identity() {_CMatrix tm = _CMatrix(); PIMM_FOR_WB(c, r) tm.m[c][r] = (c == r ? Type(1) : Type(0)); return tm;}
|
||
|
||
inline uint cols() const {return Cols;}
|
||
inline uint rows() const {return Rows;}
|
||
inline _CMCol col(uint index) {_CMCol tv; PIMM_FOR_R(i) tv[i] = m[index][i]; return tv;}
|
||
inline _CMRow row(uint index) {_CMRow tv; PIMM_FOR_C(i) tv[i] = m[i][index]; return tv;}
|
||
inline _CMatrix & setCol(uint index, const _CMCol & v) {PIMM_FOR_R(i) m[index][i] = v[i]; return *this;}
|
||
inline _CMatrix & setRow(uint index, const _CMRow & v) {PIMM_FOR_C(i) m[i][index] = v[i]; return *this;}
|
||
inline _CMatrix & swapRows(uint r0, uint r1) {Type t; PIMM_FOR_C(i) {t = m[i][r0]; m[i][r0] = m[i][r1]; m[i][r1] = t;} return *this;}
|
||
inline _CMatrix & swapCols(uint c0, uint c1) {Type t; PIMM_FOR_R(i) {t = m[c0][i]; m[c0][i] = m[c1][i]; m[c1][i] = t;} return *this;}
|
||
inline _CMatrix & fill(const Type & v) {PIMM_FOR_WB(c, r) m[c][r] = v; return *this;}
|
||
//inline _CMatrix & set(Type fval, ...) {m[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) m[i] = va_arg(vl, Type); va_end(vl); return *this;}
|
||
//inline void normalize() {Type tv = length(); if (tv == Type(1)) return; PIMV_FOR(i, 0) m[i] /= tv;}
|
||
inline bool isSquare() const {return cols() == rows();}
|
||
inline bool isIdentity() const {PIMM_FOR_WB(c, r) if ((c == r) ? m[c][r] != Type(1) : m[c][r] != Type(0)) return false; return true;}
|
||
inline bool isNull() const {PIMM_FOR_WB(c, r) if (m[c][r] != Type(0)) return false; return true;}
|
||
|
||
inline Type & at(uint col, uint row) {return m[col][row];}
|
||
inline Type at(uint col, uint row) const {return m[col][row];}
|
||
inline PIVector<Type> & operator [](uint col) {return m[col];}
|
||
inline PIVector<Type> operator [](uint col) const {return m[col];}
|
||
inline void operator =(const _CMatrix & sm) {m = sm.m;}
|
||
inline bool operator ==(const _CMatrix & sm) const {PIMM_FOR_WB(c, r) if (m[c][r] != sm.m[c][r]) return false; return true;}
|
||
inline bool operator !=(const _CMatrix & sm) const {return !(*this == sm);}
|
||
inline void operator +=(const _CMatrix & sm) {PIMM_FOR_WB(c, r) m[c][r] += sm.m[c][r];}
|
||
inline void operator -=(const _CMatrix & sm) {PIMM_FOR_WB(c, r) m[c][r] -= sm.m[c][r];}
|
||
inline void operator *=(const Type & v) {PIMM_FOR_WB(c, r) m[c][r] *= v;}
|
||
inline void operator /=(const Type & v) {PIMM_FOR_WB(c, r) m[c][r] /= v;}
|
||
inline _CMatrix operator -() {_CMatrix tm; PIMM_FOR_WB(c, r) tm.m[c][r] = -m[c][r]; return tm;}
|
||
inline _CMatrix operator +(const _CMatrix & sm) {_CMatrix tm = _CMatrix(*this); PIMM_FOR_WB(c, r) tm.m[c][r] += sm.m[c][r]; return tm;}
|
||
inline _CMatrix operator -(const _CMatrix & sm) {_CMatrix tm = _CMatrix(*this); PIMM_FOR_WB(c, r) tm.m[c][r] -= sm.m[c][r]; return tm;}
|
||
inline _CMatrix operator *(const Type & v) {_CMatrix tm = _CMatrix(*this); PIMM_FOR_WB(c, r) tm.m[c][r] *= v; return tm;}
|
||
inline _CMatrix operator /(const Type & v) {_CMatrix tm = _CMatrix(*this); PIMM_FOR_WB(c, r) tm.m[c][r] /= v; return tm;}
|
||
|
||
_CMatrix & toUpperTriangular(bool * ok = 0) {
|
||
if (Cols != Rows) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
_CMatrix smat(*this);
|
||
bool ndet;
|
||
uint crow;
|
||
Type mul;
|
||
for (uint i = 0; i < Cols; ++i) {
|
||
ndet = true;
|
||
for (uint j = 0; j < Rows; ++j) if (smat.m[i][j] != 0) ndet = false;
|
||
if (ndet) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
for (uint j = 0; j < Cols; ++j) if (smat.m[j][i] != 0) ndet = false;
|
||
if (ndet) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
}
|
||
for (uint i = 0; i < Cols; ++i) {
|
||
crow = i;
|
||
while (smat.m[i][i] == Type(0))
|
||
smat.swapRows(i, ++crow);
|
||
for (uint j = i + 1; j < Rows; ++j) {
|
||
mul = smat.m[i][j] / smat.m[i][i];
|
||
for (uint k = i; k < Cols; ++k) smat.m[k][j] -= mul * smat.m[k][i];
|
||
}
|
||
if (i < Cols - 1) {
|
||
if (fabs(smat.m[i+1][i+1]) < Type(1E-100)) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
}
|
||
}
|
||
if (ok != 0) *ok = true;
|
||
m = smat.m;
|
||
return *this;
|
||
}
|
||
|
||
_CMatrix & invert(bool * ok = 0) {
|
||
if (Cols != Rows) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
_CMatrix mtmp = _CMatrix::identity(), smat(*this);
|
||
bool ndet;
|
||
uint crow;
|
||
Type mul, iddiv;
|
||
for (uint i = 0; i < Cols; ++i) {
|
||
ndet = true;
|
||
for (uint j = 0; j < Rows; ++j) if (smat.m[i][j] != 0) ndet = false;
|
||
if (ndet) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
for (uint j = 0; j < Cols; ++j) if (smat.m[j][i] != 0) ndet = false;
|
||
if (ndet) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
}
|
||
for (uint i = 0; i < Cols; ++i) {
|
||
crow = i;
|
||
while (smat.m[i][i] == Type(0)) {
|
||
++crow;
|
||
smat.swapRows(i, crow);
|
||
mtmp.swapRows(i, crow);
|
||
}
|
||
for (uint j = i + 1; j < Rows; ++j) {
|
||
mul = smat.m[i][j] / smat.m[i][i];
|
||
for (uint k = i; k < Cols; ++k) smat.m[k][j] -= mul * smat.m[k][i];
|
||
for (uint k = 0; k < Cols; ++k) mtmp.m[k][j] -= mul * mtmp.m[k][i];
|
||
}
|
||
//cout << i << endl << smat << endl;
|
||
if (i < Cols - 1) {
|
||
if (fabs(smat.m[i+1][i+1]) < Type(1E-100)) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
}
|
||
iddiv = smat.m[i][i];
|
||
for (uint j = i; j < Cols; ++j) smat.m[j][i] /= iddiv;
|
||
for (uint j = 0; j < Cols; ++j) mtmp.m[j][i] /= iddiv;
|
||
}
|
||
for (uint i = Cols - 1; i > 0; --i) {
|
||
for (uint j = 0; j < i; ++j) {
|
||
mul = smat.m[i][j];
|
||
smat.m[i][j] -= mul;
|
||
for (uint k = 0; k < Cols; ++k) mtmp.m[k][j] -= mtmp.m[k][i] * mul;
|
||
}
|
||
}
|
||
if (ok != 0) *ok = true;
|
||
m = mtmp.m;
|
||
return *this;
|
||
}
|
||
inline _CMatrix inverted(bool * ok = 0) {_CMatrix tm(*this); tm.invert(ok); return tm;}
|
||
inline _CMatrixI transposed() {_CMatrixI tm; PIMM_FOR_WB(c, r) tm[r][c] = m[c][r]; return tm;}
|
||
|
||
private:
|
||
inline void resize(uint cols, uint rows, const Type & new_value = Type()) {m.resize(cols); PIMM_FOR_C(i) m[i].resize(rows, new_value);}
|
||
PIVector<PIVector<Type> > m;
|
||
|
||
};
|
||
|
||
template<uint Cols, uint Rows, typename Type>
|
||
inline std::ostream & operator <<(std::ostream & s, const PIMathMatrixT<Cols, Rows, Type> & m) {s << '{'; PIMM_FOR_I(c, r) s << m[c][r]; if (c < Cols - 1 || r < Rows - 1) s << ", ";} if (r < Rows - 1) s << endl << ' ';} s << '}'; return s;}
|
||
|
||
/// Multiply matrices {CR x Rows0} on {Cols1 x CR}, result is {Cols1 x Rows0}
|
||
template<uint CR, uint Rows0, uint Cols1, typename Type>
|
||
inline PIMathMatrixT<Cols1, Rows0, Type> operator *(const PIMathMatrixT<CR, Rows0, Type> & fm,
|
||
const PIMathMatrixT<Cols1, CR, Type> & sm) {
|
||
PIMathMatrixT<Cols1, Rows0, 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[k][j] * sm[i][k];
|
||
tm[i][j] = t;
|
||
}
|
||
}
|
||
return tm;
|
||
}
|
||
|
||
/// Multiply matrix {Cols x Rows} on vector {Cols}, result is vector {Rows}
|
||
template<uint Cols, uint Rows, typename Type>
|
||
inline PIMathVectorT<Rows, Type> operator *(const PIMathMatrixT<Cols, Rows, Type> & fm,
|
||
const PIMathVectorT<Cols, Type> & sv) {
|
||
PIMathVectorT<Rows, Type> tv;
|
||
Type t;
|
||
for (uint i = 0; i < Rows; ++i) {
|
||
t = Type(0);
|
||
for (uint j = 0; j < Cols; ++j)
|
||
t += fm[j][i] * sv[j];
|
||
tv[i] = t;
|
||
}
|
||
return tv;
|
||
}
|
||
|
||
typedef PIMathMatrixT<2u, 2u, int> PIMathMatrixT22i;
|
||
typedef PIMathMatrixT<3u, 3u, int> PIMathMatrixT33i;
|
||
typedef PIMathMatrixT<4u, 4u, int> PIMathMatrixT44i;
|
||
typedef PIMathMatrixT<2u, 2u, double> PIMathMatrixT22d;
|
||
typedef PIMathMatrixT<3u, 3u, double> PIMathMatrixT33d;
|
||
typedef PIMathMatrixT<4u, 4u, double> PIMathMatrixT44d;
|
||
|
||
|
||
template<typename Type>
|
||
class PIMathMatrix;
|
||
|
||
#undef PIMV_FOR
|
||
#undef PIMM_FOR
|
||
#undef PIMM_FOR_WB
|
||
#undef PIMM_FOR_I
|
||
#undef PIMM_FOR_I_WB
|
||
#undef PIMM_FOR_C
|
||
#undef PIMM_FOR_R
|
||
|
||
/// Vector
|
||
|
||
#define PIMV_FOR(v, s) for (uint v = s; v < size_; ++v)
|
||
|
||
template<typename Type>
|
||
class PIMathVector {
|
||
typedef PIMathVector<Type> _CVector;
|
||
public:
|
||
PIMathVector(const uint size = 3) {resize(size);}
|
||
PIMathVector(const uint size, Type fval, ...) {resize(size); c[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] = va_arg(vl, Type); va_end(vl);}
|
||
PIMathVector(const PIVector<Type> & val) {resize(val.size); PIMV_FOR(i, 0) c[i] = val[i];}
|
||
PIMathVector(const _CVector & st, const _CVector & fn) {resize(st.size()); PIMV_FOR(i, 0) c[i] = fn[i] - st[i];}
|
||
|
||
inline uint size() const {return size_;}
|
||
inline _CVector & resize(uint size, const Type & new_value = Type()) {size_ = size; c.resize(size, new_value); return *this;}
|
||
inline _CVector resized(uint size, const Type & new_value = Type()) {_CVector tv = _CVector(*this); tv.resize(size, new_value); return tv;}
|
||
inline _CVector & fill(const Type & v) {PIMV_FOR(i, 0) c[i] = v; return *this;}
|
||
inline _CVector & set(Type fval, ...) {c[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] = va_arg(vl, Type); va_end(vl); return *this;}
|
||
inline _CVector & move(const Type & v) {PIMV_FOR(i, 0) c[i] += v; return *this;}
|
||
inline _CVector & move(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i]; return *this;}
|
||
inline _CVector & move(Type fval, ...) {c[0] += fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) c[i] += va_arg(vl, Type); va_end(vl); return *this;}
|
||
inline _CVector & swap(uint fe, uint se) {piSwap<Type>(c[fe], c[se]); return *this;}
|
||
inline Type lengthSqr() const {Type tv(0); PIMV_FOR(i, 0) tv += (c[i] * c[i]); return tv;}
|
||
inline Type length() const {return sqrt(lengthSqr());}
|
||
inline Type manhattanLength() const {Type tv(0); PIMV_FOR(i, 0) tv += fabs(c[i]); return tv;}
|
||
inline Type angleCos(const _CVector & v) const {Type tv = v.length() * length(); return (tv == Type(0) ? Type(0) : ((*this) ^ v) / tv);}
|
||
inline Type angleSin(const _CVector & v) const {Type tv = angleCos(v); return sqrt(Type(1) - tv * tv);}
|
||
inline Type angleRad(const _CVector & v) const {return acos(angleCos(v));}
|
||
inline Type angleDeg(const _CVector & v) const {return acos(angleCos(v)) * rad2deg;}
|
||
inline _CVector projection(const _CVector & v) {Type tv = v.length(); return (tv == Type(0) ? _CVector() : v * (((*this) ^ v) / tv));}
|
||
inline _CVector & normalize() {Type tv = length(); if (tv == Type(1)) return *this; PIMV_FOR(i, 0) c[i] /= tv; return *this;}
|
||
inline _CVector normalized() {_CVector tv(*this); tv.normalize(); return tv;}
|
||
inline bool isNull() const {PIMV_FOR(i, 0) if (c[i] != Type(0)) return false; return true;}
|
||
inline bool isOrtho(const _CVector & v) const {return ((*this) ^ v) == Type(0);}
|
||
|
||
inline Type & at(uint index) {return c[index];}
|
||
inline Type at(uint index) const {return c[index];}
|
||
inline Type & operator [](uint index) {return c[index];}
|
||
inline Type operator [](uint index) const {return c[index];}
|
||
inline void operator =(const _CVector & v) {c = v.c;}
|
||
inline bool operator ==(const _CVector & v) const {PIMV_FOR(i, 0) if (c[i] != v[i]) return false; return true;}
|
||
inline bool operator !=(const _CVector & v) const {return !(*this == c);}
|
||
inline void operator +=(const _CVector & v) {PIMV_FOR(i, 0) c[i] += v[i];}
|
||
inline void operator -=(const _CVector & v) {PIMV_FOR(i, 0) c[i] -= v[i];}
|
||
inline void operator *=(const Type & v) {PIMV_FOR(i, 0) c[i] *= v;}
|
||
inline void operator *=(const _CVector & v) {PIMV_FOR(i, 0) c[i] *= v[i];}
|
||
inline void operator /=(const Type & v) {PIMV_FOR(i, 0) c[i] /= v;}
|
||
inline void operator /=(const _CVector & v) {PIMV_FOR(i, 0) c[i] /= v[i];}
|
||
inline _CVector operator -() {_CVector tv; PIMV_FOR(i, 0) tv[i] = -c[i]; return tv;}
|
||
inline _CVector operator +(const _CVector & v) {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] += v[i]; return tv;}
|
||
inline _CVector operator -(const _CVector & v) {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] -= v[i]; return tv;}
|
||
inline _CVector operator *(const Type & v) {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] *= v; return tv;}
|
||
inline _CVector operator /(const Type & v) {_CVector tv = _CVector(*this); PIMV_FOR(i, 0) tv[i] /= v; return tv;}
|
||
inline _CVector operator *(const _CVector & v) {if (size_ < 3) return _CVector(); _CVector tv; tv.fill(Type(1)); tv[0] = c[1]*v[2] - v[1]*c[2]; tv[1] = v[0]*c[2] - c[0]*v[2]; tv[2] = c[0]*v[1] - v[0]*c[1]; return tv;}
|
||
inline Type operator ^(const _CVector & v) const {Type tv(0); PIMV_FOR(i, 0) tv += c[i] * v[i]; return tv;}
|
||
|
||
//inline operator PIMathMatrix<1, Size, Type>() {return PIMathMatrix<1, Size, Type>(c);}
|
||
inline Type distToLine(const _CVector & lp0, const _CVector & lp1) {
|
||
_CVector a(lp0, lp1), b(lp0, *this), c(lp1, *this);
|
||
Type f = fabs(a[0]*b[1] - a[1]*b[0]) / a.length();//, s = b.length() + c.length() - a.length();
|
||
return f;}
|
||
|
||
template<typename Type1>
|
||
inline PIMathVector turnTo(uint size) {PIMathVector<Type1> tv; uint sz = piMin<uint>(size_, size); for (uint i = 0; i < sz; ++i) tv[i] = c[i]; return tv;}
|
||
|
||
private:
|
||
uint size_;
|
||
PIVector<Type> c;
|
||
|
||
};
|
||
|
||
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;}
|
||
|
||
typedef PIMathVector<int> PIMathVectori;
|
||
typedef PIMathVector<double> PIMathVectord;
|
||
|
||
/// Matrix
|
||
|
||
#define PIMM_FOR(c, r) for (uint c = 0; c < cols_; ++c) { for (uint r = 0; r < rows_; ++r) {
|
||
#define PIMM_FOR_WB(c, r) for (uint c = 0; c < cols_; ++c) for (uint r = 0; r < rows_; ++r) // without brakes
|
||
#define PIMM_FOR_I(c, r) for (uint r = 0; r < rows_; ++r) { for (uint c = 0; c < cols_; ++c) {
|
||
#define PIMM_FOR_I_WB(c, r) for (uint r = 0; r < rows_; ++r) for (uint c = 0; c < cols_; ++c) // without brakes
|
||
#define PIMM_FOR_C(v) for (uint v = 0; v < cols_; ++v)
|
||
#define PIMM_FOR_R(v) for (uint v = 0; v < rows_; ++v)
|
||
|
||
template<typename Type>
|
||
class PIMathMatrix {
|
||
typedef PIMathMatrix<Type> _CMatrix;
|
||
typedef PIMathVector<Type> _CMCol;
|
||
typedef PIMathVector<Type> _CMRow;
|
||
public:
|
||
PIMathMatrix(const uint cols = 3, const uint rows = 3) {resize(cols, rows);}
|
||
PIMathMatrix(const uint cols, const uint rows, Type fval, ...) {resize(cols, rows); va_list vl; va_start(vl, fval); PIMM_FOR_I_WB(c, r) m[c][r] = (r + c == 0 ? fval : va_arg(vl, Type)); va_end(vl);}
|
||
PIMathMatrix(const uint cols, const uint rows, const PIVector<Type> & val) {resize(cols, rows); int i = 0; PIMM_FOR_I_WB(c, r) m[c][r] = val[i++];}
|
||
|
||
static _CMatrix identity(const uint cols_, const uint rows_) {_CMatrix tm(cols_, rows_); PIMM_FOR_WB(c, r) tm.m[c][r] = (c == r ? Type(1) : Type(0)); return tm;}
|
||
|
||
inline uint cols() const {return cols_;}
|
||
inline uint rows() const {return rows_;}
|
||
inline _CMCol col(uint index) {_CMCol tv; PIMM_FOR_R(i) tv[i] = m[index][i]; return tv;}
|
||
inline _CMRow row(uint index) {_CMRow tv; PIMM_FOR_C(i) tv[i] = m[i][index]; return tv;}
|
||
inline _CMatrix & resize(const uint cols, const uint rows, const Type & new_value = Type()) {cols_ = cols; rows_ = rows; m.resize(cols); PIMM_FOR_C(i) m[i].resize(rows, new_value); return *this;}
|
||
inline _CMatrix & setCol(uint index, const _CMCol & v) {PIMM_FOR_R(i) m[index][i] = v[i]; return *this;}
|
||
inline _CMatrix & setRow(uint index, const _CMRow & v) {PIMM_FOR_C(i) m[i][index] = v[i]; return *this;}
|
||
inline _CMatrix & swaprows_(uint r0, uint r1) {Type t; PIMM_FOR_C(i) {t = m[i][r0]; m[i][r0] = m[i][r1]; m[i][r1] = t;} return *this;}
|
||
inline _CMatrix & swapcols_(uint c0, uint c1) {Type t; PIMM_FOR_R(i) {t = m[c0][i]; m[c0][i] = m[c1][i]; m[c1][i] = t;} return *this;}
|
||
inline _CMatrix & fill(const Type & v) {PIMM_FOR_WB(c, r) m[c][r] = v; return *this;}
|
||
//inline _CMatrix & set(Type fval, ...) {m[0] = fval; va_list vl; va_start(vl, fval); PIMV_FOR(i, 1) m[i] = va_arg(vl, Type); va_end(vl); return *this;}
|
||
//inline void normalize() {Type tv = length(); if (tv == Type(1)) return; PIMV_FOR(i, 0) m[i] /= tv;}
|
||
inline bool isSquare() const {return cols() == rows();}
|
||
inline bool isIdentity() const {PIMM_FOR_WB(c, r) if ((c == r) ? m[c][r] != Type(1) : m[c][r] != Type(0)) return false; return true;}
|
||
inline bool isNull() const {PIMM_FOR_WB(c, r) if (m[c][r] != Type(0)) return false; return true;}
|
||
|
||
inline Type & at(uint col, uint row) {return m[col][row];}
|
||
inline Type at(uint col, uint row) const {return m[col][row];}
|
||
inline PIVector<Type> & operator [](uint col) {return m[col];}
|
||
inline PIVector<Type> operator [](uint col) const {return m[col];}
|
||
inline void operator =(const _CMatrix & sm) {m = sm.m;}
|
||
inline bool operator ==(const _CMatrix & sm) const {PIMM_FOR_WB(c, r) if (m[c][r] != sm.m[c][r]) return false; return true;}
|
||
inline bool operator !=(const _CMatrix & sm) const {return !(*this == sm);}
|
||
inline void operator +=(const _CMatrix & sm) {PIMM_FOR_WB(c, r) m[c][r] += sm.m[c][r];}
|
||
inline void operator -=(const _CMatrix & sm) {PIMM_FOR_WB(c, r) m[c][r] -= sm.m[c][r];}
|
||
inline void operator *=(const Type & v) {PIMM_FOR_WB(c, r) m[c][r] *= v;}
|
||
inline void operator /=(const Type & v) {PIMM_FOR_WB(c, r) m[c][r] /= v;}
|
||
inline _CMatrix operator -() {_CMatrix tm(*this); PIMM_FOR_WB(c, r) tm.m[c][r] = -m[c][r]; return tm;}
|
||
inline _CMatrix operator +(const _CMatrix & sm) {_CMatrix tm(*this); PIMM_FOR_WB(c, r) tm.m[c][r] += sm.m[c][r]; return tm;}
|
||
inline _CMatrix operator -(const _CMatrix & sm) {_CMatrix tm(*this); PIMM_FOR_WB(c, r) tm.m[c][r] -= sm.m[c][r]; return tm;}
|
||
inline _CMatrix operator *(const Type & v) {_CMatrix tm(*this); PIMM_FOR_WB(c, r) tm.m[c][r] *= v; return tm;}
|
||
inline _CMatrix operator /(const Type & v) {_CMatrix tm(*this); PIMM_FOR_WB(c, r) tm.m[c][r] /= v; return tm;}
|
||
|
||
_CMatrix & toUpperTriangular(bool * ok = 0) {
|
||
if (cols_ != rows_) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
_CMatrix smat(*this);
|
||
bool ndet;
|
||
uint crow;
|
||
Type mul;
|
||
for (uint i = 0; i < cols_; ++i) {
|
||
ndet = true;
|
||
for (uint j = 0; j < rows_; ++j) if (smat.m[i][j] != 0) ndet = false;
|
||
if (ndet) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
for (uint j = 0; j < cols_; ++j) if (smat.m[j][i] != 0) ndet = false;
|
||
if (ndet) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
}
|
||
for (uint i = 0; i < cols_; ++i) {
|
||
crow = i;
|
||
while (smat.m[i][i] == Type(0))
|
||
smat.swaprows_(i, ++crow);
|
||
for (uint j = i + 1; j < rows_; ++j) {
|
||
mul = smat.m[i][j] / smat.m[i][i];
|
||
for (uint k = i; k < cols_; ++k) smat.m[k][j] -= mul * smat.m[k][i];
|
||
}
|
||
if (i < cols_ - 1) {
|
||
if (fabs(smat.m[i+1][i+1]) < Type(1E-100)) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
}
|
||
}
|
||
if (ok != 0) *ok = true;
|
||
m = smat.m;
|
||
return *this;
|
||
}
|
||
|
||
_CMatrix & invert(bool * ok = 0, _CMCol * sv = 0) {
|
||
if (cols_ != rows_) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
_CMatrix mtmp = _CMatrix::identity(cols_, rows_), smat(*this);
|
||
bool ndet;
|
||
uint crow;
|
||
Type mul, iddiv;
|
||
for (uint i = 0; i < cols_; ++i) {
|
||
ndet = true;
|
||
for (uint j = 0; j < rows_; ++j) if (smat.m[i][j] != 0) ndet = false;
|
||
if (ndet) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
for (uint j = 0; j < cols_; ++j) if (smat.m[j][i] != 0) ndet = false;
|
||
if (ndet) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
}
|
||
for (uint i = 0; i < cols_; ++i) {
|
||
crow = i;
|
||
while (smat.m[i][i] == Type(0)) {
|
||
++crow;
|
||
smat.swaprows_(i, crow);
|
||
mtmp.swaprows_(i, crow);
|
||
if (sv != 0) sv->swap(i, crow);
|
||
}
|
||
for (uint j = i + 1; j < rows_; ++j) {
|
||
mul = smat.m[i][j] / smat.m[i][i];
|
||
for (uint k = i; k < cols_; ++k) smat.m[k][j] -= mul * smat.m[k][i];
|
||
for (uint k = 0; k < cols_; ++k) mtmp.m[k][j] -= mul * mtmp.m[k][i];
|
||
if (sv != 0) (*sv)[j] -= mul * (*sv)[i];
|
||
}
|
||
//cout << i << endl << smat << endl;
|
||
if (i < cols_ - 1) {
|
||
if (fabs(smat.m[i+1][i+1]) < Type(1E-100)) {
|
||
if (ok != 0) *ok = false;
|
||
return *this;
|
||
}
|
||
}
|
||
iddiv = smat.m[i][i];
|
||
for (uint j = i; j < cols_; ++j) smat.m[j][i] /= iddiv;
|
||
for (uint j = 0; j < cols_; ++j) mtmp.m[j][i] /= iddiv;
|
||
if (sv != 0) (*sv)[i] /= iddiv;
|
||
}
|
||
for (uint i = cols_ - 1; i > 0; --i) {
|
||
for (uint j = 0; j < i; ++j) {
|
||
mul = smat.m[i][j];
|
||
smat.m[i][j] -= mul;
|
||
for (uint k = 0; k < cols_; ++k) mtmp.m[k][j] -= mtmp.m[k][i] * mul;
|
||
if (sv != 0) (*sv)[j] -= mul * (*sv)[i];
|
||
}
|
||
}
|
||
if (ok != 0) *ok = true;
|
||
m = mtmp.m;
|
||
return *this;
|
||
}
|
||
inline _CMatrix inverted(bool * ok = 0) {_CMatrix tm(*this); tm.invert(ok); return tm;}
|
||
inline _CMatrix transposed() {_CMatrix tm(rows_, cols_); PIMM_FOR_WB(c, r) tm[r][c] = m[c][r]; return tm;}
|
||
|
||
private:
|
||
uint cols_, rows_;
|
||
PIVector<PIVector<Type> > m;
|
||
|
||
};
|
||
|
||
template<typename Type>
|
||
inline std::ostream & operator <<(std::ostream & s, const PIMathMatrix<Type> & m) {s << '{'; for (uint r = 0; r < m.rows(); ++r) { for (uint c = 0; c < m.cols(); ++c) { s << m[c][r]; if (c < m.cols() - 1 || r < m.rows() - 1) s << ", ";} if (r < m.rows() - 1) s << endl << ' ';} s << '}'; return s;}
|
||
|
||
/// Multiply matrices {CR x Rows0} on {Cols1 x CR}, result is {Cols1 x Rows0}
|
||
template<typename Type>
|
||
inline PIMathMatrix<Type> operator *(const PIMathMatrix<Type> & fm,
|
||
const PIMathMatrix<Type> & sm) {
|
||
uint cr = fm.cols(), rows0 = fm.rows(), cols1 = sm.cols();
|
||
PIMathMatrix<Type> tm(cols1, rows0);
|
||
if (fm.cols() != sm.rows()) return tm;
|
||
Type t;
|
||
for (uint j = 0; j < rows0; ++j) {
|
||
for (uint i = 0; i < cols1; ++i) {
|
||
t = Type(0);
|
||
for (uint k = 0; k < cr; ++k)
|
||
t += fm[k][j] * sm[i][k];
|
||
tm[i][j] = t;
|
||
}
|
||
}
|
||
return tm;
|
||
}
|
||
|
||
/// Multiply matrix {Cols x Rows} on vector {Cols}, result is vector {Rows}
|
||
template<typename Type>
|
||
inline PIMathVector<Type> operator *(const PIMathMatrix<Type> & fm,
|
||
const PIMathVector<Type> & sv) {
|
||
uint c = fm.cols(), r = fm.rows();
|
||
PIMathVector<Type> tv(r);
|
||
if (c != sv.size()) return tv;
|
||
Type t;
|
||
for (uint i = 0; i < r; ++i) {
|
||
t = Type(0);
|
||
for (uint j = 0; j < c; ++j)
|
||
t += fm[j][i] * sv[j];
|
||
tv[i] = t;
|
||
}
|
||
return tv;
|
||
}
|
||
|
||
typedef PIMathMatrix<int> PIMathMatrixi;
|
||
typedef PIMathMatrix<double> PIMathMatrixd;
|
||
|
||
#undef PIMV_FOR
|
||
#undef PIMM_FOR
|
||
#undef PIMM_FOR_WB
|
||
#undef PIMM_FOR_I
|
||
#undef PIMM_FOR_I_WB
|
||
#undef PIMM_FOR_C
|
||
#undef PIMM_FOR_R
|
||
|
||
|
||
/*
|
||
Fast Fourier Transformation: direct (complement= false)
|
||
and complement (complement = true). 'x' is data source.
|
||
'x' contains 2^T items.
|
||
*/
|
||
void fft(complexd * x, int T, bool complement);
|
||
|
||
/// Differential evaluations
|
||
|
||
struct TransferFunction { // Для задания передаточной функции
|
||
PIVector<double> vector_Bm, vector_An;
|
||
};
|
||
|
||
// Класс, служащий для перевода передаточной функции в систему ОДУ первого порядка
|
||
// реализованы след. методы решения дифф. ур-ний:
|
||
// Эйлера
|
||
// Рунге-Кутта 4-го порядка
|
||
// Адамса-Бэшфортса-Моултона 2, 3, 4 порядков
|
||
class Solver
|
||
{
|
||
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
|
||
};
|
||
|
||
Solver() {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);
|
||
inline void solvePA2(double u, double h) {if (step > 0) solvePA(u, h, 2); else solvePA(u, h, 1);}
|
||
inline void solvePA3(double u, double h) {if (step > 1) solvePA(u, h, 3); else solvePA2(u, h);}
|
||
inline void solvePA4(double u, double h) {if (step > 2) solvePA(u, h, 4); else solvePA3(u, h);}
|
||
inline 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:
|
||
inline 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 // PIMATH_H
|