Share examples DF/QF

[ckormanyos]

Here is an issue related to collecting and sharing short and cool DF/QF examples.

I post the example below, which I have used in various forms. it uses an integral representation of cylindrical Bessel functions to compute 3 Bessel function values for real argument.

I tested this particular example on MSVC/GCC for both number<cpp_double_float<double>, et_off> as well as number<cpp_quad_float<double>, et_off>

It is wonderful to see DF/QF ready to roll on such cool example(s).

Feel free to try/share an example for our docs. Here in this thread!

///////////////////////////////////////////////////////////////////////////////
//  Copyright 2021 Fahad Syed.
//  Copyright 2021 Christopher Kormanyos.
//  Copyright 2021 Janek Kozicki.
//  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)
//
// Test for correctness of arithmetic operators of cpp_double_float<>

// cd /mnt/c/Users/User/Documents/Ks/PC_Software/Test
// g++ -O3 -Wall -march=native -std=c++11 -I/mnt/c/MyGitRepos/BoostGSoC21_multiprecision/include -I/mnt/c/boost/boost_1_76_0 test.cpp -o test_double_float.exe

// Handle interaction with Boost's wrap of libquadmath __float128.
// g++ -O3 -Wall -march=native -std=gnu++11 -I/mnt/c/MyGitRepos/BoostGSoC21_multiprecision/include -I/mnt/c/boost/boost_1_76_0 -DBOOST_MATH_USE_FLOAT128 test.cpp -lquadmath -o test_double_float.exe

#include <iomanip>
#include <iostream>

#include <boost/config.hpp>
#include <boost/multiprecision/number.hpp>
#ifdef BOOST_MATH_USE_FLOAT128
#include <boost/multiprecision/float128.hpp>
#endif
#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/cpp_quad_float.hpp>

namespace local
{
  template<typename RealValueType,
           typename RealFunctionType>
  RealValueType integral(const RealValueType& a,
                         const RealValueType& b,
                         const RealValueType& tol,
                         RealFunctionType real_function)
  {
    using real_value_type = RealValueType;

    std::uint_fast32_t n2(1);

    real_value_type step = ((b - a) / 2U);

    real_value_type result = (real_function(a) + real_function(b)) * step;

    const std::uint_fast8_t k_max = UINT8_C(32);

    for(std::uint_fast8_t k = UINT8_C(0); k < k_max; ++k)
    {
      real_value_type sum(0);

      for(std::uint_fast32_t j(0U); j < n2; ++j)
      {
        const std::uint_fast32_t two_j_plus_one = (j * UINT32_C(2)) + UINT32_C(1);

        sum += real_function(a + (step * two_j_plus_one));
      }

      const real_value_type tmp = result;

      result = (result / 2U) + (step * sum);

      using std::fabs;

      const real_value_type ratio = fabs(tmp / result);

      const real_value_type delta = fabs(ratio - 1U);

      if((k > UINT8_C(1)) && (delta < tol))
      {
        break;
      }

      n2 *= 2U;

      step /= 2U;
    }

    return result;
  }

  template<typename float_type>
  float_type cyl_bessel_j(const std::uint_fast8_t n,
                          const float_type& x)
  {
    const float_type epsilon = std::numeric_limits<float_type>::epsilon();

    using std::cos;
    using std::sin;
    using std::sqrt;

    const float_type tol = sqrt(epsilon);

    const float_type integration_result = integral(float_type(0),
                                                   boost::math::constants::pi<float_type>(),
                                                   tol,
                                                   [&x, &n](const float_type& t) -> float_type
                                                   {
                                                     return cos(x * sin(t) - (t * n));
                                                   });

    const float_type jn = integration_result / boost::math::constants::pi<float_type>();

    return jn;
  }
}

