diff --git a/Code/precision b/Code/precision deleted file mode 100755 index 625d64d08b3acadee85323b8dcea62712f0ce3f3..0000000000000000000000000000000000000000 Binary files a/Code/precision and /dev/null differ diff --git a/Code/precision.cpp b/Code/precision.cpp new file mode 100644 index 0000000000000000000000000000000000000000..f9bca1fb07470c2842ddb59d42dc3156851d1824 --- /dev/null +++ b/Code/precision.cpp @@ -0,0 +1,137 @@ +#include <boost/numeric/interval.hpp> +#include <boost/multiprecision/cpp_dec_float.hpp> +#include <boost/multiprecision/cpp_int.hpp> +#include <eigen3/Eigen/Dense> +#include <random> +#include <algorithm> +#include <iostream> + +//types +using interval = boost::numeric::interval<double> ; +using rational = boost::multiprecision::cpp_rational ; +using MatDX = Eigen::MatrixXd ; +using VecDX = Eigen::VectorXd ; +using MatIX = Eigen::Matrix<interval, Eigen::Dynamic, Eigen::Dynamic> ; +using VecIX = Eigen::Matrix<interval, Eigen::Dynamic, 1> ; +using MatRX = Eigen::Matrix<rational, Eigen::Dynamic, Eigen::Dynamic> ; +using VecRX = Eigen::Matrix<rational, Eigen::Dynamic, 1> ; + +namespace Eigen { + namespace internal { + template<typename X, typename S, typename P> + struct is_convertible<X,boost::numeric::interval<S,P> > { + enum { value = is_convertible<X,S>::value }; + }; + + template<typename S, typename P1, typename P2> + struct is_convertible<boost::numeric::interval<S,P1>,boost::numeric::interval<S,P2> > { + enum { value = true }; + }; + } +} + +//command line +char* get_cmd_option( char ** begin, char ** end, + const std::string & option + ) { + char ** itr = std::find(begin, end, "--" + option); + if (itr != end && ++itr != end) + { + return *itr; + } + return 0; +} + +bool cmd_option_exists( char** begin, char** end, + const std::string& option + ) { + return std::find(begin, end, "--" + option) != end; +} + +template<typename T> +bool get_cmd_value( + char ** begin, + char** end, + const std::string& option, + T* output + ) { + if(cmd_option_exists(begin, end, option)) { + std::stringstream ss ; + ss << get_cmd_option(begin, end, option) ; + T t ; + ss >> t ; + if(!ss.fail()) { + *output = t ; + return true ; + } + } + return false ; +} + +int main(int argc, char** argv) { + //test variables + unsigned int dim = 10 ; + get_cmd_value(argv, argv + argc, "dim", &dim) ; + + //alea + std::random_device seed ; + std::mt19937_64 mt(seed()) ; + std::uniform_real_distribution<double> random_double(-1, 1) ; + + //initialization + MatDX matd(dim, dim) ; + VecDX vecd(dim) ; + MatIX mati(dim, dim) ; + VecIX veci(dim) ; + MatRX matr(dim, dim) ; + VecRX vecr(dim) ; + + for(unsigned int i = 0; i < dim; ++i) { + for(unsigned int j = 0; j < dim; ++j) { + double d = random_double(mt) ; + matd(i, j) = d ; + mati(i, j) = d ; + matr(i, j) = d ; + } + double d = random_double(mt) ; + vecd(i) = d ; + veci(i) = d ; + vecr(i) = d ; + } + + rational detr = matr.fullPivLu().determinant() ; + + std::cout << "==========" << std::endl ; + { + std::cout << "solve fullPivLu" << std::endl ; + double detd = matd.fullPivLu().determinant() ; + try { + interval deti = mati.fullPivLu().determinant() ; + std::cout << detd + << " in [" << deti.lower() << " -- " << deti.upper() + << "] error " << detd - (double) detr + << " wrt. " << deti.upper() - deti.lower() << std::endl ; + } catch(std::exception& e) { + std::cout << "interval computation failed" << std::endl ; + std::cout << detd << " error " << detd - (double) detr << std::endl ; + } + std::cout << "==========" << std::endl ; + } + { + std::cout << "solve partialPivLu" << std::endl ; + double detd = matd.partialPivLu().determinant() ; + try { + interval deti = mati.partialPivLu().determinant() ; + std::cout << detd + << " in [" << deti.lower() << " -- " << deti.upper() + << "] error " << detd - (double) detr + << " wrt. " << deti.upper() - deti.lower() << std::endl ; + } catch(std::exception& e) { + std::cout << "interval computation failed" << std::endl ; + std::cout << detd << " error " << detd - (double) detr << std::endl ; + } + std::cout << "==========" << std::endl ; + } + + return 0 ; +}