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
log1p_serieslog1p_series20 log1p_series(T x)
21 : k(0), m_mult(-x), m_prod(-1) {}
22
operator ()log1p_series23 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>
log1p(T x)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
log1p_serieslog1p_series63 log1p_series(std::complex<T> x)
64 : k(0), m_mult(-x), m_prod(-1) {}
65
operator ()log1p_series66 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>
log1p(std::complex<T> x)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
main()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 }
108