Ensure function tests running

[ckormanyos]

The purpose of this PR is to check the original GSoC algos together with all improvements on formal conversoins, limits, elementary functions, string read, etc.

The actual source code is capable of switching from the original GSoC algos to the Briggs style algos with a new compiler switch. in this PR the switch is set for tunning the origninal GSoC algos.

Cc: @sinandredemption

[ckormanyos]

Hi Fahad (@sinandredemption) It seems like we now have 5 test files failing. This last push should run through but is not intended have any functional changes, as it was simply a Git merge from develop into this branch.

If it stays so good, we will merge to develop. Feel free to merge this branch and PR into develop after CI runs, please. Then I'll start removing some of my own redundant branches, as this gives us here probaly what we deed to finish cpp_double_fp<>.

In detail, one of the failing test cases that seem interesting includes the following, in which cpp_double_double< 8-byte-double > fails, but number< cpp_dec_float<51> > passes.

Let's see where it overflows or underflows. What does normal double do in this situation?

#include <iostream>

#include <boost/math/special_functions/bessel.hpp>
#include <boost/multiprecision/cpp_double_fp.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

int main()
{
  using boost::multiprecision::cpp_double_double;

  using dec_float_type = boost::multiprecision::number<boost::multiprecision::cpp_dec_float<51>, boost::multiprecision::et_off>;

  // N[BesselK[-1000, 700], 101]
  // 6.5156197914473581890355385260638331298440936198412822153940457066208936122301163064677298386173546220*10^-31

  const auto kd1 = boost::math::cyl_bessel_k(cpp_double_double(-1000), cpp_double_double(700));
  const auto kd2 = boost::math::cyl_bessel_k(dec_float_type   (-1000), dec_float_type   (700));

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

  std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<cpp_double_double>::digits10)) << kd1 << std::endl;
  std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<dec_float_type   >::digits10)) << kd2 << std::endl;

  std::cout.flags(flg);
}

[ckormanyos]

Hi Fahad, it looks like buit-in double gives a reasonable result.

I wonder if the limited exponent range of the cpp_double_double is causing trouble? Or is there more going on in the depths of the operations?

#include <iostream>

#include <boost/math/special_functions/bessel.hpp>
#include <boost/multiprecision/cpp_double_fp.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

int main()
{
  using boost::multiprecision::cpp_double_double;

  using dec_float_type = boost::multiprecision::number<boost::multiprecision::cpp_dec_float<51>, boost::multiprecision::et_off>;

  // N[BesselK[-1000, 700], 101]
  // 6.5156197914473581890355385260638331298440936198412822153940457066208936122301163064677298386173546220*10^-31

  const auto kd1 = boost::math::cyl_bessel_k(cpp_double_double(-1000), cpp_double_double(700));
  const auto kd2 = boost::math::cyl_bessel_k(dec_float_type   (-1000), dec_float_type   (700));
  const auto kd3 = boost::math::cyl_bessel_k(double           (-1000), double           (700));

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

  std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<cpp_double_double>::digits10)) << kd1 << std::endl;
  std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<dec_float_type   >::digits10)) << kd2 << std::endl;
  std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<double           >::digits10)) << kd3 << std::endl;

  std::cout.flags(flg);
}

[ckormanyos]

Hi Fahad (@sinandredemption), and finally, test cases similar to this particular test case starts to fail when the x-coordinate of the Bessel-K function increases. I seem to recall the multiplication of two extreme numbers for the Bessel-K. Could it be that x-coordinate 700 simply can't be used with cpp_double_fp< 8-byte-double >?

#include <iostream>

#include <boost/math/special_functions/bessel.hpp>
#include <boost/multiprecision/cpp_double_fp.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

