ceph/src/boost/libs/math/test/agm_test.cpp

   1 /*
   2  * Copyright Nick Thompson, 2019
   3  * Use, modification and distribution are subject to the
   4  * Boost Software License, Version 1.0. (See accompanying file
   5  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6  */
   7
   8 #include "math_unit_test.hpp"
   9 #include <numeric>
  10 #include <utility>
  11 #include <random>
  12 #include <boost/lexical_cast.hpp>
  13 #include <boost/core/demangle.hpp>
  14 #include <boost/math/tools/agm.hpp>
  15 #ifdef BOOST_HAS_FLOAT128
  16 #include <boost/multiprecision/float128.hpp>
  17 using boost::multiprecision::float128;
  18 #endif
  19
  20 using boost::math::tools::agm;
  21
  22 template<class Real>
  23 void test_gauss_constant()
  24 {
  25     // http://oeis.org/A014549/constant
  26     Real G_expected = boost::lexical_cast<Real>(".83462684167407318628142973279904680899399301349034700244982737010368199270952641186969116035127532412906785");
  27
  28     Real G_computed = 1/agm(sqrt(Real(2)), Real(1));
  29     if(!CHECK_ULP_CLOSE(G_expected, G_computed, 2)) {
  30         std::cerr << "  Gauss constant not computed correctly.\n";
  31     }
  32 }
  33
  34 template<typename Real>
  35 void test_scaling()
  36 {
  37     Real a = 2;
  38     Real g = 1;
  39     Real scale = 7;
  40     Real expected = agm(scale*a, scale*g);
  41     Real computed = scale*agm(a, g);
  42     if(!CHECK_ULP_CLOSE(expected, computed, 2)) {
  43         std::cerr << "  Scaling property agm(kx,ky) = k*agm(x, y) is violated.\n";
  44     }
  45
  46     expected = 0;
  47     computed = agm(a, Real(0));
  48     if(!CHECK_ULP_CLOSE(expected, computed, 0)) {
  49         std::cerr << "  agm(a, 0) != 0.\n";
  50     }
  51
  52     computed = agm(Real(0), Real(0));
  53     if(!CHECK_ULP_CLOSE(expected, computed, 0)) {
  54         std::cerr << "  agm(0, 0) != 0.\n";
  55     }
  56
  57     expected = 1;
  58     computed = agm(Real(1), Real(1));
  59     if(!CHECK_ULP_CLOSE(expected, computed, 0)) {
  60         std::cerr << "  agm(1, 1) != 1.\n";
  61     }
  62
  63     expected = 7;
  64     computed = agm(Real(7), Real(7));
  65     if(!CHECK_ULP_CLOSE(expected, computed, 0)) {
  66         std::cerr << "  agm(7, 7) != 1.\n";
  67     }
  68
  69     // Properties I found at: https://mathworld.wolfram.com/Arithmetic-GeometricMean.html
  70     // agm(x,y) = agm((x+y)/2, sqrt(xy))
  71     expected = agm(Real(3), Real(1));
  72     computed = agm(Real(2), sqrt(Real(3)));
  73     if(!CHECK_ULP_CLOSE(expected, computed, 0)) {
  74         std::cerr << "  agm(x, y) != agm((x+y)/2, sqrt(xy)).\n";
  75     }
  76
  77     //computed = agm(std::numeric_limits<Real>::infinity(), Real(7));
  78     //std::cout << "Computed at infinity = " << computed << "\n";
  79
  80     for (Real x = 0; x < 1; x += Real(1)/128) {
  81         expected = agm(Real(1), sqrt(1-x*x));
  82         computed = agm(1+x, 1-x);
  83         if(!CHECK_ULP_CLOSE(expected, computed, 0)) {
  84             std::cerr << "  agm(1, sqrt(1-x^2) != agm(1+x,1-x).\n";
  85         }
  86     }
  87 }
  88
  89
  90 int main()
  91 {
  92     test_gauss_constant<float>();
  93     test_gauss_constant<double>();
  94     test_gauss_constant<long double>();
  95
  96     test_scaling<float>();
  97     test_scaling<double>();
  98     test_scaling<long double>();
  99
 100     #ifdef BOOST_HAS_FLOAT128
 101     test_gauss_constant<float128>();
 102     test_scaling<float128>();
 103     #endif
 104     return boost::math::test::report_errors();
 105 }