1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/epoc32/include/stdapis/boost/math/special_functions/detail/series.hpp	Tue Mar 16 16:12:26 2010 +0000
     1.3 @@ -0,0 +1,52 @@
     1.4 +//  (C) Copyright John Maddock 2005.
     1.5 +//  Use, modification and distribution are subject to the
     1.6 +//  Boost Software License, Version 1.0. (See accompanying file
     1.7 +//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
     1.8 +
     1.9 +#ifndef BOOST_MATH_SERIES_INCLUDED
    1.10 +#define BOOST_MATH_SERIES_INCLUDED
    1.11 +
    1.12 +#include <cmath>
    1.13 +
    1.14 +#ifdef BOOST_NO_STDC_NAMESPACE
    1.15 +namespace std{ using ::pow; using ::fabs; }
    1.16 +#endif
    1.17 +
    1.18 +
    1.19 +namespace boost{ namespace math{ namespace detail{
    1.20 +
    1.21 +//
    1.22 +// Algorithm kahan_sum_series invokes Functor func until the N'th
    1.23 +// term is too small to have any effect on the total, the terms 
    1.24 +// are added using the Kahan summation method.
    1.25 +//
    1.26 +// CAUTION: Optimizing compilers combined with extended-precision
    1.27 +// machine registers conspire to render this algorithm partly broken:
    1.28 +// double rounding of intermediate terms (first to a long double machine
    1.29 +// register, and then to a double result) cause the rounding error computed
    1.30 +// by the algorithm to be off by up to 1ulp.  However this occurs rarely, and 
    1.31 +// in any case the result is still much better than a naive summation.
    1.32 +//
    1.33 +template <class Functor>
    1.34 +typename Functor::result_type kahan_sum_series(Functor& func, int bits)
    1.35 +{
    1.36 +   typedef typename Functor::result_type result_type;
    1.37 +   result_type factor = std::pow(result_type(2), bits);
    1.38 +   result_type result = func();
    1.39 +   result_type next_term, y, t;
    1.40 +   result_type carry = 0;
    1.41 +   do{
    1.42 +      next_term = func();
    1.43 +      y = next_term - carry;
    1.44 +      t = result + y;
    1.45 +      carry = t - result;
    1.46 +      carry -= y;
    1.47 +      result = t;
    1.48 +   }
    1.49 +   while(std::fabs(result) < std::fabs(factor * next_term));
    1.50 +   return result;
    1.51 +}
    1.52 +
    1.53 +} } } // namespaces
    1.54 +
    1.55 +#endif // BOOST_MATH_SERIES_INCLUDED