#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;
}
