1. 程式人生 > 實用技巧 >BZOJ #3453. tyvj 1858 XLkxc

BZOJ #3453. tyvj 1858 XLkxc

題意:

\[f(i)=1^k+2^k+\cdots + i^k\\ g(i)=f(1)+f(2)+\cdots+ f(i) \]

給定 \(k,a,i,d\) ,求

\[h(i)=g(a)+g(a+d)+g(a+2d)+\cdots+ g(a+id) \]

引理:若一個多項式差分 \(k+1\) 次後為 \(0\) ,則這是一個 \(k\) 次多項式(若 \(g(x)=f(x+1)-f(x)\),則 \(g(x)\) 為差分一次後的多項式)。

根據之前的知識,我們知道 \(f(i)\) 是一個關於 \(i\)\(k+1\) 次多項式,而 \(g(i)\) 是一個差分一次後是 \(f(i+1)\)

(也是一個 \(k+1\) 次多項式),所以 \(g(i)\) 是一個 \(k+2\) 次多項式;同理,\(h(i)\) 是一個 \(k+3\) 次多項式。

\(g(i)\) 可以用拉格朗日插值在 \(O(k)\) 的時間求出來,那麼我們可以暴力求出 \(h(1)+h(2)+\cdots+h(k+4)\) 的值,再用拉格朗日插值求出要求的 \(h(i)\)。這個故事告訴我們插值可以巢狀。

程式碼:

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#define int long long
#include<assert.h>

using namespace std;

typedef long long LL;
const int K=200,M=1234567891;
int a[K+9],p[K],cnt;
int k,A,n,D;
LL d[K+9],pre[K+9],suf[K+9],fac[K+9],inv_fac[K+9],h[K+9];

int ksm(int a,int b)
{
	int res=1;
	while(b)
	{
		if(b&1)
			res=1ll*res*a%M;
		b>>=1,a=1ll*a*a%M;
	}
	return res;
}

void prework()
{
	d[1]=1;
	for (int i=2;i<=K;i++)
		a[i]=0;
	cnt=0;
	for (int i=2;i<=K;i++)
	{
		if(!a[i])
			a[i]=i,p[++cnt]=i,d[i]=ksm(i,k);
		for (int j=1;j<=cnt;j++)
		{
			if(p[j]>a[i]||p[j]>K/i)
				break;
			a[p[j]*i]=p[j];
			d[p[j]*i]=d[p[j]]*d[i]%M;
		}
	}
	for (int i=1;i<=K;i++)
		d[i]=(d[i]+d[i-1])%M;
	for (int i=1;i<=K;i++)
		d[i]=(d[i]+d[i-1])%M;
	
}

void init()
{
	scanf("%lld %lld %lld %lld",&k,&A,&n,&D);
	prework();
}

LL Get_Sum(LL n,LL d[],int k)
{
	pre[0]=1;
	for (int i=1;i<=k+3;i++)
		pre[i]=pre[i-1]*((n-i)%M)%M;
	suf[k+4]=1;
	for (int i=k+3;i>=1;i--)
		suf[i]=suf[i+1]*((n-i)%M)%M;
	fac[0]=1;
	for (int i=1;i<=k+3;i++)
		fac[i]=fac[i-1]*i%M;
	inv_fac[k+3]=ksm(fac[k+3],M-2);
	for (int i=k+2;i>=0;i--)
		inv_fac[i]=inv_fac[i+1]*(i+1)%M;
	for (int i=1;i<=k+3;i++)
		assert(fac[i]*inv_fac[i]%M==1);
	LL ans=0;
	for (int i=1;i<=k+3;i++)
		ans=(ans+d[i]*pre[i-1]%M*suf[i+1]%M*inv_fac[i-1]%M*inv_fac[k+3-i]%M*(k+3-i&1?-1:1))%M;
	return ans;
}

void work()
{
	h[0]=Get_Sum(A,d,k);
	for (int i=1;i<=k+4;i++)
		h[i]=(h[i-1]+Get_Sum(A+1ll*i*D,d,k))%M;
	printf("%lld\n",(Get_Sum(n,h,k+1)+M)%M);
}

signed main()
{
	int T;
	scanf("%lld",&T);
	while(T--)
	{
		init();
		work();
	}
	return 0;
}