bzoj 3453数論
15437 ワード
まず,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).これにより,a配列のk+1未知数に対してGauss消元を用いてこの問題を解決できる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)を乗じた係数群遷移から得ることができる.
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).これにより,a配列のk+1未知数に対してGauss消元を用いてこの問題を解決できる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;
}