[原创题] 雾
题目：http://162.105.80.126/contest/%E3%80%90%E5%BC%B1%E7%9C%81%E8%83%A1%E7%AD%96%E3%80%91Round%20%233/%E9%9B%BE
感谢绍兴一中三位神犇前来捧(nue)场，感谢山东神犇。。。。。。
让我们来看一下题目，
$$Ans = \sum_{i=1}^{n}\sum_{j=1}^{i-1}\prod_{k=1}^{min(i-j,j)}[\sum_{l=1}^k a(l)]^{\varphi(j-k)+\varphi(i-j-k)}\prod_{i=min(i-j,j)+1}^{max(i-j,j)}f(i,j,k)$$
若$i-j>j$ , $f(i,j,k)=[\sum_{l=1}^ka(l)]^{\varphi(i-j-k)}$
否则 $f(i,j,k)=[\sum_{l=1}^ka(l)]^{\varphi(j-k)}$。
注意到$\sum_{l=1}^k a(l)$和其他的式子并没有关系，所以可以预先前缀和，下面我们所说的$a(i)$都是前缀和后的$a(i)$。
$$Ans = \sum_{i=1}^{n}\sum_{j=1}^{i-1}\prod_{k=1}^{min(i-j,j)}a(k)^{\varphi(j-k)+\varphi(i-j-k)}\prod_{i=min(i-j,j)+1}^{max(i-j,j)}f(i,j,k)$$
$10 points$： 按照定义$O(n^3)$
显然这个式子是难以优化的，我们考虑加以变形，注意到$f(i,j,k)$函数的定义恰好和前面的连乘的内容相似，我们可以将$f(i,j,k)$代入化简，然后将左侧的连乘拆开为
$$\prod_{k=1}^{min(i-j,j)}a(k)^{\varphi(j-k)}\cdot a(k)^{\varphi(i-j-k)}$$
然后做变形
$$\prod_{k=1}^{min(i-j,j)}a(k)^{\varphi(j-k)} \prod_{k=1}^{min(i-j,j)}a(k)^{\varphi(i-j-k)}$$
然后发现 $f(i,j,k)$ 其实可以接到左侧式子或者右侧式子（取决于$i-j$和$j$的大小关系）,整理得到：
$$Ans = \sum_{i=1}^{n}\sum_{j=1}^{i-1}\prod_{k=1}^{j}a(k)^{\varphi(j-k)}\prod_{i=1}^{i-j}a(k)^{\varphi(i-j-k)}$$
然后，仔细观察后面那一大坨连乘和前面的求和无关。
设$C(n) = \prod_{i=1}^{n}a(i)^{\varphi(n-i)}$
然后发现原式变成了：
$$Ans = \sum_{i=1}^{n}\sum_{j=1}^{i-1}{C(i)\cdot C(i-j)}$$
(注意这里的卷积形式，要在首项补$0$,保证$j$是从$i-1$结束)
$30 points$:
假如你不会$FFT$也不会$NTT$那么你可以$O(n^2)$预处理，计算时直接用。
总复杂度：$O(n^2)$
标准的卷积形式，可以用$NTT$ $O(nlogn)$ 解决。
然后回来看$C(n)$，左右用原根代换：
$$ gn^{c'(n)} = \prod_{i=1}^{n}gn^{a'(i) \cdot {\varphi(n-i)}}$$
求出 $C'(n)$ 之后我们就可以 $O(nlogn)$ 快速幂求出 $C(n)$ 了。
而注意到其指数恰为卷积的形式，但是根据费马小定理，应该取$mod-1$的余，模数不太对，不能用 $NTT%$，而 $FFT$ 又炸精度。
考虑别的方法，由于数据中的指数是很小的，最大的情况下数字可以是$4 \cdot 115000000000$，取一个超级大的质数，我是$1899956092796929$，原根为$7$，然后就可以了吗，朴素乘法会溢出，快速乘或者用$long \ double$的大数相乘取余函数。（然而快速乘多$logn$；$long \ double$乘有概率挂，所以不用是上策）
此外，确定一个数字是原根的多少次方时可以用$BSGS$ $O(n\sqrt{n})$注意这里的n小于等于$5\cdot 10^4$.
总效率$O(nlog^2n)$还是不能通过所有的数据呀，预计$40$~$50$分。
$50points$:
最后，我们发现我们可以选两个$10^9$级的满足条件的质数，最后$CRT$合并。
然后就可以$50$了（此处应$orz$ $yjp$）
$70points$:
继续改进我们的算法，可以发现更本不用$BSGS$，完全可以预处理出所有的值，用一个$map$或$hash$存下来，取时$O(logn)$查询，总体复杂度：$O(nlogn)$。
$100points$:
可以发现，之所以$70$分，并不是因为复杂度，而是$NTT$太慢了，注意到指数很小，所以可以$FFT$，或者卡卡常数即可得到全分。
注意精度。
#include <cstdio> #include <cmath> #include <ctime> #include <algorithm> #include <map> #define N 400010 #define LL unsigned long long #define mod 1004535809LL #define gn 3 #define mod2 998244353LL using namespace std; const double pi=acos(-1.0); struct Ex{ double r, i; Ex(){} Ex(double r, double i):r(r), i(i){} friend Ex operator+(const Ex &a,const Ex &b){return Ex(a.r+b.r,a.i+b.i);} friend Ex operator-(const Ex &a,const Ex &b){return Ex(a.r-b.r,a.i-b.i);} friend Ex operator*(const Ex &a,const Ex &b){return Ex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);} }a0[N],b0[N],c0[N]; map<LL,int> MT; inline LL qmult(LL a,LL b,LL P){ LL ans=0; for(;b;b>>=1,a=(a<<1)%P) if(b&1) ans=(ans+a)%P; return ans; } inline void swap(LL &a,LL &b){ if(a^b) a^=b,b^=a,a^=b; } inline LL qpow(LL x,LL n,LL P){ LL ans=1; for(;n;n>>=1,x=x*x%P) if(n&1) ans=ans*x%P; return ans; } int R[N]; inline LL Mod(LL x,LL P){ if(x>=P) x-=P; if(x>=P) x-=P; if(x>=P) x%=P; return x; } LL wt[N]; inline void fnt(LL *x,int n,int t,LL P){ for(int i=0;i<n;i++) if(i<R[i]) swap(x[i],x[R[i]]); for(int m=1;m<n;m<<=1){ LL wn=qpow(gn,(P-1)/(m<<1),P); if(t==-1) wn=qpow(wn,P-2,P); wt[0]=1; for(int i=1;i<m;i++) wt[i]=Mod(wt[i-1]*wn,P); for(int k=0;k<n;k+=(m<<1)) for(int i=0;i<m;i++){ LL &A=x[i+m+k]; LL &B=x[i+k],t=Mod(A*wt[i],P); A=Mod(B-t+P,P); B=Mod(B+t,P); } } if(t==-1){ LL tmp=qpow(n,P-2,P); for(int i=0;i<n;i++) x[i]=x[i]*tmp%P; } } inline void fft(Ex x[],int n,int t){ for(int i=0;i<n;i++) if(i<R[i]) swap(x[i],x[R[i]]); for(int m=1;m<n;m<<=1){ Ex wn(cos(pi/m*t),sin(pi/m*t)); for(int k=0;k<n;k+=(m<<1)){ Ex wt(1,0); for(int i=0;i<m;i++,wt=wt*wn){ Ex &A=x[k+m+i],&B=x[k+i],t=wt*A; A=B-t; B=B+t; } } } if(t==-1) for(int i=0;i<=n;i++) x[i].r/=n; } LL a1[N],b1[N],c1[N]; inline void multpoly2(LL *A,LL *B,LL *C,int n){ int nt,L=0; for(nt=1;nt<=n+n;nt<<=1) L++; for(int i=0;i<nt;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1)); for(int i=0;i<nt;i++) a0[i].r=b0[i].r=c0[i].r=0; for(int i=0;i<nt;i++) a0[i].i=b0[i].i=c0[i].i=0; for(int i=0;i<n;i++) a0[i].r=A[i],b0[i].r=B[i]; fft(a0,nt,1); fft(b0,nt,1); for(int i=0;i<nt;i++) c0[i]=a0[i]*b0[i]; fft(c0,nt,-1); for(int i=0;i<n;i++) C[i]=(c0[i].r+0.5); } inline void multpoly1(LL *A,LL *B,LL *C,int n,LL P){ int nt,L=0; for(nt=1;nt<=n+n;nt<<=1) L++; for(int i=0;i<nt;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1)); for(int i=0;i<nt;i++) a1[i]=b1[i]=c1[i]=0; for(int i=0;i<n;i++) a1[i]=A[i],b1[i]=B[i]; fnt(a1,nt,1,P); fnt(b1,nt,1,P); for(int i=0;i<nt;i++) c1[i]=a1[i]*b1[i]%P; fnt(c1,nt,-1,P); for(int i=0;i<n;i++) C[i]=c1[i]; } inline LL CRT(LL a1,LL a2){ LL ans=0,M=mod*(LL)mod2; ans+=qmult(qmult(a1,mod2,M),qpow(mod2,mod-2,mod),M); ans+=qmult(qmult(a2,mod,M),qpow(mod,mod2-2,mod2),M); return (ans%M+M)%M; } int n,phi[N]; LL ans=0,a[N],b[N]; LL c[N],ap[N]; double ti; int main(){ scanf("%d",&n); phi[1]=1; for(int i=2;i<=n;i++){ if(phi[i]) continue; for(int j=i;j<=n;j+=i){ if(!phi[j]) phi[j]=j; phi[j]=phi[j]/i*(i-1); } } for(int i=0;i<=100000;i++) MT[qpow(gn,i,mod)]=i; for(int i=0;i<n;i++) scanf("%lld",&a[i]); for(int i=1;i<n;i++) a[i]=(a[i]+a[i-1])%mod; for(int i=0;i<n;i++) ap[i]=MT[a[i]]; for(int i=0;i<n;i++) b[i]=phi[i]%23; multpoly2(ap,b,c,n); for(int i=0;i<n;i++) c[i]=c[i]%(mod-1); for(int i=0;i<n;i++) c[i]=qpow(gn,c[i],mod); for(int i=n-1;~i;i--) c[i+1]=c[i],c[i]=0; multpoly1(c,c,c,n+1,mod); for(int i=0;i<=n;i++) ans=(ans+c[i])%mod; printf("%lld\n",ans); return 0; }
