UVA 766 Sum of powers

最新推荐文章于 2020-09-09 14:47:34 发布

Killerfear

最新推荐文章于 2020-09-09 14:47:34 发布

阅读量596

点赞数

分类专栏：数学：计数

本文链接：https://blog.csdn.net/wsx1754175/article/details/42226449

版权

数学：计数专栏收录该内容

18 篇文章

订阅专栏

介绍了一种快速计算固定自然数k下从1到n的幂和的方法，避免了传统逐项累加方式随n增大而增加计算量的问题。通过将幂和表达为关于n的k+1次多项式，实现了恒定数量的运算步骤，提供了具体实现代码。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Sum of powers

A young schoolboy would like to calculate the sum

$\begin{displaymath}S_k(n) = \sum_{i=1}^n {i^k}\end{displaymath}$

for some fixed natural k and different natural n . He observed that calculating i ^k for all i ( $1 \le i \le n$ ) and summing up results is a too slow way to do it, because the number of required arithmetical operations increases as n increases. Fortunately, there is another method which takes only a constant number of operations regardless of n . It is possible to show that the sum S _k ( n ) is equal to some polynomial of degree k +1 in the variable n with rational coefficients, i.e.,

$\begin{displaymath}S_k(n) = {1 \over M} \left( a_{k+1} n^{k+1} + a_k n^k + \dots + a_1 n+ a_0 \right)\end{displaymath}$

for some integer numbers $M, a_{k+1}, a_k, \dots, a_1, a_0$ .

We require that integer M be positive and as small as possible. Under this condition the entire set of such numbers (i.e. $M, a_{k+1}, a_k, \dots, a_1, a_0$ ) will be unique for the given k. You have to write a program to find such set of coefficients to help the schoolboy make his calculations quicker.

1

2

Sample Output

6 2 3 1 0

Miguel Revilla
2001-01-22

思路：

我们把S(n,k)当成一个多项式，然后重载各种+，-，*，/ 就行了。至于M只需要取多项式的每一项的分母的最小公倍数就行了。

代码：

#include <iostream>
#include <vector>
#include <algorithm>
#include <string.h>
#include <cstring>
#include <time.h>
#include <map>
#include <set>
#include <stdio.h>
#include <cmath>
#include <string>
#include <cassert>
#include <math.h>
#define rep(i,a,b) for(int i=(a);i<(b);++i)
#define rrep(i,b,a) for(int i = (b); i >= (a); --i)
#define clr(a,x) memset(a,(x),sizeof(a))
#define LL long long
#define eps 1e-9
#define mp make_pair
using namespace std;
const int maxn = 20+5;

LL gcd(LL a,LL b)
{
    while (a && b) {
        if (a > b) a %= b;
        else b %= a;
    }
    return a + b;
}

LL lcm(LL a,LL b) { return a * b / gcd(a,b); }
struct Fact
{
    LL a,b;
    Fact() { a = 0; b = 1; }
    Fact(LL _a,LL _b)
    {
        LL g = gcd(abs(_a),abs(_b));
        a = _a / g; b = _b / g;
        if (b < 0) a = -a, b = -b;
        if (a == 0) b = 1;
    }
    Fact operator + (const Fact & f) const
    {
        return Fact(a * f.b + f.a * b, b * f.b);
    }
    Fact operator - (const Fact & f) const
    {
        return Fact(a * f.b - f.a * b, b * f.b);
    }
    Fact operator * (const Fact & f) const
    {
        return Fact(a * f.a, b * f.b);
    }
    Fact operator * (LL x) const
    {
        return Fact(a * x, b);
    }
    Fact operator / (LL x) const
    {
        return Fact(a, b * x);
    }
    Fact operator / (const Fact & f) const
    {
        return Fact(a * f.b, b * f.a);
    }
    void print() const
    {
        printf("%lld",a);
        if (b > 1) printf("/%lld",b);
    }
};

struct Poly
{
    Fact a[maxn];
    Poly() { }
    Poly operator * (LL x) const
    {
        Poly ans = *this;
        rep(i,0,maxn) ans.a[i] = ans.a[i] * x;
        return ans;
    }
    Poly operator / (LL x) const
    {
        Poly ans = *this;
        rep(i,0,maxn) ans.a[i] = ans.a[i] / x;
        return ans;
    }

    Poly operator + (const Poly & p) const
    {
        Poly ans;
        rep(i,0,maxn) ans.a[i] = a[i] + p.a[i];
        return ans;
    }
    Poly operator - (const Poly & p) const
    {
        Poly ans;
        rep(i,0,maxn) ans.a[i] = a[i] - p.a[i];
        return ans;
    }

};

Poly ans[maxn];
LL M[maxn];
LL C[maxn][maxn];

void pre_init()
{
    rep(i,0,maxn) {
        C[i][0] = 1;
        rep(j,1,i+1) C[i][j] = C[i-1][j] + C[i-1][j-1];
    }
    M[0] = 1;
    ans[0].a[1] = Fact(1,1);
    rep(k,1,maxn) {
        Poly nk;
        rep(j,0,k+2) nk.a[j] = Fact(C[k+1][k+1-j],1);
        rep(j,2,k+2)
        nk = nk - (ans[k-j+1] * C[k+1][j]);
        nk.a[0] = nk.a[0] - Fact(1,1);
        nk = nk / (k+1);
        ans[k] = nk;
        M[k] = 1;
        rep(i,0,k+2) M[k] = lcm(M[k],nk.a[i].b);
    }
    rep(k,1,maxn) {
        ans[k] = ans[k] * M[k];
    }
}


int main()
{
    #ifdef ACM
        freopen("in.txt", "r", stdin);
       // freopen("out.txt","w",stdout);
    #endif // ACM
    pre_init();
    //cout << bell[900].sz << endl;
    int T; cin >> T;
    while (T--) {
        int n; scanf("%d",&n);
        printf("%lld ",M[n]);
        rrep(i,n+1,0) {
           // printf("%d",ans[n].a[i]);
            ans[n].a[i].print();
            if (i != 0) printf(" ");
        }
        puts("");
        if (T) puts("");
    }
}