51nod 1161 Partial Sums,1172 Partial Sums V2

这两题原理是一样的,不过第二题数据量大一些。这个累加操作相当于一个矩阵乘法,然后用矩阵中的第一列数和输入的数组做卷积,比如这个样例处理2次矩阵就是这样的


然后取出第一列数和输入的数做卷积,也就是多项式乘法

$\left(1+2x+3x^2+4x^3\right) \left(1+3x+5x^2+6x^3\right)=1 + 5 x + 14 x^2 + 29 x^3 + 39 x^4 + 38 x^5 + 24 x^6$

前四项的系数1,5,14,29就是结果了

这第一列数我们可以不用矩阵乘法做,经过观察发现它们是一些有关联的组合数,而且可以递推O(n)求出第一列数

$C_0=1,C_i=\frac{C_{i-1} (i+k-1)}{i}$

 除法用逆元处理,第一题数据量小用O(n^2)模拟乘法就行了,Python写的

n,k=map(int,raw_input().split())
c=[1]*n
p=10**9+7
a=[int(input()) for i in xrange(n)]
for i in xrange(1,n):
    c[i]=c[i-1]*(k-1+i)*pow(i,p-2,p)%p
c.reverse()
for i in xrange(1,n+1):
    s=0
    for j in xrange(i):
        s=s+a[j]*c[-(i-j)]
    print(s%p)
View Code

第二题50000数据就需要FFT了,还需要模1e9+7,我们用数论变换+中国剩余定理做了

#include <iostream>
using namespace std;

#define N 340030
#define M 340020
#define LL long long
#define mod 1000000007
#define K 3

const int m[K] = {1004535809, 998244353, 104857601};
const int G = 3;

LL qpow(LL x, LL k, LL p) {
    int ret = 1;
    while(k) {
        if(k & 1) ret = 1LL * ret%p * x % p;
        k >>= 1;
        x = 1LL * x%p * x % p;
    }
    return ret;
}
struct _NTT {
    LL wn[25], p;

    void init(LL _p) {
        p = _p;
        for(int i = 1; i <= 21; ++i) {
            int t = 1 << i;
            wn[i] = qpow(G, (p - 1) / t, p);
        }
    }
    void change(LL *y, int len) {
        for(int i = 1, j = len / 2; i < len - 1; ++i) {
            if(i < j) swap(y[i], y[j]);
            int k = len / 2;
            while(j >= k) j -= k, k /= 2;
            j += k;
        }
    }
    void NTT(LL y[], int len, int on) {
        change(y, len);
        int id = 0;
        for(int h = 2; h <= len; h <<= 1) {
            ++id;
            for(int j = 0; j < len; j += h) {
                LL w = 1;
                for(int k = j; k < j + h / 2; ++k) {
                    int u = y[k];
                    int t = y[k+h/2] * w % p;
                    y[k] = u + t;
                    if(y[k] >= p) y[k] -= p;
                    y[k+h/2] = u + p - t;
                    if(y[k+h/2] >= p) y[k+h/2] -= p;
                    w = w * wn[id] % p;
                }
            }
        }
        if(on == -1) {
            for(int i = 1; i < len / 2; ++i) swap(y[i], y[len-i]);
            int inv = qpow(len, p - 2, p);
            for(int i = 0; i < len; ++i)
                y[i] = 1LL * y[i] * inv % p;
        }
    }
    void mul(LL A[], LL B[], LL len) {
        NTT(A, len, 1);
        NTT(B, len, 1);
        for(int i = 0; i < len; ++i) A[i] = 1LL * A[i] * B[i] % p;
        NTT(A, len, -1);
    }
}ntt[K];

LL tmp[N][K], t1[N], t2[N];
LL x1[N], x2[N];
LL r[K][K];

LL CRT(LL a[]) {
    LL x[K];
    for(int i = 0; i < K; ++i) {
        x[i] = a[i];
        for(int j = 0; j < i; ++j) {
            int t = (x[i] - x[j]) % m[i];
            if(t < 0) t += m[i];
            x[i] = 1LL * t * r[j][i] % m[i];
        }
    }
    LL mul = 1, ret = x[0] % mod;
    for(int i = 1; i < K; ++i) {
        mul = 1LL * mul * m[i-1] % mod;
        ret += 1LL * x[i] * mul % mod;
        if(ret >= mod) ret -= mod;
    }
    return ret;
}

void mul(LL A[], LL B[], LL len) {
    for(int id = 0; id < K; ++id) {
        for(int i = 0; i < len; ++i) {
            t1[i] = A[i], t2[i] = B[i];
        }
        ntt[id].mul(t1, t2, len);
        for(int i = 0; i < len; ++i) {
            tmp[i][id] = t1[i];
        }
    }
    for(int i = 0; i < len; ++i) A[i] = CRT(tmp[i]);
}

LL a[N],b[N];
LL inv[50005];

void init() {
    for(int i = 0; i < K; ++i) {
        for(int j = 0; j < i; ++j)
            r[j][i] = qpow(m[j], m[i] - 2, m[i]);
    }

    for(int i = 0; i < K; ++i) ntt[i].init(m[i]);
    inv[1]=1;
    for(int x=2;x<50003;x++)inv[x]=(mod-mod/x)*inv[mod%x]%mod;
}


int main()
{
    init();
    int n,k,i;
    cin>>n>>k;
    b[0]=1;
    for(i=1;i<n;i++)
    {
        b[i]=b[i-1]*(k-1+i)%mod*inv[i]%mod;
        b[i]=(b[i]+mod)%mod;
    }
    for(i=0;i<n;i++)
    {
        cin>>a[i];
    }
    mul(a,b,1<<17);
    for(i=0;i<n;i++)
    {
        cout<<(a[i]+mod)%mod<<endl;
    }
    return 0;
}
View Code

 

posted @ 2017-01-22 14:05  czp001  阅读(198)  评论(0编辑  收藏  举报