From c03a14181f958c906b46513ecd0a9c94aa060a6b Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Fri, 7 Jul 2023 15:55:58 -0700 Subject: [PATCH 1/5] Fix GaussChebyshev1 and add unit test --- .../integratorxx/quadratures/gausscheby1.hpp | 10 ++++--- test/1d_quadratures.cxx | 26 ++++++++++++++++++- 2 files changed, 31 insertions(+), 5 deletions(-) diff --git a/include/integratorxx/quadratures/gausscheby1.hpp b/include/integratorxx/quadratures/gausscheby1.hpp index c9b28b7..0054b55 100644 --- a/include/integratorxx/quadratures/gausscheby1.hpp +++ b/include/integratorxx/quadratures/gausscheby1.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include namespace IntegratorXX { @@ -64,16 +65,17 @@ struct quadrature_traits> { weight_container weights(npts); const weight_type pi_ov_npts = M_PI / npts; - const weight_type two_npts_x_pi = 2 * npts * M_PI; + const weight_type pi_ov_2npts = pi_ov_npts / 2; for (size_t i = 0; i < npts; ++i) { // The standard nodes and weights are given by - points[i] = std::cos((2.0 * (i + 1) - 1.) / two_npts_x_pi); - weights[i] = pi_ov_npts; + points[i] = -std::cos((2.0 * (i + 1) - 1.) * pi_ov_2npts); + weights[i] = pi_ov_npts; // However, as we're integrating f(x) not \frac{f(x)}{\sqrt{1 - // x^2}}, we must factor the \sqrt{1-x^2} into the weight - weights[i] *= std::sqrt(1. - std::pow(points[i],2)); + // sqrt(1 - cos^2(x)) = abs(sin(x)) + weights[i] *= std::sin((2.0 * (i + 1) - 1.) * pi_ov_2npts); //std::sqrt(1. - std::pow(points[i],2)); } return std::make_tuple(points, weights); diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 310ffbe..d0e7559 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -141,8 +141,9 @@ TEST_CASE( "Gauss-Legendre Quadratures", "[1d-quad]" ) { auto f = [=]( double x ){ return gaussian(x); }; double res = 0.; - for( auto i = 0; i < quad.npts(); ++i ) + for( auto i = 0; i < quad.npts(); ++i ) { res += wgt[i] * f(pts[i]); + } CHECK( res == Catch::Approx(ref_gaussian_int(-1.,1.)) ); } @@ -167,6 +168,29 @@ TEST_CASE( "Gauss-Legendre Quadratures", "[1d-quad]" ) { } +TEST_CASE( "Gauss-Chebyshev Quadratures", "[1d-quad]") { + + constexpr unsigned order = 200; + auto integrate = [&](auto& quad) { + const auto& pts = quad.points(); + const auto& wgt = quad.weights(); + + auto f = [=]( double x ){ return gaussian(x); }; + + double res = 0.; + for( auto i = 0; i < quad.npts(); ++i ) { + res += wgt[i] * f(pts[i]); + } + + CHECK( res == Catch::Approx(ref_gaussian_int(-1.,1.)) ); + }; + + SECTION("First Kind") { + IntegratorXX::GaussChebyshev1 quad(order, -1., 1.); + integrate(quad); + } +} + TEST_CASE( "Euler-Maclaurin Quadratures", "[1d-quad]" ) { IntegratorXX::MurrayHandyLaming quad( 150 ); From 0e44a86a56cf73f9bcffec8552b1dd47b6f15b97 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Fri, 7 Jul 2023 16:10:01 -0700 Subject: [PATCH 2/5] Cleanup and UT for GaussChebyshev2 --- include/integratorxx/quadratures/gausscheby1.hpp | 5 +++-- include/integratorxx/quadratures/gausscheby2.hpp | 12 ++++++++---- test/1d_quadratures.cxx | 4 ++++ 3 files changed, 15 insertions(+), 6 deletions(-) diff --git a/include/integratorxx/quadratures/gausscheby1.hpp b/include/integratorxx/quadratures/gausscheby1.hpp index 0054b55..1de1ccc 100644 --- a/include/integratorxx/quadratures/gausscheby1.hpp +++ b/include/integratorxx/quadratures/gausscheby1.hpp @@ -68,14 +68,15 @@ struct quadrature_traits> { const weight_type pi_ov_2npts = pi_ov_npts / 2; for (size_t i = 0; i < npts; ++i) { + const auto ti = (2.0 * (i + 1) - 1.) * pi_ov_2npts; // The standard nodes and weights are given by - points[i] = -std::cos((2.0 * (i + 1) - 1.) * pi_ov_2npts); + points[i] = -std::cos(ti); weights[i] = pi_ov_npts; // However, as we're integrating f(x) not \frac{f(x)}{\sqrt{1 - // x^2}}, we must factor the \sqrt{1-x^2} into the weight // sqrt(1 - cos^2(x)) = abs(sin(x)) - weights[i] *= std::sin((2.0 * (i + 1) - 1.) * pi_ov_2npts); //std::sqrt(1. - std::pow(points[i],2)); + weights[i] *= std::sin(ti); //std::sqrt(1. - std::pow(points[i],2)); } return std::make_tuple(points, weights); diff --git a/include/integratorxx/quadratures/gausscheby2.hpp b/include/integratorxx/quadratures/gausscheby2.hpp index aaefbe2..0f76a48 100644 --- a/include/integratorxx/quadratures/gausscheby2.hpp +++ b/include/integratorxx/quadratures/gausscheby2.hpp @@ -65,14 +65,18 @@ struct quadrature_traits> { weight_container weights(npts); point_container points(npts); for (size_t i = 0; i < npts; ++i) { + const auto ti = (i + 1) * pi_ov_npts_p_1; // The standard nodes and weights are given by - points[i] = std::cos((i + 1) * pi_ov_npts_p_1); - weights[i] = - pi_ov_npts_p_1 * std::pow(std::sin((i + 1) * pi_ov_npts_p_1), 2); + points[i] = std::cos(ti); + // Avoid float point cancellations - form modified weight directly in following statements + //weights[i] = + // pi_ov_npts_p_1 * std::pow(std::sin((i + 1) * pi_ov_npts_p_1), 2); // However, since we want the rule with a unit weight factor, we // divide the weights by sqrt(1-x^2). - weights[i] /= std::sqrt(1.0 - std::pow(points[i], 2)); + // PI / (n+1) * sin^2(x) / sqrt(1 - cos^2(x)) = PI / (n+1) * sin(x) + //weights[i] /= std::sqrt(1.0 - std::pow(points[i], 2)); + weights[i] = pi_ov_npts_p_1 * std::sin(ti); } return std::make_tuple(points, weights); diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index d0e7559..392f59c 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -189,6 +189,10 @@ TEST_CASE( "Gauss-Chebyshev Quadratures", "[1d-quad]") { IntegratorXX::GaussChebyshev1 quad(order, -1., 1.); integrate(quad); } + SECTION("Second Kind") { + IntegratorXX::GaussChebyshev2 quad(order, -1., 1.); + integrate(quad); + } } TEST_CASE( "Euler-Maclaurin Quadratures", "[1d-quad]" ) { From 1e0aa8937bb7583e244714db8469cd11cb675075 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Fri, 7 Jul 2023 16:32:43 -0700 Subject: [PATCH 3/5] Fix GaussChebyshev3 and add UT --- include/integratorxx/quadratures/gausscheby2.hpp | 1 + include/integratorxx/quadratures/gausscheby3.hpp | 8 ++++---- test/1d_quadratures.cxx | 4 ++++ 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/include/integratorxx/quadratures/gausscheby2.hpp b/include/integratorxx/quadratures/gausscheby2.hpp index 0f76a48..3c8bac4 100644 --- a/include/integratorxx/quadratures/gausscheby2.hpp +++ b/include/integratorxx/quadratures/gausscheby2.hpp @@ -68,6 +68,7 @@ struct quadrature_traits> { const auto ti = (i + 1) * pi_ov_npts_p_1; // The standard nodes and weights are given by points[i] = std::cos(ti); + // Avoid float point cancellations - form modified weight directly in following statements //weights[i] = // pi_ov_npts_p_1 * std::pow(std::sin((i + 1) * pi_ov_npts_p_1), 2); diff --git a/include/integratorxx/quadratures/gausscheby3.hpp b/include/integratorxx/quadratures/gausscheby3.hpp index 49137c7..3e7a65e 100644 --- a/include/integratorxx/quadratures/gausscheby3.hpp +++ b/include/integratorxx/quadratures/gausscheby3.hpp @@ -62,14 +62,14 @@ struct quadrature_traits> { weight_container weights(npts); point_container points(npts); - for (size_t i = 1; i <= npts; ++i) { + for (size_t i = 0; i < npts; ++i) { // The standard nodes and weights are given by - points[i] = std::cos(0.5 * (2*i - 1) * pi_ov_2n_p_1); + points[i] = std::pow(std::cos(0.5 * (2*(i+1) - 1) * pi_ov_2n_p_1),2); weights[i] = 2.0*pi_ov_2n_p_1 * points[i]; // However, since we want the rule with a unit weight factor, we - // divide the weights by sqrt((1-x)/x). - weights[i] /= std::sqrt((1.0 - points[i])/points[i]); + // divide the weights by sqrt(x/(1-x)). + weights[i] *= std::sqrt((1.0 - points[i])/points[i]); // Finally, convert from [0,1] to [-1,1] points[i] = 2*points[i] - 1.0; diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 392f59c..cb77178 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -193,6 +193,10 @@ TEST_CASE( "Gauss-Chebyshev Quadratures", "[1d-quad]") { IntegratorXX::GaussChebyshev2 quad(order, -1., 1.); integrate(quad); } + SECTION("Third Kind") { + IntegratorXX::GaussChebyshev3 quad(order, -1., 1.); + integrate(quad); + } } TEST_CASE( "Euler-Maclaurin Quadratures", "[1d-quad]" ) { From d1f4b815c8c3a24fcead99ba9082602bc09101b2 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Fri, 7 Jul 2023 16:34:29 -0700 Subject: [PATCH 4/5] Add UT for GaussChebyshev2Modified --- test/1d_quadratures.cxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index cb77178..08aa7d1 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -193,6 +193,10 @@ TEST_CASE( "Gauss-Chebyshev Quadratures", "[1d-quad]") { IntegratorXX::GaussChebyshev2 quad(order, -1., 1.); integrate(quad); } + SECTION("Second Kind (Modified)") { + IntegratorXX::GaussChebyshev2Modified quad(order, -1., 1.); + integrate(quad); + } SECTION("Third Kind") { IntegratorXX::GaussChebyshev3 quad(order, -1., 1.); integrate(quad); From ea92108426dd50a206cd2c23de890ab877d459a5 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sat, 8 Jul 2023 11:33:36 -0700 Subject: [PATCH 5/5] Keep temporaries in registers where appropraite --- include/integratorxx/quadratures/gausscheby1.hpp | 13 ++++++++----- include/integratorxx/quadratures/gausscheby2.hpp | 16 ++++++++-------- .../integratorxx/quadratures/gausscheby2_mod.hpp | 7 ++++--- include/integratorxx/quadratures/gausscheby3.hpp | 12 +++++++----- test/1d_quadratures.cxx | 1 - 5 files changed, 27 insertions(+), 22 deletions(-) diff --git a/include/integratorxx/quadratures/gausscheby1.hpp b/include/integratorxx/quadratures/gausscheby1.hpp index 1de1ccc..20963b9 100644 --- a/include/integratorxx/quadratures/gausscheby1.hpp +++ b/include/integratorxx/quadratures/gausscheby1.hpp @@ -70,13 +70,16 @@ struct quadrature_traits> { for (size_t i = 0; i < npts; ++i) { const auto ti = (2.0 * (i + 1) - 1.) * pi_ov_2npts; // The standard nodes and weights are given by - points[i] = -std::cos(ti); - weights[i] = pi_ov_npts; - + const auto xi = std::cos(ti); + auto wi = pi_ov_npts; + // However, as we're integrating f(x) not \frac{f(x)}{\sqrt{1 - // x^2}}, we must factor the \sqrt{1-x^2} into the weight - // sqrt(1 - cos^2(x)) = abs(sin(x)) - weights[i] *= std::sin(ti); //std::sqrt(1. - std::pow(points[i],2)); + wi *= std::sqrt(1. - xi*xi); + + // Store into memory + points[i] = xi; + weights[i] = wi; } return std::make_tuple(points, weights); diff --git a/include/integratorxx/quadratures/gausscheby2.hpp b/include/integratorxx/quadratures/gausscheby2.hpp index 3c8bac4..8121112 100644 --- a/include/integratorxx/quadratures/gausscheby2.hpp +++ b/include/integratorxx/quadratures/gausscheby2.hpp @@ -66,18 +66,18 @@ struct quadrature_traits> { point_container points(npts); for (size_t i = 0; i < npts; ++i) { const auto ti = (i + 1) * pi_ov_npts_p_1; - // The standard nodes and weights are given by - points[i] = std::cos(ti); + // The standard nodes are given by + const auto xi = std::cos(ti); // Avoid float point cancellations - form modified weight directly in following statements - //weights[i] = - // pi_ov_npts_p_1 * std::pow(std::sin((i + 1) * pi_ov_npts_p_1), 2); - // However, since we want the rule with a unit weight factor, we // divide the weights by sqrt(1-x^2). - // PI / (n+1) * sin^2(x) / sqrt(1 - cos^2(x)) = PI / (n+1) * sin(x) - //weights[i] /= std::sqrt(1.0 - std::pow(points[i], 2)); - weights[i] = pi_ov_npts_p_1 * std::sin(ti); + // PI / (n+1) * sin^2(x) / sqrt(1 - cos^2(x)) = PI / (n+1) * sqrt(1 - cos^2(x)) + const auto wi = pi_ov_npts_p_1 * std::sqrt(1.0 - xi*xi); + + // Store into memory + points[i] = xi; + weights[i] = wi; } return std::make_tuple(points, weights); diff --git a/include/integratorxx/quadratures/gausscheby2_mod.hpp b/include/integratorxx/quadratures/gausscheby2_mod.hpp index d793cfe..c66a895 100644 --- a/include/integratorxx/quadratures/gausscheby2_mod.hpp +++ b/include/integratorxx/quadratures/gausscheby2_mod.hpp @@ -68,9 +68,10 @@ struct quadrature_traits> { point_container points(npts); weight_container weights(npts); for (size_t i = 1; i <= npts; ++i) { - auto sine = sin(i * M_PI * oonpp); - auto sinesq = sine * sine; - auto cosine = cos(i * M_PI * oonpp); + const auto ti = i * M_PI * oonpp; + const auto sine = std::sin(ti); + const auto sinesq = sine * sine; + const auto cosine = std::cos(ti); points[i - 1] = 1.0 - 2.0 * i * oonpp + M_2_PI * (1.0 + 2.0 / 3.0 * sinesq) * cosine * sine; diff --git a/include/integratorxx/quadratures/gausscheby3.hpp b/include/integratorxx/quadratures/gausscheby3.hpp index 3e7a65e..18f053a 100644 --- a/include/integratorxx/quadratures/gausscheby3.hpp +++ b/include/integratorxx/quadratures/gausscheby3.hpp @@ -63,17 +63,19 @@ struct quadrature_traits> { weight_container weights(npts); point_container points(npts); for (size_t i = 0; i < npts; ++i) { + const auto ti = 0.5 * (2*(i+1) - 1) * pi_ov_2n_p_1; + const auto cti = std::cos(ti); // The standard nodes and weights are given by - points[i] = std::pow(std::cos(0.5 * (2*(i+1) - 1) * pi_ov_2n_p_1),2); - weights[i] = 2.0*pi_ov_2n_p_1 * points[i]; + const auto xi = cti * cti; // cos^2(t) + auto wi = 2.0*pi_ov_2n_p_1 * xi; // However, since we want the rule with a unit weight factor, we // divide the weights by sqrt(x/(1-x)). - weights[i] *= std::sqrt((1.0 - points[i])/points[i]); + wi *= std::sqrt((1.0 - xi)/xi); // Finally, convert from [0,1] to [-1,1] - points[i] = 2*points[i] - 1.0; - weights[i] *= 2.0; + points[i] = 2*xi - 1.0; + weights[i] = 2.0 * wi; } return std::make_tuple(points, weights); diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 08aa7d1..159d995 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -231,7 +231,6 @@ TEST_CASE( "Ahlrichs Quadratures", "[1d-quad]" ) { double res = 0.; for( auto i = 0; i < quad.npts(); ++i ) { - //std::cout << wgt[i] << ", " << pts[i] << std::endl; res += wgt[i] * f(pts[i]); }