C++求自然数m次幂求和通式

问题

已知m要求fm(n)的表达式。

fm(n)=1m+2m+3m+…+nm(m>=0)

例如:

f0(n) = n

f1(n) = (1/2)n^2 + (1/2)n

数据结构

要求的结果为多项式,而多项式的系数是必须是分数,于是需要抽象出两个数据类型:分数(Fraction)和多项式(Polynomial)。
我将这两个类都定义到了我自己的命名空间(MathUilt)中,方便以后生成库文件重用。

分数(Fraction)类的实现:

Fraction.h

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
#ifndef FRACTION_H
#define FRACTION_H
#include <cmath>
namespace MathUilt
{
class Fraction
{
// this is a fraction n/d
int n; // numerator
unsigned int d; // denominator

public:
Fraction(int nn,int dd)
:n(std::abs(nn)), d(std::abs(dd))
{
unsigned min = (n < d) ? n : d; //the max common divisor of n and d
if(min==0)
{
n = 0;
d = 1;
return;
}
while(n % min != 0 || d % min != 0) min--;
n /= min;
d /= min;
if(nn*dd < 0) n = -n;
}

Fraction(int nn) : n(nn), d(1) {}

~Fraction(){}

friend const Fraction operator*(const Fraction& left,const Fraction& right)
{
return Fraction(left.n * right.n, left.d * right.d);
}

friend const Fraction operator*(const Fraction& left,const int& right)
{
return Fraction(left.n * right, left.d);
}

friend const Fraction operator*(const int& left,const Fraction& right)
{
return Fraction(left * right.n, right.d);
}

const Fraction& operator*=(const Fraction& right)
{
return *this = *this * right;
}

const Fraction& operator*=(const int& right)
{
return *this = *this * right;
}

friend const Fraction operator+(const Fraction& left,const Fraction& right)
{
return Fraction(left.n * right.d + right.n * left.d, left.d * right.d);
}

friend const Fraction operator+(const Fraction& left,const int& right)
{
return Fraction(left.n + right * left.d, left.d);
}

friend const Fraction operator+(const int& left,const Fraction& right)
{
return Fraction(left * right.d + right.n, right.d);
}

const Fraction& operator+=(const Fraction& right)
{
return *this = *this + right;
}

const Fraction& operator+=(const int& right)
{
return *this = *this + right;
}

friend std::ostream& operator<<(std::ostream& os, const Fraction& right)
{
if(right.d == 1)
os << right.n;
else
os << '(' << right.n << '/' << right.d << ')';
return os;
}

friend const bool operator==(const Fraction& left,const int& right)
{
if(left.n == right * left.d) return true;
return false;
}
};
}
#endif //FRACTION_H ///:~

为了代码的可重用性,多项式类用C++模板类技术实现:

多项式类(Polynomial)类的实现:

Polynomial.h

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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H
#include <vector>
#include <iostream>
namespace MathUilt
{
template<class T> class Polynomial
{
typedef unsigned int UINT;

struct Item // data type of polynomial items
{
UINT exp; // the exponent of the item
T value; // the value fo the item
Item(UINT e, T v) : exp(e), value(v) {}
};

std::vector<Item> data;

typedef typename std::vector<Item>::const_iterator item_const_itera;
typedef typename std::vector<Item>::iterator item_itera;

public:
Polynomial();
~Polynomial();

void addItem(UINT exp,T value);

bool empty() const
{
return data.empty();
}

friend const Polynomial operator*(const Polynomial& left, const Polynomial& right)
{
Polynomial<T> result;
item_const_itera lit,rit;
for(lit=left.data.begin();lit!=left.data.end();lit++)
{
for(rit=right.data.begin();rit!=right.data.end();rit++)
{
result.addItem(lit->exp + rit->exp, lit->value * rit->value);
}
}
return result;
}

friend const Polynomial operator*(const Polynomial& left, const T& right)
{
item_itera it;
Polynomial result(left);
for(it=result.data.begin();it!=result.data.end();it++)
{
it->value *= right;
}
return result;
}

friend const Polynomial operator*(const T& left,const Polynomial& right)
{
item_itera it;
Polynomial result(right);
for(it=result.data.begin();it!=result.data.end();it++)
{
it->value *= left;
}
return result;
}

Polynomial& operator*=(const Polynomial& right)
{
return *this = *this * right;
}

Polynomial& operator*=(const T& right)
{
return *this = *this * right;
}

friend const Polynomial operator+(const Polynomial& left, const Polynomial& right)
{
item_const_itera it;
Polynomial result(left);
for(it=right.data.begin();it!=right.data.end();it++)
{
result.addItem(it->exp, it->value);
}
return result;
}

friend const Polynomial operator+(const Polynomial& left, const T& right)
{
Polynomial result(left);
result.addItem(0,right);
return result;
}

friend const Polynomial operator+(const T& left, const Polynomial& right)
{
Polynomial result(right);
result.addItem(0,left);
return result;
}

Polynomial& operator+=(const Polynomial& right)
{
return *this = *this + right;
}

Polynomial& operator+=(const T& right)
{
return *this = *this + right;
}

friend std::ostream& operator<<(std::ostream& os, const Polynomial& right)
{
if(right.empty()) return os << "NULL";
item_const_itera it;
for(it=right.data.begin();it!=right.data.end();it++)
{
if(it->exp == 0)
os << it->value;
else if(it->exp == 1)
{
if(it->value == 1)
os << 'n';
else
os << it->value << 'n';
}
else
{
if(it->value == 1)
os << "n^" << it->exp;
else
os << it->value << "n^" << it->exp;
}
if(it != right.data.end()-1) os << " + ";
}
return os;
}
};

template<class T> Polynomial<T>::Polynomial()
{

}

template<class T> Polynomial<T>::~Polynomial()
{
data.clear();
}

template<class T> void Polynomial<T>::addItem(UINT exp, T value)
{
if(value == 0) return;
item_itera it;
bool isNew = true;
for(it=data.begin();it!=data.end();it++)
{
if(it->exp==exp)
{
it->value += value;
if(it->value==0) data.erase(it);
return;
}
}
if(it==data.end())
{
for(it=data.begin();;it++)
{
if(it==data.end() || it->exp < exp)
{
data.insert(it,Item(exp,value));
return;
}
}
}
}
}
#endif //POLYNOMIAL_H ///:~

