poj 1707 Sum of powers 伯努利数系数

题目大意：

　　对于方程　　 $S_{n} = \sum_{n}^{i=1} i^{k} = \frac{1}{M}\cdot \left ( a_{k+1}\cdot n ^{k+1} \dotplus a_{k}\cdot n ^{k} \dotplus \cdot \cdot \cdot \dotplus a_{0}\cdot n ^{0} \right )$

　　输入 k ( 0 <= k <= 20 ) 输出 $\left ( M, a_{k+1},a_{k},a_{k-1}...a_{0} \right )$

解题思路：

　　由数学伯努利数公式：　　 $S_{k}^{n} = \frac{1}{k+1} \sum_{k+1}^{r} \left \binom{k+1}{r} \cdot B_{k+1-r} \cdot (n+1)^{r}$

　　带入通分化简就可以了，不过要注意，分母为正...( WA了一小时错在了这里 )

　伯努利数是18世纪瑞士数学家雅各布·伯努利引入的一个数。
　　设伯努利数为B(n)，它的定义为：
　　t/(e^t-1）=∑[B(n)*(t^n)/(n!)](n:0->；∞）
　　这里|t|<2。由计算知：
　　B(0)=1,B⑴=-1/2,
　　B⑵=1/6,B⑶=0,
　　B⑷=-1/30,B⑸=0,
　　B⑹=1/42,B⑺=0,
　　B⑻=-1/30,B⑼=0),
　　B⑽=5/66,B⑾=0,
　　B⑿=-691/2730,B⒀=0,
　　B⒁=7/6,B⒂=0,
　　B⒃=-3617/510,B⒄=0,
　　B⒅=43867/798,B⒆=0,
　　B⒇=-174611/330 ……

解题代码:

View Code

#include<stdio.h>
#include<string.h>
#include<stdlib.h>
typedef long long LL;
const int N = 25;

int B1[N],B2[N];
LL a[N], b[N], C[N][N];

void init(){//预处理伯努利系数
    B1[0] = 1; B2[0] = 1;
    B1[1] = -1;B2[1] = 2;
    B1[2] = 1; B2[2] = 6;
    B1[3] = 0; B2[3] = 0;
    B1[4] = -1; B2[4] = 30;
    B1[5] = 0; B2[5] = 0;
    B1[6] = 1; B2[6] = 42;
    B1[7] = 0; B2[7] = 0;
    B1[8] = -1; B2[8] = 30;
    B1[9] = 0; B2[9] = 0;
    B1[10] = 5,B2[10] = 66;
    B1[11] = 0; B2[11] = 0;
    B1[12] = -691; B2[12] = 2730;
    B1[13] = 0; B2[13] = 0;
    B1[14] = 7; B2[14] = 6;
    B1[15] = 0; B2[15] = 0;
    B1[16] = -3617; B2[16] = 510;
    B1[17] = 0; B2[17] = 0;
    B1[18] = 43867; B2[18] = 798;
    B1[19] = 0; B2[19] = 0;
    B1[20] = -174611; B2[20] = 330;

    // 组合数
    for(int i = 0; i <= 21; i++)
        C[i][0] = C[i][i] = 1;
    for(int i = 1; i <= 21; i++)
    {
        for(int j = 1; j < i; j++)
            C[i][j] = C[i-1][j]+C[i-1][j-1];
    }
}
LL gcd( LL a, LL b ){ return b==0?a:gcd(b,a%b);}
LL lcm( LL a, LL b ){ return (a/gcd(a,b))*b; }

int main(){
    init();
    int k;
    while( scanf("%d", &k) != EOF)
    {
        memset( a, 0, sizeof(a));
        memset( b, 0, sizeof(b));
        
        LL t1[N], t2[N];
        for(int r = 1; r <= k+1; r++)
        {
            memset( t1, 0, sizeof(t1));
            memset( t2, 0, sizeof(t2));    
            for(int i = 0; i <= r; i++)
            {
                t1[i] = B1[k+1-r]*C[k+1][r]*C[r][i];
                t2[i] = B2[k+1-r]; 
            }    
            for(int i = 0; i <= r; i++)
            {
                if( t2[i] ){ 
                    if( b[i] == 0 ){
                        a[i] = t1[i]; b[i] = t2[i];    
                    }    
                    else{//通分
                        LL z = lcm( b[i], t2[i] );
                        LL x = z/b[i], y = z/t2[i];
                        a[i] = a[i]*x + t1[i]*y;
                        b[i] = z;
                    }
                }    
            }
        }
        //获取所有分母最小公倍数 z    
        LL z = 1, p = 1;
        for(int i = 0; i <= k+1; i++)
            if( b[i] )    z = lcm( z, b[i] );
        for(int i = 0; i <= k+1; i++)
            if( b[i] ){ //统一分母z
                a[i] *= (z/b[i]);
                p = a[i]; b[i] = z;    
            }        
        z *= (k+1);    
        //获取分子分母最大公约数 p    
        for(int i = 0; i <= k+1; i++)
            if( b[i] ) p = gcd( p, a[i] );    
        p = gcd( p, z );        
        //简化形式    
        for(int i = 0; i <= k+1; i++)
            if( b[i] ) a[i] /= p;
        
        z /= p;    
        int f = ( z < 0 ) ? -1 : 1;

    //    printf("p = %lld\n", p );
        printf("%lld", f*z);
        for(int i = k+1; i >= 0; i-- )
            printf(" %lld",f*a[i]);
        puts("");
    }
    return 0;
}