int main()
{
  using big_float_type = boost::multiprecision::number<boost::multiprecision::backends::cpp_quad_float<long double>, boost::multiprecision::et_off>;

  const big_float_type j2 = local::cyl_bessel_j(2U, big_float_type(123U) / 100U);
  const big_float_type j3 = local::cyl_bessel_j(3U, big_float_type(456U) / 100U);
  const big_float_type j4 = local::cyl_bessel_j(4U, big_float_type(789U) / 100U);

  using std::fabs;

  // N[BesselJ[2, 123/100], 200]
  const big_float_type control2 = big_float_type("0.16636938378681407351267852431513159437103348245332855514956220782731992705482241194987092301624901244714163881464364951544272000398183529475122904017094866518071095290847014050736582125030318499776851");

  // N[BesselJ[3, 456/100], 200]
  const big_float_type control3 = big_float_type("0.42038820486765216162613462343078475742748171202057848576174361927529718085488573577015201639074037018491599838594935660836073991213044095508174498623210403620603617576191077285889954144530140343573653");

  // N[BesselJ[4, 789/100], 200]
  const big_float_type control4 = big_float_type("-0.078506863572127438410485520328806569617327186592090417112715332281301496921809793898580280982628142335829796109155480178162062502927315279160886474735074686005433233598424782354297994378751386383997349");

  using std::fabs;

  const big_float_type closeness2 = fabs(1 - (j2 / control2));
  const big_float_type closeness3 = fabs(1 - (j3 / control3));
  const big_float_type closeness4 = fabs(1 - (j4 / control4));

  const big_float_type tol = std::numeric_limits<big_float_type>::epsilon() * 4U;

  std::cout << closeness2 << std::endl;
  std::cout << closeness3 << std::endl;
  std::cout << closeness4 << std::endl;
  std::cout << tol        << std::endl;

  const bool result_is_ok = (   (closeness2 < tol)
                             && (closeness3 < tol)
                             && (closeness4 < tol));

  std::cout << "result_is_ok: " << std::boolalpha << result_is_ok << std::endl;

  return ((result_is_ok == true) ? 0 : -1);
}

[sinandredemption]

There is always my GSoC'21 submission, which was quite elegantly modified by @ckormanyos : sinandredemption/GSoC2021#2

[sinandredemption]

I also imagine that an example involving computing the continued fraction of a rational number via classical method (mod, sub, inverse, repeat) should intuitively provide some evidence of correctness, in the case the continued fraction came out to be finite.

[ckormanyos]

my GSoC'21 submission

Yes, in fact, I do intend to work up part of that, possibly without requesting user input from cin.

[ckormanyos]

The following example uses the tgamma() function from Boost.Math's special functions collection to verify the calculation of Gamma(1/2) using the cpp_quad_float<double> number backend.

#include <iomanip>
#include <iostream>

#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/multiprecision/cpp_quad_float.hpp>

namespace local
{
   using big_float_type =
     boost::multiprecision::number<boost::multiprecision::backends::cpp_quad_float<double>, boost::multiprecision::et_off>;
}

int main()
{
   using local::big_float_type;

   // N[Sqrt[Pi], 201]
   // 1.77245385090551602729816748334114518279754945612238712821380778985291128459103218137495065673854466541622682362428257066623615286572442260252509370960278706846203769865310512284992517302895082622893210

   std::cout << std::endl << "represent_tgamma_half" << std::endl;

   const big_float_type g    = boost::math::tgamma(big_float_type(0.5F));
   const big_float_type ctrl = sqrt(boost::math::constants::pi<big_float_type>());

   const big_float_type delta = fabs(1 - (g / ctrl));

   const std::streamsize original_streamsize = std::cout.precision();
   std::cout << std::setprecision(std::numeric_limits<big_float_type>::digits10) << g    << std::endl;
   std::cout << std::setprecision(std::numeric_limits<big_float_type>::digits10) << ctrl << std::endl;
   std::cout << std::scientific << std::setprecision(4) << delta << std::endl;
   std::cout.precision(original_streamsize);
   std::cout.unsetf(std::ios::scientific);

   const bool result_is_ok = (delta < std::numeric_limits<big_float_type>::epsilon() * 40U);

   std::cout << "result_is_ok: " << std::boolalpha << result_is_ok << std::endl;

   return (result_is_ok ? 0 : -1);
}

[ckormanyos]

As the cpp_double_fp_backend class gets into better and better form, keep an eye on nice examples for Boost's docs.

[ckormanyos]

This will be handled in #124. Get any notes from here if needed.