int main()
{
  using boost::multiprecision::cpp_double_double;

  using dec_float_type = boost::multiprecision::number<boost::multiprecision::cpp_dec_float<51>, boost::multiprecision::et_off>;

  for(auto   index = static_cast<unsigned>(UINT16_C(660));
             index < static_cast<unsigned>(UINT16_C(680));
           ++index)
  {
    std::cout << "calculation for index: " << index << std::endl;

    const auto kd1 = boost::math::cyl_bessel_k(cpp_double_double(-1000), cpp_double_double(index));
    const auto kd2 = boost::math::cyl_bessel_k(dec_float_type   (-1000), dec_float_type   (index));
    const auto kd3 = boost::math::cyl_bessel_k(double           (-1000), double           (index));

    std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<cpp_double_double>::digits10)) << kd1 << std::endl;
    std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<dec_float_type   >::digits10)) << kd2 << std::endl;
    std::cout << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<double           >::digits10)) << kd3 << std::endl;
  }

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

  std::cout.flags(flg);
}

[sinandredemption]

Thank you Chris (@ckormanyos) for your detailed and helpful analysis. I am investigating the issues, and will let you know my thoughts soon.

[sinandredemption]

I wonder if the limited exponent range of the cpp_double_double is causing trouble?

Yes Chris (@ckormanyos), it seems like that is the case. I investigated the first case you sent, and it seems while computing BesselK[-1000, 700], at one point exp(-700) is computed, which is approximately 9.86 × 10^-305. An exponent of -305 is out-of-range for cpp_double_fp<double>, hence eval_exp() returns 0 which cascades to the final output being 0 as well.

It makes me wonder, is this acceptable? A casual user of this library might not expect such quirks and has to be very careful in advance if the resulting operation would underflow, given the extreme asymmetry of the exponent range. Does that mean we should optionally provide special underflow flags for this data type?

EDIT: Note that double_fp constructions have a asymmetrical exponent range, i.e. valid exponent ranges for cpp_double_fp<double> and built-in 8-byte double are (-272, 308) and (-308, 308) respectively.

[sinandredemption]

More details: the underflow occurs here in function bessel_k0_imp:

return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x)); with x == 700.

The part that underflows is exp(-x) which evaluates to exp(-700) == 0 for cpp_double_fp<double> because of the out-of-range exponent as taken care of by this check inside eval_exp.

This causes the function to return 0, which finally causes all the iterations to evaluate to 0 while computing the bessel function with function bessel_kn.

CC: @ckormanyos

[ckormanyos]

double_fp constructions have a asymmetrical exponent range, i.e. valid exponent ranges for cpp_double_fp and built-in 8-byte double are (-272, 308) and (-308, 308) respectively.

Oh yes. Absolutely. I believe we have to find the range(s) where construct and/or possibly even mul/div will lead to overflow/underflow. I suspect these might be asymmetrically reduced ranges compared to the constituent type.

Then when doing add/sub/mul/div, we usually detect overflow/underflow/infinitiy/NaN and zero purposefully and set these accorgingly to avoid nonsense such as inf/inf leading to NaN.

As it turns out, I believe our ranges and edeg-case detection maight, in fact, be the most significant blocking points on specfun at the moment.

We should do a bit more thinking and work on range checking and might even need to adjust the selected min/max values for the composite type.

We might hash out this discussion here as well as in offline chats and maybe future issues. But I think this is kind of the last big hurdle for the cpp_double_fp<> backend.

Cc: @sinandredemption

[ckormanyos]

Anyway, since this is really our, let's say, best state so far with the good CI and just a few failing files in spcefun, I'll go ahead and merge.

I'll be cleaning up some of my redundant branches, since the new develop hass kind of all the right stuff to probably finish...

Cc: @sinandredemption

[ckormanyos]

Hi Fahad (@sinandredemption), thanks for investigating these matters.

I wrote a long text on non-finite values in our offline chat, as I believe we have quite a lot more depth on this issue to discuss before formally finishing edge cases.

Regarding the exact BesselK test case, I'm not actually sure if the cpp_double_fp< 8-byte-float > type could ever be capable of computing cyl_bessel_k(-1000, 700) with the current implementation of the modified Bessel function of type K. We might have to reduce the x-value of the test case or remove this case entirely when this data type is being tested.

We might re-visit this exact test case when we get a solid grip on infinities, NaNs, zeros and possible subnormals.