Re: [boost] [Review Request] Multiprecision Arithmetic Library
Sorry, John my file was too big and admin rejected post. The candidate file is now 7-zip and should make it through. <snip>
I experimented with FFTs. We need your inputs if we *really*
want to make very high-precision available at this stage.
<snip>
maybe about 30 lines of code.
Nice, had no idea it could be implemented in such a small space!
Ooops, my bad. It was 100 lines. I guess I was in a daze by then. <snip>
So my gut feeling is to hold off for now, and maybe do things properly later (dynamic allocation, better std lib support etc)?
I agree. We should pull the reigns on it for now. We won't get the reliability we need if we try to do something fast. I could go on indefinitely with big-numbs, but we have to get that first version potentially reviewed. I should really get back to my book project soon. But I'm definitely in on any relevant review issues.
I think the existing ones work out quite nicely now thanks, I could use some ideas for integer and maybe rational number examples, but that's a whole other ball game ;-)
I'm not an integer guy. Maybe someone could help us out here. But what springs to my mind is: * Prime factorization * Greatest common denominator (GCD) * Primality testing for Mersenne Primes (has an exceptionally simple algo) * MD5, SHA-1, 2, etc. (but cool std::uint32_t reference apps already exist)<snip>
I resolved the excessive guard digit issue and reduced the guard digits to the necessary amount (and provide rationale in a code comment).
<snip>
Sure, or mail me a patch and I'll run the full tests (or you can - just cd into libs/multiprecision/test and do a "bjam --enable-specfun" and then wait for a long time!) Thanks, John.
OK. I attached my proposed version in /big_number/multiprecision/cpp_dec_float.hpp. The diffs include the removal of some non-needed constants, the trimming of max_digits10 and a minor improvement to the multiplication. Thanks, Chris.
OK. I attached my proposed version in /big_number/multiprecision/cpp_dec_float.hpp.
The diffs include the removal of some non-needed constants, the trimming of max_digits10 and a minor improvement to the multiplication.
Hmm, that fails a number of tests with error rates > 10^10 - see for example test_exp.cpp. Can you investigate? Thanks, John.
The diffs include the removal of some non-needed constants, the trimming of max_digits10 and a minor improvement to the multiplication.
Hmm, that fails a number of tests with error rates > 10^10 - see for example test_exp.cpp. Can you investigate?
Never mind - I think I have it - needed to replace uses of cpp_dec_float_max_digits10 with cpp_dec_float_total_digits10 which does what it says, where as cpp_dec_float_max_digits10 doesn't anymore. With that change the tests I've run pass so far.... but they're still running. Maybe we should rename cpp_dec_float_max_digits10 to something else? BTW, what's the logic behind the +1 in defining that one? Not that it makes much difference either which way. Cheers, John.
The diffs include the removal of some non-needed constants, the trimming of
max_digits10 and a minor improvement to the multiplication.
Hmm, that fails a number of tests with error rates > 10^10 - see for example test_exp.cpp. Can you investigate?
Never mind - I think I have it - needed to replace uses of cpp_dec_float_max_digits10 with cpp_dec_float_total_digits10 which does what it says, where as cpp_dec_float_max_digits10 doesn't anymore. With that change the tests I've run pass so far.... but they're still running.
Maybe we should rename cpp_dec_float_max_digits10 to something else? BTW, what's the logic behind the +1 in defining that one? Not that it makes much difference either which way.
Cheers, John.
Thanks John, I'm starting to lose it. My brain is getting sloppy. Sorry about all this confusion that I have caused. I still can't build the test suite with bjam. I will learn how to build the tests. There are still too many constants. I still need to eliminate a few. Best regards, Chris.
Hmm, that fails a number of tests with error rates > 10^10 - see for
example test_exp.cpp. Can you investigate?
Never mind - I think I have it - needed to replace uses of cpp_dec_float_max_digits10 with cpp_dec_float_total_digits10 which does what it says, where as cpp_dec_float_max_digits10 doesn't anymore. With that change the tests I've run pass so far.... but they're still running.
Maybe we should rename cpp_dec_float_max_digits10 to something else? Or eliminate it because it is redundant. I will report soon.
BTW, what's the logic behind the +1 in defining that one? Not that it makes much difference either which way. It's redundant and not needed, as addition of +1 should be done in std::numeric_limits anyway.
Cheers, John.
Thanks again, John. I got the tests running. What a nice test suite! I need to run the tests myself and know if the tests are OK before suggesting any changes. Is there a collected summary of the output files? Is this target report at the end of the test run of concern? ...failed updating 7 targets... ...skipped 21 targets... ...updated 574 targets... Best regards, Chris.
Is this target report at the end of the test run of concern?
...failed updating 7 targets... ...skipped 21 targets... ...updated 574 targets...
Hi John, You really have an impressive test bench worked out! It looks like some tests such as test_cyl_bessel_i, test_cyl_bessel_j, etc., some elliptic-stuff have problems with global divide. In other words, * operator / (char*, mp_type)I found it difficult to interpret the compiler's message because many macros are involved. It's probably a simple thing. If I can get your input on the correct interpretation of the test results, then I can verify test-passed and fix cpp_dec_float, thereby supplying you with the patch for your approval. Best regards, Chris.
Is this target report at the end of the test run of concern?
...failed updating 7 targets... ...skipped 21 targets... ...updated 574 targets...
Yes it means 7 tests failed, if you directed output to file then you can scroll up to find the failures - search for "..failed" and then scroll up to see the actual error messages.
You really have an impressive test bench worked out!
It looks like some tests such as test_cyl_bessel_i, test_cyl_bessel_j, etc., some elliptic-stuff have problems with global divide. In other words,
* operator / (char*, mp_type)I found it difficult to interpret the compiler's message because many macros are involved. It's probably a simple thing.
Hmm, I don't have any failures like that here - ah wait, you may (will actually) need to test against latest SVN Trunk to stand any chance of those building - sorry I should have made that clear. HTH, John.
If I can get your input on the correct interpretation of the test results, then I can verify test-passed and fix cpp_dec_float, thereby supplying you with the patch for your approval.
Thanks. I've just committed your polynomial example, plus some simple integer examples. Also tightened up some of the tests relating to cpp_dec_float - current SVN does not pass those tests - but the version you sent me with the modifications I mentioned (changing all *uses* of cpp_dec_float_max_digits10 to cpp_dec_float_total_digits10) does. I've also found a bug in my: template <unsigned D> cpp_dec_float& operator=(const cpp_dec_float<D>& f) and template <unsigned D> cpp_dec_float(const cpp_dec_float<D>& f) which needs fixing, but I'll let you complete whatever refactoring needs doing on those constants first. Cheers, John.
I've just committed your polynomial example, plus some simple integer examples.
Thanks for all that work, John. I had a good day of refactoring and testing and there is a lot to report. Cool! Thanks for properly scaling the argument to mysin(). I believe the expected output in the comment should be 8.4147098480789650665250232163029899962256306079837106567275170999e-1
Also tightened up some of the tests relating to cpp_dec_float - current SVN does not pass those tests - but the version you sent me with the modifications I mentioned (changing all *uses* of cpp_dec_float_max_digits10 to cpp_dec_float_total_digits10) does. I've also found a bug in my:
template <unsigned D> cpp_dec_float& operator=(const cpp_dec_float<D>& f) and template <unsigned D> cpp_dec_float(const cpp_dec_float<D>& f) which needs fixing, but I'll let you complete whatever refactoring needs doing on those constants first. Cheers, John.
I refactored the constants and got the tests running WRT trunk. Two tests did not compile at first: test_arithmetic.cpp (line 468) and test_float_io.cpp (line 204). The culprit is the compiler's failure to resolve frexp(Real, int64). The compiler required a plain integer argument. I believe that the compiler is correct (Visual Studio 2010, SP1, in this case). Do you also experience this problem? All tests except one pass with my refactored code. Unfortunately, test_float_io.cpp fails when round tripping from cpp_dec_float to a string with a large exponent back to cpp_dec_float. I believe that this problem can be solved because I identified a potentially incorrect rounding and/or digit inclusion when creating an exponent string in my original code. I have candidate changes in: * cpp_dec_float.hpp *trig.hpp * floating_point_examples.cpp And I would request that you please take a look at that frexp() issue mentioned above at your convenience. Would you like me to commit my changes or would you prefer to receive the files and commit them yourself? I will continue investigating the failure in test_float_io.cpp. It seems to be related to my code. Best regards, Chris
Cool! Thanks for properly scaling the argument to mysin(). I believe the expected output in the comment should be 8.4147098480789650665250232163029899962256306079837106567275170999e-1
Both sin and mysin(pi/4) output 7.0710678118654752440084436210484903928483593768847403658833986900e-01 for me - and my calculator agrees :-)
I refactored the constants and got the tests running WRT trunk. Two tests did not compile at first: test_arithmetic.cpp (line 468) and test_float_io.cpp (line 204). The culprit is the compiler's failure to resolve frexp(Real, int64). The compiler required a plain integer argument. I believe that the compiler is correct (Visual Studio 2010, SP1, in this case). Do you also experience this problem?
My mistake - I forgot to commit changes to default_ops.hpp - update again and that should get it. Sorry about that.
All tests except one pass with my refactored code. Unfortunately, test_float_io.cpp fails when round tripping from cpp_dec_float to a string with a large exponent back to cpp_dec_float. I believe that this problem can be solved because I identified a potentially incorrect rounding and/or digit inclusion when creating an exponent string in my original code.
OK good, I'd like that to work if possible, it seems like something users would want and/or expect?
I have candidate changes in: * cpp_dec_float.hpp *trig.hpp * floating_point_examples.cpp
Would you like me to commit my changes or would you prefer to receive the files and commit them yourself?
Other than the expected value for sin(pi/4) please do go ahead and commit! Many thanks, John.
I believe the expected output in the comment should be
8.4147098480789650665250232163029899962256306079837106567275170999e-1
Both sin and mysin(pi/4) output 7.0710678118654752440084436210484903928483593768847403658833986900e-01 for me - and my calculator agrees :-)
I guess I'm merely talking about the HTML documentation shown in the link below. But maybe it's already been fixed elsewhere: https://svn.boost.org/svn/boost/sandbox/big_number/libs/multiprecision/doc/h... The example code in the HTML documentation uses mysin(1). *Emphasize*: The argument is 1. Whereas the "expected output" in the code comment in the HTML documentation uses the number for the result of mysin(pi/4). *Emphasize*: The argument is pi/4. I apologize for giving you that code in an unscaled and unfinished form in the first place. Sorry. By the way, I looked at the compiler-generated assembly for that polynomial expansion using template expressions under x64 with VisualStudio2010. My jaw dropped to the floor! It's baffling what the compiler does with expression templates in combination with inlining, loop unrolling, constant folding, etc. I would like to to take a closer look at this and really comprehend the optimization potential. The state of C++ compiler technology today is truly impressive! I can hardly wait for VS2012 support for some of the new stuff.
I refactored the constants and got the tests running WRT trunk. Two tests did not compile at first: test_arithmetic.cpp (line 468) and test_float_io.cpp (line 204). The culprit is the compiler's failure to resolve frexp(Real, int64). The compiler required a plain integer argument. I believe that the compiler is correct (Visual Studio 2010, SP1, in this case). Do you also experience this problem?
My mistake - I forgot to commit changes to default_ops.hpp - update again and that should get it. Sorry about that.
Great! I'll check it out.
All tests except one pass with my refactored code. Unfortunately, test_float_io.cpp fails when round tripping from cpp_dec_float to a string with a large exponent back to cpp_dec_float. I believe that this problem can be solved because I identified a potentially incorrect rounding and/or digit inclusion when creating an exponent string in my original code.
OK good, I'd like that to work if possible, it seems like something users would want and/or expect?
I have candidate changes in: * cpp_dec_float.hpp *trig.hpp * floating_point_examples.cpp
OK. I committed: * cpp_dec_float.hpp with 23 changes---all straightforward, * trig.hpp with 2 changes I also changed the way that the digits of precision progress in Newton iteration for sqrt and atan. I needed this for high-precision. But it's maybe still not as good as a real "order of the delta" convergence check. But it doesn't matter at this stage in the game.
Other than the expected value for sin(pi/4) please do go ahead and commit! We just need to make the sample code agree with the comment in the documentation. But I didn't check anything in for this.
Many thanks, John. Thanks to you, too. I believe that we're really moving along well, now.
Best regards, Chris.
Both sin and mysin(pi/4) output 7.0710678118654752440084436210484903928483593768847403658833986900e-01 for me - and my calculator agrees :-)
I guess I'm merely talking about the HTML documentation shown in the link below. But maybe it's already been fixed elsewhere: https://svn.boost.org/svn/boost/sandbox/big_number/libs/multiprecision/doc/h...
The example code in the HTML documentation uses mysin(1). *Emphasize*: The argument is 1.
Doh! My bad, looking in the wrong place and forgot to change the argument to pi/4, will correct shortly.
By the way, I looked at the compiler-generated assembly for that polynomial expansion using template expressions under x64 with VisualStudio2010. My jaw dropped to the floor! It's baffling what the compiler does with expression templates in combination with inlining, loop unrolling, constant folding, etc. I would like to to take a closer look at this and really comprehend the optimization potential. The state of C++ compiler technology today is truly impressive! I can hardly wait for VS2012 support for some of the new stuff.
Nod, getting all those one-line expression templates expanded inline is key here, and yes, compilers are getting pretty impressive at handling this stuff! Cheers, John.
All tests except one pass with my refactored code.
Unfortunately, test_float_io.cpp fails when round tripping from cpp_dec_float to a string with a large exponent back to cpp_dec_float.
OK good, I'd like that to work if possible, it seems like something users would want and/or expect?
What I can't figure out, though, is if a floating-point type actually *should* round-trip back from a string in scientific format with (max_digits10 + 1) precision. Remember, that kind of string has (max_digits10 + 1) after the decimal point plus one digit in the mantissa. That would be like asking 8-byte double to be exact when round-tripping with 17 decimal digits of precision since the scientific notation would have 1+16 decimal digits. Would a *relative-to-epsilon* kind of test, maybe, be more correct (correct-er) in test_float_io.cpp? Is *correcter* a word other than the one who corrects? template <class T> void do_round_trip(const T& val, std::ios_base::fmtflags f) { std::stringstream ss; ss.flags(f); ss << std::setprecision(std::numeric_limits<T>::digits10) << val; T new_val = ss.str(); BOOST_CHECK_CLOSE(new_val, val, std::numeric_limits<T>::epsilon()); new_val = val.str(0, f); BOOST_CHECK_CLOSE(new_val, val, std::numeric_limits<T>::epsilon()); } Best regards, Chris.
OK good, I'd like that to work if possible, it seems like something users would want and/or expect?
What I can't figure out, though, is if a floating-point type actually *should* round-trip back from a string in scientific format with (max_digits10 + 1) precision. Remember, that kind of string has (max_digits10 + 1) after the decimal point plus one digit in the mantissa. That would be like asking 8-byte double to be exact when round-tripping with 17 decimal digits of precision since the scientific notation would have 1+16 decimal digits.
I'm not sure I follow you - what the std says is: static constexpr int digits10; 11 Number of base 10 digits that can be represented without change. 12 Meaningful for all specializations in which is_bounded != false. static constexpr int max_digits10; 13 Number of base 10 digits required to ensure that values which differ are always differentiated. Which I take as meaning that if you print at least max_digits10 digits then you should always be able to get back to the same unique value. In other words we get to choose the value of max_digits10 so that we can round trip when exactly that many digits are printed out. In the case of cpp_dec_float, given that there is no binary to decimal conversion going on, I don't see the problem (but of course that doesn't mean there isn't one)? Oh wait.... just looked at the code and I see the issue - you set max_digits10 to digits10+1 which isn't enough, should be set to cpp_dec_float_total_digits10 which is the largest number of digits possible in the number right? Making that change causes the test to pass, so committing, but please advise if that's upping max_digits10 too far ;-) Cheers, John.
<snip>
What I can't figure out, though, is if a floating-point
type actually *should* round-trip back from a string in scientific format with (max_digits10 + 1) precision.
Wait, you are right. And I finally have to solve this problem correctly! Paul Bristow and I wiere beginning to identify this problem last year. But I did not identify what may be the right solution. (See below.) <snip>
I'm not sure I follow you - what the std says is: static constexpr int digits10; 11 Number of base 10 digits that can be represented without change. 12 Meaningful for all specializations in which is_bounded != false.
Could you please tell me where this is in ISO 14882:2011? I can't find it. <snip>
Which I take as meaning that if you print at least max_digits10 digits then you should always be able to get back to the same unique value.
I agree. <snip>
Oh wait.... just looked at the code and I see the issue - you set max_digits10 to digits10+1 which isn't enough, should be set to cpp_dec_float_total_digits10 which is the largest number of digits possible in the number right? Making that change causes the test to pass, so committing, but please advise if that's upping max_digits10 too far ;-) Cheers, John.
Wait John. It took me a long time to finally figure this out. And it gets back to some of our original work with Paul Bristow. I believe that we actually need to test for equality to within 1ULP, whereby the ULP in this case is the least significant, potentially rounded and carried base-10 "digit" in the limb containing the ULP of least precision in cpp_dec_float. I need to actually add more intelligence to template <unsigned Digits10> int cpp_dec_float<Digits10>::cmp_data(const array_type& vd) const. The current equality case includes *all* the guard limbs, even though some are insignificant. I need to analyze the real number of identical base-10 digits in LHS and RHS, take rounding into account and test for equality not until the end of the limbs, but until one max_digits10, potentially rounded and carried from (max_digits10 + 1). Can you follow me? Do you agree with this? It shouldn't be too difficult to implement and I don't anticipate significant performance loss due to analyzing two limbs in the array in base-10, possibly in conjunction with lexical conversion. And I would really like to give this a try. With your permission, may I try this? Am I on the right track? I believe that we are getting cpp_dec_float to behave very much like a built-in type. I know I keep adding more time to the development cycle. But I feel that we are still discovering improvements. In my opinion, we are getting close. Best regards, Chris.
I'm not sure I follow you - what the std says is: static constexpr int digits10; 11 Number of base 10 digits that can be represented without change. 12 Meaningful for all specializations in which is_bounded != false.
Could you please tell me where this is in ISO 14882:2011? I can't find it.
Section 18.3.2.4 para 13 (I'm looking at the final draft though so there may be detail changes to the final std I guess).
<snip>
Which I take as meaning that if you print at least max_digits10 digits then you should always be able to get back to the same unique value.
I agree.
<snip>
Oh wait.... just looked at the code and I see the issue - you set max_digits10 to digits10+1 which isn't enough, should be set to cpp_dec_float_total_digits10 which is the largest number of digits possible in the number right? Making that change causes the test to pass, so committing, but please advise if that's upping max_digits10 too far ;-) Cheers, John.
Wait John.
It took me a long time to finally figure this out. And it gets back to some of our original work with Paul Bristow. I believe that we actually need to test for equality to within 1ULP, whereby the ULP in this case is the least significant, potentially rounded and carried base-10 "digit" in the limb containing the ULP of least precision in cpp_dec_float.
I need to actually add more intelligence to
template <unsigned Digits10> int cpp_dec_float<Digits10>::cmp_data(const array_type& vd) const.
The current equality case includes *all* the guard limbs, even though some are insignificant. I need to analyze the real number of identical base-10 digits in LHS and RHS, take rounding into account and test for equality not until the end of the limbs, but until one max_digits10, potentially rounded and carried from (max_digits10 + 1).
Can you follow me? Do you agree with this? It shouldn't be too difficult to implement and I don't anticipate significant performance loss due to analyzing two limbs in the array in base-10, possibly in conjunction with lexical conversion. And I would really like to give this a try. With your permission, may I try this? Am I on the right track?
I don't know ;-) Here's what I think - firstly lets "Keep It Simple" - and if the digits are present in the internal representation, then why *shouldn't* they be used for comparisons and printing out? Let me give you an example - if there are guard digits which aren't used for comparisons - could you not get into a situation where two numbers which compare equal lead to different results after some sequence of arithmetic operations? IMO that would not be good, and at the very least could lead to all kinds of confusions. Consider what GMP's mpf_t does - there are hidden guard digits that you don't see, which do take part in comparisons and all arithmetic operations, on the other hand they don't provide a way to print those out which makes round tripping impossible - which IMO is very annoying. On the other hand, MPFR takes a much less loose approach, and rounds everything to *exactly* the requested precision after each step, round tripping is also possible in this case. The equivalent in cpp_dec_float would be to have no guard digits at all, and round after each operation to exactly Digits10 decimal places. While there may be use cases for that, I think we want to Keep It Simple for now and not go down that road. So.... the situation we have now seems to be path of least resistance. It's somewhat inelegant having max_digits10 so much larger than digits10, but it seems to me that nothing really bad can happen here - just maybe that serialized numbers will occupy a touch more space on disk than expected? Cheers, John.
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of John Maddock Sent: Monday, April 16, 2012 11:38 AM To: boost@lists.boost.org Subject: Re: [boost] [Review Request] Multiprecision Arithmetic Library
I'm not sure I follow you - what the std says is: static constexpr int digits10; 11 Number of base 10 digits that can be represented without change. 12 Meaningful for all specializations in which is_bounded != false.
Could you please tell me where this is in ISO 14882:2011? I can't find it.
Section 18.3.2.4 para 13 (I'm looking at the final draft though so there may be detail changes to
the final
std I guess).
<snip>
Which I take as meaning that if you print at least max_digits10 digits then you should always be able to get back to the same unique value.
I agree.
<snip>
Oh wait.... just looked at the code and I see the issue - you set max_digits10 to digits10+1 which isn't enough, should be set to cpp_dec_float_total_digits10 which is the largest number of digits possible in the number right? Making that change causes the test to pass, so committing, but please advise if that's upping max_digits10 too far ;-) Cheers, John.
Wait John.
It took me a long time to finally figure this out. And it gets back to some of our original work with Paul Bristow. I believe that we actually need to test for equality to within 1ULP, whereby the ULP in this case is the least significant, potentially rounded and carried base-10 "digit" in the limb containing the ULP of least precision in cpp_dec_float.
I need to actually add more intelligence to
template <unsigned Digits10> int cpp_dec_float<Digits10>::cmp_data(const array_type& vd) const.
The current equality case includes *all* the guard limbs, even though some are insignificant. I need to analyze the real number of identical base-10 digits in LHS and RHS, take rounding into account and test for equality not until the end of the limbs, but until one max_digits10, potentially rounded and carried from (max_digits10 + 1).
Can you follow me? Do you agree with this? It shouldn't be too difficult to implement and I don't anticipate significant performance loss due to analyzing two limbs in the array in base-10, possibly in conjunction with lexical conversion. And I would really like to give this a try. With your permission, may I try this? Am I on the right track?
I don't know ;-)
Here's what I think - firstly lets "Keep It Simple" - and if the digits are present in the internal representation, then why *shouldn't* they be used for comparisons and printing out?
Let me give you an example - if there are guard digits which aren't used for comparisons - could you not get into a situation where two numbers which compare equal lead to different results after some sequence of arithmetic operations? IMO that would not be good, and at the very least could lead to all kinds of confusions.
Consider what GMP's mpf_t does - there are hidden guard digits that you don't see, which do take
comparisons and all arithmetic operations, on the other hand they don't provide a way to print
Like the many confusions caused by digits10 before max_digits10 was introduced. Two numbers that could be distinguished displayed the same, which was very confusing, especially to the majority whose grip on floating-point was less than secure ;-) part in those out
which makes round tripping impossible - which IMO is very annoying.
On the other hand, MPFR takes a much less loose approach, and rounds everything to *exactly* the requested precision after each step, round tripping is also possible in this case. The equivalent in cpp_dec_float would be to have no guard digits at all, and round after each operation to exactly Digits10 decimal places. While there may be use cases for that, I think we want to Keep It Simple for now and not go down that road.
Agreed.
So.... the situation we have now seems to be path of least resistance. It's somewhat inelegant having max_digits10 so much larger than digits10, but it seems to me that nothing really bad can happen here - just maybe that serialized numbers will occupy a touch more space on disk than expected?
And that any display using max_digits10 will look a bit odd, having so many potentially insignificant digits. I can't quite see why more than two or three guard digits should be needed? If we compare with the built-in floating-point, then only two are needed from Kahan's formula even at 128-bit. So lots more sort of feels untidy? But John is right in saying that nothing too evil can happen from having too many. Paul --- Paul A. Bristow, Prizet Farmhouse, Kendal LA8 8AB UK +44 1539 561830 07714330204 pbristow@hetp.u-net.com
I can't quite see why more than two or three guard digits should be needed? If we compare with thebuilt-in floating-point, then only two are needed from Kahan's formula even at 128-bit. So lots more sort of feels untidy?
Kahan is still the master of us all in float.
But John is right in saying that nothing too evil can happen from having too many. Paul
Hi Paul, Thanks for your insight. The cpp_dec_float back-end is based on my original design which did have some design trade-offs. The limbs of cpp_dec_float are 32-bit unsigned integers, in each one of which 8 decimal digits are present. However, the most significant limb and the second least significant limb can "slide" anywhere from 1...99,999,999. In addition, cpp_dec_float loses yet another limb during multiply because I only multiply half of the triangle and the carry from below simply gets lost. The digits in cpp_dec_float inherit these characteristics. They slide around within 3 extra limbs of slop, but do adhere to the requirements of the back-end and have fair performance at low precision. (We have discussed this at length.) I don't like the internal shifting and, if I ever write another big-float back-end, I will use base-2 and tight digit control such as that in MPFR. At times, I regret using base-10. Other times, when debugging, I'm glad to have it visible in recognizable form in the data array of limbs. Still, my next one (if there ever is such a thing) will be base-2. Thanks for keeping in touch, Paul. Best regards, Chris.
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Christopher Kormanyos Sent: Monday, April 16, 2012 8:52 PM To: boost@lists.boost.org Subject: Re: [boost] [Review Request] Multiprecision Arithmetic Library
I can't quite see why more than two or three guard digits should be needed? If we compare with thebuilt-in floating-point, then only two are needed from Kahan's formula even at 128-bit. So lots more sort of feels untidy?
Kahan is still the master of us all in float.
But John is right in saying that nothing too evil can happen from having too many.
Thanks for your insight. The cpp_dec_float back-end is based on my original design which did have some design trade-offs.
The limbs of cpp_dec_float are 32-bit unsigned integers, in each one of which 8 decimal digits are
present.
However, the most significant limb and the second least significant limb can "slide" anywhere from 1...99,999,999. In addition, cpp_dec_float loses yet another limb during multiply because I only multiply half of the triangle and the carry from below simply gets lost.
The digits in cpp_dec_float inherit these characteristics. They slide around within 3 extra limbs of slop, but do adhere to the requirements of the back-end and have fair performance at low precision. (We have discussed this at length.)
OK - thanks - so there is a valid reason why it 'feels untidy'. So we should stick with that and get the code out to the users. (Perhaps document in case anyone is curious/puzzled by the number of guard digits?) And anyone who wants less guard digits can always ask for digits10+3, say? Paul --- Paul A. Bristow, Prizet Farmhouse, Kendal LA8 8AB UK +44 1539 561830 07714330204 pbristow@hetp.u-net.com
<snip>
So we should stick with that and get the code out to the users.
Absolutely! Bring the numbers to the people. <snip>
(Perhaps document in case anyone is curious/puzzled by the number of guard digits?) I agree.
And anyone who wants less guard digits can always ask for digits10+3, say? Yes.
<snip> Thanks, Paul. Best regards, Chris.
Here's what I think - firstly lets "Keep It Simple" - and if the digits
are present in the internal representation, then why *shouldn't* they be used for comparisons and printing out?
Nod. I looked at a few possibilities today. But, as usual, I had underestimated the difficulty.
Let me give you an example - if there are guard digits which aren't used for comparisons - could you not get into a situation where two numbers> which compare equal lead to different results after some sequence of arithmetic operations? IMO that would not be good, and at the very least could lead to all kinds of confusions.
Consider what GMP's mpf_t does - there are hidden guard digits that you don't see, which do take part in comparisons and all arithmetic operations, on the other hand they don't provide a way to print those out which makes round tripping impossible - which IMO is very annoying.
On the other hand, MPFR takes a much less loose approach, and rounds everything to *exactly* the requested precision after each step, round tripping is also possible in this case.
The equivalent in cpp_dec_float would be to have no guard digits at all, and round after each operation to exactly Digits10 decimal places. While there may be use cases for that, I think we want to Keep It Simple for now and not go down that road.
This is too costly in base-10. If I (we) ever write another back-end, base-2 might be a better way to go.
So.... the situation we have now seems to be path of least resistance. It's somewhat inelegant having max_digits10 so much larger than digits10, but it seems to me that nothing really bad can happen here - just maybe that serialized numbers will occupy a touch more space on disk than expected? Cheers, John.
You are right again. And if I keep fiddling around with it, then I'm just wasting everyone's time. I designed the e_float back-end (the basis of cpp_dec_float) and tested it for years with its present strategy. I should just stop fiddling around and simply rely on the test suite and adherence to the back-end requirements. And we are getting very close to this, if not already there. I think we're done with the digits10, max_digits10 thing now. Thanks for your help. Best regards, Chris.
Oh wait.... just looked at the code and I see the issue - you set max_digits10 to digits10+1 which isn't enough, should be set to cpp_dec_float_total_digits10 which is the largest number of digits possible in the number right? Making that change causes the test to pass, so committing, but please advise if that's upping max_digits10 too far ;-)
Cheers, John.
Thanks, John. I have been doing quite a bit of testing and reviewing, concentrating on the implementation of cpp_dec_float. Everything is going quite well! It's a lot of fun playing with this impressive analysis system! I have several change requests which I have implemented in my SVN checkout and thoroughly tested with the test suite. John, if you would like me to commit the changes described below, I can. If you would prefer to wait until a time in the future OK. I can also provide you with source for your approval. I would like the changes now, even though we are rapidly approaching a potential review. After this, I promise to stop fiddling around with cpp_dec_float because the quality of the implementation and testing of is, in my opinion, more than adequate for a first release. 1) The version of cpp_dec_float in the sandbox has a lower-limit on the number of decimal digits of 30. I would like to reduce the lower limit to 9, thereby being sure that the internal data array retains at least two limbs. Our previous tests at 17 and 30 digits both simply tested at 30 digits. This change will smoothen the transition from built-in types to arbitrary precision. 2) Thereby, add test cases at 9 and 18 digits to the following: sqrt, sin, asin, atan, exp, log, sinh, pow, cyl_bessel_i, zeta. 3) Include one additional test at 1000 digits to sqrt, thereby extending the test data in the test case accordingly. 4) I could not find a particular dependency in the Jamfile. I need to add it for testing square root with cpp_dec_float. 5) Remove the unnecessary temporary storage array in the internal multiplication routine of cpp_dec_float. (Unfortunately operator+= still has one that I must remove at a much later time due to higher code complexity). 6) Cleaned up all the constants in cpp_dec_float. There were constants with prefixes of "ef_", "mp_" and "cpp_dec_float_", originating from 3 different development epochs. I have unified them all to "cpp_dec_float_". I am just finishing up the final testing of 564 targets. All the tests pass. Whould you like the changes committed? Best regards, Chris.
I am just finishing up the final testing of 564 targets.
All the tests pass. Would you like the changes committed?
Sure, go for it! Many thanks, John.
My changes have been committed. Thanks for providing time for this refactoring. If you experience problems with testing or would like any additional modifications, please let me know. I will send bibliographic references on the weekend to provide a literary context for the documentation. Barring any unforeseen problems, then, I don't anticipate significant changes from my side before potential review. Thanks as always and best regards, Chris.
I will send bibliographic references on the weekend
to provide a literary context for the documentation.
Hi John, Here is the Bibliography for much of my R&D on big numbs. Don't know if you want it for the docs. But here it is anyway. Best regards, Chris J. Arndt and C. Haenel "Pi Unleashed" Springer, Heidelberg, 2001 ISBN 3-540-66572-2 F. Bornemann, D. Laurie, S. Wagon and J. Waldvogel, "The SIAM 100-Digit Challenge: A Study in High-Accuracy Numerical Computing" Society for Industrial and Applied Mathematics SIAM (TM), Philadelphia, 2004. ISBN 0-898-71561-X J. M. Borwein and P. Borwein "Pi and the AGM: A Study in Analytic Number Theory and Computational Complexity" John Wiley, New York, 1998 ISBN 0-471-83138-7, Paperback R. Brent and P. Zimmermann "Modern Computer Arithmetic" Cambridge University Press, Cambridge, 2011 ISBN 978-0-521-19469-3 L. Fousse, G. Hanrot, V. Lefevre, P. Pelissier and P. Zimmermann "MPFR: A Multiple-Precision Binary Floating-Point Library With Correct Rounding" ACM ACM Trans. Math. Soft., Vol. 33, Issue 2 ACM, 2007 A. Gil, J. Segura and N. M. Temme, "Numerical Methods for Special Functions" Society for Industrial and Applied Mathematics SIAM (TM), Philadelphia, 2007 ISBN 978-0-898-71634-4 J. M. Muller "Elementary Functions: Algorithms and Implementation", 2nd Edn. Birkhaeuser, Boston, 2006 ISBN 978-0-817-64372-0 J. M. Muller, N. Brisebarre, F. de Dinechin, C. M. Jeannerod, V. Lefevre, G. Melquiond, N. Revol, D. Stehle and T. Torres "Handbook of Floating-Point Arithmetic" Birkhaeuser, Boston, 2010 ISBN 978-0-817-64704-9 C. M. Kormanyos "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations" ACM ACM Trans. Math. Soft., Vol. 37, Issue 4 ACM, 2011 D. E. Knuth "The Art of Computer Programming, Volume 3 Seminumerical Algorithms", 3rd edn. Addison Wesley, 1998
I am just finishing up the final testing of 564 targets. All the tests pass. Would you like the changes committed?
Sure, go for it! Many thanks, John.
My changes have been committed. Thanks for providing time for this refactoring.
Hello John. Matthieu's work with <complex.hpp> revealed two trivial errors in my refactor of cpp_dec_float.hpp. I also performed two or three additional cosmetic modifications. The tests all pass. I committed cpp_dec_float.hpp to the sandbox on Sunday. Sorry, I should have asked permission first. My bad. Best regards, Chris.
Matthieu's work with <complex.hpp> revealed two trivial errors in my refactor of cpp_dec_float.hpp. I also performed two or three additional cosmetic modifications.
The tests all pass. I committed cpp_dec_float.hpp to the sandbox on Sunday.
Sorry, I should have asked permission first. My bad.
No problems, cheers, John.
The tests all pass. I committed cpp_dec_float.hpp
to the sandbox on Sunday.
Sorry, I should have asked permission first. My bad.
No problems, cheers, John.
Hi John. Thanks. I would like to report another success with the proposed Boost.Multiprecision. I added partial adherence to the backend requirements for an mp_number to my fixed-point class and verified the computation of the area of a circle of radius (123 / 100) with it. The beautiful set of numbers below shows the results of the area for: 1. fixed-point Q15.16 2. built-in float of 4 bytes 3. built-in double of 8 bytes 4. cpp_dec_float_100 4.7528 4.7529159 4.7529155256159976 4.752915525615998190470133174563559913501897584314597596552993673698486406 Your wunderful creation provides unprecedented potential for generic mathematical calculations! Best regards, Chris.
The tests all pass. I committed cpp_dec_float.hpp
to the sandbox on Sunday. Sorry, I should have asked permission first. My bad.
No problems, cheers, John.
I added partial adherence to the backend requirements for an mp_number to my fixed-point class and verified the computation of the area of a circle of radius (123 / 100) with it.
... And thereby plugged this augmented fixed-point type into mp_number, of course. Best regards, Chris.
Hi John. Thanks.
I would like to report another success with the proposed Boost.Multiprecision.
I added partial adherence to the backend requirements for an mp_number to my fixed-point class and verified the computation of the area of a circle of radius (123 / 100) with it.
Your wunderful creation provides unprecedented potential for generic mathematical calculations!
Thanks, IMO it works quite well for that kind of mash-up - as the minimal backend requirements are fairly small - of course an "optimal" solution may well require more work, but then it always does I guess! John.
participants (4)
-
Christopher Kormanyos -
John Maddock -
JOHN MADDOCK
-
Paul A. Bristow