Add cpp_double_fp and exercise arithmetic_tests

[ckormanyos]

The purpose of this pull request is to add the cpp_double_fp_backend<> template class.

In the initial commit, we:

• Run the arithmetic tests.
• Modify elementary function tests to prepare for TEST_CPP_DOUBLE_FLOAT, but do not run them yet.
• Try some further small adaptions to arithmetic_tests.hpp that improve the testing logic for mixed-int128/float128 interaction.

If this runs through, I can add the function tests (stuff like sin(), exp(), etc.) as a next step. These run well and properly in the GSoC fork, but are not yet exercised in this initial commit.

[ckormanyos]

The first CI run (even with only the tiny step of arithmetic tests and modified but not running function tests) has failed due to a rookie typo involving a missing semicolon.

Running again now...

[ckormanyos]

OK the last CI run passed on GHA but failed a silly typo on performance_test on Circle CI.

The latest push goes a lot farther with an attempt for:

• correcting the above and also
• adding most of functions-and-limits.

[ckormanyos]

OK.

• GHA is running green with arithmetic tests and most of functions and limits.
• In drone I think I forgot some needed quadmath inclusions in the JAMfile.
• Circle CI failed on specfun_cpp_dec_float, which might be unrelated to the double-double work.

I will press on trying to correct the Jamfile for the drone runs.

After the next cycle, I can collect the know weaknesses - and there are a few, especially with 10-byte long double mysteries.

And although this is going well, I'd like to resolve any known mysteries (maybe review needed) before pressing on to specfun and maybe more deep conversoin tests.

[ckormanyos]

OK. I'll jump the CI a bit and think this latest commit might bring this work to a small juncture.

At this juncture, coull you John (@jzmaddock) and Matt (@mborland) please take a cursory glance at the header/tests and give me a sanity check? If you see something waaayyy wrong or out of line, please give me a heads up?

Most of the tests and things are handled by cpp_double_fp_backend<> now at the moment.

Known problems include:

1. Round-tripping and banker's rounding (which may be related) and seem to involve the junction between the low/high parts of the double-fp.
2. Problems with edge-cases of 10-byte double-long-double on platforms having this 10-byte long double type. For instance, I get non-finite log() of double-long-double's (max)(), but only for 10-byte long double. I suppose there might be a classic, Chris, dont't you remember, ... something, or a few L suffixes, or similar.

In light of this, I'll retreat to the fork a bit and try to look into at least the 10-byte long double mysteries prior to going in full for spce_fun.

After solving the mysteries, going for spec_fun and setting up some rudimentary docs, this thing will actually be close to ready.

Cc: @sinandredemption and @cosurgi

[jzmaddock]

@ckormanyos my apologies for not being able to help out much right now: but a quick scan over the sources shows nothing obviously wrong.

[ckormanyos]

a quick scan over the sources shows nothing obviously wrong

Thanks John (@jzmaddock) I'm actually doing alright. I am always amazed how powerful and deep the Multiprecision test suite(s) are. These are wonderful tests.

I've found a few dependencies on libquadmath.a in the final drone CI run that I will eliminate in the limits. Other than that (and the aforementioned problems with rounding and 10.byte long double) , I'm comfortable to simply keep going.

[mborland]

By and large this looks good. I added some minor comments that are applicable in a number of places.

[ckormanyos]

OK. This seems like it's about to go green, which means that the formal correctness of cpp_double_fp_backend is nearing completion.

There are, nonetheless, a whole bunch of steps to do.

• We are actually now farther ahead in the Boost.Multiprecision PR than on our Fork. So I will need to merge from Boost's branch "back to" our GSoC fork. I anticipate doing this today or tomorrow.
• There are a bunch of open items needing attention.

These include:

• Banker's rounding is not fully working. Why? Should it work? Why or Why not?
• The selection of cpp_double_fp_backend::my_value_max() is unclear. Can we justify this value?
• eval_exp() and eval_log() need further optimization.
• Optimize (or prove that is already optimized) the routine for division.
• Optimize (or prove that is already optimized) the routine for eval_sqrt().

The only known problem is with banker's rounding. Otherwise, it's just optimize, tests and docs.

