题目分析
莫比乌斯反演
所谓的sgcd(i,j)sgcd(i,j)sgcd(i,j),就是gcd(i,j)gcd(i,j)gcd(i,j)除以其最小的一个质因子。我们记g(x)=(xminpri(x))Kg(x)=(\frac{x}{minpri(x)})^Kg(x)=(minpri(x)x)K,答案就是求∑i=1n∑j=1ng(gcd(i,j))\sum_{i=1}^n \sum_{j=1}^n g(gcd(i,j))∑i=1n∑j=1ng(gcd(i,j))
既然涉及到gcd,就用莫比乌斯反演搞一搞:
∑d=1ng(d)∑i=1n∑j=1n[gcd(i,j)=d]\sum_{d=1}^n g(d)\sum_{i=1}^n\sum_{j=1}^n [gcd(i,j)=d]∑d=1ng(d)∑i=1n∑j=1n[gcd(i,j)=d]
∑d=1ng(d)∑i=1⌊nd⌋∑j=1⌊nd⌋∑t∣i,t∣jμ(t)\sum_{d=1}^n g(d)\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{t|i,t|j} \mu(t)∑d=1ng(d)∑i=1⌊dn⌋∑j=1⌊dn⌋∑t∣i,t∣jμ(t)
∑d=1ng(d)∑t=1⌊nd⌋μ(t)(⌊ndt⌋)2\sum_{d=1}^n g(d) \sum_{t=1}^{\lfloor \frac{n}{d} \rfloor} \mu(t) (\lfloor \frac{n}{dt} \rfloor)^2∑d=1ng(d)∑t=1⌊dn⌋μ(t)(⌊dtn⌋)2
∑T=1n(⌊nT⌋)2∑d∣Tg(d)μ(Td)\sum_{T=1}^n (\lfloor \frac{n}{T} \rfloor)^2 \sum_{d|T} g(d) \mu(\frac{T}{d})∑T=1n(⌊Tn⌋)2∑d∣Tg(d)μ(dT)
前面的TTT部分可以数论分块,所以关键是求后面那一块的前缀和。
杜教筛
求(g∗μ)(i)(g * \mu)(i)(g∗μ)(i)的前缀和,想到杜教筛。杜教筛中使用的辅助函数,就用μ\muμ常用的辅助函数III即可(这个函数是每一位都为1)
所以设S(n)S(n)S(n)表示那个卷积函数前nnn项的前缀和,有:
I(1)S(n)=∑i=1n(g∗μ∗I)(i)−∑i=2nI(i)S(⌊ni⌋)I(1)S(n)=\sum_{i=1}^n (g* \mu * I)(i) - \sum_{i=2}^n I(i)S(\lfloor \frac{n}{i} \rfloor)I(1)S(n)=∑i=1n(g∗μ∗I)(i)−∑i=2nI(i)S(⌊in⌋)
也就是:
S(n)=∑i=1ng(i)−∑i=2nS(⌊ni⌋)S(n) = \sum_{i=1}^n g(i) - \sum_{i=2}^n S(\lfloor \frac{n}{i} \rfloor)S(n)=∑i=1ng(i)−∑i=2nS(⌊in⌋)
又可以数论分块,关键问题又变成了求g(i)g(i)g(i)的前缀和
min_25筛
嘛,其实我还是喜欢叫这个东西扩展埃氏筛(因为那个下划线在用markdown的时候有点不方便,雾)。
设s(i)s(i)s(i)表示前iii个数中的质数个数。
f(i)f(i)f(i)表示前iii个数中质数的kkk次方和。
g(i)g(i)g(i)表示前iii个数中合数都去掉自己的最小质因子后留下的数字的kkk次方和(和前面的ggg意义不同了)。
前两个很简单,转移式子分别是对于每个质因子ppp去掉以其为最小质因子的合数的贡献:s(i)−=s(⌊ip⌋)−s(p−1)s(i)-=s(\lfloor \frac{i}{p} \rfloor)-s(p-1)s(i)−=s(⌊pi⌋)−s(p−1)和f(i)−=pk(f(⌊ip⌋)−f(p−1))f(i)-=p^k(f(\lfloor \frac{i}{p} \rfloor)-f(p-1))f(i)−=pk(f(⌊pi⌋)−f(p−1))
发现每次f(i)f(i)f(i)丢掉的f(⌊ip⌋)−f(p−1)f(\lfloor \frac{i}{p} \rfloor)-f(p-1)f(⌊pi⌋)−f(p−1)是一些应该加进g(i)g(i)g(i)中的贡献,所以每次把g(i)g(i)g(i)加上这一串。
那么原式子中求的ggg,就是现在求出来的这些东西中的g(i)+s(i)g(i)+s(i)g(i)+s(i)(sss是算质数的贡献)
现在又双叒叕有个问题,扩埃筛执行时,需要处理前iii个数的函数贡献和。也就是我们要求前iii个数的KKK次方和。
第二类斯特林数
这个问题我知道!用伯努利数就可以了!ヾ(≧∇≦*)ゝ
(╯°Д°)╯
……个鬼啦!模数没逆元!
那就用第二类斯特林数吧,∑i=1niK=∑i=1KS2(K,i)(n+1)i+1‾i+1\sum_{i=1}^n i^K=\sum_{i=1}^K S_2(K,i)\frac{(n+1)^{\underline {i+1}}}{i+1}∑i=1niK=∑i=1KS2(K,i)i+1(n+1)i+1
(如果你以前没听说过这个式子,请别把下降幂看成幂了)
代码
#include<bits/stdc++.h>
using namespace std;
typedef unsigned int uint;
const int N=31630;
int K;uint n,ans,sqn;int vis1[N],vis2[N];
uint s1[N],s2[N],f1[N],f2[N],g1[N],g2[N],S2[55][55],mp1[N],mp2[N];
uint ksm(uint x,uint y) {
uint re=1;
for(;y;y>>=1,x=x*x) if(y&1) re=re*x;
return re;
}
void prework() {
S2[0][0]=1;
for(uint i=1;i<=K;++i)
for(uint j=1;j<=i;++j)
S2[i][j]=S2[i-1][j-1]+j*S2[i-1][j];
}
uint power_sum(uint x) {
uint re=0;
for(uint i=0;i<=K;++i) {
uint kl=S2[K][i];
for(uint j=0;j<=i;++j)
if((x+1-j)%(i+1)) kl*=x+1-j;
else kl*=(x+1-j)/(i+1);
re+=kl;
}
return re;
}
void min25() {
if(n<=1) return;
sqn=sqrt(n);
for(uint i=1;i<=sqn;++i) {
s1[i]=i-1,s2[i]=n/i-1;
f1[i]=power_sum(i)-1,f2[i]=power_sum(n/i)-1;
}
for(uint p=2;p<=sqn;++p) {
if(s1[p]==s1[p-1]) continue;
uint tmps=s1[p-1],tmpf=f1[p-1],pk=ksm(p,K);
for(uint i=1;i<=sqn/p;++i) {
s2[i]-=s2[i*p]-tmps;
f2[i]-=pk*(f2[i*p]-tmpf);
g2[i]+=f2[i*p]-tmpf;
}
for(uint i=sqn/p+1;i<=n/p/p&&i<=sqn;++i) {
s2[i]-=s1[n/i/p]-tmps;
f2[i]-=pk*(f1[n/i/p]-tmpf);
g2[i]+=f1[n/i/p]-tmpf;
}
for(uint i=sqn;i>=p*p;--i) {
s1[i]-=s1[i/p]-tmps;
f1[i]-=pk*(f1[i/p]-tmpf);
g1[i]+=f1[i/p]-tmpf;
}
}
}
uint du_sieve(uint x) {
if(x<=sqn?vis1[x]:vis2[n/x]) return x<=sqn?mp1[x]:mp2[n/x];
uint re=(x<=sqn?g1[x]+s1[x]:g2[n/x]+s2[n/x]);
for(uint i=2,j;i<=x;i=j+1) j=x/(x/i),re-=(j-i+1)*du_sieve(x/i);
(x<=sqn?vis1[x]:vis2[n/x])=1,(x<=sqn?mp1[x]:mp2[n/x])=re;
return re;
}
int main()
{
scanf("%u%d",&n,&K);
prework(),min25();
for(uint i=1,j,las=0,now;i<=n;i=j+1,las=now)
j=n/(n/i),now=du_sieve(j),ans+=(n/i)*(n/i)*(now-las);
printf("%u\n",ans);
return 0;
}