diff --git a/.kdev4/pip.kdev4 b/.kdev4/pip.kdev4 index 0de0e122..66546cad 100644 --- a/.kdev4/pip.kdev4 +++ b/.kdev4/pip.kdev4 @@ -28,5 +28,8 @@ Project Target=pip,pip Working Directory= isExecutable=false +[MakeBuilder] +Number Of Jobs=5 + [Project] VersionControlSupport=kdevgit diff --git a/Makefile b/Makefile deleted file mode 100644 index eeafe64f..00000000 --- a/Makefile +++ /dev/null @@ -1,460 +0,0 @@ -# CMAKE generated file: DO NOT EDIT! -# Generated by "Unix Makefiles" Generator, CMake Version 2.8 - -# Default target executed when no arguments are given to make. -default_target: all -.PHONY : default_target - -#============================================================================= -# Special targets provided by cmake. - -# Disable implicit rules so canoncical targets will work. -.SUFFIXES: - -# Remove some rules from gmake that .SUFFIXES does not remove. -SUFFIXES = - -.SUFFIXES: .hpux_make_needs_suffix_list - -# Suppress display of executed commands. -$(VERBOSE).SILENT: - -# A target that is always out of date. -cmake_force: -.PHONY : cmake_force - -#============================================================================= -# Set environment variables for the build. - -# The shell in which to execute make rules. -SHELL = /bin/sh - -# The CMake executable. -CMAKE_COMMAND = /usr/bin/cmake - -# The command to remove a file. -RM = /usr/bin/cmake -E remove -f - -# The program to use to edit the cache. -CMAKE_EDIT_COMMAND = /usr/bin/ccmake - -# The top-level source directory on which CMake was run. -CMAKE_SOURCE_DIR = /home/peri4/pprojects/pip - -# The top-level build directory on which CMake was run. -CMAKE_BINARY_DIR = /home/peri4/pprojects/pip - -#============================================================================= -# Targets provided globally by CMake. - -# Special rule for the target edit_cache -edit_cache: - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake cache editor..." - /usr/bin/ccmake -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) -.PHONY : edit_cache - -# Special rule for the target edit_cache -edit_cache/fast: edit_cache -.PHONY : edit_cache/fast - -# Special rule for the target rebuild_cache -rebuild_cache: - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake to regenerate build system..." - /usr/bin/cmake -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) -.PHONY : rebuild_cache - -# Special rule for the target rebuild_cache -rebuild_cache/fast: rebuild_cache -.PHONY : rebuild_cache/fast - -# The main all target -all: cmake_check_build_system - $(CMAKE_COMMAND) -E cmake_progress_start /home/peri4/pprojects/pip/CMakeFiles /home/peri4/pprojects/pip/CMakeFiles/progress.marks - $(MAKE) -f CMakeFiles/Makefile2 all - $(CMAKE_COMMAND) -E cmake_progress_start /home/peri4/pprojects/pip/CMakeFiles 0 -.PHONY : all - -# The main clean target -clean: - $(MAKE) -f CMakeFiles/Makefile2 clean -.PHONY : clean - -# The main clean target -clean/fast: clean -.PHONY : clean/fast - -# Prepare targets for installation. -preinstall: all - $(MAKE) -f CMakeFiles/Makefile2 preinstall -.PHONY : preinstall - -# Prepare targets for installation. -preinstall/fast: - $(MAKE) -f CMakeFiles/Makefile2 preinstall -.PHONY : preinstall/fast - -# clear depends -depend: - $(CMAKE_COMMAND) -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) --check-build-system CMakeFiles/Makefile.cmake 1 -.PHONY : depend - -#============================================================================= -# Target rules for targets named pip - -# Build rule for target. -pip: cmake_check_build_system - $(MAKE) -f CMakeFiles/Makefile2 pip -.PHONY : pip - -# fast build rule for target. -pip/fast: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/build -.PHONY : pip/fast - -main.o: main.cpp.o -.PHONY : main.o - -# target to build an object file -main.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/main.cpp.o -.PHONY : main.cpp.o - -main.i: main.cpp.i -.PHONY : main.i - -# target to preprocess a source file -main.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/main.cpp.i -.PHONY : main.cpp.i - -main.s: main.cpp.s -.PHONY : main.s - -# target to generate assembly for a file -main.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/main.cpp.s -.PHONY : main.cpp.s - -piconfig.o: piconfig.cpp.o -.PHONY : piconfig.o - -# target to build an object file -piconfig.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piconfig.cpp.o -.PHONY : piconfig.cpp.o - -piconfig.i: piconfig.cpp.i -.PHONY : piconfig.i - -# target to preprocess a source file -piconfig.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piconfig.cpp.i -.PHONY : piconfig.cpp.i - -piconfig.s: piconfig.cpp.s -.PHONY : piconfig.s - -# target to generate assembly for a file -piconfig.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piconfig.cpp.s -.PHONY : piconfig.cpp.s - -piconsole.o: piconsole.cpp.o -.PHONY : piconsole.o - -# target to build an object file -piconsole.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piconsole.cpp.o -.PHONY : piconsole.cpp.o - -piconsole.i: piconsole.cpp.i -.PHONY : piconsole.i - -# target to preprocess a source file -piconsole.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piconsole.cpp.i -.PHONY : piconsole.cpp.i - -piconsole.s: piconsole.cpp.s -.PHONY : piconsole.s - -# target to generate assembly for a file -piconsole.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piconsole.cpp.s -.PHONY : piconsole.cpp.s - -piethernet.o: piethernet.cpp.o -.PHONY : piethernet.o - -# target to build an object file -piethernet.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piethernet.cpp.o -.PHONY : piethernet.cpp.o - -piethernet.i: piethernet.cpp.i -.PHONY : piethernet.i - -# target to preprocess a source file -piethernet.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piethernet.cpp.i -.PHONY : piethernet.cpp.i - -piethernet.s: piethernet.cpp.s -.PHONY : piethernet.s - -# target to generate assembly for a file -piethernet.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piethernet.cpp.s -.PHONY : piethernet.cpp.s - -pifile.o: pifile.cpp.o -.PHONY : pifile.o - -# target to build an object file -pifile.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pifile.cpp.o -.PHONY : pifile.cpp.o - -pifile.i: pifile.cpp.i -.PHONY : pifile.i - -# target to preprocess a source file -pifile.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pifile.cpp.i -.PHONY : pifile.cpp.i - -pifile.s: pifile.cpp.s -.PHONY : pifile.s - -# target to generate assembly for a file -pifile.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pifile.cpp.s -.PHONY : pifile.cpp.s - -pikbdlistener.o: pikbdlistener.cpp.o -.PHONY : pikbdlistener.o - -# target to build an object file -pikbdlistener.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pikbdlistener.cpp.o -.PHONY : pikbdlistener.cpp.o - -pikbdlistener.i: pikbdlistener.cpp.i -.PHONY : pikbdlistener.i - -# target to preprocess a source file -pikbdlistener.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pikbdlistener.cpp.i -.PHONY : pikbdlistener.cpp.i - -pikbdlistener.s: pikbdlistener.cpp.s -.PHONY : pikbdlistener.s - -# target to generate assembly for a file -pikbdlistener.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pikbdlistener.cpp.s -.PHONY : pikbdlistener.cpp.s - -piprotocol.o: piprotocol.cpp.o -.PHONY : piprotocol.o - -# target to build an object file -piprotocol.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piprotocol.cpp.o -.PHONY : piprotocol.cpp.o - -piprotocol.i: piprotocol.cpp.i -.PHONY : piprotocol.i - -# target to preprocess a source file -piprotocol.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piprotocol.cpp.i -.PHONY : piprotocol.cpp.i - -piprotocol.s: piprotocol.cpp.s -.PHONY : piprotocol.s - -# target to generate assembly for a file -piprotocol.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piprotocol.cpp.s -.PHONY : piprotocol.cpp.s - -piserial.o: piserial.cpp.o -.PHONY : piserial.o - -# target to build an object file -piserial.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piserial.cpp.o -.PHONY : piserial.cpp.o - -piserial.i: piserial.cpp.i -.PHONY : piserial.i - -# target to preprocess a source file -piserial.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piserial.cpp.i -.PHONY : piserial.cpp.i - -piserial.s: piserial.cpp.s -.PHONY : piserial.s - -# target to generate assembly for a file -piserial.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/piserial.cpp.s -.PHONY : piserial.cpp.s - -pistring.o: pistring.cpp.o -.PHONY : pistring.o - -# target to build an object file -pistring.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pistring.cpp.o -.PHONY : pistring.cpp.o - -pistring.i: pistring.cpp.i -.PHONY : pistring.i - -# target to preprocess a source file -pistring.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pistring.cpp.i -.PHONY : pistring.cpp.i - -pistring.s: pistring.cpp.s -.PHONY : pistring.s - -# target to generate assembly for a file -pistring.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pistring.cpp.s -.PHONY : pistring.cpp.s - -pithread.o: pithread.cpp.o -.PHONY : pithread.o - -# target to build an object file -pithread.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pithread.cpp.o -.PHONY : pithread.cpp.o - -pithread.i: pithread.cpp.i -.PHONY : pithread.i - -# target to preprocess a source file -pithread.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pithread.cpp.i -.PHONY : pithread.cpp.i - -pithread.s: pithread.cpp.s -.PHONY : pithread.s - -# target to generate assembly for a file -pithread.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pithread.cpp.s -.PHONY : pithread.cpp.s - -pitimer.o: pitimer.cpp.o -.PHONY : pitimer.o - -# target to build an object file -pitimer.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pitimer.cpp.o -.PHONY : pitimer.cpp.o - -pitimer.i: pitimer.cpp.i -.PHONY : pitimer.i - -# target to preprocess a source file -pitimer.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pitimer.cpp.i -.PHONY : pitimer.cpp.i - -pitimer.s: pitimer.cpp.s -.PHONY : pitimer.s - -# target to generate assembly for a file -pitimer.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pitimer.cpp.s -.PHONY : pitimer.cpp.s - -pivariable.o: pivariable.cpp.o -.PHONY : pivariable.o - -# target to build an object file -pivariable.cpp.o: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pivariable.cpp.o -.PHONY : pivariable.cpp.o - -pivariable.i: pivariable.cpp.i -.PHONY : pivariable.i - -# target to preprocess a source file -pivariable.cpp.i: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pivariable.cpp.i -.PHONY : pivariable.cpp.i - -pivariable.s: pivariable.cpp.s -.PHONY : pivariable.s - -# target to generate assembly for a file -pivariable.cpp.s: - $(MAKE) -f CMakeFiles/pip.dir/build.make CMakeFiles/pip.dir/pivariable.cpp.s -.PHONY : pivariable.cpp.s - -# Help Target -help: - @echo "The following are some of the valid targets for this Makefile:" - @echo "... all (the default if no target is provided)" - @echo "... clean" - @echo "... depend" - @echo "... edit_cache" - @echo "... pip" - @echo "... rebuild_cache" - @echo "... main.o" - @echo "... main.i" - @echo "... main.s" - @echo "... piconfig.o" - @echo "... piconfig.i" - @echo "... piconfig.s" - @echo "... piconsole.o" - @echo "... piconsole.i" - @echo "... piconsole.s" - @echo "... piethernet.o" - @echo "... piethernet.i" - @echo "... piethernet.s" - @echo "... pifile.o" - @echo "... pifile.i" - @echo "... pifile.s" - @echo "... pikbdlistener.o" - @echo "... pikbdlistener.i" - @echo "... pikbdlistener.s" - @echo "... piprotocol.o" - @echo "... piprotocol.i" - @echo "... piprotocol.s" - @echo "... piserial.o" - @echo "... piserial.i" - @echo "... piserial.s" - @echo "... pistring.o" - @echo "... pistring.i" - @echo "... pistring.s" - @echo "... pithread.o" - @echo "... pithread.i" - @echo "... pithread.s" - @echo "... pitimer.o" - @echo "... pitimer.i" - @echo "... pitimer.s" - @echo "... pivariable.o" - @echo "... pivariable.i" - @echo "... pivariable.s" -.PHONY : help - - - -#============================================================================= -# Special targets to cleanup operation of make. - -# Special rule to run CMake to check the build system integrity. -# No rule that depends on this can have commands that come from listfiles -# because they might be regenerated. -cmake_check_build_system: - $(CMAKE_COMMAND) -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) --check-build-system CMakeFiles/Makefile.cmake 0 -.PHONY : cmake_check_build_system - diff --git a/clean b/clean new file mode 100755 index 00000000..ac8674c1 --- /dev/null +++ b/clean @@ -0,0 +1,19 @@ +#! /bin/bash +VERBOSE=1 make clean +rm -vf ./lib/* +#for i in $( ls -1 ); do +# if [ "`ls -1 --file-type | grep $i | grep -o /`" = "/" ]; then +# cd $i +# rm -rvf ./CMakeFiles +# rm -vf ./CMakeCache.txt ./Makefile ./cmake_install.cmake ./install_manifest.txt ./*~ +# cd ../ +# fi +#done +rm -rvf ./CMakeFiles +rm -vf ./CMakeCache.txt ./Makefile ./cmake_install.cmake ./install_manifest.txt ./*~ ./*cxx ./moc_* ./*.o +#cd ./include +#for i in $( ls -1 ); do +# if [ "`ls -1 --file-type $i | grep -v @`" = "" ]; then +# rm -v $i +# fi +#done diff --git a/main.cpp b/main.cpp index d0c7119f..9eccf944 100644 --- a/main.cpp +++ b/main.cpp @@ -1,10 +1,26 @@ #include "pip.h" int main(int argc, char * argv[]) { - PIMathMatrix<3, 2> m0(1., 2., 3., - 4., 5., 6.); - PIMathMatrix<3, 2> m1(1., -2., 3., - -4., 5., -6.); - m0 += m1; - cout << m0 << endl; + /*PIString s("I love Kusia! <3"); + PIByteArray c(s.data(), s.size()); + c.compressHuffman(); + cout << s << endl << c << endl; + //c.decompressRLE(128); + cout << c << endl;*/ + PIBitArray a(0xFFF00Fu); + //cout << a << endl; + a.setBit(31); + cout << a << endl; + a.push_back(true); + cout << a << endl; + a.push_front(true); + cout << a << endl; + a.push_front(false); + cout << a << endl; + a.pop_front(); + cout << a << endl; + a.pop_front(); + cout << a << endl; + a.pop_front(); + cout << a << endl; }; diff --git a/make.sh b/make.sh new file mode 100755 index 00000000..d218016a --- /dev/null +++ b/make.sh @@ -0,0 +1,3 @@ +#! /bin/bash +cmake ./ +make $@ diff --git a/pibitarray.h b/pibitarray.h new file mode 100644 index 00000000..7446e049 --- /dev/null +++ b/pibitarray.h @@ -0,0 +1,84 @@ +#ifndef PIBITARRAY_H +#define PIBITARRAY_H + +#include "piincludes.h" + +class PIBitArray { +public: + inline PIBitArray(const int & size = 0) {resize(size);} + inline PIBitArray(uchar val) {resize(sizeof(val) * 8); data_[0] = val;} + inline PIBitArray(ushort val) {resize(sizeof(val) * 8); memcpy(data(), &val, sizeof(val));} + inline PIBitArray(uint val) {resize(sizeof(val) * 8); memcpy(data(), &val, sizeof(val));} + inline PIBitArray(ulong val) {resize(sizeof(val) * 8); memcpy(data(), &val, sizeof(val));} + inline PIBitArray(ullong val) {resize(sizeof(val) * 8); memcpy(data(), &val, sizeof(val));} + inline PIBitArray(uchar * bytes, uint size) {resize(size * 8); memcpy(data(), bytes, size);} + + inline uint bitSize() const {return size_;} + inline uint byteSize() const {return bytesInBits(size_);} + inline PIBitArray & resize(const uint & size) {size_ = size; data_.resize(bytesInBits(size_)); return *this;} + + inline PIBitArray & clearBit(const uint & index) {data_[index / 8] &= ~(1 << (index % 8)); return *this;} + inline PIBitArray & setBit(const uint & index) {data_[index / 8] |= (1 << (index % 8)); return *this;} + inline PIBitArray & writeBit(const uint & index, const bool & value) {if (value) setBit(index); else clearBit(index); return *this;} + inline PIBitArray & writeBit(const uint & index, const uchar & value) {return writeBit(index, static_cast(value));} + + inline PIBitArray & push_back(const bool & value) {resize(size_ + 1); writeBit(size_ - 1, value); return *this;} + inline PIBitArray & push_back(const uchar & value) {return push_back(static_cast(value));} + inline PIBitArray & insert(const uint & index, const bool & value) { + resize(size_ + 1); + uint fi = byteSize() - 1, si = index / 8, ti = index % 8; + uchar c = data_[si]; + for (uint i = fi; i > si; --i) { + data_[i] <<= 1; + if ((0x80 & data_[i - 1]) == 0x80) data_[i] |= 1; + else data_[i] &= 0xFE;} + data_[si] &= (0xFF >> (7 - ti)); + data_[si] |= ((c << 1) & (0xFF << (ti))); + if (value) data_[si] |= (1 << ti); + else data_[si] &= ~(1 << ti); + return *this;} + inline PIBitArray & insert(const uint & index, const uchar & value) {return push_back(static_cast(value));} + inline PIBitArray & push_front(const bool & value) {return insert(0, value);} + inline PIBitArray & push_front(const uchar & value) {return push_front(static_cast(value));} + inline PIBitArray & pop_back() {return resize(size_ - 1);} + inline PIBitArray & pop_front() { + if (size_ == 0) return *this; + uint fi = byteSize() - 1; + for (uint i = 0; i < fi; ++i) { + data_[i] >>= 1; + if ((1 & data_[i + 1]) == 1) data_[i] |= 0x80; + else data_[i] &= 0x7F;} + data_[fi] >>= 1; + resize(size_ - 1); + return *this;} + inline PIBitArray & append(const PIBitArray & ba) {for (uint i = 0; i < ba.bitSize(); ++i) push_back(ba[i]); return *this;} + + inline uchar * data() {return data_.data();} + inline uchar toUChar() {if (size_ == 0) return 0; return data_[0];} + inline ushort toUShort() {ushort t = 0; memcpy(&t, data(), piMin(byteSize(), sizeof(t))); return t;} + inline uint toUInt() {uint t = 0; memcpy(&t, data(), piMin(byteSize(), sizeof(t))); return t;} + inline ulong toULong() {ulong t = 0; memcpy(&t, data(), piMin(byteSize(), sizeof(t))); return t;} + inline ullong toULLong() {ullong t = 0; memcpy(&t, data(), piMin(byteSize(), sizeof(t))); return t;} + + inline bool at(const uint & index) const {return (1 & (data_[index / 8] >> (index % 8))) == 1 ? true : false;} + inline bool operator [](const uint & index) const {return at(index);} + inline void operator +=(const PIBitArray & ba) {append(ba);} + inline bool operator ==(const PIBitArray & ba) const {if (bitSize() != ba.bitSize()) return false; for (uint i = 0; i < bitSize(); ++i) if (at(i) != ba[i]) return false; return true;} + inline bool operator !=(const PIBitArray & ba) const {return !(*this == ba);} + inline void operator =(const uchar & val) {resize(sizeof(val) * 8); data_[0] = val;} + inline void operator =(const ushort & val) {resize(sizeof(val) * 8); memcpy(data(), &val, sizeof(val));} + inline void operator =(const uint & val) {resize(sizeof(val) * 8); memcpy(data(), &val, sizeof(val));} + inline void operator =(const ulong & val) {resize(sizeof(val) * 8); memcpy(data(), &val, sizeof(val));} + inline void operator =(const ullong & val) {resize(sizeof(val) * 8); memcpy(data(), &val, sizeof(val));} + +private: + inline uint bytesInBits(const uint & bits) const {return (bits + 7) / 8;} + + PIVector data_; + uint size_; + +}; + +inline std::ostream & operator <<(std::ostream & s, const PIBitArray & ba) {for (uint i = 0; i < ba.bitSize(); ++i) {s << ba[i]; if (i % 8 == 7) s << ' ';} return s;} + +#endif // PIBITARRAY_H diff --git a/pibytearray.cpp b/pibytearray.cpp new file mode 100644 index 00000000..3bac7ed2 --- /dev/null +++ b/pibytearray.cpp @@ -0,0 +1,140 @@ +#include "pibytearray.h" + + +int PIHuffman::nodeCompare(const void * f, const void * s) { + return (reinterpret_cast(const_cast(s))->freq - + reinterpret_cast(const_cast(f))->freq); +} + + +PIVector PIHuffman::compress(const PIVector & src) { + calcFrequencies(src); + return src; +} + + +void PIHuffman::calcFrequencies(const PIVector & src) { + nodes.resize(256); + for (uint i = 0; i < 256; ++i) { + nodes[i].parent = nodes[i].right = nodes[i].left = 0; + nodes[i].freq = 0; + nodes[i].word.resize(1); + nodes[i].word[0] = static_cast(i); + } + for (uint i = 0; i < src.size(); ++i) + nodes[src[i]].freq++; + std::qsort(nodes.data(), 256, sizeof(node), nodeCompare); + for (uint i = 255; i >= 0; --i) + if (nodes[i].freq > 0 && i < 255) + {nodes.remove(i + 1, 255 - i); break;} + for (uint i = 0; i < nodes.size(); ++i) + cout << string((char*)nodes[i].word.data(), 1) << ": " << nodes[i].freq << endl; +} + + +PIHuffman PIByteArray::huffman; + +PIByteArray & PIByteArray::convertToBase64() { + base64HelpStruct hs; + PIByteArray t; + if (size() == 0) return *this; + uint sz = (size() / 3) * 3; + for (uint i = 0; i < sz; ++i) { + hs.byte0 = hs.byte1 = hs.byte2 = 0; + hs.byte0 = at(i); + hs.byte1 = at(++i); + hs.byte2 = at(++i); + t.push_back(base64Table[hs.ascii0]); + t.push_back(base64Table[hs.ascii1]); + t.push_back(base64Table[hs.ascii2]); + t.push_back(base64Table[hs.ascii3]); + } + hs.byte0 = hs.byte1 = hs.byte2 = 0; sz = size() % 3; + switch (sz) { + case 1: + hs.byte0 = back(); + t.push_back(base64Table[hs.ascii0]); + t.push_back(base64Table[hs.ascii1]); + t.push_back('='); + t.push_back('='); + break; + case 2: + hs.byte0 = at(size() - 2); hs.byte1 = back(); + t.push_back(base64Table[hs.ascii0]); + t.push_back(base64Table[hs.ascii1]); + t.push_back(base64Table[hs.ascii2]); + t.push_back('='); + break; + default: break; + } + *this = t; + return *this; +} + + +PIByteArray & PIByteArray::convertFromBase64() { + base64HelpStruct hs; + PIByteArray t; + uint sz = size(); + if (sz == 0) return *this; + for (uint i = 0; i < sz; ++i) { + hs.byte0 = hs.byte1 = hs.byte2 = 0; + hs.ascii0 = base64InvTable[at(i)]; + hs.ascii1 = base64InvTable[at(++i)]; + hs.ascii2 = base64InvTable[at(++i)]; + hs.ascii3 = base64InvTable[at(++i)]; + t.push_back(hs.byte0); + t.push_back(hs.byte1); + t.push_back(hs.byte2); + } + if (back() == '=') t.pop_back(); + if (sz > 1) if (at(sz - 2) == '=') t.pop_back(); + *this = t; + return *this; +} + + +PIByteArray & PIByteArray::compressRLE(uchar threshold) { + PIByteArray t; + uchar fb, clen, mlen = 255 - threshold; + for (uint i = 0; i < size();) { + fb = at(i); + clen = 1; + while (at(++i) == fb) { + ++clen; + if (clen == mlen) + break; + } + if (clen > 1) { + t.push_back(threshold + clen); + t.push_back(fb); + continue; + } + if (fb >= threshold) { + t.push_back(threshold + 1); + t.push_back(fb); + } else + t.push_back(fb); + } + *this = t; + return *this; +} + + +PIByteArray & PIByteArray::decompressRLE(uchar threshold) { + PIByteArray t; + uchar fb, clen; + for (uint i = 0; i < size(); ++i) { + fb = at(i); + if (fb >= threshold) { + clen = fb - threshold; + fb = at(++i); + for (uint j = 0; j < clen; ++j) + t.push_back(fb); + continue; + } else + t.push_back(fb); + } + *this = t; + return *this; +} diff --git a/pibytearray.h b/pibytearray.h new file mode 100644 index 00000000..435ca46f --- /dev/null +++ b/pibytearray.h @@ -0,0 +1,114 @@ +#ifndef PIBYTEARRAY_H +#define PIBYTEARRAY_H + +#include "pibitarray.h" + +#pragma pack(push, 1) + +static const char base64Table[64] = { +0x41, 0x42, 0x43, 0x44, 0x45, 0x46, 0x47, 0x48, +0x49, 0x4a, 0x4b, 0x4c, 0x4d, 0x4e, 0x4f, 0x50, +0x51, 0x52, 0x53, 0x54, 0x55, 0x56, 0x57, 0x58, +0x59, 0x5a, 0x61, 0x62, 0x63, 0x64, 0x65, 0x66, +0x67, 0x68, 0x69, 0x6a, 0x6b, 0x6c, 0x6d, 0x6e, +0x6f, 0x70, 0x71, 0x72, 0x73, 0x74, 0x75, 0x76, +0x77, 0x78, 0x79, 0x7a, 0x30, 0x31, 0x32, 0x33, +0x34, 0x35, 0x36, 0x37, 0x38, 0x39, 0x2b, 0x2f}; + +static const char base64InvTable[256] = { +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x3E, 0x0, 0x0, 0x0, 0x3F, +0x34, 0x35, 0x36, 0x37, 0x38, 0x39, 0x3A, 0x3B, +0x3C, 0x3D, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x1, 0x2, 0x3, 0x4, 0x5, 0x6, +0x7, 0x8, 0x9, 0xA, 0xB, 0xC, 0xD, 0xE, +0xF, 0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, +0x17, 0x18, 0x19, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x1A, 0x1B, 0x1C, 0x1D, 0x1E, 0x1F, 0x20, +0x21, 0x22, 0x23, 0x24, 0x25, 0x26, 0x27, 0x28, +0x29, 0x2A, 0x2B, 0x2C, 0x2D, 0x2E, 0x2F, 0x30, +0x31, 0x32, 0x33, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, +0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0}; + +#pragma pack(pop) + +class PIByteArray; + +class PIHuffman { +public: + PIVector compress(const PIVector & src); + +private: + struct node { + int freq; + PIVector word; + PIBitArray path; + node * parent; + node * right; + node * left; + }; + static int nodeCompare(const void * f, const void * s); + void calcFrequencies(const PIVector & src); + PIVector nodes; +}; + +class PIByteArray: public PIVector { +public: + PIByteArray() {;} + PIByteArray(const char * data, uint size) {for (uint i = 0; i < size; ++i) push_back(data[i]);} + PIByteArray(const uchar * data, uint size) {for (uint i = 0; i < size; ++i) push_back(data[i]);} + + PIByteArray & convertToBase64(); + PIByteArray & convertFromBase64(); + PIByteArray & compressRLE(uchar threshold = 192); + PIByteArray & decompressRLE(uchar threshold = 192); + inline PIByteArray & compressHuffman() {*this = huffman.compress(*this); return *this;} + + inline PIByteArray toBase64() {PIByteArray ba(*this); ba.convertToBase64(); return ba;} + inline PIByteArray fromBase64() {PIByteArray ba(*this); ba.convertFromBase64(); return ba;} + inline PIByteArray compressedRLE(uchar threshold = 192) {PIByteArray ba(*this); ba.compressedRLE(threshold); return ba;} + inline PIByteArray decompressedRLE(uchar threshold = 192) {PIByteArray ba(*this); ba.decompressedRLE(threshold); return ba;} + + inline void operator =(const PIVector & d) {resize(d.size()); for (uint i = 0; i < size(); ++i) (*this)[i] = d[i];} + +private: + union base64HelpStruct { + struct { + uchar ascii0: 6; + uchar ascii1: 6; + uchar ascii2: 6; + uchar ascii3: 6; + }; + struct { + uchar byte0; + uchar byte1; + uchar byte2; + }; + }; + + static PIHuffman huffman; + +}; + +inline std::ostream & operator <<(std::ostream & s, const PIByteArray & ba) {for (uint i = 0; i < ba.size(); ++i) s << ba[i]; return s;} + +#endif // PIBYTEARRAY_H diff --git a/piconfig.cpp b/piconfig.cpp index 742e1404..8983907d 100644 --- a/piconfig.cpp +++ b/piconfig.cpp @@ -8,6 +8,8 @@ PIConfig::PIConfig(const PIString & path, Flags mode): PIFile(path, mode) } +PIStringList PIConfig::getValue(const PIString & vname, const PIStringList & def, bool * exist) const { + return getValue(vname, def.join("%|%"), exist).split("%|%");} bool PIConfig::getValue(const PIString & vname, const bool def, bool * exist) const { return atob(getValue(vname, btos(def), exist));} char PIConfig::getValue(const PIString & vname, const char def, bool * exist) const { @@ -32,6 +34,8 @@ double PIConfig::getValue(const PIString & vname, const double def, bool * exist return getValue(vname, dtos(def), exist).toDouble();} +void PIConfig::setValue(const PIString & name, const PIStringList & value, bool write) { + setValue(name, value.join("%|%"), "l", write);} void PIConfig::setValue(const PIString & name, const char * value, bool write) { setValue(name, PIString(value), "s", write);} void PIConfig::setValue(const PIString & name, const bool value, bool write) { diff --git a/piconfig.h b/piconfig.h index d465a719..6e738895 100644 --- a/piconfig.h +++ b/piconfig.h @@ -10,7 +10,8 @@ public: ~PIConfig() {;} PIString getValue(const PIString & vname, const PIString & def = "", bool * exist = 0) const; - PIString getValue(const PIString & vname, const char * def = "", bool * exist = 0) const; + PIString getValue(const PIString & vname, const char * def = "", bool * exist = 0) const {return getValue(vname, PIString(def), exist);} + PIStringList getValue(const PIString & vname, const PIStringList & def = PIStringList(), bool * exist = 0) const; bool getValue(const PIString & vname, const bool def = false, bool * exist = 0) const; char getValue(const PIString & vname, const char def = 0, bool * exist = 0) const; short getValue(const PIString & vname, const short def = 0, bool * exist = 0) const; @@ -24,6 +25,7 @@ public: double getValue(const PIString & vname, const double def = 0., bool * exist = 0) const; void setValue(const PIString & name, const PIString & value, const PIString & type = "s", bool write = true); + void setValue(const PIString & name, const PIStringList & value, bool write = true); void setValue(const PIString & name, const char * value, bool write = true); void setValue(const PIString & name, const bool value, bool write = true); void setValue(const PIString & name, const char value, bool write = true); diff --git a/pidir.cpp b/pidir.cpp new file mode 100644 index 00000000..23eae4af --- /dev/null +++ b/pidir.cpp @@ -0,0 +1,76 @@ +#include "pidir.h" + + +PIDir::PIDir(const PIString & dir) { + path_ = dir.stdString(); + dir_ = 0; + open(); +} + + +bool PIDir::operator==(const PIDir & d) const { + return true; +} + + +bool PIDir::open() { + if (dir_ != 0) + if (!close()) + return false; + dir_ = opendir(path_.c_str()); + return (dir_ != 0); +} + + +bool PIDir::close() { + if (dir_ == 0) return true; + if (closedir(dir_) == 0) { + dir_ = 0; + return true; + } + return false; +} + + +PIString PIDir::absolutePath() { + //di; + return PIString(); +} + + +PIDir PIDir::current() { + char rc[1024]; + if (getcwd(rc, 1024) == 0) return PIString(); + return PIDir(rc); +} + + +PIDir PIDir::home() { + return PIDir(getenv("HOME")); +} + + +PIDir PIDir::temporary() { + char * rc; + rc = tmpnam(0); + PIString s(rc); + return PIDir(s.left(s.findLast(PIDir::separator))); +} + + +bool PIDir::mkDir(const PIString & path) { + string s = path.stdString(); + int ret = mkdir(s.c_str(), 16877); + if (ret == 0) return true; + printf("[PIDir] mkDir(\"%s\") error: %s\n", s.c_str(), strerror(errno)); + return false; +} + + +bool PIDir::rmDir(const PIString & path) { + string s = path.stdString(); + int ret = rmdir(s.c_str()); + if (ret == 0) return true; + printf("[PIDir] rmDir(\"%s\") error: %s\n", s.c_str(), strerror(errno)); + return false; +} diff --git a/pidir.h b/pidir.h new file mode 100644 index 00000000..35b53867 --- /dev/null +++ b/pidir.h @@ -0,0 +1,42 @@ +#ifndef PIDIR_H +#define PIDIR_H + +#include "pifile.h" +#include + +class PIDir +{ +public: + PIDir() {dir_ = 0;} + PIDir(const PIString & dir); + ~PIDir() {close();} + + inline const bool isExists() {return (dir_ != 0);} + inline PIString path() {return PIString(path_);} + PIString absolutePath(); + + bool operator ==(const PIDir & d) const; + +#ifdef WINDOWS + static const char separator = '\\'; +#else + static const char separator = '/'; +#endif + + static PIDir current(); + static PIDir home(); + static PIDir temporary(); + static bool mkDir(const PIString & path); + static bool rmDir(const PIString & path); + //static bool rmDir(const PIString & path); + +private: + bool open(); + bool close(); + + string path_; + DIR * dir_; + +}; + +#endif // PIDIR_H diff --git a/pifile.h b/pifile.h index dfb828b4..7899bb73 100644 --- a/pifile.h +++ b/pifile.h @@ -1,9 +1,9 @@ #ifndef PIFILE_H #define PIFILE_H -#include #include "piincludes.h" #include "pistring.h" +#include using std::fstream; @@ -52,6 +52,7 @@ public: PIFile & operator <<(const char & v) {stream.write(&v, 1); return *this;} //PIFile & operator <<(const string & v) {stream.write(v.c_str(), v.size()); return *this;} PIFile & operator <<(const PIString & v) {stream.write(v.data(), v.size()); return *this;} + PIFile & operator <<(const PIByteArray & v) {stream.write((char * )v.data(), v.size()); return *this;} PIFile & operator <<(const short & v) {stream << v; return *this;} PIFile & operator <<(const int & v) {stream << v; return *this;} PIFile & operator <<(const long & v) {stream << v; return *this;} diff --git a/pigeometry.h b/pigeometry.h index ceead339..5ebdd898 100644 --- a/pigeometry.h +++ b/pigeometry.h @@ -12,7 +12,9 @@ public: PIPoint() {x = y = 0;}; PIPoint(Type x_, Type y_) {set(x_, y_);} - inline void set(Type x_, Type y_) {x = x_; y = y_;} + inline PIPoint & set(Type x_, Type y_) {x = x_; y = y_; return *this;} + inline PIPoint & move(Type x_, Type y_) {x += x_; y += y_; return *this;} + inline PIPoint & move(const PIPoint & p) {x += p.x; y += p.y; return *this;} inline double angleRad() const {return atan2(y, x);} inline int angleDeg() const {return round(atan2(y, x) * 180. / M_PI);} @@ -27,6 +29,9 @@ public: inline bool operator !=(const PIPoint & p) const {return (x != p.x || y != p.y);} }; +template +inline std::ostream & operator <<(std::ostream & s, const PIPoint & v) {s << '{' << v.x << ", " << v.y << '}'; return s;} + template class PIRect { public: @@ -38,8 +43,10 @@ public: PIRect() {x0 = y0 = x1 = y1 = 0;}; PIRect(Type x, Type y, Type w, Type h) {set(x, y, w, h);} PIRect(const PIPoint & tl, const PIPoint & br) {set(tl.x, tl.y, br.x, br.y);} + PIRect(const PIPoint & p0, const PIPoint & p1, const PIPoint & p2) {set(piMin(p0.x, p1.x, p2.x), piMin(p0.y, p1.y, p2.y), + piMax(p0.x, p1.x, p2.x), piMax(p0.y, p1.y, p2.y));} - inline void set(Type x, Type y, Type w, Type h) {x0 = x; y0 = y; x1 = x + w; y1 = y + h;} + inline PIRect & set(Type x, Type y, Type w, Type h) {x0 = x; y0 = y; x1 = x + w; y1 = y + h; return *this;} inline bool pointIn(Type x, Type y) const {return (x <= x1 && x >= x0 && y <= y1 && y >= y0);} inline bool pointIn(const PIPoint & p) const {return pointIn(p.x, p.y);} inline bool isEmpty() const {return (x1 - x0 == 0 && y1 - y0 == 0);} @@ -51,8 +58,12 @@ public: inline PIRect & scale(const PIPoint & p) {setWidth(width() * p.x); setHeight(height() * p.y); return *this;} inline PIRect scaled(Type x, Type y) {PIRect r(*this); r.scale(x, y); return r;} inline PIRect scaled(const PIPoint & p) {PIRect r(*this); r.scale(p); return r;} - //inline void unite(const PIRect & r) {;} - //inline PIRect & united(const PIRect & r) {unite(r); return *this;} + inline PIRect & normalize() {if (x0 > x1) piSwap(x0, x1); if (y0 > y1) piSwap(y0, y1); return *this;} + inline PIRect normalized() {PIRect r(*this); r.normalize(); return r;} + inline PIRect & unite(const PIRect & r) {x0 = piMin(x0, r.x0); y0 = piMin(y0, r.y0); x1 = piMax(x1, r.x1); y1 = piMax(y1, r.y1); return *this;} + inline PIRect united(const PIRect & rect) {PIRect r(*this); r.unite(rect); return r;} + inline PIRect & intersect(const PIRect & r) {x0 = piMax(x0, r.x0); y0 = piMax(y0, r.y0); x1 = piMin(x1, r.x1); y1 = piMin(y1, r.y1); if (x0 > x1 || y0 > y1) x0 = x1 = y0 = y1 = Type(0); return *this;} + inline PIRect intersected(const PIRect & rect) {PIRect r(*this); r.intersect(rect); return r;} inline Type top() const {return y0;} inline Type left() const {return x0;} inline Type right() const {return x1;} @@ -77,13 +88,24 @@ public: inline void operator -=(const PIPoint & p) {translate(-p);} inline void operator *=(Type p) {x0 *= p; x1 *= p; y0 *= p; y1 *= p;} inline void operator /=(Type p) {x0 /= p; x1 /= p; y0 /= p; y1 /= p;} + inline void operator |=(const PIRect & r) {unite(r);} + inline void operator &=(const PIRect & r) {intersect(r);} + inline PIRect operator +(const PIPoint & p) {return PIRect(*this).translated(p);} + inline PIRect operator -(const PIPoint & p) {return PIRect(*this).translated(-p);} + inline PIRect operator |(const PIRect & r) {return PIRect(*this).united(r);} + inline PIRect operator &(const PIRect & r) {return PIRect(*this).intersected(r);} }; +template +inline std::ostream & operator <<(std::ostream & s, const PIRect & v) {s << '{' << v.x0 << ", " << v.y0 << "; " << v.x1 - v.x0 << ", " << v.y1 - v.y0 << '}'; return s;} + typedef PIPoint PIPointi; +typedef PIPoint PIPointu; typedef PIPoint PIPointf; typedef PIPoint PIPointd; typedef PIRect PIRecti; +typedef PIRect PIRectu; typedef PIRect PIRectf; typedef PIRect PIRectd; diff --git a/piincludes.h b/piincludes.h index 4229160a..86de29c8 100644 --- a/piincludes.h +++ b/piincludes.h @@ -10,8 +10,16 @@ #include #include -#include -#include +#ifndef QNX +# include +# include +#else +# include +# include +#endif +#include +#include +#include #include #include #include @@ -51,6 +59,7 @@ public: inline Type & at(uint index) {return (*this)[index];} inline const Type * data(uint index = 0) const {return &(*this)[index];} inline Type * data(uint index = 0) {return &(*this)[index];} + inline void fill(const Type & t) {vector::assign(vector::size(), t);} inline void pop_front() {vector::erase(vector::begin());} inline void push_front(const Type & t) {vector::insert(vector::begin(), t);} inline void remove(uint num) {vector::erase(vector::begin() + num);} @@ -58,6 +67,17 @@ public: inline void insert(uint pos, const Type & t) {vector::insert(vector::begin() + pos, t);} }; +template > +class PIList: public list { +public: + inline const Type * data(uint index = 0) const {return &(*this)[index];} + inline Type * data(uint index = 0) {return &(*this)[index];} + inline void fill(const Type & t) {list::assign(list::size(), t);} + inline void remove(uint num) {list::erase(list::begin() + num);} + inline void remove(uint num, uint count) {list::erase(list::begin() + num, list::begin() + num + count);} + inline void insert(uint pos, const Type & t) {list::insert(list::begin() + pos, t);} +}; + template class Flags { private: @@ -88,6 +108,11 @@ public: inline int random() {return rand();} #endif +template inline void piSwap(Type & f, Type & s) {Type t = f; f = s; s = t;} +template inline Type piMin(const Type & f, const Type & s) {return (f > s) ? s : f;} +template inline Type piMax(const Type & f, const Type & s) {return (f < s) ? s : f;} +template inline Type piMin(const Type & f, const Type & s, const Type & t) {return (f < s && f < t) ? f : ((s < t) ? s : t);} +template inline Type piMax(const Type & f, const Type & s, const Type & t) {return (f > s && f > t) ? f : ((s > t) ? s : t);} inline bool atob(const string & str) { return str == "1" ? true : false;}; inline string btos(const bool num) { return num ? "0" : "1";}; inline string itos(const int num) { diff --git a/pimath.cpp b/pimath.cpp new file mode 100644 index 00000000..d83ffaa7 --- /dev/null +++ b/pimath.cpp @@ -0,0 +1,382 @@ +#include "pimath.h" + +/* +* Fast Fourier Transformation +* ==================================================== +* Coded by Miroslav Voinarovsky, 2002 +* This source is freeware. +*/ + +// This array contains values from 0 to 255 with reverse bit order +static uchar reverse256[]= { + 0x00, 0x80, 0x40, 0xC0, 0x20, 0xA0, 0x60, 0xE0, + 0x10, 0x90, 0x50, 0xD0, 0x30, 0xB0, 0x70, 0xF0, + 0x08, 0x88, 0x48, 0xC8, 0x28, 0xA8, 0x68, 0xE8, + 0x18, 0x98, 0x58, 0xD8, 0x38, 0xB8, 0x78, 0xF8, + 0x04, 0x84, 0x44, 0xC4, 0x24, 0xA4, 0x64, 0xE4, + 0x14, 0x94, 0x54, 0xD4, 0x34, 0xB4, 0x74, 0xF4, + 0x0C, 0x8C, 0x4C, 0xCC, 0x2C, 0xAC, 0x6C, 0xEC, + 0x1C, 0x9C, 0x5C, 0xDC, 0x3C, 0xBC, 0x7C, 0xFC, + 0x02, 0x82, 0x42, 0xC2, 0x22, 0xA2, 0x62, 0xE2, + 0x12, 0x92, 0x52, 0xD2, 0x32, 0xB2, 0x72, 0xF2, + 0x0A, 0x8A, 0x4A, 0xCA, 0x2A, 0xAA, 0x6A, 0xEA, + 0x1A, 0x9A, 0x5A, 0xDA, 0x3A, 0xBA, 0x7A, 0xFA, + 0x06, 0x86, 0x46, 0xC6, 0x26, 0xA6, 0x66, 0xE6, + 0x16, 0x96, 0x56, 0xD6, 0x36, 0xB6, 0x76, 0xF6, + 0x0E, 0x8E, 0x4E, 0xCE, 0x2E, 0xAE, 0x6E, 0xEE, + 0x1E, 0x9E, 0x5E, 0xDE, 0x3E, 0xBE, 0x7E, 0xFE, + 0x01, 0x81, 0x41, 0xC1, 0x21, 0xA1, 0x61, 0xE1, + 0x11, 0x91, 0x51, 0xD1, 0x31, 0xB1, 0x71, 0xF1, + 0x09, 0x89, 0x49, 0xC9, 0x29, 0xA9, 0x69, 0xE9, + 0x19, 0x99, 0x59, 0xD9, 0x39, 0xB9, 0x79, 0xF9, + 0x05, 0x85, 0x45, 0xC5, 0x25, 0xA5, 0x65, 0xE5, + 0x15, 0x95, 0x55, 0xD5, 0x35, 0xB5, 0x75, 0xF5, + 0x0D, 0x8D, 0x4D, 0xCD, 0x2D, 0xAD, 0x6D, 0xED, + 0x1D, 0x9D, 0x5D, 0xDD, 0x3D, 0xBD, 0x7D, 0xFD, + 0x03, 0x83, 0x43, 0xC3, 0x23, 0xA3, 0x63, 0xE3, + 0x13, 0x93, 0x53, 0xD3, 0x33, 0xB3, 0x73, 0xF3, + 0x0B, 0x8B, 0x4B, 0xCB, 0x2B, 0xAB, 0x6B, 0xEB, + 0x1B, 0x9B, 0x5B, 0xDB, 0x3B, 0xBB, 0x7B, 0xFB, + 0x07, 0x87, 0x47, 0xC7, 0x27, 0xA7, 0x67, 0xE7, + 0x17, 0x97, 0x57, 0xD7, 0x37, 0xB7, 0x77, 0xF7, + 0x0F, 0x8F, 0x4F, 0xCF, 0x2F, 0xAF, 0x6F, 0xEF, + 0x1F, 0x9F, 0x5F, 0xDF, 0x3F, 0xBF, 0x7F, 0xFF, +}; + +//This is array exp(-2*pi*j/2^n) for n= 1,...,32 +//exp(-2*pi*j/2^n) = complexd( cos(2*pi/2^n), -sin(2*pi/2^n) ) +static complexd W2n[32] = { + complexd(-1.00000000000000000000000000000000, 0.00000000000000000000000000000000), // W2 calculator (copy/paste) : po, ps + complexd( 0.00000000000000000000000000000000, -1.00000000000000000000000000000000), // W4: p/2=o, p/2=s + complexd( 0.70710678118654752440084436210485, -0.70710678118654752440084436210485), // W8: p/4=o, p/4=s + complexd( 0.92387953251128675612818318939679, -0.38268343236508977172845998403040), // p/8=o, p/8=s + complexd( 0.98078528040323044912618223613424, -0.19509032201612826784828486847702), // p/16= + complexd( 0.99518472667219688624483695310948, -9.80171403295606019941955638886e-2), // p/32= + complexd( 0.99879545620517239271477160475910, -4.90676743274180142549549769426e-2), // p/64= + complexd( 0.99969881869620422011576564966617, -2.45412285229122880317345294592e-2), // p/128= + complexd( 0.99992470183914454092164649119638, -1.22715382857199260794082619510e-2), // p/256= + complexd( 0.99998117528260114265699043772857, -6.13588464915447535964023459037e-3), // p/(2y9)= + complexd( 0.99999529380957617151158012570012, -3.06795676296597627014536549091e-3), // p/(2y10)= + complexd( 0.99999882345170190992902571017153, -1.53398018628476561230369715026e-3), // p/(2y11)= + complexd( 0.99999970586288221916022821773877, -7.66990318742704526938568357948e-4), // p/(2y12)= + complexd( 0.99999992646571785114473148070739, -3.83495187571395589072461681181e-4), // p/(2y13)= + complexd( 0.99999998161642929380834691540291, -1.91747597310703307439909561989e-4), // p/(2y14)= + complexd( 0.99999999540410731289097193313961, -9.58737990959773458705172109764e-5), // p/(2y15)= + complexd( 0.99999999885102682756267330779455, -4.79368996030668845490039904946e-5), // p/(2y16)= + complexd( 0.99999999971275670684941397221864, -2.39684498084182187291865771650e-5), // p/(2y17)= + complexd( 0.99999999992818917670977509588385, -1.19842249050697064215215615969e-5), // p/(2y18)= + complexd( 0.99999999998204729417728262414778, -5.99211245264242784287971180889e-6), // p/(2y19)= + complexd( 0.99999999999551182354431058417300, -2.99605622633466075045481280835e-6), // p/(2y20)= + complexd( 0.99999999999887795588607701655175, -1.49802811316901122885427884615e-6), // p/(2y21)= + complexd( 0.99999999999971948897151921479472, -7.49014056584715721130498566730e-7), // p/(2y22)= + complexd( 0.99999999999992987224287980123973, -3.74507028292384123903169179084e-7), // p/(2y23)= + complexd( 0.99999999999998246806071995015625, -1.87253514146195344868824576593e-7), // p/(2y24)= + complexd( 0.99999999999999561701517998752946, -9.36267570730980827990672866808e-8), // p/(2y25)= + complexd( 0.99999999999999890425379499688176, -4.68133785365490926951155181385e-8), // p/(2y26)= + complexd( 0.99999999999999972606344874922040, -2.34066892682745527595054934190e-8), // p/(2y27)= + complexd( 0.99999999999999993151586218730510, -1.17033446341372771812462135032e-8), // p/(2y28)= + complexd( 0.99999999999999998287896554682627, -5.85167231706863869080979010083e-9), // p/(2y29)= + complexd( 0.99999999999999999571974138670657, -2.92583615853431935792823046906e-9), // p/(2y30)= + complexd( 0.99999999999999999892993534667664, -1.46291807926715968052953216186e-9), // p/(2y31)= +}; + +/* + * x: x - array of items + * T: 1 << T = 2 power T - number of items in array + * complement: false - normal (direct) transformation, true - reverse transformation + */ +void fft(complexd * x, int T, bool complement) +{ + uint I, J, Nmax, N, Nd2, k, m, mpNd2, Skew; + uchar *Ic = (uchar*) &I; + uchar *Jc = (uchar*) &J; + complexd S; + complexd * Wstore, * Warray; + complexd WN, W, Temp, *pWN; + + Nmax = 1 << T; + + //first interchanging + for(I = 1; I < Nmax - 1; I++) + { + Jc[0] = reverse256[Ic[3]]; + Jc[1] = reverse256[Ic[2]]; + Jc[2] = reverse256[Ic[1]]; + Jc[3] = reverse256[Ic[0]]; + J >>= (32 - T); + if (I < J) + { + S = x[I]; + x[I] = x[J]; + x[J] = S; + } + } + + //rotation multiplier array allocation + Wstore = new complexd[Nmax / 2]; + Wstore[0] = complexd(1., 0.); + + //main loop + for(N = 2, Nd2 = 1, pWN = W2n, Skew = Nmax >> 1; N <= Nmax; Nd2 = N, N += N, pWN++, Skew >>= 1) + { + //WN = W(1, N) = exp(-2*pi*j/N) + WN= *pWN; + if (complement) + WN = complexd(WN.real(), -WN.imag()); + for(Warray = Wstore, k = 0; k < Nd2; k++, Warray += Skew) + { + if (k & 1) + { + W *= WN; + *Warray = W; + } + else + W = *Warray; + + for(m = k; m < Nmax; m += N) + { + mpNd2 = m + Nd2; + Temp = W; + Temp *= x[mpNd2]; + x[mpNd2] = x[m]; + x[mpNd2] -= Temp; + x[m] += Temp; + } + } + } + + delete[] Wstore; + + if (complement) + { + for( I = 0; I < Nmax; I++ ) + x[I] /= Nmax; + } +} + + +const char Solver::methods_desc[] = "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"; + +Solver::Method Solver::method_global = Solver::AdamsBashfortMoulton_4; + + +void Solver::solve(double u, double h) { + switch (method) { + case Global: + switch (method_global) { + case Eyler1: solveEyler1(u, h); break; + case Eyler2: 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 Eyler1: solveEyler1(u, h); break; + case Eyler2: 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 Solver::fromTF(const TransferFunction & TF) { + if (TF.vector_An.size() >= TF.vector_Bm.size()) + size = TF.vector_An.size()-1; + else { + cout << "Solver error: {A} should be greater than {B}" << endl; + 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]; // Рассчитываем вектор d + for (uint i = 1; i < size + 1; ++i) { + sum = 0.; + for (uint m = 0; m < i; ++m) + sum += a1[i - m] * d[m]; + d[i] = b1[i] - sum; + } + + for (uint i = 0; i < size - 1; ++i) // Заполняем матрицу А + 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 Solver::solveEyler1(double u, double h) { + /*for (uint i = 0; i < size; ++i) { + * sum = 0.; + * for (uint j = 0; j < size; ++j) + * sum += A[j][i] * X[j]; + * xx[i] = sum + d[i] * u; + }*/ + F[0] = A * X + d * u; + X += F[0] * h; + moveF(); +} + + +void Solver::solveEyler2(double u, double h) { + F[0] = A * X + d * u; + X += (F[0] + F[1]) * h / 2.; + moveF(); +} + + +void Solver::solveRK4(double u, double h) { + k1 = A * X + d * u; xx = k1 * h / 2.; XX = X + xx; + k2 = A * (XX + k1 * h / 2.) + d * (u + xx[0]); xx = k2 * h / 2.; XX += xx; + k3 = A * (XX + k2 * h / 2.) + d * (u + xx[0]); xx = k3 * h; XX += xx; + k4 = A * (XX + k3 * h) + d * (u + xx[0]); + //cout << k1 << k2 << k3 << k4 << endl; + X += (k1 + k2 * 2. + k3 * 2. + k4) * h / 6.; +} + + +void Solver::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 Solver::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 Solver::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 Solver::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(); +} diff --git a/pimath.h b/pimath.h index d38587ad..0f7215c5 100644 --- a/pimath.h +++ b/pimath.h @@ -1,10 +1,14 @@ -#ifndef PIMathATH_H -#define PIMathATH_H +#ifndef PIMATH_H +#define PIMATH_H -#include #include #include #include "piincludes.h" +#ifndef QNX +# include +#else +# include +#endif using std::complex; @@ -13,29 +17,35 @@ typedef complex complexf; typedef complex complexd; typedef complex complexld; -const double deg2rad = atan(1) / 45.; -const double rad2deg = 45. / atan(1); +const double deg2rad = atan(1.) / 45.; +const double rad2deg = 45. / atan(1.); inline int pow2(int p) {return (int)1 << p;} -template inline int pimin(const Type & f, const Type & s) {return (f > s) ? s : f;} -template inline int pimax(const Type & f, const Type & s) {return (f < s) ? s : f;} +inline double sqr(const double & v) {return v * v;} -/// Vector +template +class PIMathMatrixT; + +/// Vector templated #define PIMV_FOR(v, s) for (uint v = s; v < Size; ++v) template -class PIMathVector { - typedef PIMathVector _CVector; +class PIMathVectorT { + typedef PIMathVectorT _CVector; public: - PIMathVector() {resize(Size);} - PIMathVector(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() {resize(Size);} + 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 & val) {resize(Size); PIMV_FOR(i, 0) c[i] = val[i];} + PIMathVectorT(const _CVector & st, const _CVector & fn) {resize(Size); PIMV_FOR(i, 0) c[i] = fn[i] - st[i];} + 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 & move(const Type & v) {PIMV_FOR(i, 0) c[i] += v; return *this;} - inline Type lengthSqr() const {Type tv(0); PIMV_FOR(i, 0) tv += c[i] * c[i]; return tv;} + 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);} @@ -43,10 +53,11 @@ public: 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 void normalize() {Type tv = length(); if (tv == Type(1)) return; PIMV_FOR(i, 0) c[i] /= 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];} @@ -57,31 +68,46 @@ public: 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 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 /// vector {Size, Type} to vector {Size1, Type1} + inline PIMathVectorT turnTo() {PIMathVectorT tv; uint sz = piMin(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 c; - + }; template -inline std::ostream & operator <<(std::ostream & s, const PIMathVector & v) {s << '{'; PIMV_FOR(i, 0) {s << v[i]; if (i < Size - 1) s << ", ";} s << '}'; return s;} +inline std::ostream & operator <<(std::ostream & s, const PIMathVectorT & v) {s << '{'; PIMV_FOR(i, 0) {s << v[i]; if (i < Size - 1) s << ", ";} s << '}'; return s;} -typedef PIMathVector<2u, int> PIMathVector2i; -typedef PIMathVector<3u, int> PIMathVector3i; -typedef PIMathVector<4u, int> PIMathVector4i; -typedef PIMathVector<2u, double> PIMathVector2d; -typedef PIMathVector<3u, double> PIMathVector3d; -typedef PIMathVector<4u, double> PIMathVector4d; +//template /// vector {Size0, Type0} to vector {Size1, Type1} +//inline operator PIMathVectorT(const PIMathVectorT & v) {PIMathVectorT tv; uint sz = piMin(Size0, Size1); for (uint i = 0; i < sz; ++i) tv[i] = v[i]; return tv;} -/// Matrix +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 @@ -91,28 +117,33 @@ typedef PIMathVector<4u, double> PIMathVector4d; #define PIMM_FOR_R(v) for (uint v = 0; v < Rows; ++v) template -class PIMathMatrix { - typedef PIMathMatrix _CMatrix; - typedef PIMathMatrix _CMatrixI; - typedef PIMathVector _CMCol; - typedef PIMathVector _CMRow; +class PIMathMatrixT { + typedef PIMathMatrixT _CMatrix; + typedef PIMathMatrixT _CMatrixI; + typedef PIMathVectorT _CMCol; + typedef PIMathVectorT _CMRow; public: - PIMathMatrix() {resize(Cols, Rows);} - PIMathMatrix(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() {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 & 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 & operator [](uint col) {return m[col];} @@ -129,25 +160,534 @@ public: 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 > m; - + }; template -inline std::ostream & operator <<(std::ostream & s, const PIMathMatrix & 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;} +inline std::ostream & operator <<(std::ostream & s, const PIMathMatrixT & 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;} -typedef PIMathMatrix<2u, 2u, int> PIMathMatrix22i; -typedef PIMathMatrix<3u, 3u, int> PIMathMatrix33i; -typedef PIMathMatrix<4u, 4u, int> PIMathMatrix44i; -typedef PIMathMatrix<2u, 2u, double> PIMathMatrix22d; -typedef PIMathMatrix<3u, 3u, double> PIMathMatrix33d; -typedef PIMathMatrix<4u, 4u, double> PIMathMatrix44d; +/// Multiply matrices {CR x Rows0} on {Cols1 x CR}, result is {Cols1 x Rows0} +template +inline PIMathMatrixT operator *(const PIMathMatrixT & fm, + const PIMathMatrixT & sm) { + PIMathMatrixT 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 +inline PIMathVectorT operator *(const PIMathMatrixT & fm, + const PIMathVectorT & sv) { + PIMathVectorT 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 +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 -#endif // PIMathATH_H +/// Vector + +#define PIMV_FOR(v, s) for (uint v = s; v < size_; ++v) + +template +class PIMathVector { + typedef PIMathVector _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 & 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(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 + inline PIMathVector turnTo(uint size) {PIMathVector tv; uint sz = piMin(size_, size); for (uint i = 0; i < sz; ++i) tv[i] = c[i]; return tv;} + +private: + uint size_; + PIVector c; + +}; + +template +inline std::ostream & operator <<(std::ostream & s, const PIMathVector & v) {s << '{'; for (uint i = 0; i < v.size(); ++i) {s << v[i]; if (i < v.size() - 1) s << ", ";} s << '}'; return s;} + +typedef PIMathVector PIMathVectori; +typedef PIMathVector 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 +class PIMathMatrix { + typedef PIMathMatrix _CMatrix; + typedef PIMathVector _CMCol; + typedef PIMathVector _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 & 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 & operator [](uint col) {return m[col];} + inline PIVector 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 > m; + +}; + +template +inline std::ostream & operator <<(std::ostream & s, const PIMathMatrix & 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 +inline PIMathMatrix operator *(const PIMathMatrix & fm, + const PIMathMatrix & sm) { + uint cr = fm.cols(), rows0 = fm.rows(), cols1 = sm.cols(); + PIMathMatrix 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 +inline PIMathVector operator *(const PIMathMatrix & fm, + const PIMathVector & sv) { + uint c = fm.cols(), r = fm.rows(); + PIMathVector 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 PIMathMatrixi; +typedef PIMathMatrix 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 vector_Bm, vector_An; +}; + +// Класс, служащий для перевода передаточной функции в систему ОДУ первого порядка +// реализованы след. методы решения дифф. ур-ний: +// Эйлера +// Рунге-Кутта 4-го порядка +// Адамса-Бэшфортса-Моултона 2, 3, 4 порядков +class Solver +{ +public: + enum Method {Global = -1, + Eyler1 = 01, + Eyler2 = 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 F; + PIVector times; + uint size, step; + Method method; + double sum, td, ct, lp, dh, t, x1, x0; + bool ok; + +}; + + +#endif // PIMATH_H diff --git a/pip.h b/pip.h index beeb7f18..da53d67c 100644 --- a/pip.h +++ b/pip.h @@ -1,6 +1,6 @@ #include "pitimer.h" #include "pivariable.h" #include "piconsole.h" -#include "pigeometry.h" #include "pimath.h" +#include "pidir.h" #include "piprotocol.h" diff --git a/pistring.cpp b/pistring.cpp index c468b01c..e43cf32a 100644 --- a/pistring.cpp +++ b/pistring.cpp @@ -22,6 +22,13 @@ PIString & PIString::operator +=(const PIString & str) { } +PIString & PIString::operator +=(const PIByteArray & ba) { + uint l = ba.size(); + for (uint i = 0; i < l; ++i) push_back(ba[i]); + return *this; +} + + bool PIString::operator ==(const PIString & str) const { uint l = str.size(); if (size() != l) return false; @@ -106,6 +113,20 @@ PIString & PIString::trim() { } +PIStringList PIString::split(const PIString & delim) { + PIString ts(*this); + PIStringList sl; + int ci = ts.find(delim); + while (ci >= 0) { + sl << ts.left(ci); + ts.cutLeft(ci + delim.length()); + ci = ts.find(delim); + } + if (ts.length() > 0) sl << ts; + return sl; +} + + int PIString::find(const char str, const int start) { for (int i = start; i < length(); ++i) if (at(i) == str) diff --git a/pistring.h b/pistring.h index 90a09dd3..87828d21 100644 --- a/pistring.h +++ b/pistring.h @@ -1,7 +1,9 @@ #ifndef PISTRING_H #define PISTRING_H -#include "piincludes.h" +#include "pibytearray.h" + +class PIStringList; class PIString: public PIVector { @@ -11,6 +13,7 @@ public: PIString(const char * str) {*this += str;} PIString(const string & str) {*this += str;} PIString(const PIString & str) {*this += str;} + PIString(const PIByteArray & ba) {*this += ba;} PIString(const int len, const char c = ' ') {for (int i = 0; i < len; ++i) push_back(c);} operator const char*() {return data();} @@ -19,21 +22,25 @@ public: inline const char operator [](const int pos) const {return at(pos);} PIString & operator +=(const PIString & str); + PIString & operator +=(const PIByteArray & ba); inline PIString & operator +=(const char & c) {push_back(c); return *this;} PIString & operator +=(const char * str); PIString & operator +=(const string & str); bool operator ==(const PIString & str) const; + inline bool operator ==(const PIByteArray & ba) const {return *this == PIString(ba);} inline bool operator ==(const char & c) const {return *this == PIString(c);} inline bool operator ==(const char * str) const {return *this == PIString(str);} inline bool operator ==(const string & str) const {return *this == PIString(str);} bool operator !=(const PIString & str) const; + inline bool operator !=(const PIByteArray & ba) const {return *this != PIString(ba);} inline bool operator !=(const char & c) const {return *this != PIString(c);} inline bool operator !=(const char * str) const {return *this != PIString(str);} inline bool operator !=(const string & str) const {return *this != PIString(str);} - PIString & operator <<(const PIString & str) {*this += str; return *this;} + inline PIString & operator <<(const PIString & str) {*this += str; return *this;} + inline PIString & operator <<(const PIByteArray & ba) {*this += ba; return *this;} inline PIString & operator <<(const char & c) {*this += c; return *this;} inline PIString & operator <<(const char * str) {*this += str; return *this;} inline PIString & operator <<(const string & str) {*this += str; return *this;} @@ -53,6 +60,7 @@ public: PIString & trim(); const char * data() const {return stdString().c_str();} string stdString() const {return (size() == 0) ? string() : string(&at(0), size());} + PIStringList split(const PIString & delim); PIString toUpperCase() const; PIString toLowerCase() const; @@ -74,6 +82,7 @@ public: long toLong() const {return atol(data());} float toFloat() const {return (float)atof(data());} double toDouble() const {return atof(data());} + inline PIByteArray toBase64() {return PIByteArray(data(), size()).toBase64();} inline PIString & setNumber(const int value) {clear(); *this += itos(value); return *this;} inline PIString & setNumber(const short value) {clear(); *this += itos(value); return *this;} @@ -82,6 +91,7 @@ public: inline PIString & setNumber(const double value) {clear(); *this += dtos(value); return *this;} inline static PIString fromNumber(const int value) {return PIString(itos(value));} + inline static PIString fromNumber(const uint value) {return PIString(itos(value));} inline static PIString fromNumber(const short value) {return PIString(itos(value));} inline static PIString fromNumber(const long value) {return PIString(ltos(value));} inline static PIString fromNumber(const float value) {return PIString(ftos(value));} @@ -91,10 +101,12 @@ public: inline PIString operator +(const PIString & str, const PIString & f) {PIString s(str); s += f; return s;} +inline PIString operator +(const PIString & f, const PIByteArray & ba) {PIString s(f); s += ba; return s;} inline PIString operator +(const PIString & f, const char & c) {PIString s(f); s.push_back(c); return s;} inline PIString operator +(const PIString & f, const char * str) {PIString s(f); s += str; return s;} inline PIString operator +(const PIString & f, const string & str) {PIString s(f); s += str; return s;} +inline PIString operator +(const PIByteArray & ba, const PIString & f) {return PIString(ba) + f;} inline PIString operator +(const char & c, const PIString & f) {return PIString(c) + f;} inline PIString operator +(const char * str, const PIString & f) {return PIString(str) + f;} inline PIString operator +(const string & str, const PIString & f) {return PIString(str) + f;} @@ -102,4 +114,25 @@ inline PIString operator +(const string & str, const PIString & f) {return PIStr char chrUpr(char c); char chrLwr(char c); +class PIStringList: public PIVector +{ +public: + PIStringList() {;} + PIStringList(const PIString & str) {push_back(str);} + + inline PIString join(const PIString & delim) const {PIString s; for (uint i = 0; i < size(); ++i) {s += at(i); if (i < size() - 1) s += delim;} return s;} + + inline PIStringList & operator <<(const PIString & str) {push_back(str); return *this;} + inline PIStringList & operator <<(const PIByteArray & ba) {push_back(ba); return *this;} + inline PIStringList & operator <<(const char & c) {push_back(PIString(c)); return *this;} + inline PIStringList & operator <<(const char * str) {push_back(PIString(str)); return *this;} + inline PIStringList & operator <<(const string & str) {push_back(str); return *this;} + inline PIStringList & operator <<(const int & num) {push_back(PIString::fromNumber(num)); return *this;} + inline PIStringList & operator <<(const short & num) {push_back(PIString::fromNumber(num)); return *this;} + inline PIStringList & operator <<(const long & num) {push_back(PIString::fromNumber(num)); return *this;} + inline PIStringList & operator <<(const float & num) {push_back(PIString::fromNumber(num)); return *this;} + inline PIStringList & operator <<(const double & num) {push_back(PIString::fromNumber(num)); return *this;} + +}; + #endif // PISTRING_H