bzoj 3453 数论

　　首先我们知道对于f(x)来说，它是一个k次的多项式，那么f(x)的通项公式可以表示成一个k+1次的式子，且因为f(x)没有常数项，所以我们设这个式子为

　　　　f(x)=Σ(a[i]*x^i) (1<=i<=k+1)

那么比较显然的是f(x+1)-f(x)=(x+1)^k，因为(x+1)^k=Σc(k,i)*x^i (0<=i<=k)，所以我们可以将这个式子的左右展开，可以得到

　　　　f(x+1)-f(x)=(x+1)^k    Σ(a[i]*(x+1)^i)-Σ(a[i]*x^i)=(x+1)^k   Σa[i](Σc(i,j)*x^j (0<=j<=i)-x^i) (1<=i<=k+1)=Σc(k,i)*x^i (0<=i<=k)

　　　　Σa[i]Σc(i,j)*x^j (0<=j<=i-1) (1<=i<=k+1)=Σc(k,i)*x^i (0<=i<=k)

那么我们发现，式子的左右两边x的指数都是属于[0,k]的，那么对应项系数相等，可以得出k+1个式子，则对于任意一个系数j，我们有Σa[i]*c(i,j) (j+1<=i<=k+1)=c(k,j)。这样我们就有了k+1个式子，对于a数组的k+1个未知数，我们可以用高斯消元来解决这个问题。

　　下面对于g(x)，我们可以表示为

　　　　g(x)=Σ(x-i+1)*i^k (1<=i<=x)

　　我们可以进一步整理为

　　　　g(x)=(x+1)Σi^k-Σi^(k+1) (1<=i<=x)

　　那么我们发现，这个式子的前面就是(x+1)乘上f(x)，后面就是一个类似f(x)的东西，只是系数变成了k+1，那么我们设这个表达式为w(x)且系数为b[i]，那么类似于求f(x)的过程我们可以将b[i]求出来，那么g(x)=(x+1)*f(x)-w(x)，那么对于相同指数的x我们可以直接合并，那么我们就可以得到g(x)的系数组。

　　那么我们构造矩阵A[i](1,i^1,i^2....,i^(k+2))，那么对于(i+1)^j=Σc(j,l)*i^l (0<=l<=j)，那么我们显然可以由上一层累加得到，系数为c(j,l)。

　　那么我们构造矩阵B[j](1,i^1,i^2....,i^(k+2),Σg(s+j*d)) i=s+j*d，初始的时候最后一项可以由A[i]矩阵的s次幂转移乘上g的系数组得到，那么对于这个矩阵的转移矩阵，前半部分显然可以由为A[i]的转移矩阵的d次自乘得到，对于最后一项的转移则为上一次的Σ值累加上g(s+i*d)的值，我们可以由前面的i的若干次幂乘上g(x)的系数组转移得到。

　　

/**************************************************************
    Problem: 3453
    User: BLADEVIL
    Language: C++
    Result: Accepted
    Time:11180 ms
    Memory:1444 kb
****************************************************************/
 
//By BLADEVIL
#include <cstdio>
#include <cstring>
#include <algorithm>
#define LL long long
#define maxk 150
#define P 1234567891
 
using namespace std;
 
struct mat{
    int n,m;
    int a[maxk][maxk];
    mat (int n,int m):n(n),m(m){memset(a,0,sizeof a);};
};
 
int K,S,N,D;
int C[maxk][maxk];
int a[maxk][maxk],fac[3][maxk];
 
 
mat operator *(const mat &a,const mat &b) {
    mat c(a.n,b.m);
    for (int i=0;i<a.n;i++)
        for (int j=0;j<b.m;j++) 
            for (int k=0;k<a.m;k++) 
                c.a[i][j]=(c.a[i][j]+(LL)a.a[i][k]*b.a[k][j])%P;
    return c;
}
 
mat pwr(mat a,int x) {
    mat ans(a.n,a.n);
    for (int i=0;i<a.n;i++) ans.a[i][i]=1;
    while (x) {
        if (x&1) ans=ans*a;
        a=a*a;
        x>>=1;
    }
    return ans;
}
 
void pre(){
    C[0][0]=1;
    for (int i=1;i<maxk;i++) {
        C[i][0]=C[i][i]=1;
        for (int j=1;j<i;j++) 
            C[i][j]=((LL)C[i-1][j-1]+C[i-1][j])%P;
    }
}
 
int pwr(int x,int y) {
    int ans=1;
    while (y) {
        if (y&1) ans=((LL)ans*x)%P;
        x=((LL)x*x)%P;
        y>>=1;
    }
    return ans;
}
 
void gauss(int a[][maxk],int n) {
    for (int i=0;i<n;i++) {
        int k=0,t;
        for (k=i;(k<n)&&(!a[k][i]);k++);
        for (int j=0;j<=n;j++) swap(a[i][j],a[k][j]);
        t=pwr(a[i][i],P-2);
        for (int j=i;j<=n;j++) a[i][j]=((LL)a[i][j]*t)%P;
        for (int j=0;j<n;j++) if (j!=i)
            for (k=n;k>=i;k--) a[j][k]=(a[j][k]+(LL)a[j][i]*(P-a[i][k]))%P;
    }
}
 
void solve() {
    memset(a,0,sizeof a); memset(fac,0,sizeof fac);
    for (int i=0;i<=K;i++) {
        for (int j=i+1;j<=K+1;j++) a[i][j-1]=C[j][i];
        a[i][K+1]=C[K][i];
    }
    gauss(a,K+1);
    for (int i=0;i<=K;i++) fac[0][i+1]=a[i][K+1]; K++;
    memset(a,0,sizeof a);
    for (int i=0;i<=K;i++) {
        for (int j=i+1;j<=K+1;j++) a[i][j-1]=C[j][i];
        a[i][K+1]=C[K][i];
    }
    gauss(a,K+1);
    for (int i=0;i<=K;i++) fac[1][i+1]=a[i][K+1]; K--;
    for (int i=2;i<=K+2;i++) fac[2][i]=fac[0][i-1];
    for (int i=1;i<=K+1;i++) fac[2][i]=((LL)fac[2][i]+fac[0][i])%P;
    for (int i=1;i<=K+2;i++) fac[2][i]=((LL)fac[2][i]-fac[1][i]+P)%P;
    mat DP(K+3,K+3);
    for (int i=0;i<=K+2;i++) 
        for (int j=0;j<=i;j++) DP.a[j][i]=C[i][j];
    DP=pwr(DP,D); DP.n++; DP.m++;
    for (int i=0;i<=K+2;i++)
        for (int j=0;j<=K+2;j++)
            DP.a[j][K+3]=(DP.a[j][K+3]+(LL)DP.a[j][i]*fac[2][i])%P;
    DP.a[K+3][K+3]=1;   
    mat A(1,K+4); A.a[0][0]=1;
    for (int i=1;i<=K+2;i++) A.a[0][i]=((LL)A.a[0][i-1]*S)%P;
    for (int i=1;i<=K+2;i++) A.a[0][K+3]=(A.a[0][K+3]+(LL)fac[2][i]*A.a[0][i])%P;
    A=A*pwr(DP,N);
    printf("%d
",A.a[0][K+3]);
}
 
int main() {
    pre();
    int task;
    scanf("%d",&task);
    while (task--) {
        scanf("%d%d%d%d",&K,&S,&N,&D);
        solve();
    }
    return 0;   
}