The constexpr-ification is moving along very well. In fact, large parts of the double-double backend are already now constexpr. For this, I have written quite a bit of <cmath> functionality in a true constexpr-fashion. Here I borrowed inspiration from Matt (@mborland) in <ccmath.hpp>. But there is more work to go. We need to constexpr-ify ::floor() and maybe even ::log() and/or ::exp() for float, double, long double and __float128. This is a super-cool area to march forward in.

The next step for Boost testing will be to run the cpp_double_fp_backend through the gallery of special-function tests, known as the test-suite "specfun". We will probably want to do this locally at first. Also I honestly don't know exactly how to hook up a new Multiprecision backend to the specfun tests.

So I'll ask John (@jzmaddock), very briefly in a few keywords, how do I hook up a new backend to specfun testing in CI? briefly.

Cc: @sinandredemption

[ckormanyos]

briefly in a few keywords, how do I hook up a new backend to specfun testing in CI?

Hi @jzmaddock I think I'm seeing the pattern.

I need to:

• Pick up the test files in Jamfile for my particular group.
• Define the setup in the appropriate location(s) in source code.
• Define the test suite in Jamfile. Let's call itspecfun_cpp_double_double.
• Select my backend of choice.
• And the run a CI job (or command-line job) for the new specfun_cpp_double_double-suite.

I'll try in my fork first.

[jzmaddock]

Yup, duplicate this rule:

multiprecision/test/Jamfile.v2

Line 1472 in 8de4503

rule get_specfun_tests_float128

and rename the tests that are generated to match your new types name, then add the rule invocation below to actually create the testing target like this:

multiprecision/test/Jamfile.v2

Line 1521 in 8de4503

test-suite specfun_mpfr : [ get_specfun_tests_mpfr ] ;

Then modify https://github.com/boostorg/multiprecision/blob/develop/test/math/setup.hpp to set up testing of the new type, again copy the float128 case is the easiest way for what you have.

Finally, duplicate and modify this line:

multiprecision/test/math/instances/Jamfile.v2

Line 73 in 8de4503

