loj2409. 「THUPC 2017」小 L 的计算题 / Sum

题意:

f k = ∑ i = 1 n a i k f_k=\sum_{i=1}^n a_i^k fk=i=1naik
求出 f 1 f_1 f1~ f n f_n fn

题解:

考虑构造 f f f的生成函数 F F F
A ( x ) = ∑ j a i j x j A(x)=\sum_j a_i^jx^j A(x)=jaijxj
F ( x ) = ∑ i n A ( i ) F(x)=\sum_i^n A(i) F(x)=inA(i)
F F F的系数即是 ∑ i = 1 n a i k \sum_{i=1}^n a_i^k i=1naik
所以有 F ( x ) = ∑ i = 1 n 1 1 − a i x F(x)=\sum_{i=1}^n \frac{1}{1-a_ix} F(x)=i=1n1aix1
这样就可以分治fft+多项式求逆 O ( n l o g 2 n ) O(nlog^2n) O(nlog2n)
然而常数巨大无法AC
接下来的操作就很奇奇妙妙
= ∑ i = 1 n 1 + a i x 1 − a i x =\sum_{i=1}^n 1+\frac{a_ix}{1-a_ix} =i=1n1+1aixaix

= n − x ∑ i = 1 n − a i 1 − a i x =n-x\sum_{i=1}^n \frac{-a_i}{1-a_ix} =nxi=1n1aixai

= n − x ∑ i = 1 n [ l n ( 1 − a i x ) ] ′ =n-x\sum_{i=1}^n [ln(1-a_ix)]' =nxi=1n[ln(1aix)]

= n − x [ ∑ i = 1 n l n ( 1 − a i x ) ] ′ =n-x \left[ \sum_{i=1}^n ln(1-a_ix) \right]' =nx[i=1nln(1aix)]

= n − x { l n [ ∏ i = 1 n ( 1 − a i x ) ] } ′ =n-x \left \{ ln \left[ \prod_{i=1}^n (1-a_ix) \right] \right \}' =nx{ln[i=1n(1aix)]}
用分治求这个式子就快很多,然后再求 l n ln ln,求导即可
听说是套路,自闭
code:

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#define mem(x) memset(x,0,sizeof(x))
#define LL long long
using namespace std;
const LL p=998244353,yg=3;
LL bin[1600010],a[1600010],g[1600010],h[1600010],tmp[1600010],b[1600010];
LL Pow(LL a,LL b)
{
	LL ans=1;
	while(b)
	{
		if(b&1) ans=ans*a%p;
		a=a*a%p;b>>=1;
	}
	return ans;
}
void ntt(LL *a,LL n,LL op)
{
	for(LL i=0;i<n;i++) bin[i]=(bin[i>>1]>>1)|((i&1)*(n>>1));
	for(LL i=0;i<n;i++) if(i<bin[i]) swap(a[i],a[bin[i]]);
	for(LL i=1;i<n;i<<=1)
	{
		LL wn=Pow(yg,op==1?(p-1)/(i<<1):p-1-(p-1)/(i<<1)),t,w;
		for(LL j=0;j<n;j+=i<<1)
		{
			w=1;
			for(LL k=0;k<i;k++)
			{
				t=a[i+j+k]*w%p;w=wn*w%p;
				a[i+j+k]=(a[j+k]-t)%p;a[j+k]=(a[j+k]+t)%p;
			}
		}
	}
	if(op==-1)
	{
		LL Inv=Pow(n,p-2);
		for(LL i=0;i<n;i++) a[i]=a[i]*Inv%p;
	}
}
void pol_inv(LL n,LL *a,LL *g)
{
	if(n==1) {g[0]=Pow(a[0],p-2);return;}
	pol_inv((n+1)>>1,a,g);
	LL N=1;while(N<n<<1) N<<=1;
	copy(a,a+n,tmp);fill(tmp+n,tmp+N,0);
	ntt(g,N,1);ntt(tmp,N,1);
	for(LL i=0;i<N;i++) g[i]=(2-tmp[i]*g[i]%p)*g[i]%p;
	ntt(g,N,-1);fill(g+n,g+N,0);
}
void pol_ln(LL n,LL *a,LL *h)
{
	LL N=1;while(N<n<<1) N<<=1;
	fill(g,g+N,0);pol_inv(n,a,g);
	for(LL i=0;i<n;i++) h[i]=a[i+1]*(i+1)%p;
	fill(h+n,h+N,0);ntt(g,N,1);ntt(h,N,1);
	for(LL i=0;i<N;i++) h[i]=h[i]*g[i]%p;
	ntt(h,N,-1);
	for(LL i=n-1;i>=1;i--) h[i]=h[i-1]*Pow(i,p-2)%p;
	fill(h+n,h+N,0);
}
LL n,c[400010];
void solve(LL l,LL r,LL N)
{
	if(l==r) {c[l]=-c[l];return;}
	LL mid=(l+r)/2;
	solve(l,mid,N>>1);solve(mid+1,r,N>>1);
	a[0]=b[0]=1;
	for(int i=1;i<N;i++) a[i]=b[i]=0;
	for(LL i=l;i<=mid;i++) a[i-l+1]=c[i];
	for(LL i=mid+1;i<=r;i++) b[i-mid]=c[i];
	ntt(a,N,1);ntt(b,N,1);
	for(int i=0;i<N;i++) a[i]=a[i]*b[i]%p;
	ntt(a,N,-1);
	for(LL i=l;i<=r;i++) c[i]=a[i-l+1];
}
LL read()
{
	char ch;
	while(!isdigit(ch=getchar()));
	LL x=ch-48;
	while(isdigit(ch=getchar())) x=x*10+ch-48;
	return x;
}
int main()
{
	LL T=read();
	while(T--)
	{
		n=read();
		for(LL i=1;i<=n;i++) c[i]=read();
		LL N=1;while(N<(n+1)<<1) N<<=1;
		solve(1,n,N);
		mem(a);a[0]=1;for(int i=1;i<=n;i++) a[i]=c[i];
		pol_ln(n+1,a,h);for(LL i=0;i<=n;i++) h[i]=h[i+1]*(i+1)%p;
		for(int i=n+1;i>=1;i--) h[i]=-h[i-1];h[0]=n;
		LL ans=0;for(LL i=1;i<=n;i++) ans^=(h[i]+p)%p;
		printf("%lld\n",ans);
	}
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值