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)类的实现:
#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++模板类技术实现:
#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中实现:
#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 #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)
下面就是实现这个算法了,不算太难:
#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再大时候基本上意义也不大,于是就没有再费劲优化了。
欢迎讨论。