Re: [boost] [mpl-cf] [RFC] MPL Real Numbers using Continued Fractions

Subject: Re: [boost] [mpl-cf] [RFC] MPL Real Numbers using Continued Fractions
From: Hal Finkel (half_at_[hidden])
Date: 2010-12-01 09:25:23

On Wed, 2010-12-01 at 09:16 +0000, Paul A. Bristow wrote:
> > -----Original Message-----
> > From: boost-bounces_at_[hidden] [mailto:boost-bounces_at_[hidden]]
> > On Behalf Of Hal Finkel
> > Sent: Tuesday, November 30, 2010 7:12 PM
> > To: boost_at_[hidden]
> > 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
> [phone]  07714330204
> pbristow_at_[hidden]
> 
> 
> 
> 
> _______________________________________________
> Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost