【题解】【原创题目】第〇类循环数·行
传送门：第〇类循环数·行
【题目描述】
给出 (n,m,l,k,x)，求 (frac{1}{nml}sumlimits_{a=1}^{n}sumlimits_{b=1}^{m}sumlimits_{c=1}^{l}f(k,(ax^2+bx+c)mod k) pmod{1004535809})，其中 (f(i,j)) 定义如下：
【分析】
【Subtask #1】
抄题目描述总会吧。
本测试包轻微卡常，需要提前把 (f) 全部预处理出来。
复杂度 (O(k^2+n^3T))，期望得分：(1pts)。
【Code #1】
in(T),inv[1]=1;
for(Re i=2;i<=M-3;++i)inv[i]=(LL)inv[P%i]*(P-P/i)%P;
for(Re i=1;i<=M-3;++i){
f[i][0]=i*i+i;
for(Re j=1;j<=i-1;++j)
f[i][j]=((j?(
(LL)inv[i]*inv[j]%P*f[i-1][j-1]%P+
(LL)inv[i]*f[i-1][j-1]%P+
(LL)inv[j]*f[i-1][j-1]%P+
f[i-1][j-1]
)%P:0)+
(j<i-1?((LL)inv[i]*f[i-1][j]%P+f[i-1][j])%P:0))%P;
}
while(T--){
in(n),in(m),in(l),in(K),in(x),x%=K,x2=(LL)x*x%K,ans=0;
for(Re i=1;i<=n;++i)
for(Re j=1;j<=m;++j)
for(Re k=1;k<=l;++k)
(ans+=f[K][((LL)i*x2%K+(LL)j*x%K+k)%K])%=P;
printf("%lld
",(LL)ans*Inv(n)%P*Inv(m)%P*Inv(l)%P);
}
【Subtask #2】
存在特殊性质 (k|x)，说明 (q) 与 (a,b) 无关，即求 (frac{1}{l}sum_{c=1}^{l}f(k,cmod k))，发现暴力枚举有很多重复计算，让 (l) 对 (k) 取个模即可降到 (O(k)) 。
本测试包轻微卡常。
复杂度 (O(k^2+kT))，期望得分：(1+9=10pts)。
【Code #2】
while(T--){
in(n),in(m),in(l),in(K),in(x),x%=K,x2=(LL)x*x%K,ans=0;
Re tmp=0;
for(Re k=0;k<K;++k)(tmp+=f[K][k])%=P;
(ans+=(LL)l/K%P*tmp%P)%=P;
for(Re k=1;k<=(l%K);++k)(ans+=f[K][k])%=P;
printf("%lld
",(LL)ans*Inv(l)%P);
}
【Subtask #3】
用类似 ( ext{Subtask 2}) 的降级做法，依次处理两个东西：(g_1(a)=sum_{i=1}^{l}f_{k}(a+i),) (g_2(a)=sum_{i=1}^{m}g_1(a+ix))，最后答案为 (sum_{i=1}^{n}g_{2}(ix^2)) 。
本做法较卡常，需要把加法取模优化掉。
复杂度 (O(k^2+k^2T))，期望得分：(1+9+11=21pts)。
【Code #3】
inline void mo(Re &x,Re mod=P){x=x>=mod?x-mod:x;}
inline int Mo(Re x,Re mod=P){return x>=mod?x-mod:x;}
int main(){
while(T--){
in(n),in(m),in(l),in(K),in(x),x%=K,x2=(LL)x*x%K,ans=0;
for(Re a=0;a<K;++a){
Re tmp=0;
if(l>=K)for(Re k=0,Tmp=a;k<K;++k)mo(tmp+=f[K][Tmp]),mo(++Tmp,K);
g1[a]=(LL)l/K%P*tmp%P;
for(Re k=1,Tmp=Mo(a+1,K),ed=l%K;k<=ed;++k)mo(g1[a]+=f[K][Tmp]),mo(++Tmp,K);
}
for(Re a=0;a<K;++a){
Re tmp=0;
if(m>=K)for(Re k=0,Tmp=a;k<K;++k)mo(tmp+=g1[Tmp]),mo(Tmp+=x,K);
g2[a]=(LL)m/K%P*tmp%P;
for(Re k=1,Tmp=Mo(a+x,K),ed=m%K;k<=ed;++k)mo(g2[a]+=g1[Tmp]),mo(Tmp+=x,K);
}
Re tmp=0;
if(n>=K)for(Re k=0,Tmp=0;k<K;++k)mo(tmp+=g2[Tmp]),mo(Tmp+=x2,K);
ans=(LL)n/K%P*tmp%P;
for(Re k=1,Tmp=x2,ed=n%K;k<=ed;++k)mo(ans+=g2[Tmp]),mo(Tmp+=x2,K);
printf("%lld
",(LL)ans*Inv(n)%P*Inv(m)%P*Inv(l)%P);
}
}
【Subtask #4】
观察 (f) 递推式，首先通分合并同类项：
两边同时除以 ((i+1)(j+1))：
发现是组合数的杨辉三角递推式。设 (g(i,j)=frac{f(i,j)}{(i+1)(j+1)})，则有 (g(i,j)=g(i-1,j-1)+[j eq i-1]g(i-1,j)) 。
注意此处的细节：通常组合数递推为 (dbinom{n}{m}=dbinom{n-1}{m-1}+dbinom{n-1}{m})，当 (m=n) 时 (dbinom{n-1}{m}) 才为 (0)，但 (g(i-1,j)) 前面的系数为 ([j eq i-1])，因此 (g(i,j)) 应该等于(dbinom{i}{j+1}) 而不是 (dbinom{i}{j}) ，否则就与递推式不相符。
得出 (f(i,j)=(i+1)(j+1)dbinom{i}{j+1})，这个东西算出来的 (f(i,0)) 和前面定义是符合的。
注意到 (n) 比较小，预处理出组合数后直接沿用 ( ext{Subtask 3}) 的做法即可。
复杂度 (O(min(n,k)kT))，期望得分：(1+9+11+19=40pts)。
【Code #4】
代码略。
【Subtask #5】
推出 (f) 通项式后直接沿用 ( ext{Subtask 2}) 的做法。
此做法良心出题人没有卡常。
复杂度 (O(kT))，期望得分：(9+18=27pts)，结合前面算法可得 (58pts)。
【Code #5】
代码略。
【Subtask #6】
前置知识：单位根反演。【常见公式汇总】
按照套路，枚举 (f) 的参数并交换和式：
([(ax^2+bx+c)mod k=d]) 可化为 ([k|(ax^2+bx+c-d)])，套用单位根反演的柿子，得到：
提一个 (w_{k}^{-id}) 出来，并交换和式：
右边那一坨显然可以简化为三个等比数列和的乘积，记一个 (calc) 函数：
用组合数吸收公式去掉 (f) 柿子里的 ((j+1))，即 (f(i,j)=(i+1)(j+1)dbinom{i}{j+1}=i(i+1)dbinom{i-1}{j})，然后代进去用二项式定理合并：
本测试包不卡常，可以直接暴算。
复杂度 (O(Tklog n)))，期望得分：(91pts) 。
【Code #6】
inline int calc(Re x,Re n){//注意等比数列求和要特判公比等于1的情况
return x==1?n:(LL)x*(1-mi(x,n)+P)%P*Inv((1-x+P)%P)%P;
}
inline int calc(Re w){
return (LL)calc(mi(w,x2),n)*calc(mi(w,x),m)%P*calc(w,l)%P;
}
int main(){
in(T);
while(T--){
in(n),in(m),in(l),in(K),in(x),x%=K,x2=(LL)x*x%K,ans=0;
omega[0]=1,omega[1]=mi(G,(P-1)/K);
for(Re i=2;i<=K;++i)omega[i]=(LL)omega[i-1]*omega[1]%P;
for(Re i=0;i<=K-1;++i)
(ans+=(LL)mi(omega[K-i]+1,K-1)*calc(omega[i])%P)%=P;
ans=(LL)ans*(K+1)%P;
printf("%lld
",(LL)ans*Inv(n)%P*Inv(m)%P*Inv(l)%P);
}
}
【Subtask #7】
最后两个测试点需要大力卡常。
大常数主要来源于快速幂，枚举一个 (i) 要调用 (9) 次，极限数据跑了约 (5.8s)（洛谷上要快 (0.6s)）。
实际上其中有 (5) 次都在算单位根 (w_{k}^{1}) 的若干次幂，由于 (w_{k}^{i}=w_{k}^{imod k})，只要预处理出 (0sim k-1) 次幂就能 (O(1)) 查询了。
然后是 (3) 次计算逆元，同样可以预处理，具体做法见 乘法逆元 (2) ( ext{[P5431]}) 。
其实单位根还可以放到询问前面预处理，对于不同的 (k) 根据 (w_{k}^{a}=w_{kn}^{an}) 这一性质直接得到。
但瓶颈不在于此，常数影响不大。
复杂度 (O(Tklog k)))，期望得分：(100pts) 。
【Code #7】
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<cmath>
#define LD double
#define LL long long
#define Re register int
using namespace std;
const int N=1e9+3,M=1048576+3,P=1004535809,G=3;
int n,m,l,K,T,x2,ans,Sm[M],invo[M],omega[M];LL x;
template<typename T>inline void in(T &x){
int f=0;x=0;char ch=getchar();
while(ch<'0'||ch>'9')f|=ch=='-',ch=getchar();
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
x=f?-x:x;
}
inline int mi(Re x,Re k){
Re s=1;
while(k){
if(k&1)s=(LL)s*x%P;
x=(LL)x*x%P,k>>=1;
}
return s;
}
inline int Inv(Re x){return mi(x,P-2);}
inline void mo(Re &x,Re mod=P){x=x>=mod?x-mod:x;}
inline int Mo(Re x,Re mod=P){return x>=mod?x-mod:x;}
inline int calc(Re i,Re n){//注意等比数列求和要特判公比等于1的情况
return (omega[i]==1?n:(LL)omega[i]*(omega[(LL)i*n%K]-1+P)%P*invo[i]%P);
}
int main(){
in(T);
while(T--){
in(n),in(m),in(l),in(K),in(x),x%=K,x2=(LL)x*x%K,ans=0;
omega[0]=Sm[0]=1,Sm[1]=(omega[1]=mi(G,(P-1)/K))-1;
for(Re i=2;i<=K;++i)omega[i]=(LL)omega[i-1]*omega[1]%P,Sm[i]=(LL)Sm[i-1]*(omega[i]-1)%P;
Re tmp=Inv(Sm[K-1]);
for(Re i=K-1;i>=1;--i)invo[i]=(LL)tmp*Sm[i-1]%P,tmp=(LL)tmp*(omega[i]-1)%P;
for(Re i=0,mp1=0,mp2=0;i<K;++i)
mo(ans+=(LL)mi(omega[K-i]+1,K-1)*calc(mp1,n)%P*calc(mp2,m)%P*calc(i,l)%P),mo(mp1+=x2,K),mo(mp2+=x,K);
ans=(LL)ans*(K+1)%P;
printf("%lld
",(LL)ans*Inv(n)%P*Inv(m)%P*Inv(l)%P);
}
}
【题外话—命题报告】
出题人：辰星凌
验题人：鏡音リン
（感谢神仙 ( ext{karry}) 解惑，( ext{ezoixx130}) 划水验题）
(7) 月 (28) 晚，翻看自己整理的数学知识点汇总时，发现单位根反演这个东西题不多，就想搞一个出来玩玩儿。但数学题并不好出，直接放柿子的话会被喷惨。但我又比较菜，斗不出组合意义，顿时感觉希望渺茫。
抱着试一试的心态，随便选了个方向：求 (sum_{i,j}f(imod j))，其中 (f) 先待定成普通多项式和下降幂多项式两种情况。推推推，最后弄出来一长串柿子，用了各种各样的技巧，结果最后还是 (O(n^2)) ！甚至常数还更大了！
人傻了....
后来调整数据范围，另外搞了个柿子，并斗出一个简单的 (f) 函数，消去了很多东西，似乎能出成一道题了。
但这还不够，因为只是给出了一个柿子让做题者去推，这样是没有意义的。而且单反就那两种变化，由于白兔之舞的存在，似乎无论怎么出都会感觉很套路，多半会被狂 ( ext{D}) 。
那就....加上构造！
把 (f) 函数里的组合数拆开，弄出一个递推式，在题面上加点迷惑元素，让出题者去尝试构造 (f)，并在边界处设置一个坑人的细节问题。
这样看起来似乎就不那么板了~(≥▽≤)/~
（其实还是很套路，只要对组合数够熟悉就能一眼秒）。
另外，关于用 (f) 递推式求通项表示式，一开始是想用生成函数暴力解的，但貌似要求导积分，小蒟蒻对这一类的东西不太熟悉，如果有直接硬推成功解出来的巨佬，请接受蒟蒻真诚的膜拜stO
emm....从巨佬那里得知 (f) 可以在 ( ext{oeis}) 上找到，自闭....
还是不加入到公开赛了，免得被骂搬原式。
( ext{updata:}) 经 ( ext{jklover}) 巨佬提醒，( ext{Subtask 3,4}) 里的做法是可以直接用 ( ext{DFT}) 优化的，总复杂度也为 (O(nlog n))，虽然常数大一些，但实际运行效果只比前面的无常数单 (log) 做法慢 (0.3s)，过于柯啪.....
懒得写新 std 和重新造 Subtask 卡常了，就酱吧。