cpp_double_fp_backend<4-byte float> division fails accuracy threshold

[sinandredemption]

TLDR: cpp_double_fp_backend<4-byte float> is expected to preserve 44 bits of accuracy based on it's epsilon value, but only preserves 43 bits.

Hello @ckormanyos @cosurgi

I ran a modification of a old test I wrote during the GSoC days to assess how many bits of accuracy are preserved while performing add/sub/mul/div operations by both the David Bailey's implementation and our Briggs/Shoup derived implementation. The only "special" thing about that test is the random values generated are log-normally distributed, which has helped me catch quite a few bugs in the past.

Based on the value of my_value_eps(), we expect cpp_double_fp_backend< 4-byte float > to preserve 44 bits of accuracy, and cpp_double_fp_backend< 8-byte double > to preserve 102 bits. Note that I could not run these tests for float128.
After running a million rounds of each operation here are some observations:

• For cpp_double_fp_backend< 4-byte float >, division preserves only 43 bits of accuracy (vs. 44 expected) for both the algorithms.
• Both the algorithms preserve 102 bits of accuracy for cpp_double_fp_backend< 8-byte double >
• The accuracy is pretty much bottle-necked by division: add/sub preserve 104-bits easily, mul works fine for 103 bits, but division fails above 102 bits for both the algorithms.
• The only difference is in terms of accuracy between the two implementations seems to be in mul, where David Bailey's implementation preserves one extra bit (104 bits vs 103 bits) at the cost of speed, presumably.

Overall, it seems like the new algos are fast and reliable. However, as we have previously decided that my_value_eps() will reflect the accuracy of the type, this issue needs to be addressed. Possible ways forward include specializing numeric_limits<> for float but want to hear your thoughts before that.

The test I used is a somewhat of a kludge, but here is the link to it anyway.

[ckormanyos]

Hi Fahad (@sinandredemption) I am very glad to read your post and see the real dedicated effort to assess our precision at the bit-expected level.

could not run these tests for float128

Quick note here... If you are using GCC, then compile with -std=gnu++14 and the __float128 instantiation of the type should become available. Detection can use the compiler preprocessor querying on BOOST_HAS_FLOAT128 if you activate the -std=gnu++XX standard(s).

[ckormanyos]

it seems like the new algos are fast and reliable. However, as we have previously decided that my_value_eps() will reflect the accuracy of the type, this issue needs to be addressed.

Ummm... Just thinking out loud here, we could use traditional Briggs/Shoup for add/sub/mul. Maybe mix these with 1 round of Newton-Raphson for div/sqrt?

Maybe a combination will give us the best compromise on speed/accuracy/precision?

Cc: @cosurgi

[sinandredemption]

we could use traditional Briggs/Shoup for add/sub/mul. Maybe mix these with 1 round of Newton-Raphson for div/sqrt?

I should've mentioned it earlier: division fails 102 bits for cpp_double_fp_backend<8-byte double> for both the algorithms.

The only difference is in terms of accuracy between the two implementations seems to be in mul, where David Bailey's implementation preserves one extra bit (104 bits vs 103 bits) at the cost of speed, presumably.

[ckormanyos]

should've mentioned it earlier

Sorry Fahad (@sinandredemption) you had actually mentioned that. But I was not being a careful reader.

only difference is in terms of accuracy between the two implementations seems to be in mul, where David Bailey's implementation preserves one extra bit (104 bits vs 103 bits) at the cost of speed, presumably.

To me, that one is a real toss-up. In my experience, either of the choice of sets of these algorithms are still all several factors faster than our own cpp_bin_float at comparable bit counts. So if that one bit is the question, I would actually mix-and-match the algorithms to actually get that single bit of improvement back at the potential cost of speed.

Tough choice, but I'd go for the bit.

I would encourage you to run some cpp_double_fp-versus-cpp_bin_float performance tests, such as those here and see what you tend to think. I'd go for the bit back to 104.

[cosurgi]

I'd go for the bit back to 104.

I agree with this.

[ckormanyos]

So in light of 1.82 I've fired up this effort again.

I watched with great interest as an extreme value case was handled in Math, as i suspect this was hampering extreme values on cpp_double_double in the specfun tests.

I forgot a few details where we actually were, but i did just hammer out a specfun test locally. As it turns out now with the Math improvements for 1.82, the cpp_double_double<T> backend Multiprecision specfun tests now fail only 2 cases.

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

So I'll try to re-orient to 1.82 and finally get cpp_double_double<T> closer to boost readiness in this cycle prior to 1.83.

[ckormanyos]

P.S. Those failing cpp_double_double<T> tests prior to the fixes for 1.82 in Math were in regions of specfuntests involving ibeta() and friends.