Trigonometric functions for overly large arguments

[ckormanyos]

Based on the detailed discussion in #135, the following test case (shown below) leads to partly unusual behavior for the sin() function.

The bahaviors are:

1. Correct answer with expected bit loss for cpp_double_fp< 8-byte-double >. OK.
2. Wrong answer, but no crash for original GSoC algorithms (add/sub/mul/div) for cpp_double_fp< 4-byte-float >.
3. Crash in eval_sin() for new (Briggs/Shoup-style) algorithms (add/sub/mul/div) for cpp_double_fp< 4-byte-float >.

Getting the wrong answer in 2 does not really bother me due to the low precision of the composite type.

But crashing in 3 or seg-faulting, etc. does bother me. It would be nice if possible to see what's going on there.

Cc: @sinandredemption and @cosurgi

///////////////////////////////////////////////////////////////
//  Copyright Christopher Kormanyos 2023.
//  Distributed under the Boost
//  Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
//

#include <iomanip>
#include <iostream>

#include <boost/multiprecision/cpp_double_fp.hpp>

namespace local
{
  //using double_fp_type = boost::multiprecision::cpp_double_float;
  using double_fp_type = boost::multiprecision::cpp_double_double;

  const auto control =
    static_cast<double_fp_type>
    (
      "-0.59543111331725610500680489068465901094029367239081751915820526461172514714208535369834909216413765236"
    );

  const auto euler_gamma =
    static_cast<double_fp_type>
    (
      "0.577215664901532860606512090082402431042159335939923598805767234884867726777664670936947063291746749514631447249807082480960504014486542836224173997644923536253500333742937337737673942792595258247094916008735203948165670853233151776611528621199501507984793745085705740029921354786146694029604325421519"
    );
}


int main()
{
  const auto flg = std::cout.flags();

  using local::double_fp_type;
  using local::euler_gamma;
  using local::control;

  // N[Sin[EulerGamma 4294962396], 101]
  // -0.59543111331725610500680489068465901094029367239081751915820526461172514714208535369834909216413765236

  constexpr auto crazy_k
  {
    std::uint32_t { UINT32_C(4294962396) }
  };

  const double_fp_type arg = double_fp_type(euler_gamma * crazy_k);

  const double_fp_type val = sin(arg);

  std::cout << std::setprecision(std::streamsize { std::numeric_limits<double_fp_type>::max_digits10 } )
            << val
            << std::endl;

  std::cout.flags(flg);
}

[ckormanyos]

At the moment, the new algos are in develop and the original algos are less-commonly available, but are still present in the get_specfun_running_in_ci branch.