[mpl-cf] [RFC] MPL Real Numbers using Continued Fractions
Hello, A few months ago, I created an MPL extension type for representing real numbers. The indented use is for compile-time generation of floating-point constants, although it could also be used for MPL-based graph algorithms, etc. The primary design requirement was that almost any compiler would be able to completely inline the conversion to floating-point (which means no functional recursion can be used), and I think I've achieved that. The current implementation uses a modified version of MPL's vector type (modified to add the conversion to floating point, add "overloads", and to remove the typeof-based implementation variant (which is incompatible with the floating-point conversion)). I've posted the code to the Boost vault under the "Template Metaprogramming" section; the file name is boost-mpl-cf-20100531.zip Continued fractions were chosen for the representation because the associated algorithms are simple, and simple proved important for keeping down compilation times. Even so, the compilation times might be too long for practical use. You can decrease the compilation time by using a shorter maximum cf-length by defining BOOST_MPL_LIMIT_CF_SIZE to be 10 instead of 20. The current default is 20 to allow an effective precision near that of a double. Also, continued fractions, being "based on" rational numbers, can exactly represent any rational (for example, 1/3) so the computations can be more accurate than those using a floating-point representation. If you don't know what a continued fraction is, and would like to know, you can check out the wikipedia page: http://en.wikipedia.org/wiki/Continued_fraction The current implementation supports addition, subtraction, multiplication and division, comparison ops, some constants like pi, and I've included sin()/cos() as well. I've not yet written any real documentation, but the archive has a couple of test programs; the usage is pretty simple. Each test program has many individual tests, so its compilation time will be representative of a more-complicated computation, if you just need to do a few things, it would be much faster. If there is interest, I'll add support for generating numbers using metafunctions, write some documentation, jam files, etc. Let me know what you think, and I welcome suggestions, criticisms, etc. Also, I've currently tested using only the GNU and Intel compilers, so let me know whether it does or does not work when using something else. Thanks again, Hal P.S. I know that there is a mpl_math archive on the vault from back in 2007. I did not discover this until after I had almost finished writing my mpl-cf implementation, but, the mailing-list messages from the time indicate that it was too slow to be used practically, since I think the continued fractions algorithms are simpler, mpl-cf might be better, but I've not tested it to find out.
On 28/11/2010 22:30, Hal Finkel wrote:
Hello,
A few months ago, I created an MPL extension type for representing real numbers.
I find that idea very interesting.
The current default is 20 to allow an effective precision near that of a double. Also, continued fractions, being "based on" rational numbers, can exactly represent any rational (for example, 1/3) so the computations can be more accurate than those using a floating-point representation.
It would be nice to know what minimum value is necessary to provide correctly rounded results for IEEE754 floats and doubles.
On Mon, 2010-11-29 at 11:07 +0100, Mathias Gaunard wrote:
On 28/11/2010 22:30, Hal Finkel wrote:
Hello,
A few months ago, I created an MPL extension type for representing real numbers.
I find that idea very interesting.
The current default is 20 to allow an effective precision near that of a double. Also, continued fractions, being "based on" rational numbers, can exactly represent any rational (for example, 1/3) so the computations can be more accurate than those using a floating-point representation.
It would be nice to know what minimum value is necessary to provide correctly rounded results for IEEE754 floats and doubles.
This depends on the number: For a rational, you need only two integers ;) For a general real number, it depends on the rate of convergence. In my experience, you get about one decimal digit per integer coefficient, so 10 is sufficient for single precision, 20 for double precision (and long double too on an x86, since the extended type there is only 18 digits). The MPL containers generally come in max-size-increments of 10, but it would not be too hard to generate something else. -Hal
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
On Mon, 2010-11-29 at 10:24 -0500, Hal Finkel wrote:
On Mon, 2010-11-29 at 11:07 +0100, Mathias Gaunard wrote:
On 28/11/2010 22:30, Hal Finkel wrote:
Hello,
A few months ago, I created an MPL extension type for representing real numbers.
I find that idea very interesting.
The current default is 20 to allow an effective precision near that of a double. Also, continued fractions, being "based on" rational numbers, can exactly represent any rational (for example, 1/3) so the computations can be more accurate than those using a floating-point representation.
It would be nice to know what minimum value is necessary to provide correctly rounded results for IEEE754 floats and doubles.
This depends on the number: For a rational, you need only two integers ;)
Actually, ignore that statement. You only need two integers for a0 + 1/a1, in general you need more. For example, if you look at the test program (basic.cpp in the archive), it provides some examples (where the floating-point output is float, double, long double): 3/2 + 1/4 = [ 1 1 3] 1.75 1.75 1.75 1/7 + 1/15 = [ 0 4 1 3 2 2] 0.209524 0.20952380952381 0.209523809523809524 There are other continued fraction representations where any rational can be represented with a fixed number of integers (a0 + b0/(a1 + ...)), but the algorithms are a little more complicated, and I fear the compile times would increase. -Hal
For a general real number, it depends on the rate of convergence. In my experience, you get about one decimal digit per integer coefficient, so 10 is sufficient for single precision, 20 for double precision (and long double too on an x86, since the extended type there is only 18 digits). The MPL containers generally come in max-size-increments of 10, but it would not be too hard to generate something else.
-Hal
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
Actually, ignore that statement. You only need two integers for a0 + 1/a1, in general you need more. For example, if you look at the test program (basic.cpp in the archive), it provides some examples (where the floating-point output is float, double, long double):
3/2 + 1/4 = [ 1 1 3] 1.75 1.75 1.75
1/7 + 1/15 = [ 0 4 1 3 2 2] 0.209524 0.20952380952381 0.209523809523809524
There are other continued fraction representations where any rational can be represented with a fixed number of integers (a0 + b0/(a1 + ...)), but the algorithms are a little more complicated, and I fear the compile times would increase.
I was looking for something very similar to that, altough I believe compile times are to long for realistic problems. Especially the example with sine and cosine took approximately 10 minutes on my PC. But who knows? Maybe compilers are fast enough in some years... I have also two little questions: 1. Can one express for example 1.25e-10 in an easy manner, for example by simply specifying the exponent? 2. Can one use different representations of rational number, maybe boost::ratio? I guess this could be difficult to work with the elementary operations.
On Tue, 2010-11-30 at 16:11 +0100, Karsten Ahnert wrote:
Actually, ignore that statement. You only need two integers for a0 + 1/a1, in general you need more. For example, if you look at the test program (basic.cpp in the archive), it provides some examples (where the floating-point output is float, double, long double):
3/2 + 1/4 = [ 1 1 3] 1.75 1.75 1.75
1/7 + 1/15 = [ 0 4 1 3 2 2] 0.209524 0.20952380952381 0.209523809523809524
There are other continued fraction representations where any rational can be represented with a fixed number of integers (a0 + b0/(a1 + ...)), but the algorithms are a little more complicated, and I fear the compile times would increase.
I was looking for something very similar to that, altough I believe compile times are to long for realistic problems. Especially the example with sine and cosine took approximately 10 minutes on my PC. But who knows? Maybe compilers are fast enough in some years...
What compiler were you using (just curious)? I think it takes g++ about that long on my netbook too. The trig functions are slow, probably too slow for (most) practical application, but having real numbers themselves might be useful. Maybe someone with more experience could speed up the entire implementation ;)
I have also two little questions: 1. Can one express for example 1.25e-10 in an easy manner, for example by simply specifying the exponent?
1.25e-10 = [0, 8000000000], but you'd have to do the calculation beforehand and write: cf_c<0, 8000000000ll> or, for some extra compile time, write: boost::mpl::times< cf_c<1, 4> /* 1 + 1/4 */ , cf_c<0, 10000000000ll> /* 0 + 1/10^10 */ >::type
2. Can one use different representations of rational number, maybe boost::ratio? I guess this could be difficult to work with the elementary operations.
Not currently, but it is the same as writing boost::mpl::divides< cf_c<a>, cf_c<b> >::type, so it would be simple to add. I know that boost units has mpl rationals internally, so interfacing with them should also be possible. Or do you mean converting directly to rationals? That might work now (I've not tried it), but the value() function is templated on the return type. The advantage of working with continued fractions instead of rationals is that the individual integer elements of the continued fractions tend to stay small, even when the rational numerator and denominator get large (and would overflow your integer type). -Hal
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
On Mon, Nov 29, 2010 at 07:24, Hal Finkel <half@halssoftware.com> wrote:
On Mon, 2010-11-29 at 11:07 +0100, Mathias Gaunard wrote:
It would be nice to know what minimum value is necessary to provide correctly rounded results for IEEE754 floats and doubles.
This depends on the number: For a rational, you need only two integers ;) For a general real number, it depends on the rate of convergence.
Wouldn't it be bounded by the number of terms needed for a sufficiently-precise Phi, since it's the slowest-converging CF? ~ Scott
On Mon, Nov 29, 2010 at 9:11 PM, Scott McMurray <me22.ca+boost@gmail.com> wrote:
On Mon, Nov 29, 2010 at 07:24, Hal Finkel <half@halssoftware.com> wrote:
On Mon, 2010-11-29 at 11:07 +0100, Mathias Gaunard wrote:
It would be nice to know what minimum value is necessary to provide correctly rounded results for IEEE754 floats and doubles.
This depends on the number: For a rational, you need only two integers ;) For a general real number, it depends on the rate of convergence.
Wouldn't it be bounded by the number of terms needed for a sufficiently-precise Phi, since it's the slowest-converging CF?
For some inexplicable reason, the very fact that there's someone on this list who knows enough about math to ask that question makes me really, really happy. -- Dave Abrahams BoostPro Computing http://www.boostpro.com
On Mon, Nov 29, 2010 at 18:19, Dave Abrahams <dave@boostpro.com> wrote:
For some inexplicable reason, the very fact that there's someone on this list who knows enough about math to ask that question makes me really, really happy.
I'm glad :) For anyone else who wants to know just about as much as me: <http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/cfINTRO.html> :) I did a quick experiment, and it seems for phi, you need about 20 ones for a float and 40 ones for a double: typedef float T; // or double cout << setprecision(numeric_limits<T>::digits10+2); cout << (sqrt((T)5) + 1) / 2 << "\n\n"; int steps = 1; T value = 1; while (1 + 1/(1 + 1/value) != value) { ++steps; value = 1 + 1/value; cout << steps << ": " << value << "\n"; } I don't think that's practical to specify, as such :P However, perhaps you could limit it to some smaller number, but allow the last term to be a repeating specifier of some kind. It would be much more palatable to specify (ignoring that the constants would probably have to be int_c): typedef cf<1, repeat<1>> phi; typedef cf<0, 2> half; typedef cf<1, repeat<2>> root_2; typedef cf<1, repeat<8, 2>> half_of_root_5; This still isn't a full solution, though, since the one everyone will ask about can't be nicely shortened: typedef cf<3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3> pi; // hopefully enough... And it obviously makes implementing math and comparisons much harder. BTW, how do you avoid floating point annoyances when calculating the runtime value? I'd assume that cascading inversions would be a terrible idea. ~ Scott
On Tue, 2010-11-30 at 00:20 -0800, Scott McMurray wrote:
On Mon, Nov 29, 2010 at 18:19, Dave Abrahams <dave@boostpro.com> wrote:
For some inexplicable reason, the very fact that there's someone on this list who knows enough about math to ask that question makes me really, really happy.
I'm glad :)
For anyone else who wants to know just about as much as me: <http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/cfINTRO.html> :)
I did a quick experiment, and it seems for phi, you need about 20 ones for a float and 40 ones for a double:
typedef float T; // or double
cout << setprecision(numeric_limits<T>::digits10+2); cout << (sqrt((T)5) + 1) / 2 << "\n\n";
int steps = 1; T value = 1; while (1 + 1/(1 + 1/value) != value) { ++steps; value = 1 + 1/value; cout << steps << ": " << value << "\n"; }
I don't think that's practical to specify, as such :P
I had not thought of that, but yes, in the convention used by the implementation, that probably is the bound (+-1).
However, perhaps you could limit it to some smaller number, but allow the last term to be a repeating specifier of some kind. It would be much more palatable to specify (ignoring that the constants would probably have to be int_c):
In the current implementation, I provide cf_c<...> for that purpose.
typedef cf<1, repeat<1>> phi; typedef cf<0, 2> half; typedef cf<1, repeat<2>> root_2; typedef cf<1, repeat<8, 2>> half_of_root_5;
In case you did not notice, I provide these constants (except alf_of_root_5) to length 20 in boost/mpl/cf/constants.hpp I had thought it may be useful to allow for metafunctions which act as infinite stream generators, since then numbers like sqrt(2) could be represented "exactly." Care would be necessary to make sure that operations on such objects did not infinitely recurse. Also, it would be important to keep the performance impact in mind; When I added the integer-overflow checks into the arithmetic routines the compile times increased by a factor of 2, and those are relatively simple checks.
This still isn't a full solution, though, since the one everyone will ask about can't be nicely shortened:
typedef cf<3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3> pi; // hopefully enough...
There are other (more-complicated) representations in which pi has a simple form, but the everything else becomes more complicated... See, for example: http://en.wikipedia.org/wiki/Pi
And it obviously makes implementing math and comparisons much harder.
BTW, how do you avoid floating point annoyances when calculating the runtime value? I'd assume that cascading inversions would be a terrible idea.
That was also my original assumption, but I found it worked pretty well in practice. For one thing, I've chosen a convention such that all of the (non-zero) integers have the same sign, This means that there is no "catastrophic cancellation," and also makes implementing the comparison ops relatively simple (comparison is cheap compared to arithmetic). To be fair, my original use case was the generation of ratios of pi for use in FFT kernels, and as you've noted above, all of the integers in the pi expansion are fairly small (and there are a lot of 1s), so maybe it is not entirely representative, but I've yet to encounter a case where this has caused a problem. That having been said, I've not really tried to find such a case either. I can certainly think about this more carefully. -Hal
~ Scott _______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Hal Finkel Sent: Tuesday, November 30, 2010 7:12 PM To: boost@lists.boost.org Subject: Re: [boost] [mpl-cf] [RFC] MPL Real Numbers using Continued Fractions
That was also my original assumption, but I found it worked pretty well in practice. For one thing, I've chosen a convention such that all of the (non-zero) integers have the same sign, This means that there is no "catastrophic cancellation," and also makes implementing the comparison ops relatively simple (comparison is cheap compared to arithmetic).
To be fair, my original use case was the generation of ratios of pi for use in FFT kernels, and as you've noted above, all of the integers in the pi expansion are fairly small (and there are a lot of 1s), so maybe it is not entirely representative, but I've yet to encounter a case where this has caused a problem. That having been said, I've not really tried to find such a case either. I can certainly think about this more carefully.
I'm not sure what the OP really wanted to calculate, (if he really needed MPL calculation or MPL was just convenient) but I thought I would just add that if the need is something as simple as ratios of pi, or indeed fractions, one might like to consider using NTL www.shoup.net or GMP which have far higher precision reals (100+ decimal digits), enough to allow calculation of values to be pasted or cast into C++ builtin types, ensuring the accuracy of the float/double/long double value in C++. This method was used for the Boost.Math tables (many using continued fractions) and for the small number of constants, including pi, and its fractions and multiples. (The necessary patches and info to use NTL are in the Boost.Math package). Finally, John Maddock and I are currently giving some thought how to make these Boost.Math constants (and perhaps others - like a bigger collection of pi fractions for FFT for example) more 'user friendly' for wider use (and perhaps providing more of them). We will expose these ideas to discussion fairly soon. Paul --- Paul A. Bristow, Prizet Farmhouse, Kendal LA8 8AB UK +44 1539 561830 07714330204 pbristow@hetp.u-net.com
On Wed, 2010-12-01 at 09:16 +0000, Paul A. Bristow wrote:
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Hal Finkel Sent: Tuesday, November 30, 2010 7:12 PM To: boost@lists.boost.org Subject: Re: [boost] [mpl-cf] [RFC] MPL Real Numbers using Continued Fractions
That was also my original assumption, but I found it worked pretty well in practice. For one thing, I've chosen a convention such that all of the (non-zero) integers have the same sign, This means that there is no "catastrophic cancellation," and also makes implementing the comparison ops relatively simple (comparison is cheap compared to arithmetic).
To be fair, my original use case was the generation of ratios of pi for use in FFT kernels, and as you've noted above, all of the integers in the pi expansion are fairly small (and there are a lot of 1s), so maybe it is not entirely representative, but I've yet to encounter a case where this has caused a problem. That having been said, I've not really tried to find such a case either. I can certainly think about this more carefully.
I'm not sure what the OP really wanted to calculate, (if he really needed MPL calculation or MPL was just convenient) but I thought I would just add that if the need is something as simple as ratios of pi, or indeed fractions, one might like to consider using NTL www.shoup.net or GMP which have far higher precision reals (100+ decimal digits), enough to allow calculation of values to be pasted or cast into C++ builtin types, ensuring the accuracy of the float/double/long double value in C++.
Here's the backstory: I maintain a code which does a lot of 3-D FFTs (Fast Fourier Transforms), currently using either FFTW or Intel's MKL, and I saw the pair of articles by Vlodymyr Myrnyy (http://www.drdobbs.com/cpp/199500857 and http://www.drdobbs.com/cpp/199702312) where he claims that a technique whereby an FFT kernel generated using inlining with C++ template recursion outperforms FFTW by a significant margin (the size of the transform is known at compile time). To realize the performance benefit, the compiler must completely inline the floating-point constants used for the computation, and these are numbers of the form: sin(A*pi/B) and cos(A*pi/B) where A and B are integers. Vlodymyr uses a scheme to compute these at compile time using the Taylor expansion(s) of sin and cos based on having the compiler completely inline static template functions like (copied from the article): template<unsigned M, unsigned N, unsigned B, unsigned A> struct SinCosSeries { static double value() { return 1-(A*M_PI/B)*(A*M_PI/B)/M/(M+1) *SinCosSeries<M+2,N,B,A>::value(); } }; template<unsigned N, unsigned B, unsigned A> struct SinCosSeries<N,N,B,A> { static double value() { return 1.; } }; And the article shows benchmark results using g++ as the compiler, which does indeed always inline these calls completely, leaving only the floating-point constants in the assembly output. *but*, I discovered that several of the commercial compilers which I also need to have my code support, will not do the same: they do not inline those calls and leave the computation of the associated floating-point constants to runtime (thus killing any performance benefit associated with the technique). Not all compilers have a force_inline keyword (or similar extension), and arbitrarily increasing the "inlining depth" had a negative impact on other parts of the code (including the FFT code), and that is not really "portable." To have the technique work in a portable way, I needed a portable way to "force" the compiler to inline the computation of the floating-point constants (meaning function to inline needed to have a maximum call depth of 1), thus the motivation for a real-number MPL type; and it works (although evaluating the sin/cos functions is *very* expensive in terms of compile time). I don't think using MPL-CF for sin/cos computation for everyone, since the compiling of sin/cos is so slow, but I thought others may yet derive some benefit from having compile-time reals (perhaps for some optimization routine built into a Boost.Proto-based DSL evaluator, just a thought). If you just need to evaluate a few values using simple arithmetic (and comparisons are cheap), then MPL-CF seems to work quite well. Could you just use rationals? Sure, in theory, but continued fractions are much less susceptible to integer overflow.
This method was used for the Boost.Math tables (many using continued fractions) and for the small number of constants, including pi, and its fractions and multiples. (The necessary patches and info to use NTL are in the Boost.Math package).
Finally, John Maddock and I are currently giving some thought how to make these Boost.Math constants (and perhaps others - like a bigger collection of pi fractions for FFT for example) more 'user friendly' for wider use (and perhaps providing more of them). We will expose these ideas to discussion fairly soon.
I'm looking forward to hearing about it; I use Boost.Math constants on a regular basis. -Hal
Paul
--- Paul A. Bristow, Prizet Farmhouse, Kendal LA8 8AB UK +44 1539 561830 07714330204 pbristow@hetp.u-net.com
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
participants (6)
-
Dave Abrahams -
Hal Finkel -
Karsten Ahnert -
Mathias Gaunard -
Paul A. Bristow -
Scott McMurray