Gym - 101480F[FFT或者推公式]

7 篇文章 0 订阅

题目链接:https://vjudge.net/contest/298078#problem/F

 

解题思路:

很显然是去考虑每个点的贡献:

对于每个(1,i)上的值ti,它的贡献是ti*C(2*n-2-i,n-i)*b^{n-1}a^{n-i},为什么是2*n-2-i呢,因为第一步肯定是往上走的。

对于每个(i,1)上的值li,它的贡献是li*C(2*n-2-i,n-i)*a^{n-1}b^{n-i},还是因为要先向右走。

对于这两种贡献可以直接用O(n)求出,剩下就是考虑每个点(i,j)的c对答案的贡献。

对于c的贡献次数是:

                                 \sum_{i=2}^{n}\sum_{j=2}^{n}C(2*n-i-j,n-i)a^{n-i}b^{n-j}

可以再变换成:    

                              \sum_{i=0}^{n-2}\sum_{j=0}^{n-2}C(i+j,i)*a^{i}b^{j}

把组合数变为阶乘的形式,最后变为:

                            \sum_{i=0}^{2*n-4}i!\sum_{j=0}^{i}\frac{b^{i-j}}{(i-j)!}\frac{a^{j}}{j!}[i-j<=n-2 and j<n-2]

很明显的一个FFT卷积式子, 由于精度问题所以我们可以把系数x拆成(a*N+b)的形式,上下都拆做四次。

时间复杂度O(大常数*n*log(n))

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef complex<double> comp;
const double pi = acos(-1);
const int mx = 4e5 + 10;
const int mod = 1e6 + 3;
const int N = 1e3;
int n,rev[1<<20];
ll inv[mx],fac[mx],w1[1<<20],w2[1<<20];
ll ret[mx];
int a,b,c,row[mx],col[mx];
comp a1[1<<20],b1[1<<20];
void init()
{
	inv[0] = fac[0] = 1;
	inv[1] = fac[1] = 1;
	for(int i=2;i<mx;i++){
		fac[i] = fac[i-1]*i%mod;
		inv[i] = (mod-mod/i)*inv[mod%i]%mod;
	}
	for(int i=2;i<mx;i++) inv[i] = inv[i]*inv[i-1]%mod;
}
ll qpow(ll x,ll y){
	ll ans = 1;
	while(y){
		if(y&1) ans = ans*x%mod;
		y >>= 1;
		x = x*x%mod;
	}
	return ans;
}
ll C(int N,int M){
	return fac[N]*inv[M]%mod*inv[N-M]%mod;
}
void get_rev(int len)
{
	for(int i=1;i<(1<<len);i++)
	rev[i] = (rev[i>>1]>>1)|((i&1)<<(len-1));
}
void fft(comp *p,int len,int v)
{
	for(int i=0;i<len;i++)
	if(i<rev[i]) swap(p[i],p[rev[i]]);
	for(int i=1;i<len;i<<=1){
		comp tep = exp(comp(0,v*pi/i));
		for(int j=0;j<len;j+=(i<<1)){
			comp mul(1,0);
			for(int k=j;k<i+j;k++){
				comp x = p[k];
				comp y = mul*p[k+i];
				p[k] = x + y,p[k+i] = x - y;
				mul *= tep;
			}
		}
	}
	if(v==-1) for(int i=0;i<len;i++) p[i] /= len;
}
void mul(int len,ll d){
	fft(a1,len,1);
	fft(b1,len,1);
	for(int i=0;i<len;i++) a1[i] *= b1[i];
	fft(a1,len,-1);
	for(int i=0;i<2*n;i++)
	ret[i] = (ret[i]+ll(a1[i].real()+0.5)*d)%mod;
}
int main()
{ 	
	init();
	scanf("%d%d%d%d",&n,&a,&b,&c);
	ll ans = 0;
	int bit = 0,ma = 2*n,len;
	for(int i=1;i<=n;i++) scanf("%d",col+i);
	for(int i=1;i<=n;i++) scanf("%d",row+i);
	for(int i=2;i<=n;i++){
		ans += row[i]*C(2*n-i-2,n-i)%mod*qpow(b,n-1)%mod*qpow(a,n-i)%mod;
		ans %= mod;
		ans += col[i]*C(2*n-i-2,n-i)%mod*qpow(a,n-1)%mod*qpow(b,n-i)%mod;
		ans %= mod;
	}
	for(int i=0;i<=n-2;i++){
		w1[i] = qpow(b,i)*qpow(fac[i],mod-2)%mod;
		w2[i] = qpow(a,i)*qpow(fac[i],mod-2)%mod;
	}
	while((1<<bit)<ma) bit++;
	len = (1<<bit);get_rev(bit);
	for(int i=0;i<len;i++){
		int c1 = w1[i]/N,c2 = w2[i]/N;
		a1[i] = comp(c1,0);
		b1[i] = comp(c2,0);
	}
	mul(len,N*N);
	for(int i=0;i<len;i++){
		int c1 = w1[i]/N,c2 = w2[i]%N;
		a1[i] = comp(c1,0);
		b1[i] = comp(c2,0);
	}
	mul(len,N);
	for(int i=0;i<len;i++){
		int c1 = w1[i]%N,c2 = w2[i]/N;
		a1[i] = comp(c1,0);
		b1[i] = comp(c2,0);
	}
	mul(len,N);
	for(int i=0;i<len;i++){
		int c1 = w1[i]%N,c2 = w2[i]%N;
		a1[i] = comp(c1,0);
		b1[i] = comp(c2,0);
	}
	mul(len,1);
	for(int i=0;i<=2*n-4;i++)
	ans = (ans+fac[i]*ret[i]%mod*c%mod)%mod;
	printf("%lld\n",ans);
	return 0;
}