lib test_instances_float128 : [ generate_objs "_float128" : 1 : <define>TEST_FLOAT128 ] : [ check-target-builds ../../../config//has_float128 : : <build>no ] ;

to create the template instantiations of the special functions. This shouldn't be necessary, but it speeds up compilation to do things this way.

And I think/hope that's the lot.

[ckormanyos]

Yup, duplicate this rule: ...

Copy that John (@jzmaddock). It looks like the struggle with VC 14.0 continues a tiny bit. And maybe some other lint. But I'm ready to give specfun for cpp_double_double a try. I expect some tolerances will fail and some compiler issues may crop up.

I keep looking ahead to the future, but after all that (and I do expect success), then the adaption of quad-float (from the original feature-request) should get a lot easier...

Cc: @sinandredemption

[ckormanyos]

And what does all this mean in the practical sense.

I actually think we have a fair shake (and i'd like to give it a fair try) to get cpp_double_fp_backend<T> into 1.82 (with docs, funcs and tests and CI).

Then I think later in the year we can take a more coherent run at quad-float, maybe even after getting some of the field experience with double-float.

Cc: @jzmaddock and @sinandredemption and @mborland

[jzmaddock]

Tolerances should be broadly in agreement with similar backends - so say cpp_bin_float_quad. The special fun tests are actually very good at really hammering basic arithmetic and have uncovered more than a few bugs in the past, so if you see anything wildly out, I would suspect a bug somewhere.

[ckormanyos]

The special fun tests are actually very good at really hammering basic arithmetic and have uncovered more than a few bugs in the past, so if you see anything wildly out, I would suspect a bug somewhere.

Indeed that's for sure John (@jzmaddock). I'm psyched on finding out how it goes. I guess there are still a few remaining problems in CI so I'll work through these prior to going for specfun.

This is the first time I've tried to add an entire backend on my own. It's challenging but a great learning experience.

Also to be clear, I'm only (we're only) going for cpp_double_fp_backend in this PR and (maybe) for 1.82. The quead stuff will come later.

Cc: @sinandredemption

[ckormanyos]

This first step took a while and has now gone green.
But we are not yet ready to merge.

We still need the following:

1. Test and pass specfun (might take a while)
2. Write and build small docs on cpp_double_fp_backend<T>.
3. Review.

I am happy with this first preliminary green milestone and will now turn to points 1 and 2, beginning with specfun.

At first I'll prepare this in our GSoC fork.

Cc: @jzmaddock and @sinandredemption and @cosurgi and @mborland

[ckormanyos]

Hi John (@jzmaddock),

I wonder if you could help just a tiny bit with guidance? Much to my delight my local tests of specfun running on g++-9, -std=c++17 were remarkably successful on a first run. I received the following last lines of a build log...

I got ...updated 337 targets... but also ...failed updating 12 targets... (see second text block below).

Now I need to see the failures.

From these failures on Zeta, do you have any guidance on how I can actually run the relevant test case and put it under a microscopic debug session? I feel like I do not know what the row or mean are all about? See first text block below:

For Zeta

====== BEGIN OUTPUT ======
Running 1 test case...
Tests run with GNU C++ version 9.3.0, GNU libstdc++ version 20200808, linux
Zeta: Random values greater than 1 with type cpp_double_double
zeta<cpp_double_double> Max = 0.1958 RMS Mean=0.0706
    worst case at row: 12
    { 4.988127231597900390625, 1.037268683773212337925891820825986369163 }


Zeta: Random values less than 1 with type cpp_double_double
zeta<cpp_double_double> Max = 147.7 RMS Mean=13.23
    worst case at row: 69
    { -17.99382781982421875, -0.08425805633923885570999830270532256204479 }


Zeta: Values close to and greater than 1 with type cpp_double_double
zeta<cpp_double_double> Max = 0.19 RMS Mean=0.08633
    worst case at row: 22
    { 1.49527740478515625, 2.631125733364186076058384433696635480314 }


Zeta: Values close to and less than 1 with type cpp_double_double
zeta<cpp_double_double> Max = 1.003e+04 RMS Mean=2735
    worst case at row: 16
    { 0.999987542629241943359375, -80273.18355079243243774478653269383556504 }
Peak error greater than expected value of 1000
../../../libs/math/test/handle_test_result.hpp(165): �[1;31;49merror: in "test_main": check bounds.first >= max_error_found has failed�[0;39;49m
Mean error greater than expected value of 200
../../../libs/math/test/handle_test_result.hpp(170): �[1;31;49merror: in "test_main": check bounds.second >= mean_error_found has failed�[0;39;49m


Zeta: Integer arguments with type cpp_double_double
zeta<cpp_double_double> Max = 8.453 RMS Mean=2.34
    worst case at row: 54
    { 56.0, 1.0000000000000000138777878097252328 }



�[1;31;49m*** 2 failures are detected in the test module "Master Test Suite"
�[0;39;49m
EXIT STATUS: 201
====== END OUTPUT ======

Build result:

...failed updating 12 targets...
...skipped 12 targets...
...updated 337 targets...

Now I'm having difficulty identifying exactly which lines/cases/tests are needing special attention.

[jzmaddock]

So if we look at:

Zeta: Values close to and less than 1 with type cpp_double_double
zeta<cpp_double_double> Max = 1.003e+04 RMS Mean=2735
    worst case at row: 16
    { 0.999987542629241943359375, -80273.18355079243243774478653269383556504 }
Peak error greater than expected value of 1000

The case causing the failure is zeta(0.999987542629241943359375) which should yield -80273.18355079243243774478653269383556504 but is 10^4eps out compared to an expected tolerance of 10^3. Looking at the expected results for other types, this does indeed look like a slight outlier.

That case should go through:

template <class T, class Policy>
T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant<int, 113>&)

which is a rational approximation and should be super accurate - maybe check that the constant coefficients are being initialized correctly? For comparison, maybe check the error rate for cpp_bin_float_quad.

BTW I've been meaning to create a "dual precision" backend for multiprecision that would let you run the same calculation on two backends simultaneously and see where they they start to differ.

[jzmaddock]

Just tried:

   using mp = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<106, boost::multiprecision::digit_base_2> >;
   mp q("0.999987542629241943359375");
   auto z = boost::math::zeta(q);

   std::cout << boost::math::epsilon_difference(z, mp("-80273.18355079243243774478653269383556504")) << std::endl;

I get an error of 0, which looks good to me, so the difference with the double double looks like it might be real?

[ckormanyos]

Just tried:
get an error of 0, which looks good to me, so the difference with the double double looks like it might be real?

Indeed. I stepped through the impl. Your exact result makes sense to me. It especially makes sense, now that I study the implementation of zeta() in that region, which is solely based on algebraic polynomial expansion.

I fear this might take us all the way back to some fundamental bit-losses in the basic operations add/sub/mul/div in the low-level implementation of cpp_double_fp_backend<>.

We have an open issue (one of the first issues from July 2021) in the GSoC about this.

I fear that this type (its low-level eval_ops might not be fully cooked yet.

Cc: @sinandredemption This is something we still need to fully understand.

[ckormanyos]

get an error of 0, which looks good to me

Ahhh. OK, John (@jzmaddock), well, I figured out at least partly what's going on, but I don't know the best way to fix it entirely in the backend/test/etc. world.

• The input argument is an exact double. However the initialization of the argument is carried out with string-conversion, since we are above 64 digits, yet below 113 digits in the type.
• This causes imprecision in the digits immediately following the ...375 in the input.
• From that point on, the calculatoin carries this fuzz through.

If I simply use the double representation of the number to construct the argument x to zeta(x), then the input and following calculation retain exactness.

I'll need to see if the backend (on its own) can handle the case decision to initialize with the native floating-point type OR string. I think I might need more constructors. Or if the backend can't handle this alone, then we need to think harder.

[jzmaddock]

Is that not an error in string conversion then: I would assume it should convert to an exact binary value?

[ckormanyos]

Is that not an error in string conversion then: I would assume it should convert to an exact binary value?

I suppose maybe. That's what I was thinking and I do hope to find a solution there.

I went (lazily) directly to boost::multiprecision::detail::convert_from_string(). But there are many places where I suppose I could have an error. Will investigate.

[ckormanyos]

Is that not an error in string conversion then: I would assume it should convert to an exact binary value?

hope to find a solution

Hi John (@jzmaddock) regarding string input for cpp_double_fp_backend, after a while of thinking it's clear that you already solved this problem with cpp_bin_float<>. I racked my brain and found that for a binary backend, that's the working model. We need to (basically just like a compiler) support another multiprecision to fully interpret multiprecision strings.

This will pull cpp_int into the game. But I'm OK with that.

Along the lines of getting it working and then optimizing it, I find a usage of the cpp_bin_float-way for string-input is what is needed. Anything else does not capture the digits. And... Why should I write a different one when your code (adapted) will work for both of us?

I do not know if my ramblings make sense. But If so, I'd duplicate your code, test it, then with you together, find a way to make it publicly available for binary backends.

Thoughts? Am I crazy?

[jzmaddock]

Thoughts? Am I crazy?

Probably crazy to want to take on that much coding, but yes, I fear that is the correct way to go. The hard part is probably packing the bits in the resultant integers into the float type.

[sinandredemption]

Probably crazy to want to take on that much coding

@ckormanyos Please do let me know in chats if you need an extra hand. I've been trying to start contributing again.

[cosurgi]

It is great to see you going forward with this. Unfortunately I cannot join right now: too much work at the university at the moment :( But I will keep monitoring your changes, hoping to jump in when it becomes possible.

[ckormanyos]

you already solved this problem with cpp_bin_float<>

want to take on that much coding, but yes, I fear that is the correct way to go

Thanks for your input John. Yes. I will try - based on Your code that is already fully tested, well-established and cleanly available in cpp_bin_float.

This will have a bit of coding depth to it, so best to hammer this particular detail in our fork.

[ckormanyos]

let me know in chats if you need an extra hand.

Hi Fahad (@sinandredemption) we'll be moving back and forth to/from our GSoC fork for a couple of the more in-depth corrections.

This PR has already shown that:

• cpp_double_fp_backend seems robust enough to become a reliable badkend.
• We still need to work on exact string read.
• I'd like to study if we do/do-not have optimal algorithms for add/sub/mul/div/sqrt.
• Passing all of specfun still needs a bunch of work locally like in our fork.

I am going to start the first ponit back in our fork. I'll be posting more details on the above-mentioned issues of the fork. We can work out the heavy details there, and we can keep this PR open as we progress.