Re: [boost] [MPL.Math] updates for double

From: Peder Holt (peder.holt_at_[hidden])
Date: 2005-11-03 13:23:57

On 11/3/05, Cromwell Enage <sponage_at_[hidden]> wrote:
> --- Peder Holt wrote:
> > I tried implementing a power using fractions, but it
> > turned out that for larger exponents its convergence
> > rate was very poor.  (I needed ~100 recursions for
> > something like 10^20).
>
> Yeah, that stinks.
>
> > What we should do, is add a specialization for power
> > with integral exponent.  This is a very fast
> > algorithm (also in respect to compile-time
> > efficiency).
>
> The current implementation calls integral_power if it
> detects that exponent::tag is the same as
> integral_c_tag, unsigned_big_integral_tag, or
> big_integral_tag.

Great! Most of the job is done already, then.

>
> > What we could do, is to split power(z,a) into two:
> > power(z,floor(a))*power(z,a-floor(a));
>
> Ah, a specialization for fractional exponents.  We
> could use integral_part for the floor function, but if
> a itself is very large (>= 2^30), then
> minus<a,integral_part<a> > will yield a mixed_number,
> not a double.  fractional_part<a> returns a rational
> regardless of magnitude.
>
> What I'll do is:
> 1) Create an internal metafunction in the double_::aux
> namespace that does what fractional_part does now.
> 2) Change fractional_part so that it returns a double.
> 3) Implement numerator and denominator specializations
> that use the metafunction in 1).
>
> Then, with fractional_power as a template nested
> within power_impl, we can return:
>
>     times<
>         integral_power<z,integral_part<a> >
>       , fractional_power<z, fractional_part<a> >
>     >

Hmm. Looking closer at the issue, it seems that the continued fraction
approximation of power is poor no matter what :(
Evaluating the expression
pow(20.5,0.98)
A continuing fractions implementation needs 73recursions to get the
same accuracy as the runtime equivalent for double.
The series:
z^a == Sum[(Log[z]^k/k!) a^k, {k, 0, Infinity}]
requires 27 recursions, but the code is very simple:

template <typename Tag1,typename Tag2>
struct power_impl {
   template<typename LogZ,typename A,typename I,typename SeriesCount>
   struct series {
   private:
            typedef typename eval_if<
                        less<typename I,SeriesCount>
                      , series<LogZ,A,typename I::next,SeriesCount>
                      , mpl::int_<1>
                    >::type
                    next_term;
   public:
           typedef typename
plus<1,times<divides<LogZ,I>,a,next_term>::type type;
   };
};

Regards,
Peder

>
>                               Cromwell D. Enage
>
>
> __________________________________________________
> Do You Yahoo!?
> Tired of spam?  Yahoo! Mail has the best spam protection around
> http://mail.yahoo.com
> _______________________________________________
> Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
>