题意简述
求
∑
i
=
1
n
∑
j
=
1
n
i
j
gcd
(
i
,
j
)
\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\dfrac{ij}{\gcd (i,j)}
i=1∑nj=1∑ngcd(i,j)ij。
对 998244353 998244353 998244353 取模。 n < = 1 0 7 , m < = 1 0 14 n<=10^7,m<=10^{14} n<=107,m<=1014。
注:空间限制只有64MB,只够开
1
0
7
10^7
107 的int
型数组。
思路
基础操作
设
n
≤
m
n\le m
n≤m,因为
n
,
m
n,m
n,m 对称。
设
g
=
gcd
(
i
,
j
)
g=\gcd(i,j)
g=gcd(i,j) (简写)
∑
i
=
1
n
∑
j
=
1
m
i
j
g
2
\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} \dfrac{ij}{g^2}
i=1∑nj=1∑mg2ij
∑
d
=
1
n
∑
i
=
1
n
/
d
∑
j
=
1
m
/
d
[
g
=
1
]
i
j
\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{m/d}[g=1]ij
d=1∑ni=1∑n/dj=1∑m/d[g=1]ij (枚举
gcd
\gcd
gcd)
∑
q
=
1
n
μ
(
q
)
∑
d
=
1
n
∑
i
=
1
n
/
d
∑
j
=
1
m
/
d
[
q
∣
g
]
i
j
\sum\limits_{q=1}^{n}\mu(q)\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{m/d}[q|g]ij
q=1∑nμ(q)d=1∑ni=1∑n/dj=1∑m/d[q∣g]ij
∑
q
=
1
n
μ
(
q
)
∑
d
=
1
n
∑
i
=
1
n
/
d
q
∑
j
=
1
m
/
d
q
i
q
×
j
q
\sum\limits_{q=1}^{n}\mu(q)\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{n/dq}\sum\limits_{j=1}^{m/dq}iq\times jq
q=1∑nμ(q)d=1∑ni=1∑n/dqj=1∑m/dqiq×jq
∑
q
=
1
n
μ
(
q
)
q
2
∑
d
=
1
n
∑
i
=
1
n
/
d
q
∑
j
=
1
m
/
d
q
i
j
\sum\limits_{q=1}^{n}\mu(q)q^2\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{n/dq}\sum\limits_{j=1}^{m/dq}ij
q=1∑nμ(q)q2d=1∑ni=1∑n/dqj=1∑m/dqij
设
s
(
x
)
=
∑
i
=
1
x
i
s(x)=\sum\limits_{i=1}^x i
s(x)=i=1∑xi
原式化为:
∑
q
=
1
n
μ
(
q
)
q
2
∑
d
=
1
n
s
(
n
/
d
q
)
s
(
m
/
d
q
)
\sum\limits_{q=1}^n\mu(q)q^2\sum\limits_{d=1}^{n}s(n/dq)s(m/dq)
q=1∑nμ(q)q2d=1∑ns(n/dq)s(m/dq)
我们发现,
d
q
>
n
dq>n
dq>n 时,
s
(
n
/
d
q
)
=
0
s(n/dq)=0
s(n/dq)=0,也就不会有贡献了。那么我们珂以缩小
d
d
d 的范围:
∑
q
=
1
n
μ
(
q
)
q
2
∑
d
=
1
n
/
q
s
(
n
/
d
q
)
s
(
m
/
d
q
)
\sum\limits_{q=1}^{n}\mu(q)q^2\sum\limits_{d=1}^{n/q}s(n/dq)s(m/dq)
q=1∑nμ(q)q2d=1∑n/qs(n/dq)s(m/dq)
整除分块两次,是 O ( n ) O(n) O(n) 的(?)。
但是似乎这个复杂度还和 m m m 有点关系,我不是很清楚…总之,它TLE了,只有 80 80 80 分。(亲测)
那么我们考虑枚举 q × d q\times d q×d,设为 T T T。它的贡献就是 s ( n / T ) s ( m / T ) s(n/T)s(m/T) s(n/T)s(m/T) 。而它被算到的此时,显然, T = d q T=dq T=dq 是 q q q 的倍数,所以 q q q 是 T = d q T=dq T=dq 的因数。那 d q = T dq=T dq=T 被算到的次数就是 ∑ q ∣ T μ ( q ) q 2 \sum\limits_{q|T} \mu(q)q^2 q∣T∑μ(q)q2。将这个值设为 f ( T ) f(T) f(T)。
式子化一下,变成: ∑ T = 1 n f ( T ) × s ( n / T ) × s ( m / T ) \sum\limits_{T=1}^{n}f(T)\times s(n/T)\times s(m/T) T=1∑nf(T)×s(n/T)×s(m/T)。
那么问题就在于我们怎么筛这个 f ( T ) f(T) f(T) ,还要严格线性,因为 n < = 1 0 7 n<=10^7 n<=107。
如何筛f函数
显然, f ( 1 ) = 1 f(1)=1 f(1)=1,而且当 p p p 和 q q q 互质时, f ( p × q ) = f ( p ) × f ( q ) f(p\times q)=f(p)\times f(q) f(p×q)=f(p)×f(q)。这说明它也许能线性筛。
考虑质数情况。当 p p p 为质数时, f ( p ) = 1 − p 2 f(p)=1-p^2 f(p)=1−p2。
接着, f ( p k ) = ∑ d ∣ p k μ ( d ) d 2 f(p^k)=\sum\limits_{d|p^k}\mu(d)d^2 f(pk)=d∣pk∑μ(d)d2。显然, p k p^k pk 的因数只有 p 0 , p 1 ⋯ p k p^0,p^1\cdots p^k p0,p1⋯pk。而对于 p j p^j pj ( 2 ≤ j ≤ k 2\le j \le k 2≤j≤k),显然, μ ( p j ) = 0 \mu(p^j)=0 μ(pj)=0,也就不会有贡献了。因此,有贡献的只有 p 0 p^0 p0 和 p 1 p^1 p1 。
那也就是说, f ( p k ) = f ( p ) f(p^k)=f(p) f(pk)=f(p)!!!(换句话说,同一个质因子乘多少遍都不会改变 f f f 的值)
那么我们在线性筛的时候有这样一步:枚举
i
i
i,找一个质数
u
u
u。如果
u
u
u 不是
i
i
i 的因数,那么显然
i
i
i 和
u
u
u 互质,
f
(
i
×
u
)
=
f
(
i
)
×
f
(
u
)
f(i\times u)=f(i)\times f(u)
f(i×u)=f(i)×f(u)。否则
u
u
u 是
i
i
i 的质因子,然后我们要计算
f
(
i
×
u
)
f(i\times u)
f(i×u) 的值,并 break
。
由上面那个性质得, u u u 是 i i i 的质因子,而我们又把它乘了一遍,并不会改变 f f f 的值。也就是 f ( i × u ) = f ( i ) × f ( u ) f(i\times u)=f(i)\times f(u) f(i×u)=f(i)×f(u)。
知道了这些,就能线性筛这个 f f f 函数了!!!
线性筛出来了 f f f 之后,整除分块都不用,直接 O ( n ) O(n) O(n) 暴力就能求出原式的值了。而且是很稳的 O ( n ) O(n) O(n) 哦!
极 限 卡 常
上面说了,空间限制只有 64MB ,只能开的下一个
1
0
7
10^7
107 个int
数组,两个就不行了。
说的具体点,能开的下
16
777
216
16\ 777\ 216
16 777 216个 int
。
但是我们常见的线性筛,需要一个int primes[1e7]
保存质数,还有一个bool notp[1e7]
标记是否不是素数,然后才是我们的f
数组。这可怎么办呢?
首先,
[
1
,
n
]
[1,n]
[1,n] 中的质数个数和
n
n
n 绝对不是同一个级别。具体的说,大概是
n
ln
(
n
)
\dfrac{n}{\ln(n)}
ln(n)n 级别的。
[
1
,
1
e
7
]
[1,1e7]
[1,1e7] 中的质数,经过测试,还不足
1
0
6
10^6
106 个。那么我们的 primes
数组就珂以开的小一点了。
接着,bool notp[1e7]
珂以用一个 bitset
存储,空间直接少掉
32
32
32 倍,你说爽不爽。
然后我们就能卡进空间常数了!nice!
还有,记得不要全开long long
啊!
代码
#include <bits/stdc++.h>
using namespace std;
namespace Flandre_Scarlet
{
#define N 10000007
#define mod 998244353ll
#define i6 166374059ll
#define i2 499122177ll
// 模数,2的逆元,6的逆元
// 只是打个板子,并不是所有的数都有用上
#define ll long long
#define F(i,l,r) for(int i=l;i<=r;++i)
#define D(i,r,l) for(int i=r;i>=l;--i)
#define Fs(i,l,r,c) for(int i=l;i<=r;c)
#define Ds(i,r,l,c) for(int i=r;i>=l;c)
#define MEM(x,a) memset(x,a,sizeof(x))
#define FK(x) MEM(x,0)
#define Tra(i,u) for(int i=G.Start(u),__v=G.To(i);~i;i=G.Next(i),__v=G.To(i))
#define p_b push_back
#define sz(a) ((int)a.size())
#define iter(a,p) (a.begin()+p)
//并不想用快读,数据小
ll n,m;
void Input()
{
scanf("%lld%lld",&n,&m);
}
int f[N];
int primes[N/10]; //N/10足够了
bitset<N> notp; //STL万岁
void Init()
{
int &cnt=primes[0];
int n=1e7;
f[1]=notp[1]=1;
F(i,2,n)
{
if (!notp[i])
{
f[i]=(mod+1-1ll*i*i%mod)%mod;
// 上面说了,i为质数的时候,f[i]=1-i*i
// 别忘了取模啊
primes[++cnt]=i;
}
for(int j=1;j<=cnt and i*primes[j]<=n;++j)
{
int u=primes[j];
notp[i*u]=1;
if (i%u!=0)
{
f[i*u]=1ll*f[i]*f[u]%mod;
//利用积性性质,f[i*u]=f[i]*f[u]
}
else
{
f[i*u]=f[i];
//上面特意讲了,此时f[i*u]=f[i]
break;
}
}
}
}
ll s(ll x){x%=mod;return x*(x+1)%mod*i2%mod;}
//计算1+2+3...+x的和,也就是x*(x+1)/2
void Soviet()
{
if (n>m) swap(n,m);
//令n<=m
ll ans=0;
F(i,1,n)
{
ans+=1ll*f[i]*s(n/i)%mod*s(m/i)%mod;
//每次ans+=f[i]*s(n/i)*s(m/i)
//别忘了取模啊!
ans%=mod;
}
printf("%lld\n",ans);
}
#define Flan void
Flan IsMyWife()
{
Init();
Input();
Soviet();
}
}
int main()
{
Flandre_Scarlet::IsMyWife();
getchar();getchar();
return 0;
}