
Boost : 
From: Peder Holt (peder.holt_at_[hidden])
Date: 20051103 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 compiletime
> > 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,afloor(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
>
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk