First cut at new rational_adaptor. #366
Conversation
|
If anyone (@ckormanyos ?) can take a look and review/double check this I'd really appreciate it. This removes the dependency to Boost.Rational. I still need to run some performance comparisons, but other than that it should be ready to go now. |
Thank you , John for diving into this topic and trying to remove the rational dependency --- a great step toward potential standalone.
|
[CI SKIP]
Will fix.
The definition in here: https://github.com/boostorg/multiprecision/blob/1b31fbea21bcf5b5b1e7d5ef4a8870aabed2161e/include/boost/multiprecision/traits/std_integer_traits.hpp is almost the same as
Not good, good catch, will fix shortly. |
Tidy up whitespace, and use nullptr where appropriate.
|
At the end of arithmetic operations, instead of simplifying by dividing with the gcd, the gcd function is called several times to simplify intermediate expressions. Does anybody know if that pays off for all sizes of integers or if there is a threshold. |
Good question, I'll experiment. In the mean time, here's the preliminary performance comparisons, mostly the new code does OK, but there are a couple of strange outliers I need to look into: Table 1.72. Operator *
Table 1.73. Operator *(int)
Table 1.74. Operator *(unsigned long long)
Table 1.75. Operator *=(unsigned long long)
Table 1.76. Operator +
Table 1.77. Operator +(int)
Table 1.78. Operator +(unsigned long long)
Table 1.79. Operator +=(unsigned long long)
Table 1.80. Operator -
Table 1.81. Operator -(int)
Table 1.82. Operator -(unsigned long long)
Table 1.83. Operator -=(unsigned long long)
Table 1.84. Operator /
Table 1.85. Operator /(int)
Table 1.86. Operator /(unsigned long long)
Table 1.87. Operator /=(unsigned long long)
Table 1.88. Operator construct
Table 1.89. Operator construct(unsigned long long)
Table 1.90. Operator construct(unsigned)
Table 1.91. Operator str
|
|
Update: there's a bug in the add-scalar routine that does an unnecessary gcd. Testing fixes now. |
Leaving aside experiments for a moment, there is a good theoretical reason for this: if we assume that gcd is O(N^2) then it is better to do an O(N) operation twice, than an O(2N) operation once, ie 2N^2 is better than (2N)^2. |
|
Much better for +- a scalar now: Operator +(int)
|
I coded this up quickly for multiply/divide, and while the code is noticeably simpler, it's also slower at all the bit-counts I tested:
|
|
Thank you for the benchmarks. Could you put them in this PR or on gist.github.com? |
The benchmark program is libs/multiprecision/performance/performance_test.cpp (though this does need updating and documenting a bit better). For the benefit of future googlers, the alternative multiply divide code is: |
What I've been keeping my eye out for, but haven't found yet, is a nice example of a rational number sequence - not too crazy hard to compute, but hard enough to properly exercise the arithmetic operators. Any one any ideas? |
Bernoulli numbers. But they diverge rapidly. There are some really neat results here. I wonder if something like |
|
Consider 10 terms of the expansion of The table is, for instance, here The result of the product is, I believe: |
What we have in our project (www.cgal,org) are things like orientation tests for plane/point, or sphere/point , both computing essentially signs of determinants of small matrices. |
Fix generic_interconvert.hpp for unsigned types. Update performance testing code to include testing with ::value_type. Start testing unsigned and checked integer types with rational_adaptor. Update arithmetic tests to test mixed arithmetic with ::value_type.
[CI SKIP]
|
OK, there are a couple of new Google benchmark programs added, one calculates Bernoulli(m) using the Louis Saalschütz explicit definition, the other calculates zeta(18) using Chris' equivalence. These results were run on Windows, with the unoptimized "generic" MPIR library, I'd expect gmp to perform much better once architecture specific optimizations are enabled: if someone would like to run these on other platforms that would be cool. zeta[18] calculation:
This is a relatively trivial calculation where the numbers don't grow too large, I suspect cpp_int's better memory usage for small values is the main driver here, though mpq_class is doing rather well - I need to check to see if it's eluding some temporaries somewhere that rational_adaptor isn't. Bernoulli calculation is a rather sterner test of arithmetic performance, here's what I'm seeing, all the result except those for cpp_rational are relative to cpp_rational's performance, so smaller is better:
|
In the expression-template case these do not return expression templates as a result of #175. However, the previous fix was inefficient and created unnecessary temporaries.
Speeds up construction from 2 components (rationals etc).
Correct ET operators for pre-C++17. Add test case for allocator usage.
|
OK big update to reduce temporary usage, memory allocations look much better now, here for Bernoulli[m] compiled with gcc-11 mingw:
Runtime performance looks rather better too, again mingw, and again performance figures are relative to cpp_rational, so smaller is better:
|
|
Wouldn't it make sense to cherry-pick your fix of the gcd into this branch. It may have an impact on performance, and without it computations might be plain wrong. |
|
@jzmaddock do you have a (vague) schedule for an integration into the develop branch of boost or a boost release. It's a question coming from @lrineau, the release manager of the CGAL project. |
Yep.
It's basically good to go, I just need to work on documentation a bit, and do a quick matrix determinant performance test. I just got a bit distracted by bug reports elsewhere ;) |
Done, and 3x3 determinant calculation added as a benchmark - it's a surprisingly tough calculation for rational values isn't it? Determinant calculation shows mpq significantly ahead of cpp_rational, this calculation mainly just hammers gcd I suspect, but if I have time, I'll investigate some more. As before, column for cpp_rational shows the actual time, the others are relative to that:
Bernoulli calculations are largely the same as before:
|
Includes updating all the performance results.
|
Barring CI SNAFU's, this is now ready to go. |
| operator*(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& a, number<B, ET>&& b) | ||
| { | ||
| b *= a; | ||
| return static_cast<number<B, et_on>&&>(b); |
There was a problem hiding this comment.
Seeing the CI results at https://drone.cpp.al/boostorg/multiprecision/292/179/2, can you comment on that static_cast? The type of b is number<B, ET>&&. Why should it cast nicely into number<B, et_on>&&?
| operator*(number<B, ET>&& a, const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& b) | ||
| { | ||
| a *= b; | ||
| return static_cast<number<B, et_on>&&>(a); |
There was a problem hiding this comment.
Same here: why et_on and not ET?
There was a problem hiding this comment.
Good catch: too much blind refactoring in one go!
I'm a bit concerned that didn't show up in local testing, also thinking I should be using std::foward here.
There was a problem hiding this comment.
std::forward and the static_cast (where right) are completely equivalent. The advantage of using std::forward is a better understanding of the intent of the developer.
Edit: fix the typo indent/intent
There was a problem hiding this comment.
Brilliant! Now that I see your commit e4827ce... of course it was not std::forward, but std::move, and now the code is a lot clearer to read! Thanks.
This is a first cut at a new rational_adaptor that doesn't use Boost.Rational and is better optimised for arbitrary precision integers. It passes our current tests, but more tests are required.