LGP6825题解
科技的力量！！！！！！我德意志科技天下第一！！！
这是一篇需要一点儿科技的题解，但实际上这个科技我认为甚至算不上科技，太 simple 了。
首先是推柿子：
\[\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)^{i+j}
\]
\[\sum_{d=1}^n\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=d]d^{i+j}
\]
\[\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac n d \rfloor}\sum_{j=1}^{\lfloor \frac n d \rfloor}[\gcd(i,j)=1](d^d)^{i+j}
\]
\[\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac n d \rfloor}\sum_{j=1}^{\lfloor \frac n d \rfloor}\sum_{k|i,k|j}\mu(k)(d^d)^{i+j}
\]
\[\sum_{d=1}^n\sum_{k=1}^{\lfloor \frac n d \rfloor}\sum_{i=1}^{\lfloor \frac n {dk} \rfloor}\sum_{j=1}^{\lfloor \frac n {dk} \rfloor
}\mu(k)(d^{dk})^{i+j}\]
设
\[f(d,n)=\sum_{i=1}^n\sum_{j=1}^nd^{i+j}
\]
特别的，当 \(d=1\) 时，\(f(1,n)=n^2\)
否则上面那玩意儿就是：
\[(\sum_{i=1}^nd^i)^2
\]
等比数列求和，在 \(O(\log n)\) 的时间内算出来。
回到原柿：
\[\sum_{d=1}^n\sum_{k=1}^{\lfloor \frac n d \rfloor}\mu(k)f(d^{dk},\lfloor \frac n {dk} \rfloor)
\]
考虑将 \((d^d)^k\) 带入 \(f\)。可以发现 \(f\) 只与 \(d^d\) 的若干次方有关，并且这个幂的上界是 \(O(k \times \lfloor \frac n {dk} \rfloor)=O(n)\)。
考虑对于所有的 \(d^d\) 预处理光速幂数组。可以参考这里，复杂度是 \(O(\sqrt [4] {lim})-O(1)\)。本题的 \(lim\) 为 \(1.5 \times 10^6\)，不过实际上当 \(d=1\) 时不需要光速幂数组，所以实际上上界不到 \(10^6\)。别小看这一点，\(2^{20}>10^6\)，这导致了空间和时间复杂度得到了改善。
于是我们发现唯一的问题是求 \((d^d)^k-1\) 的逆元。我们使用离线求逆元，这里使用的是封装好的离线求逆元。
以及，这个屑题卡空间，所以我们对 \([1,n]\) 分成若干段，对每一段使用离线求逆元和光速幂。离线求逆元当缓存满了的时候需要清空。
我们发现，每段的长度最好是 \(8\)，因为离线求逆元最多只会执行 \(10^5\) 次，且光速幂的预处理不会炸空间。还有一个好处就是，这使得离线求逆元的缓存区很短，能够卡进执行文件的缓存，大大降低了常数。
总复杂度为 \(O(n\sqrt [4] {n}+n\ln n)\)。
优化1：只枚举 \(\mu(k) \neq 0\) 的 \(k\)，常数降低。
优化2：当 \(d\) 较大时，能够与其匹配的 \(k\) 过少，此时使用快速幂比光速幂更优。
不过不知道为啥只跑了 1.16s，跑得飞快。。。
code:
#include<cstdio>
#include<vector>
typedef unsigned uint;
typedef __uint128_t L;
typedef unsigned long long ull;
const uint M=1.5e6+5;
uint P,top,pri[M],mu[M];bool zhi[M];
struct FastMod{
ull b,m;
FastMod(ull b):b(b),m(ull((L(1)<<64)/b)){}
friend inline ull operator%(const ull&a,const FastMod&mod){
ull r=a-(L(mod.m)*a>>64)*mod.b;
return r>=mod.b?r-mod.b:r;
}
}mod(2);
struct Pow{
std::vector<uint>p1,p2,p3,p4;
inline void clear(){
std::vector<uint>().swap(p1);std::vector<uint>().swap(p2);
std::vector<uint>().swap(p3);std::vector<uint>().swap(p4);
}
inline void init(register uint x){
register uint i;
p1.reserve(33);p2.reserve(33);p3.reserve(33);p4.reserve(33);
p1[0]=p2[0]=p3[0]=p4[0]=1;
for(p1[1]=x,i=2;i<=32;++i)p1[i]=1ull*p1[i-1]*x%mod;
for(x=p2[1]=p1[32],i=2;i<=32;++i)p2[i]=1ull*p2[i-1]*x%mod;
for(x=p3[1]=p2[32],i=2;i<=32;++i)p3[i]=1ull*p3[i-1]*x%mod;
for(x=p4[1]=p3[32],i=2;i<=32;++i)p4[i]=1ull*p4[i-1]*x%mod;
}
inline uint pow(const uint&x){
return 1ull*p1[x&31]*p2[x>>5&31]%mod*p3[x>>10&31]%mod*p4[x>>15&31]%mod;
}
Pow(){
clear();
}
}p[M>>1];
inline uint pow(uint a,uint b){
uint ans=1;
for(;b;b>>=1,a=1ull*a*a%mod)if(b&1)ans=1ull*ans*a%mod;
return ans;
}
inline void sieve(const uint&M){
register uint i,j,x;mu[1]=1;
for(i=2;i<=M;++i){
if(!zhi[i])mu[pri[++top]=i]=-1;
for(j=1;j<=top&&(x=i*pri[j])<=M;++j){
zhi[x]=true;if(!(i%pri[j]))break;mu[x]=-mu[i];
}
}
}
inline uint calc(const uint&d,const uint&D,const uint&len,const uint&k,const uint&inv){
uint S=(1ull*D*p[d].pow(len*k)%mod-p[d].pow(k)+P)*inv%mod;
return 1ull*S*S%mod;
}
struct GetAns{
uint ans,len,d[131073],k[131073],l[131073],f[131073],S[131073];bool o[131073];
inline void flush(){
uint i,s;if(!len)return;
for(i=1;i<=len;++i)S[i]=1ull*S[i-1]*f[i]%mod;s=pow(S[len],P-2);
do{
if(o[len])ans=(ans+calc(d[len],f[len]+1,l[len],k[len],1ull*s*S[len-1]%mod))%mod;
else ans=(ans+P-calc(d[len],f[len]+1,l[len],k[len],1ull*s*S[len-1]%mod))%mod;
s=1ull*s*f[len]%mod;
}while(--len);
}
inline void Insert(const uint&cd,const uint&ck,const uint&cl,const uint&cf,const bool&co){
++len;d[len]=cd;k[len]=ck;l[len]=cl;f[len]=cf;o[len]=co;
if(!(len&131071))flush();
}
}q;
inline uint Solve(const uint&n){
const uint&lim4=n/4,&lim3=n/3,lim2=n/2;
register uint i,k,x,ans=0,lst=2;q.ans=q.len=0;
for(i=1;i<=n;++i){
if(mu[i]==1)ans=(ans+1ull*(n/i)*(n/i))%mod;
else if(mu[i])ans=(ans+P-1ull*(n/i)*(n/i)%mod)%mod;
}
for(i=2;(i<<1)<=n;++i){
if(!(i&7)){
q.flush();while(lst^i)p[lst++].clear();
}
p[i].init(pow(i,i));
for(k=1;i*k<=n;++k){
if(mu[k]){
x=p[i].pow(k)-1;
if(!x){
if(mu[k]==1)ans=(ans+1ull*(n/(i*k))*(n/(i*k)))%mod;
else if(mu[k])ans=(ans+P-1ull*(n/(i*k))*(n/(i*k))%mod)%mod;
}
else q.Insert(i,k,n/(i*k),x,mu[k]==1);
}
}
}
q.flush();while(lst^i)p[lst++].clear();
for(ans=(ans+q.ans)%mod;i<=n;++i)k=pow(i,i),ans=(ans+1ull*k*k)%mod;
return ans;
}
signed main(){
uint T,n1,n2,m1,m2;q.S[0]=1;
scanf("%u",&T);
if(T==1)scanf("%u%u",&n1,&m1),sieve(n1),mod=FastMod(P=m1),printf("%u",Solve(n1));
else{
scanf("%u%u%u%u",&n1,&m1,&n2,&m2);
sieve(n1>n2?n1:n2);
mod=FastMod(P=m1);printf("%u\n",Solve(n1));
mod=FastMod(P=m2);printf("%u",Solve(n2));
}
}
闲话：这个 \(O(n\log n)\) 可比官方题解的 \(O(n\log n)\) 靠谱多了。
