[51nod1690][FFT]区间求和2

Description

给出一个长度为n的数组a。区间[L,R]的值为 ∑i=0R−La[L+i]∗a[R−i]

求所有长度为质数的区间的值的总和。

Input

第一行一个数n(1<=n<=100000)
第二行n个数，表示数组a(0<=a[i]<=1000)

Output

一个数表示答案，答案对10^9+7取模。

Sample Input

4
1 2 3 4

Sample Output

75

题解

考虑两个数$a_i,a_j(i<=j)$的贡献
有两种情况
第一种是$i+j<=n+1$的
这个时候有长度为$[j-i+1,j-i+1+2\ast (i-1)]$的区间能让他们做出贡献
化简完就是$[j-i+1,i+j-1]$，显然是这些长度区间中的质数个数
写出来就是$a[i]\ast a[j]\ast (s[i+j-1]-s[j-i])$，$s$表示质数的前缀和
第二种是$i+j>n+1$的
这个时候是长度为$[j-i+1,2n-(i+j-1)]$的区间
贡献写作$a[i]\ast a[j]\ast (s[2n-(i+j-1)]-s[j-i])$
这两个柿子都可以FFT做了

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<queue>
#include<vector>
#include<ctime>
#include<map>
#include<bitset>
#define LL long long
#define mp(x,y) make_pair(x,y)
#define lb long double
#define mod 1000000007
using namespace std;
const lb PI=acos(-1.0);
inline int read()
{
	int f=1,x=0;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}
inline void write(LL x)
{
	if(x<0)putchar('-'),x=-x;
	if(x>9)write(x/10);
	putchar(x%10+'0');
}
inline void pr1(int x){write(x);printf(" ");}
inline void pr2(LL x){write(x);puts("");}
const LL inv=500000004;
const int MAXN=210000;
bool vis[MAXN];int pr[MAXN],plen;
void getpr()
{
	memset(vis,true,sizeof(vis));
	for(int i=2;i<=200000;i++)
	{
		if(vis[i])pr[++plen]=i;
		for(int j=1;j<=plen&&i*pr[j]<=200000;j++)
		{
			vis[i*pr[j]]=false;
			if(!(i%pr[j]))break;
		}
	}
}
LL s[MAXN];
struct Complex
{
	long double r,i;
	Complex(){}
	Complex(long double _r,long double _i){r=_r;i=_i;}
	friend Complex operator +(Complex u,Complex v){return Complex(u.r+v.r,u.i+v.i);}
	friend Complex operator -(Complex u,Complex v){return Complex(u.r-v.r,u.i-v.i);}
	friend Complex operator *(Complex u,Complex v){return Complex(u.r*v.r-u.i*v.i,u.r*v.i+u.i*v.r);}
}A[MAXN*4],B[MAXN*4];int R[MAXN*4],L,ln;
void fft(Complex *y,LL len,int on)
{
	for(int i=0;i<len;i++)if(i<R[i])swap(y[i],y[R[i]]);
	for(int i=1;i<len;i<<=1)
	{
		Complex wn(cos(PI/i),sin(on*PI/i));
		for(int j=0;j<len;j+=(i<<1))
		{
			Complex w(1,0);
			for(int k=0;k<i;k++)
			{
				Complex u=y[j+k];
				Complex v=y[j+k+i]*w;
				y[j+k]=u+v;
				y[j+k+i]=u-v;
				w=w*wn;
			}
		}
	}
	if(on==-1)for(int i=0;i<len;i++)y[i].r/=(long double)len;
}
LL s1[MAXN*4],s2[MAXN*4],ans;
int n,a[MAXN*4];
void ad(LL &x,LL y){x+=y;if(x>=mod)x-=mod;}
void dl(LL &x,LL y){x-=y;if(x<0)x+=mod;}
int main()
{
//	freopen("a.in","r",stdin);
	getpr();vis[2]=false;
	for(int i=2;i<=200000;i++)s[i]=s[i-1]+vis[i];
	n=read();
	for(int i=1;i<=n;i++)scanf("%d",&a[i]);
	for(ln=1;ln<=2*n;ln<<=1)L++;
	for(int i=0;i<ln;i++)R[i]=(R[i>>1]>>1)|(i&1)<<(L-1);
	for(int i=0;i<ln;i++)A[i]=B[i]=Complex(a[i+1],0);
	fft(A,ln,1);fft(B,ln,1);
	for(int i=0;i<ln;i++)A[i]=A[i]*B[i];
	fft(A,ln,-1);
	for(int i=0;i<ln;i++)
	{
		s1[i]=(LL)(A[i].r+0.5)%mod;
		if(i==195752)
		{
			int tee;
			tee++;
		}
//		if(s1[i]<0)printf("GG %d %lld\n",i,s1[i]);
	}
//	for(int i=0;i<ln;i++)pr2(s1[i]);
	for(int i=0;i<n;i++)if(!(i&1))
	{
		LL temp=(s1[i]%mod);
		ad(ans,(LL)temp*s[i+1]%mod);
	}
	for(int i=n;i<2*n;i++)if(!(i&1))
	{
		LL temp=(s1[i]%mod);
		ad(ans,(LL)temp*s[2*n-(i+1)]%mod);
//		if(ans<0)printf("NO %d %lld\n ",i,ans);
	}
	for(int i=0;i<ln;i++)A[i]=Complex(a[i+1],0),B[i]=Complex(0,0);
	for(int i=0;i<n;i++)B[i]=Complex(a[n-i],0);
	fft(A,ln,1);fft(B,ln,1);
	for(int i=0;i<ln;i++)A[i]=A[i]*B[i];
	fft(A,ln,-1);
	for(int i=0;i<ln;i++)s2[i]=(LL)(A[i].r+0.5)%mod;
	for(int i=0;i<n;i++)if(!((n-i-1)&1))
	{
		int ll=n-i-1;
		dl(ans,2*s2[i]%mod*s[ll]%mod);
	}
//	for(int i=2;i<=n/2;i+=2)ad(ans,s1[i]*s[i+3]%mod);
//	for(int i=0;i<n/2;i++)
//	{
//		int temp=n/2-1-i;
//		dl(ans,s2[i]*s[temp]%mod);
//	}
	for(int i=1;i<n;i++)
	{
		ad(ans,(LL)a[i]*a[i+1]%mod);
		ad(ans,(LL)a[i+1]*a[i]%mod);
	}
	pr2(ans);
	return 0;
}