On 11/3/05, Cromwell Enage <sponage@yahoo.com> 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