Blame view

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