
Boost : 
Subject: Re: [boost] [mplcf] [RFC] MPL Real Numbers using Continued Fractions
From: Hal Finkel (half_at_[hidden])
Date: 20101201 09:25:23
On Wed, 20101201 at 09:16 +0000, Paul A. Bristow wrote:
> > Original Message
> > From: boostbounces_at_[hidden] [mailto:boostbounces_at_[hidden]]
> > On Behalf Of Hal Finkel
> > Sent: Tuesday, November 30, 2010 7:12 PM
> > To: boost_at_[hidden]
> > Subject: Re: [boost] [mplcf] [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
> (nonzero)
> > 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 3D 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 floatingpoint 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
floatingpoint 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 floatingpoint 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 floatingpoint constants (meaning function to inline
needed to have a maximum call depth of 1), thus the motivation for a
realnumber MPL type; and it works (although evaluating the sin/cos
functions is *very* expensive in terms of compile time).
I don't think using MPLCF 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 compiletime reals (perhaps for some
optimization routine built into a Boost.Protobased DSL evaluator, just
a thought). If you just need to evaluate a few values using simple
arithmetic (and comparisons are cheap), then MPLCF 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_at_[hidden]
>
>
>
>
> _______________________________________________
> Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk