diff --git a/doc/equations/agm_bound.svg b/doc/equations/agm_bound.svg new file mode 100644 index 000000000..3981d3224 --- /dev/null +++ b/doc/equations/agm_bound.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/equations/agm_inequality.svg b/doc/equations/agm_inequality.svg new file mode 100644 index 000000000..7787bcd66 --- /dev/null +++ b/doc/equations/agm_inequality.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/equations/agm_sequence.svg b/doc/equations/agm_sequence.svg new file mode 100644 index 000000000..441e37c68 --- /dev/null +++ b/doc/equations/agm_sequence.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/graphs/agm_ulps_plot.svg b/doc/graphs/agm_ulps_plot.svg new file mode 100644 index 000000000..987e3b77e --- /dev/null +++ b/doc/graphs/agm_ulps_plot.svg @@ -0,0 +1,2544 @@ + + + + + + + +-1.5 + +-1 + +-0.5 + +0.5 + +1 + +1.5 + +1001 + +2001 + +3001 + +4001 + +5000 + +6000 + +7000 + +8000 + +9000 + +1e+04 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/internals/agm.qbk b/doc/internals/agm.qbk new file mode 100644 index 000000000..28f585dad --- /dev/null +++ b/doc/internals/agm.qbk @@ -0,0 +1,77 @@ +[/ + Copyright Nick Thompson, John Maddock, 2020 + Distributed under 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). +] + +[section:agm Arithmetic-Geometric Mean] + +[h4 Synopsis] + + #include + + namespace boost::math::tools { + + template + Real agm(Real a0, Real g0); + + } // namespaces + + +The function `agm` produces the limiting value of the sequence + +[$../equations/agm_sequence.svg] + + +A basic usage is + + double G = boost::math::tools::agm(sqrt(2.0), 1.0); + +The AGM inequality + +[$../equations/agm_sequence.svg] + +shows that + +[$../equations/agm_bound.svg] + +We use this condition internally to measure convergence; however, there is no need to worry about putting arguments in the correct order since we extend `agm` to a symmetric function by definition. +Both arguments must be non-negative, as the sequence becomes complex for negative arguments. +(We have not implemented the complex-valued AGM sequence.) +The function `agm` is "essentially" one-dimensional, as the scale invariance `agm(k*x, k*y) == k*agm(x,y)` always allows us to take one argument to be unity. +The following ULP plot has been generated with the function `agm(x, Real(1))`: + +[$../graphs/agm_ulps_plot.svg] + +The graph above shows an ulps plot of the Boost implementation of `agm(x, Real(1))`. +An ~2 ULP bound is to be expected. + +A google benchmark for various types is available in `boost/libs/math/reporting/performance/test_agm.cpp`; some results on a consumer laptop are provided for convenience: + +``` +Run on (16 X 2300 MHz CPU s) +CPU Caches: + L1 Data 32K (x8) + L1 Instruction 32K (x8) + L2 Unified 262K (x8) + L3 Unified 16777K (x1) +Load Average: 2.02, 2.14, 2.00 +------------------------------------------------------------------------------- +Benchmark Time CPU Iterations +------------------------------------------------------------------------------- +AGM 8.52 ns 8.51 ns 59654685 +AGM 13.5 ns 13.5 ns 51709746 +AGM 30.6 ns 30.6 ns 18745247 +AGM 2332 ns 2332 ns 299303 +``` + +If any inputs are NaNs, the result is a NaN. +If any inputs are +∞, the result is +∞, unless the other argument fails NaN or negative validation. + +[heading References] + +* Steven R. Finch. ['Mathematical Constants] Cambridge, 2003. + + +[endsect] diff --git a/doc/math.qbk b/doc/math.qbk index 89053b0ad..06569ad05 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -748,6 +748,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include internals/internals_overview.qbk] [section:internals Internal tools] [include internals/series.qbk] +[include internals/agm.qbk] [include internals/fraction.qbk] [include internals/recurrence.qbk] [/include internals/rational.qbk] [/moved to tools] diff --git a/example/Jamfile.v2 b/example/Jamfile.v2 index feaf77faa..e7076a2f9 100644 --- a/example/Jamfile.v2 +++ b/example/Jamfile.v2 @@ -74,6 +74,7 @@ test-suite examples : [ run inverse_chi_squared_example.cpp ] [ run legendre_stieltjes_example.cpp : : : [ requires cxx11_auto_declarations cxx11_defaulted_functions cxx11_lambdas ] ] [ run airy_ulps_plot.cpp : : : [ requires cxx17_std_apply cxx17_if_constexpr ] ] + [ run agm_example.cpp : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx17_std_apply cxx17_if_constexpr ] ] #[ # run inverse_chi_squared_find_df_example.cpp ] #[ run lambert_w_basic_example.cpp ] [ run lambert_w_basic_example.cpp : : : [ requires cxx11_numeric_limits ] ] diff --git a/example/agm_example.cpp b/example/agm_example.cpp new file mode 100644 index 000000000..13208d575 --- /dev/null +++ b/example/agm_example.cpp @@ -0,0 +1,21 @@ +// (C) Copyright Nick Thompson 2020. +// 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 +#include +#include +#include + +// This example computes the lemniscate constant to high precision using the agm: +using boost::math::tools::agm; +using boost::math::constants::pi; +using boost::multiprecision::float128; + +int main() { + using Real = float128; + Real G = agm(sqrt(Real(2)), Real(1)); + std::cout << std::setprecision(std::numeric_limits::max_digits10); + std::cout << " Gauss's lemniscate constant = " << pi()/G << "\n"; + std::cout << "Expected lemniscate constant = " << "2.62205755429211981046483958989111941368275495143162316281682170380079058707041425023029553296142909344613\n"; +} diff --git a/include/boost/math/tools/agm.hpp b/include/boost/math/tools/agm.hpp new file mode 100644 index 000000000..157f389e1 --- /dev/null +++ b/include/boost/math/tools/agm.hpp @@ -0,0 +1,44 @@ +// (C) Copyright Nick Thompson 2020. +// 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) + +#ifndef BOOST_MATH_TOOLS_AGM_HPP +#define BOOST_MATH_TOOLS_AGM_HPP +#include + +namespace boost::math::tools { + +template +Real agm(Real a, Real g) +{ + if (a < g) + { + // Mathematica, mpfr, and mpmath are all symmetric functions: + return agm(g, a); + } + // Use: M(rx, ry) = rM(x,y) + if (a <= 0 || g <= 0) { + if (a < 0 || g < 0) { + return std::numeric_limits::quiet_NaN(); + } + return Real(0); + } + + // The number of correct digits doubles on each iteration. + // Divide by 512 for some leeway: + const Real scale = sqrt(std::numeric_limits::epsilon())/512; + while (a-g > scale*g) + { + Real anp1 = (a + g)/2; + g = sqrt(a*g); + a = anp1; + } + + // Final cleanup iteration recovers down to ~2ULPs: + return (a + g)/2; +} + + +} +#endif diff --git a/reporting/accuracy/test_agm.cpp b/reporting/accuracy/test_agm.cpp new file mode 100644 index 000000000..ca28feaba --- /dev/null +++ b/reporting/accuracy/test_agm.cpp @@ -0,0 +1,36 @@ +// (C) Copyright Nick Thompson 2020. +// 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 +#include +#include +#include +#include + +using boost::math::tools::ulps_plot; +using boost::math::tools::agm; + +int main() { + using PreciseReal = boost::multiprecision::float128; + using CoarseReal = float; + + auto agm_coarse = [](CoarseReal x) { + return agm(x, CoarseReal(1)); + }; + auto agm_precise = [](PreciseReal x) { + return agm(x, PreciseReal(1)); + }; + + std::string filename = "agm_" + boost::core::demangle(typeid(CoarseReal).name()) + ".svg"; + int samples = 2500; + int width = 1100; + PreciseReal clip = 100; + auto plot = ulps_plot(agm_precise, CoarseReal(0), CoarseReal(10000), samples); + plot.clip(clip).width(width); + std::string title = "AGM ULP plot at " + boost::core::demangle(typeid(CoarseReal).name()) + " precision"; + //plot.title(title); + plot.vertical_lines(10); + plot.add_fn(agm_coarse); + plot.write(filename); +} diff --git a/reporting/performance/test_agm.cpp b/reporting/performance/test_agm.cpp new file mode 100644 index 000000000..53ee4b358 --- /dev/null +++ b/reporting/performance/test_agm.cpp @@ -0,0 +1,32 @@ +// (C) Copyright Nick Thompson 2020. +// 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 +#include +#include +#include + +using boost::math::tools::agm; +template +void AGM(benchmark::State& state) +{ + std::random_device rd; + std::mt19937_64 mt(rd()); + std::uniform_real_distribution unif(1,100); + + Real x = static_cast(unif(mt)); + for (auto _ : state) + { + benchmark::DoNotOptimize(agm(x,Real(1))); + } +} + +BENCHMARK_TEMPLATE(AGM, float); +BENCHMARK_TEMPLATE(AGM, double); +BENCHMARK_TEMPLATE(AGM, long double); +BENCHMARK_TEMPLATE(AGM, boost::multiprecision::float128); + + +BENCHMARK_MAIN(); diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index d3f7e68c0..7eb9b5520 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -945,6 +945,7 @@ test-suite misc : [ run daubechies_scaling_test.cpp : : : msvc:/bigobj [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run daubechies_wavelet_test.cpp : : : msvc:/bigobj [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run wavelet_transform_test.cpp : : : msvc:/bigobj [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] + [ run agm_test.cpp : : : msvc:/bigobj [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ compile compile_test/daubechies_filters_incl_test.cpp : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ compile compile_test/daubechies_scaling_incl_test.cpp : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run whittaker_shannon_test.cpp : : : [ requires cxx11_auto_declarations cxx11_constexpr cxx11_smart_ptr cxx11_defaulted_functions ] ] diff --git a/test/agm_test.cpp b/test/agm_test.cpp new file mode 100644 index 000000000..863d3c5d2 --- /dev/null +++ b/test/agm_test.cpp @@ -0,0 +1,105 @@ +/* + * Copyright Nick Thompson, 2019 + * 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 +#include +#include +#include +#include +#include +#ifdef BOOST_HAS_FLOAT128 +#include +using boost::multiprecision::float128; +#endif + +using boost::math::tools::agm; + +template +void test_gauss_constant() +{ + // http://oeis.org/A014549/constant + Real G_expected = boost::lexical_cast(".83462684167407318628142973279904680899399301349034700244982737010368199270952641186969116035127532412906785"); + + Real G_computed = 1/agm(sqrt(Real(2)), Real(1)); + if(!CHECK_ULP_CLOSE(G_expected, G_computed, 2)) { + std::cerr << " Gauss constant not computed correctly.\n"; + } +} + +template +void test_scaling() +{ + Real a = 2; + Real g = 1; + Real scale = 7; + Real expected = agm(scale*a, scale*g); + Real computed = scale*agm(a, g); + if(!CHECK_ULP_CLOSE(expected, computed, 2)) { + std::cerr << " Scaling property agm(kx,ky) = k*agm(x, y) is violated.\n"; + } + + expected = 0; + computed = agm(a, Real(0)); + if(!CHECK_ULP_CLOSE(expected, computed, 0)) { + std::cerr << " agm(a, 0) != 0.\n"; + } + + computed = agm(Real(0), Real(0)); + if(!CHECK_ULP_CLOSE(expected, computed, 0)) { + std::cerr << " agm(0, 0) != 0.\n"; + } + + expected = 1; + computed = agm(Real(1), Real(1)); + if(!CHECK_ULP_CLOSE(expected, computed, 0)) { + std::cerr << " agm(1, 1) != 1.\n"; + } + + expected = 7; + computed = agm(Real(7), Real(7)); + if(!CHECK_ULP_CLOSE(expected, computed, 0)) { + std::cerr << " agm(7, 7) != 1.\n"; + } + + // Properties I found at: https://mathworld.wolfram.com/Arithmetic-GeometricMean.html + // agm(x,y) = agm((x+y)/2, sqrt(xy)) + expected = agm(Real(3), Real(1)); + computed = agm(Real(2), sqrt(Real(3))); + if(!CHECK_ULP_CLOSE(expected, computed, 0)) { + std::cerr << " agm(x, y) != agm((x+y)/2, sqrt(xy)).\n"; + } + + //computed = agm(std::numeric_limits::infinity(), Real(7)); + //std::cout << "Computed at infinity = " << computed << "\n"; + + for (Real x = 0; x < 1; x += Real(1)/128) { + expected = agm(Real(1), sqrt(1-x*x)); + computed = agm(1+x, 1-x); + if(!CHECK_ULP_CLOSE(expected, computed, 0)) { + std::cerr << " agm(1, sqrt(1-x^2) != agm(1+x,1-x).\n"; + } + } +} + + +int main() +{ + test_gauss_constant(); + test_gauss_constant(); + test_gauss_constant(); + + test_scaling(); + test_scaling(); + test_scaling(); + + #ifdef BOOST_HAS_FLOAT128 + test_gauss_constant(); + test_scaling(); + #endif + return boost::math::test::report_errors(); +}