直接推公式的常数小的O(n*log(n))的链接点击这里,这里我就不推了= =(太懒了),不过要注意a或者b等于1时式子不适用。

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int mx = 4e5 + 10;
const int mod = 1e6 + 3;
int n,a,b,c,col[mx],row[mx];
ll inv[mx],fac[mx],ga[mx],f[mx],gb[mx];
void init()
{
	inv[0] = fac[0] = 1;
	inv[1] = fac[1] = 1;
	for(int i=2;i<mx;i++){
		fac[i] = fac[i-1]*i%mod;
		inv[i] = (mod-mod/i)*inv[mod%i]%mod;
	}
	for(int i=2;i<mx;i++) inv[i] = inv[i]*inv[i-1]%mod;
}
ll C(int N,int M){
	return fac[N]*inv[M]%mod*inv[N-M]%mod;
}
ll qpow(ll x,ll y){
	ll ans = 1;
	while(y){
		if(y&1) ans = ans*x%mod;
		y >>= 1;
		x = x*x%mod;
	}
	return ans;
}
int main()
{
	init();
	scanf("%d%d%d%d",&n,&a,&b,&c);
	ll ans = 0;
	for(int i=1;i<=n;i++) scanf("%d",col+i);
	for(int i=1;i<=n;i++) scanf("%d",row+i);
	for(int i=2;i<=n;i++){
		ans += row[i]*C(2*n-i-2,n-i)%mod*qpow(b,n-1)%mod*qpow(a,n-i)%mod;
		ans %= mod;
		ans += col[i]*C(2*n-i-2,n-i)%mod*qpow(a,n-1)%mod*qpow(b,n-i)%mod;
		ans %= mod;
	}
	f[0] = 1,ga[1] = gb[1] = 1;
	for(int i=2;i<=n;i++){
		if(a==1) ga[i] = C(2*i+1,i)-C(2*i,i);
		else{
			ga[i] = (ga[i-1]+C(2*i-2,i-1)*qpow(a,i-1)%mod-C(2*i-1,i-1)*qpow(a,i)%mod);
			ga[i] = ga[i]*qpow(1-a+mod,mod-2)%mod;
		}
		if(b==1) gb[i] = C(2*i+1,i)-C(2*i,i);
		else{
			gb[i] = (gb[i-1]+C(2*i-2,i-1)*qpow(b,i-1)%mod-C(2*i-1,i-1)*qpow(b,i)%mod);
			gb[i] = gb[i]*qpow(1-b+mod,mod-2)%mod;
		}
	}  
	for(int i=1;i<=n;i++){
		f[i] = f[i-1] + qpow(a,i)*gb[i]%mod + qpow(b,i)*ga[i]%mod;
		f[i] = (f[i]+C(2*i,i)*qpow(a,i)%mod*qpow(b,i)%mod)%mod; 
	}
	printf("%lld\n",(ans+mod+f[n-2]*c%mod)%mod);
	return 0;
}

 

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值