online compiler and debugger for c/c++

#include <boost/numeric/ublas/vector.hpp> #include <boost/numeric/ublas/matrix.hpp> #include <boost/numeric/ublas/lu.hpp> #include <boost/numeric/ublas/io.hpp> #include <cmath> #include <vector> #include <functional> #include <iostream> #include <chrono> template <typename F> class SuperLinearInterpolation { private: F m_f; std::vector<double> m_xs; std::vector<double> m_ys; double m_delta_x; public: SuperLinearInterpolation(F f, std::function<double(double)> integral_f, std::function<double(double)> integral_xf, double x0, double xN, int N) : m_f(f), m_xs(N+1), m_ys(N+1), m_delta_x((xN - x0) / N) { double y0 = f(x0); double yN = f(xN); for (int i = 0; i <= N; ++i) { m_xs[i] = x0 + i * m_delta_x; } std::vector<double> F0(N); std::vector<double> F1(N); for (int i = 0; i < N; ++i) { F0[i] = integral_f(m_xs[i+1]) - integral_f(m_xs[i]); F1[i] = integral_xf(m_xs[i+1]) - integral_xf(m_xs[i]); } boost::numeric::ublas::matrix<double> M = boost::numeric::ublas::zero_matrix<double>(N-1, N-1); boost::numeric::ublas::vector<double> C = boost::numeric::ublas::zero_vector<double>(N-1); for (int i = 0; i < N-1; ++i) { C(i) = -F1[i] + m_xs[i]*F0[i] + F1[i+1] - m_xs[i+2]*F0[i+1]; M(i,i) = -2*m_delta_x*m_delta_x/3; if (i == 0) { C(i) += y0 * m_delta_x*m_delta_x/6; M(i,i+1) = -m_delta_x*m_delta_x/6; } else if (i == N-2) { C(i) += yN * m_delta_x*m_delta_x/6; M(i,i-1) = -m_delta_x*m_delta_x/6; } else { M(i,i-1) = -m_delta_x*m_delta_x/6; M(i,i+1) = -m_delta_x*m_delta_x/6; } } // create a working copy of the input boost::numeric::ublas::matrix<double> Mcopy(M); // create a permutation matrix for the LU-factorization boost::numeric::ublas::permutation_matrix<std::size_t> pm(M.size1()); // perform LU-factorization int res = boost::numeric::ublas::lu_factorize(Mcopy, pm); if (res != 0) throw std::runtime_error("Couldn't LU factorize matrix in SuperLinearInterpolation"); // create identity matrix of "inverse" boost::numeric::ublas::matrix<double> inverse = boost::numeric::ublas::identity_matrix<double>(Mcopy.size1()); // backsubstitute to get the inverse boost::numeric::ublas::lu_substitute(Mcopy, pm, inverse); boost::numeric::ublas::vector<double> y = boost::numeric::ublas::prod(inverse, C); m_ys[0] = y0; for (int i = 0; i < N-1; ++i) { m_ys[i+1] = y[i]; } m_ys[N] = yN; } const std::vector<double>& xs() const { return m_xs; } const std::vector<double>& ys() const { return m_ys; } double operator()(double x) const { if (x < m_xs.front() || x >= m_xs.back()) return m_f(x); int i = (int)((x - m_xs.front()) / m_delta_x); return m_ys[i] + (x - m_xs[i]) * (m_ys[i+1] - m_ys[i]) / m_delta_x; } }; SuperLinearInterpolation<std::function<double(double)>> super_sin( [](double x) -> double { return std::sin(x); }, [](double x) -> double { return -std::cos(x); }, [](double x) -> double { return std::sin(x) - x*std::cos(x); }, 0, 2*M_PI, 90 ); double fast_sin(double theta) { while (theta < 0) theta += 2*M_PI; while (theta >= 2*M_PI) theta -= 2*M_PI; return super_sin(theta); } int main() { int N = 10000000; double sum = 0; auto start = std::chrono::system_clock::now(); for (double x = -M_PI; x <= 3*M_PI; x += (4*M_PI)/N) { sum += std::sin(x); } auto end = std::chrono::system_clock::now(); auto diff = std::chrono::duration_cast<std::chrono::nanoseconds>(end - start); std::cout << "std::sin, sum " << sum << ", time per call (ns): " << 1.0*diff.count()/N << std::endl; sum = 0; start = std::chrono::system_clock::now(); for (double x = -M_PI; x <= 3*M_PI; x += (4*M_PI)/N) { sum += fast_sin(x); } end = std::chrono::system_clock::now(); diff = std::chrono::duration_cast<std::chrono::nanoseconds>(end - start); std::cout << "fast_sin, sum " << sum << ", time per call (ns): " << 1.0*diff.count()/N << std::endl; double max_error = 0; double error; for (double x = 0; x <= 2*M_PI; x += 2*M_PI/100) { error = std::abs(fast_sin(x) - std::sin(x)); max_error = std::max(error, max_error); } std::cout << "max absolute error: " << max_error << std::endl; return 0; }