这两个头文件Fraction.h 和 Polynomial.h都包含在MathUilt.h中,该头文件也定义了一些数学工具函数,在MathUilt.cpp中实现:

MathUilt.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#ifndef MATHUILT_H
#define MATHUILT_H
#include <ostream>
#include "Fraction.h"
#include "Polynomial.h"
namespace MathUilt
{
class Fraction;

template<class T> class Polynomial;

typedef unsigned int UINT;

UINT factorial(UINT n, UINT k = 1);
// requires n >= 0
// return the factorial from n to 1 if k is not given
// else return the factorial from n to k

UINT cmb(UINT n, UINT m);
// requires n <= m
// return the combination of n in m
}

#endif //MATHUILT_H ///:~

MathUilt.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// MathUilt.cpp
#include "MathUilt.h"
using namespace MathUilt;

UINT MathUilt::factorial(UINT n, UINT k /*= 1*/)
{
if(n==0 || n==1) return 1;
if(n==k) return k;
else return n*factorial(n-1,k);
}

UINT MathUilt::cmb(UINT n, UINT m)
{
// if n > m throw exception
return factorial(m,m-n+1) / factorial(n);
}
///:~

算法

实现算法的公式就是递推,具体的数学推导就不说了,就是把(n+1)m+1用二项式展开,最后得到下面这个式子:

fm(n) = {(n+1)m+1 - 1 - [C(2,m+1)fm-1(n) + C(3,m+1)fm-2(n) + … + C(m+1,m+1)f0(n)]} / C(1,m+1)

下面就是实现这个算法了,不算太难:

main.cpp

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
#include <iostream>
#include "MathUilt.h"
using namespace std;
using namespace MathUilt;
#define PRINT(X) cout << #X " = " << X << endl

Polynomial<Fraction>* summation_p(unsigned int m, Polynomial<Fraction>* result_p[]);

Polynomial<Fraction> summation(unsigned int m)
// this function return a polynomial which is
// the generial formula of the summation of the k-th
// power of natural numbers from 1 to n
{
Polynomial<Fraction> **result_p,result;
result_p = new Polynomial<Fraction>*[m+1];
for(int i=0;i<m+1;i++) result_p[i] = NULL;
summation_p(m,result_p);
result = *result_p[m];
for(int i=0;i<m+1;i++) delete result_p[i];
delete[] result_p;
return result;
}

int main()
{
int m;
cin >> m;
cout << summation(m) << endl;
return 0;
}

Polynomial<Fraction>* summation_p(unsigned int m, Polynomial<Fraction>* result_p[])
{
Polynomial<Fraction> *result = new Polynomial<Fraction>();
if(m==0) result->addItem(1,Fraction(1));
else
{
Polynomial<Fraction> temp;
temp.addItem(0,Fraction(1));
temp.addItem(1,Fraction(1));
*result += temp;
for(int i=0;i<m;i++)
{
*result *= temp;
}

result->addItem(0,Fraction(-1));

for(unsigned int i=2;i<=m+1;i++)
{
if(result_p[m+1-i]==NULL) summation_p(m+1-i,result_p);
*result += *result_p[m+1-i] * Fraction(-cmb(i,m+1));
}

*result *= Fraction(1,cmb(1,m+1));
}
result_p[m] = result;
return result;
}

用了指针做返回值,相比直接用堆栈中的对象来做递归的返回值(需要调用大量的拷贝构造和析构)提高了调度的效率;而用指针数组保存堆中已经计算完成的fk(n)的对象也可以大幅度的减少计算时间。

虽然做了这些优化,但是计算到m=19的时候已经耗时比较长了(优化前只有一个函数返回堆栈中的对象来递归,m=15的时候计算时间就已经相当长了)。根据计算的结果看出m=19时分数的分母部分已经是相当大的整数了,m再大时候基本上意义也不大,于是就没有再费劲优化了。

欢迎讨论。