UVA 766 Sum of powers



  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.

Input 

The input file contains several datasets, each of them containing a single integer  k  (  $0 \le k \le 20$ ).

The first line of the input contains the number of datasets, and it's followed by a blank line. There's also a blank line between datasets.

Output 

For each dataset, write integer numbers  $M, a_{k+1}, a_k, \dots, a_1, a_0$  to the output file in the given order. Numbers should be separated by one space. Remember that you should write the answer with the smallest positive  M possible.

Print a blank line between datasets.

Sample input 

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("");
    }
}





评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值