mirror of
https://github.com/boostorg/math.git
synced 2025-05-11 21:33:52 +00:00
126 lines
3.5 KiB
C++
126 lines
3.5 KiB
C++
/*
|
|
* Copyright Nick Thompson, 2021
|
|
* Use, modification and distribution are subject to the
|
|
* Boost Software License, Version 1.0. (See accompanying file
|
|
* LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
|
*/
|
|
|
|
#include "math_unit_test.hpp"
|
|
#include <numeric>
|
|
#include <random>
|
|
#include <array>
|
|
#include <boost/core/demangle.hpp>
|
|
#include <boost/math/interpolators/bilinear_uniform.hpp>
|
|
#ifdef BOOST_HAS_FLOAT128
|
|
#include <boost/multiprecision/float128.hpp>
|
|
using boost::multiprecision::float128;
|
|
#endif
|
|
|
|
#if __has_include(<stdfloat>)
|
|
# include <stdfloat>
|
|
#endif
|
|
|
|
using boost::math::interpolators::bilinear_uniform;
|
|
|
|
template<class Real>
|
|
void test_four_values()
|
|
{
|
|
Real x0 = 0;
|
|
Real y0 = 0;
|
|
Real dx = 1;
|
|
Real dy = 1;
|
|
Real value = Real(1.5);
|
|
std::vector<Real> v(2*2, value);
|
|
auto v_copy = v;
|
|
auto ub = bilinear_uniform<decltype(v)>(std::move(v_copy), 2, 2, dx, dy, x0, y0);
|
|
for (Real x = x0; x <= x0 + dx; x += dx/8) {
|
|
for (Real y = y0; y <= y0 + dx; y += dy/8) {
|
|
CHECK_ULP_CLOSE(value, ub(x, y), 1);
|
|
}
|
|
}
|
|
|
|
// Now we test the unit square:
|
|
std::random_device rd;
|
|
std::uniform_real_distribution<Real> dis(1,2);
|
|
|
|
int i = 0;
|
|
while (i++ < 300) {
|
|
v[0] = dis(rd);
|
|
v[1] = dis(rd);
|
|
v[2] = dis(rd);
|
|
v[3] = dis(rd);
|
|
|
|
// See https://en.wikipedia.org/wiki/Bilinear_interpolation, section: Unit square
|
|
auto f = [&v](Real x, Real y) {
|
|
return v[0]*(1-x)*(1-y) + v[1]*x*(1-y) + v[2]*(1-x)*y + v[3]*x*y;
|
|
};
|
|
|
|
v_copy = v;
|
|
ub = bilinear_uniform<decltype(v_copy)>(std::move(v_copy), 2, 2, dx, dy, x0, y0);
|
|
for (Real x = x0; x <= x0 + dx; x += dx/16) {
|
|
for (Real y = y0; y <= y0 + dx; y += dy/16) {
|
|
CHECK_ULP_CLOSE(f(x,y), ub(x, y), 3);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void test_linear()
|
|
{
|
|
std::random_device rd;
|
|
std::uniform_real_distribution<Real> dis(1,2);
|
|
std::array<Real, 4> a{dis(rd), dis(rd), dis(rd), dis(rd)};
|
|
auto f = [&a](Real x, Real y) {
|
|
return a[0] + a[1]*x + a[2]*y + a[3]*x*y;
|
|
};
|
|
|
|
for (int rows = 2; rows < 20; ++rows) {
|
|
for (int cols = 2; cols < 20; ++cols) {
|
|
Real dx = dis(rd);
|
|
Real dy = dis(rd);
|
|
Real x0 = dis(rd);
|
|
Real y0 = dis(rd);
|
|
std::vector<Real> v(rows*cols, std::numeric_limits<Real>::quiet_NaN());
|
|
for (int i = 0; i < cols; ++i) {
|
|
for (int j = 0; j < rows; ++j) {
|
|
v[j*cols + i] = f(x0 + i*dx, y0 + j*dy);
|
|
}
|
|
}
|
|
auto ub = bilinear_uniform<decltype(v)>(std::move(v), rows, cols, dx, dy, x0, y0);
|
|
|
|
for (Real x = x0; x < x0 + (cols-1)*dx; x += dx/8) {
|
|
for (Real y = y0; y < y0 + (rows-1)*dy; y += dy/8) {
|
|
if (!CHECK_ULP_CLOSE(f(x,y), ub(x, y), 13)) {
|
|
std::cerr << " f(" << x << ", " << y << ") = " << f(x,y) << "\n";
|
|
std::cerr << "ub(" << x << ", " << y << ") = " << ub(x,y) << "\n";
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
int main()
|
|
{
|
|
#ifdef __STDCPP_FLOAT32_T__
|
|
test_four_values<std::float32_t>();
|
|
#else
|
|
test_four_values<float>();
|
|
#endif
|
|
|
|
#ifdef __STDCPP_FLOAT64_T__
|
|
test_four_values<std::float64_t>();
|
|
test_linear<std::float64_t>();
|
|
#else
|
|
test_four_values<double>();
|
|
test_linear<double>();
|
|
#endif
|
|
|
|
test_four_values<long double>();
|
|
|
|
return boost::math::test::report_errors();
|
|
}
|