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

问题：已知m要求f_m(n)的表达式。

f_m(n) = 1^m+2^m+3^m+...+n^m(m>=0)

例如：

f₀(n) = n

f₁(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用二项式展开，最后得到下面这个式子：

f_m(n) = {(n+1)^m+1- 1 - [C(2,m+1)f_m-1(n) + C(3,m+1)f_m-2(n) + ... + C(m+1,m+1)f₀(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;
}

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

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

欢迎讨论。
原文链接: https://www.cnblogs.com/zorro-x/archive/2013/04/02/2994440.html

欢迎关注

微信关注下方公众号，第一时间获取干货硬货；公众号内回复【pdf】免费获取数百本计算机经典书籍