libs/math/example/series.cpp

   1 //  (C) Copyright John Maddock 2018.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #include <boost/math/tools/series.hpp>
   7 #include <boost/assert.hpp>
   8
   9 #include <iostream>
  10 #include <complex>
  11 #include <cassert>
  12
  13 //[series_log1p
  14 template <class T>
  15 struct log1p_series
  16 {
  17    // we must define a result_type typedef:
  18    typedef T result_type;
  19
  20    log1p_series(T x)
  21       : k(0), m_mult(-x), m_prod(-1) {}
  22
  23    T operator()()
  24    {
  25       // This is the function operator invoked by the summation
  26       // algorithm, the first call to this operator should return
  27       // the first term of the series, the second call the second
  28       // term and so on.
  29       m_prod *= m_mult;
  30       return m_prod / ++k;
  31    }
  32
  33 private:
  34    int k;
  35    const T m_mult;
  36    T m_prod;
  37 };
  38 //]
  39
  40 //[series_log1p_func
  41 template <class T>
  42 T log1p(T x)
  43 {
  44    // We really should add some error checking on x here!
  45    BOOST_ASSERT(std::fabs(x) < 1);
  46
  47    // Construct the series functor:
  48    log1p_series<T> s(x);
  49    // Set a limit on how many iterations we permit:
  50    boost::uintmax_t max_iter = 1000;
  51    // Add it up, with enough precision for full machine precision:
  52    return boost::math::tools::sum_series(s, std::numeric_limits<T>::epsilon(), max_iter);
  53 }
  54 //]
  55
  56 //[series_clog1p_func
  57 template <class T>
  58 struct log1p_series<std::complex<T> >
  59 {
  60    // we must define a result_type typedef:
  61    typedef std::complex<T> result_type;
  62
  63    log1p_series(std::complex<T> x)
  64       : k(0), m_mult(-x), m_prod(-1) {}
  65
  66    std::complex<T> operator()()
  67    {
  68       // This is the function operator invoked by the summation
  69       // algorithm, the first call to this operator should return
  70       // the first term of the series, the second call the second
  71       // term and so on.
  72       m_prod *= m_mult;
  73       return m_prod / T(++k);
  74    }
  75
  76 private:
  77    int k;
  78    const std::complex<T> m_mult;
  79    std::complex<T> m_prod;
  80 };
  81
  82
  83 template <class T>
  84 std::complex<T> log1p(std::complex<T> x)
  85 {
  86    // We really should add some error checking on x here!
  87    BOOST_ASSERT(abs(x) < 1);
  88
  89    // Construct the series functor:
  90    log1p_series<std::complex<T> > s(x);
  91    // Set a limit on how many iterations we permit:
  92    boost::uintmax_t max_iter = 1000;
  93    // Add it up, with enough precision for full machine precision:
  94    return boost::math::tools::sum_series(s, std::complex<T>(std::numeric_limits<T>::epsilon()), max_iter);
  95 }
  96 //]
  97
  98 int main()
  99 {
 100    using namespace boost::math::tools;
 101
 102    std::cout << log1p(0.25) << std::endl;
 103
 104    std::cout << log1p(std::complex<double>(0.25, 0.25)) << std::endl;
 105
 106    return 0;
 107 }