题目描述
记(sgcd(i,j))为(i,j)的次大公约数。
给你(n)，求
[sum_{i=1}^nsum_{j=1}^n{sgcd(i,j)}^k
]
对(2^{32})取模。
(nleq {10}^9,kleq 50)
题解
记(f(n))为(n)的次大因数
显然(sgcd(i,j)=f(gcd(i,j)))
先推一波式子。
[egin{align}
&sum_{i=1}^nsum_{j=1}^n{sgcd(i,j)}^k\
=&sum_{i=1}^n{f(i)}^k(2sum_{j=1}^frac{n}{i}varphi(j)-1)
end{align}
]
后面那个(varphi(i))的前缀和可以杜教筛/min_25筛+数论分块解决，所以只用关心前面这部分。
观察min_25筛求质数的(k)次方的前缀和的过程，发现在求(f_{n,j}=sum_{i=2}^n[i ext{是质数或}i ext{的每个质因子都}>p_j]i^k)的时候，减掉的那部分就是（最小质因子(=p_j)的数除以(p_j)后的值）的(k)次方和。所以直接把这些东西加起来就是我们要求的答案了。
还要加上质数的答案，即区间质数个数。
自然数幂求和要用第二类斯特林数那个做法。
[egin{align}
S_m(n)&=sum_{i=1}^ni^m\
&=sum_{i=1}^nsum_{j=1}^megin{Bmatrix}m\jend{Bmatrix}i^underline{j}\
&=sum_{j=1}^megin{Bmatrix}m\jend{Bmatrix}sum_{i=1}^ni^underline{j}\
&=sum_{j=1}^megin{Bmatrix}m\jend{Bmatrix}frac{{(n+1)}^underline{j+1}}{j+1}
end{align}
]
时间复杂度：(O(sqrt{n}k^2+frac{n^frac{3}{4}}{log n}))或(O(sqrt{n}k^2+frac{n^frac{3}{4}}{log n}+n^frac{2}{3}))
代码
1
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<utility>
#include<cmath>
#include<functional>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
void sort(int &a,int &b)
{
if(a>b)
swap(a,b);
}
void open(const char *s)
{
#ifndef ONLINE_JUDGE
char str[100];
sprintf(str,"%s.in",s);
freopen(str,"r",stdin);
sprintf(str,"%s.out",s);
freopen(str,"w",stdout);
#endif
}
int rd()
{
int s=0,c,b=0;
while(((c=getchar())<'0'||c>'9')&&c!='-');
if(c=='-')
{
c=getchar();
b=1;
}
do
{
s=s*10+c-'0';
}
while((c=getchar())>='0'&&c<='9');
return b?-s:s;
}
void put(int x)
{
if(!x)
{
putchar('0');
return;
}
static int c[20];
int t=0;
while(x)
{
c[++t]=x%10;
x/=10;
}
while(t)
putchar(c[t--]+'0');
}
int upmin(int &a,int b)
{
if(b<a)
{
a=b;
return 1;
}
return 0;
}
int upmax(int &a,int b)
{
if(b>a)
{
a=b;
return 1;
}
return 0;
}
ll n;
int k;
int fp(int a,int b)
{
int s=1;
for(;b;b>>=1,a*=a)
if(b&1)
s*=a;
return s;
}
int S[60][60];
int s[60];
int calc(ll x)
{
s[0]=x;
for(int i=1;i<=k;i++)
{
s[i]=1;
ll v=(x+1)/(i+1)*(i+1);
s[i]*=(x+1)/(i+1);
for(ll j=x+1;j>=x-i+1;j--)
s[i]*=(j==v?1:j);
for(int j=0;j<i;j++)
s[i]-=S[i][j]*s[j];
}
return s[k];
}
void init()
{
S[0][0]=1;
for(int i=1;i<=k;i++)
for(int j=1;j<=i;j++)
S[i][j]=S[i-1][j-1]-(i-1)*S[i-1][j];
}
namespace gao1
{
const int m=100000;
const int M=100010;
int b[M],pri[M];
int cnt;
int f1[M],f2[M];
int g1[M],g2[M];
int ans1[M],ans2[M];
void init()
{
for(int i=2;i<=m;i++)
{
if(!b[i])
pri[++cnt]=i;
for(int j=1;j<=cnt&&i*pri[j]<=m;j++)
{
b[i*pri[j]]=1;
if(i%pri[j]==0)
break;
}
}
}
void gao()
{
init();
for(int i=1;i<=m;i++)
{
f1[i]=calc(i)-1;
g1[i]=i-1;
}
for(int i=1;n/i>m;i++)
{
f2[i]=calc(n/i)-1;
g2[i]=n/i-1;
}
for(int i=1;i<=cnt;i++)
{
int j;
int v=fp(pri[i],k);
for(j=1;n/pri[i]/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
{
ans2[j]+=f2[pri[i]*j]-f1[pri[i]-1];
f2[j]-=v*(f2[pri[i]*j]-f1[pri[i]-1]);
g2[j]-=g2[pri[i]*j]-g1[pri[i]-1];
}
for(;n/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
{
ans2[j]+=f1[n/pri[i]/j]-f1[pri[i]-1];
f2[j]-=v*(f1[n/pri[i]/j]-f1[pri[i]-1]);
g2[j]-=g1[n/pri[i]/j]-g1[pri[i]-1];
}
for(j=m;j>=2&&j>=(ll)pri[i]*pri[i];j--)
{
ans1[j]+=f1[j/pri[i]]-f1[pri[i]-1];
f1[j]-=v*(f1[j/pri[i]]-f1[pri[i]-1]);
g1[j]-=g1[j/pri[i]]-g1[pri[i]-1];
}
}
for(int i=1;i<=m;i++)
ans1[i]+=g1[i];
for(int i=1;n/i>m;i++)
ans2[i]+=g2[i];
}
int query(ll x)
{
return x<=m?ans1[x]:ans2[n/x];
}
}
namespace gao2
{
const int m=1000000;
const int M=1000010;
int pri[M],b[M],cnt;
int phi[M],s[M];
int b2[M],s2[M];
void init()
{
phi[1]=1;
for(int i=2;i<=m;i++)
{
if(!b[i])
{
pri[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt&&i*pri[j]<=m;j++)
{
b[i*pri[j]]=1;
if(i%pri[j]==0)
{
phi[i*pri[j]]=phi[i]*pri[j];
break;
}
phi[i*pri[j]]=phi[i]*phi[pri[j]];
}
}
for(int i=1;i<=m;i++)
s[i]=s[i-1]+phi[i];
}
void gao()
{
init();
}
int query(ll x)
{
if(x<=m)
return s[x];
if(b2[n/x])
return s2[n/x];
b2[n/x]=1;
int &res=s2[n/x];
res=x*(x+1)/2;
for(ll i=2,j;i<=x;i=j+1)
{
j=x/(x/i);
res-=query(x/i)*(j-i+1);
}
return res;
}
}
int main()
{
open("51nod1847");
scanf("%lld%d",&n,&k);
init();
gao1::gao();
gao2::gao();
int ans=0;
for(ll i=2,j;i<=n;i=j+1)
{
j=n/(n/i);
ans+=(gao1::query(j)-gao1::query(i-1))*(2*gao2::query(n/i)-1);
}
printf("%u
",(unsigned)ans);
return 0;
}
2
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<utility>
#include<cmath>
#include<functional>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
void sort(int &a,int &b)
{
if(a>b)
swap(a,b);
}
void open(const char *s)
{
#ifndef ONLINE_JUDGE
char str[100];
sprintf(str,"%s.in",s);
freopen(str,"r",stdin);
sprintf(str,"%s.out",s);
freopen(str,"w",stdout);
#endif
}
int rd()
{
int s=0,c,b=0;
while(((c=getchar())<'0'||c>'9')&&c!='-');
if(c=='-')
{
c=getchar();
b=1;
}
do
{
s=s*10+c-'0';
}
while((c=getchar())>='0'&&c<='9');
return b?-s:s;
}
void put(int x)
{
if(!x)
{
putchar('0');
return;
}
static int c[20];
int t=0;
while(x)
{
c[++t]=x%10;
x/=10;
}
while(t)
putchar(c[t--]+'0');
}
int upmin(int &a,int b)
{
if(b<a)
{
a=b;
return 1;
}
return 0;
}
int upmax(int &a,int b)
{
if(b>a)
{
a=b;
return 1;
}
return 0;
}
ll n;
int k;
int fp(int a,int b)
{
int s=1;
for(;b;b>>=1,a*=a)
if(b&1)
s*=a;
return s;
}
int S[60][60];
int s[60];
int calc(ll x)
{
s[0]=x;
for(int i=1;i<=k;i++)
{
s[i]=1;
ll v=(x+1)/(i+1)*(i+1);
s[i]*=(x+1)/(i+1);
for(ll j=x+1;j>=x-i+1;j--)
s[i]*=(j==v?1:j);
for(int j=0;j<i;j++)
s[i]-=S[i][j]*s[j];
}
return s[k];
}
void init()
{
S[0][0]=1;
for(int i=1;i<=k;i++)
for(int j=1;j<=i;j++)
S[i][j]=(S[i-1][j-1]-(i-1)*S[i-1][j]);
}
namespace gao1
{
const int m=100000;
const int M=100010;
int b[M],pri[M];
int cnt;
int f1[M],f2[M];
int ans1[M],ans2[M];
// vector<int> s1[M],s2[M];
// vector<short> b1[M],b2[M];
void init()
{
for(int i=2;i<=m;i++)
{
if(!b[i])
pri[++cnt]=i;
for(int j=1;j<=cnt&&i*pri[j]<=m;j++)
{
b[i*pri[j]]=1;
if(i%pri[j]==0)
break;
}
}
}
void gao()
{
init();
for(int i=1;i<=m;i++)
f1[i]=calc(i)-1;
for(int i=1;n/i>m;i++)
f2[i]=calc(n/i)-1;
for(int i=1;i<=cnt;i++)
{
int j;
int v=fp(pri[i],k);
for(j=1;n/pri[i]/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
{
int x=f2[pri[i]*j]-f1[pri[i]-1];
ans2[j]+=x;
f2[j]-=v*x;
}
for(;n/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
{
int x=f1[n/pri[i]/j]-f1[pri[i]-1];;
ans2[j]+=x;
f2[j]-=v*x;
}
for(j=m;j>=(ll)pri[i]*pri[i];j--)
{
int x=f1[j/pri[i]]-f1[pri[i]-1];
ans1[j]+=x;
f1[j]-=v*x;
}
}
}
int query(ll x)
{
return x<=m?ans1[x]:ans2[n/x];
}
}
namespace gao2
{
const int m=100000;
const int M=100010;
int pri[M],b[M],cnt;
int f1[M],f2[M],g1[M],g2[M];
int s1[M],s2[M];
int b1[M],b2[M];
void init()
{
for(int i=2;i<=m;i++)
{
if(!b[i])
pri[++cnt]=i;
for(int j=1;j<=cnt&&i*pri[j]<=m;j++)
{
b[i*pri[j]]=1;
if(i%pri[j]==0)
break;
}
}
pri[cnt+1]=m+1;
pri[0]=1;
}
void gao()
{
init();
for(int i=1;i<=m;i++)
{
f1[i]=(ll)(i+2)*(i-1)/2;
g1[i]=i-1;
}
for(int i=1;n/i>m;i++)
{
f2[i]=(n/i&1?(n/i-1)/2*(n/i+2):(n/i+2)/2*(n/i-1));
g2[i]=n/i-1;
}
for(int i=1;i<=cnt;i++)
{
int j;
int x=f1[pri[i]-1];
int y=g1[pri[i]-1];
for(j=1;n/pri[i]/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
{
f2[j]-=pri[i]*(f2[pri[i]*j]-x);
g2[j]-=g2[pri[i]*j]-y;
}
for(;n/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
{
f2[j]-=pri[i]*(f1[n/pri[i]/j]-x);
g2[j]-=g1[n/pri[i]/j]-y;
}
for(j=m;j>=2&&j>=(ll)pri[i]*pri[i];j--)
{
f1[j]-=pri[i]*(f1[j/pri[i]]-x);
g1[j]-=g1[j/pri[i]]-y;
}
}
for(int i=1;i<=m;i++)
gao1::ans1[i]+=g1[i];
for(int i=1;n/i>m;i++)
gao1::ans2[i]+=g2[i];
for(int i=1;i<=m;i++)
f1[i]-=g1[i];
for(int i=1;n/i>m;i++)
f2[i]-=g2[i];
for(int i=1;i<=m;i++)
s1[i]=f1[i];
for(int i=1;n/i>m;i++)
s2[i]=f2[i];
for(int i=cnt;i>=1;i--)
{
int x=pri[i];
for(int j=1;n/j>m&&n/j>=(ll)x*x;j++)
{
s2[j]-=x-1;
ll x2=n/j/x;
int s=x-1;
for(;x2;x2/=x,s*=x)
s2[j]+=((x2<=m?s1[x2]-f1[min((int)x2,pri[i])]:s2[n/x2]-f1[pri[i]])+1)*s;
}
for(int j=m;j>=(ll)x*x;j--)
{
s1[j]-=x-1;
int x2=j/x;
int s=x-1;
for(;x2;x2/=x,s*=x)
s1[j]+=(s1[x2]-f1[min(x2,pri[i])]+1)*s;
}
}
}
int query(ll x)
{
return (x<=m?s1[x]:s2[n/x])+1;
}
}
int main()
{
open("51nod1847");
scanf("%lld%d",&n,&k);
init();
gao1::gao();
gao2::gao();
int ans=0;
for(ll i=2,j;i<=n;i=j+1)
{
j=n/(n/i);
ans+=(gao1::query(j)-gao1::query(i-1))*(2*gao2::query(n/i)-1);
}
printf("%u
",ans);
return 0;
}