题目大意
求:
∑
i
=
1
n
∑
j
=
1
n
gcd
(
i
,
j
)
k
l
c
m
(
i
,
j
)
[
gcd
i
s
p
r
i
m
e
]
\sum_{i=1}^n \sum_{j=1}^n \gcd(i,j)^k lcm(i,j) ~[\gcd~is~prime]
i=1∑nj=1∑ngcd(i,j)klcm(i,j) [gcd is prime]
1
≤
n
≤
1
0
10
,
1
≤
k
≤
100
1 \leq n \leq 10^{10},~1 \leq k \leq 100
1≤n≤1010, 1≤k≤100
多测,
T
≤
5
T \leq 5
T≤5,时限 10s
\\
\\
\\
题解
∑
i
=
1
n
∑
j
=
1
n
gcd
(
i
,
j
)
k
l
c
m
(
i
,
j
)
[
gcd
i
s
p
r
i
m
e
]
=
∑
d
=
1
n
d
k
+
1
[
d
i
s
p
r
i
m
e
]
∑
i
′
=
1
⌊
n
d
⌋
∑
j
′
=
1
⌊
n
d
⌋
i
′
j
′
[
gcd
(
i
,
j
)
=
1
]
\sum_{i=1}^n \sum_{j=1}^n \gcd(i,j)^k lcm(i,j) ~[\gcd~is~prime] \\ =\sum_{d=1}^n d^{k+1}[d~is~prime] \sum_{i'=1}^{\lfloor \frac nd \rfloor} \sum_{j'=1}^{\lfloor \frac nd \rfloor} i'j'[\gcd(i,j)=1] \\
i=1∑nj=1∑ngcd(i,j)klcm(i,j) [gcd is prime]=d=1∑ndk+1[d is prime]i′=1∑⌊dn⌋j′=1∑⌊dn⌋i′j′[gcd(i,j)=1]
后面那玩意儿是个经典转换:
∑
i
=
1
m
∑
j
=
1
m
i
j
[
gcd
(
i
,
j
)
=
1
]
=
(
2
×
∑
i
=
2
m
i
∑
j
=
1
i
−
1
j
[
gcd
(
i
,
j
)
=
1
]
)
+
(
∑
i
=
1
m
i
2
[
gcd
(
i
,
i
)
=
1
]
)
=
(
2
×
∑
i
=
2
m
i
i
ϕ
(
i
)
2
)
+
1
=
∑
i
=
1
m
i
2
ϕ
(
i
)
\begin{aligned} & \sum_{i=1}^m \sum_{j=1}^m ij[\gcd(i,j)=1] \\ =& \bigg( 2 \times \sum_{i=2}^m i \sum_{j=1}^{i-1} j[\gcd(i,j)=1] \bigg) + \bigg( \sum_{i=1}^m i^2[\gcd(i,i)=1] \bigg)\\ =& \bigg( 2 \times \sum_{i=2}^m i \frac{i \phi(i)}2 \bigg) +1 \\ =& \sum_{i=1}^m i^2 \phi(i) \end{aligned}
===i=1∑mj=1∑mij[gcd(i,j)=1](2×i=2∑mij=1∑i−1j[gcd(i,j)=1])+(i=1∑mi2[gcd(i,i)=1])(2×i=2∑mi2iϕ(i))+1i=1∑mi2ϕ(i)
因此原式等于
∑
d
=
1
n
d
k
+
1
[
d
i
s
p
r
i
m
e
]
∑
i
=
1
⌊
n
d
⌋
i
2
ϕ
(
i
)
\sum_{d=1}^n d^{k+1} [d~is~prime] \sum_{i=1}^{\lfloor \frac nd \rfloor} i^2 \phi(i)
d=1∑ndk+1[d is prime]i=1∑⌊dn⌋i2ϕ(i)
于是对
⌊
n
d
⌋
\lfloor \frac nd \rfloor
⌊dn⌋ 分块,前面求质数幂和用洲阁筛或 min25,后面求
∑
i
2
ϕ
(
i
)
\sum i^2\phi(i)
∑i2ϕ(i) 用杜教筛或 min25。
扯淡
symbol 推法(可参考 WC2016 第二课堂课件)居然不是万能的。。。
∑
i
=
1
n
∑
j
=
1
n
gcd
(
i
,
j
)
k
l
c
m
(
i
,
j
)
[
gcd
i
s
p
r
i
m
e
]
=
∑
d
=
1
n
d
k
+
1
[
d
i
s
p
r
i
m
e
]
∑
i
′
=
1
⌊
n
d
⌋
∑
j
′
=
1
⌊
n
d
⌋
i
′
j
′
[
gcd
(
i
,
j
)
=
1
]
=
∑
d
=
1
n
d
k
+
1
[
d
i
s
p
r
i
m
e
]
∑
D
=
1
⌊
n
d
⌋
μ
(
D
)
D
2
(
1
+
⌊
n
d
D
⌋
)
2
(
⌊
n
d
D
⌋
)
2
4
=
∑
T
=
1
n
(
1
+
⌊
n
T
⌋
)
2
(
⌊
n
T
⌋
)
2
4
T
∑
d
∣
T
,
d
i
s
p
r
i
m
e
d
k
−
1
μ
(
T
d
)
\begin{aligned} &\sum_{i=1}^n \sum_{j=1}^n \gcd(i,j)^k lcm(i,j) ~[\gcd~is~prime] \\ =&\sum_{d=1}^n d^{k+1}[d~is~prime] \sum_{i'=1}^{\lfloor \frac nd \rfloor} \sum_{j'=1}^{\lfloor \frac nd \rfloor} i'j'[\gcd(i,j)=1] \\ =&\sum_{d=1}^n d^{k+1}[d~is~prime] \sum_{D=1}^{\lfloor \frac nd \rfloor} \mu(D)D^2\frac{(1+\lfloor \frac n{dD} \rfloor)^2 (\lfloor \frac n{dD} \rfloor)^2}4 \\ =&\sum_{T=1}^n \frac{(1+\lfloor \frac nT \rfloor)^2 (\lfloor \frac nT \rfloor)^2}4 T \sum_{d|T,~d~is~prime} d^{k-1}\mu(\frac Td) \end{aligned}
===i=1∑nj=1∑ngcd(i,j)klcm(i,j) [gcd is prime]d=1∑ndk+1[d is prime]i′=1∑⌊dn⌋j′=1∑⌊dn⌋i′j′[gcd(i,j)=1]d=1∑ndk+1[d is prime]D=1∑⌊dn⌋μ(D)D24(1+⌊dDn⌋)2(⌊dDn⌋)2T=1∑n4(1+⌊Tn⌋)2(⌊Tn⌋)2Td∣T, d is prime∑dk−1μ(dT)
以前是觉得,这种变量前移的推法能通杀一切,即便跟标准做法不一样,但最后会推出来相同的东西,或者我这样推也是能做的(甚至更快)。做这题之前,这个想法都是对的。
直到遇见了这题。。。后面那坨居然不能预处理又不能筛??
代码
#include<bits/stdc++.h>
#define fo(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
typedef long long LL;
const int maxsqrtn=2e6+5, maxk=105;
const LL mo=1e9+7;
LL n;
int sqrtn,k;
LL mi(LL x,LL y)
{
LL re=1;
for(; y; y>>=1, x=x*x%mo) if (y&1) re=re*x%mo;
return re;
}
LL phi[maxsqrtn],Sp[maxsqrtn],inv[maxsqrtn];
int p0,p[maxsqrtn],Np[maxsqrtn];
bool bz[maxsqrtn];
void Pre(int n)
{
p0=0;
memset(bz,0,sizeof(bz));
memset(Np,0,sizeof(Np));
phi[1]=1;
fo(i,2,n)
{
if (!bz[i])
{
p[++p0]=i;
phi[i]=i-1;
Np[i]=1;
Sp[p0]=(Sp[p0-1]+mi(i,k))%mo;
}
fo(j,1,p0)
{
if ((LL)i*p[j]>n) break;
bz[i*p[j]]=1;
if (i%p[j]==0)
{
phi[i*p[j]]=phi[i]*p[j];
break;
} else phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
inv[1]=1;
fo(i,2,n)
{
phi[i]=(phi[i-1]+phi[i]*i%mo*i%mo)%mo;
Np[i]+=Np[i-1];
inv[i]=(-(mo/i)*inv[mo%i])%mo;
}
}
LL x[maxk],y[maxk],w[maxk],invw[maxk];
void lagrange_add(int n)
{
fo(i,1,n-1) (w[i]*=(x[i]-x[n]))%=mo;
w[n]=1;
fo(i,1,n-1) (w[n]*=(x[n]-x[i]))%=mo;
}
LL lagrange_get(int n,LL nx)
{
if (nx<=n) return y[nx];
nx%=mo;
LL sum=0, l=1;
fo(i,1,n)
{
(l*=(nx-x[i]))%=mo;
(sum+=y[i]*invw[i]%mo*(nx-x[i]<=maxsqrtn-5 ?inv[nx-x[i]] :mi(nx-x[i],mo-2)))%=mo;
}
return l*sum%mo;
}
LL mw[maxsqrtn],g[maxsqrtn],id1[maxsqrtn],id2[maxsqrtn];
int w0;
LL min25_g(LL n)
{
w0=0;
for(LL i=1, j; i<=n; i=j+1)
{
j=n/(n/i);
mw[++w0]=n/i;
if (mw[w0]<=sqrtn) id1[mw[w0]]=w0; else id2[j]=w0;
g[w0]=lagrange_get(k+2,mw[w0])-1;
}
fo(j,1,Np[sqrtn])
for(int i=1; i<=w0 && (LL)p[j]*p[j]<=mw[i]; i++)
{
int id=(mw[i]/p[j]<=sqrtn) ?id1[mw[i]/p[j]] :id2[n/(mw[i]/p[j])];
(g[i]-=(Sp[j]-Sp[j-1])*(g[id]-Sp[j-1]))%=mo;
}
}
unordered_map<LL,LL> f;
LL inv2,inv6;
LL sum(LL n,int id)
{
n%=mo;
if (id==2) return n*(n+1)%mo*inv6%mo*(2*n+1)%mo;
else if (id==3)
{
LL t=n*(n+1)%mo*inv2%mo;
return t*t%mo;
}
}
LL Sphi(LL x,int id)
{
if (x<=maxsqrtn-5) return phi[x];
if (f.count(x)) return f[x];
LL re=sum(x,id+1);
for(LL i=2, j; i<=x; i=j+1)
{
j=x/(x/i);
(re-=(sum(j,id)-sum(i-1,id)+mo)%mo*Sphi(x/i,id)%mo)%=mo;
}
(re+=mo)%=mo;
return f[x]=re;
}
int T;
int main()
{
inv2=mi(2,mo-2);
inv6=mi(6,mo-2);
scanf("%d",&T);
while (T--)
{
scanf("%lld %d",&n,&k); k++;
sqrtn=sqrt(n);
Pre(maxsqrtn-5);
fo(i,1,k+2)
{
x[i]=i, y[i]=(y[i-1]+mi(i,k))%mo;
lagrange_add(i);
}
fo(i,1,k+2) invw[i]=mi(w[i],mo-2);
min25_g(n);
LL ans=0;
for(LL i=1, j; i<=n; i=j+1)
{
j=n/(n/i);
int idi=((i-1)<=sqrtn) ?id1[(i-1)] :id2[n/(i-1)] ;
int idj=(j<=sqrtn) ?id1[j] :id2[n/j] ;
(ans+=(g[idj]-g[idi]+mo)%mo*Sphi(n/i,2))%=mo;
}
printf("%lld\n",(ans+mo)%mo);
}
}