11.23A T4方（思维好题----杜教筛（总复习））
Description
给一个 n × m n\times m n×m的网格（一共 ( n + 1 ) × ( m + 1 ) (n+1)\times(m+1) (n+1)×(m+1)个格点），定义格点正方形的权值时它完全包含的格子数，求网格中所有正方形（可以是斜的）的权值和。
Solution
容易发现一个正方形的权值就是其面积减去两条斜边穿过的格子。
面积易求。穿过的格子怎么办？让人没有头绪。
考虑穿过的格子坐标满足的条件？那么线段应该与格子有两条边相交，限制比较复杂。
但在思考的时候我们可以发现 格子是否被穿过 似乎可以与 边的相交 联系起来。
于是有这样的思路：当线段穿过一条网格边时，它就进入了一个新格子，且它不会再次回到它。而经过一条横边则意味着线段在纵轴上升高了1，纵边同理。
不过上述方式有一个前提—经过的是网格边；那么经过格点时呢，也是进入了一个新格子，为了方便，我们在上一种情况中一起计算，再减去重复的就好。
如果设斜边在坐标轴上的投影分别为
a
,
b
a,b
a,b ，那么有：
c
n
t
(
a
,
b
)
=
a
+
b
−
gcd
(
a
,
b
)
\large cnt(a,b)=a+b-\gcd (a,b)
cnt(a,b)=a+b−gcd(a,b)
则可以通过枚举
a
,
b
a,b
a,b 得到答案：
a
n
s
(
n
,
m
)
=
∑
a
=
0
n
−
1
∑
b
=
1
min
(
n
,
m
)
−
a
(
a
2
+
b
2
−
2
c
n
t
(
a
,
b
)
)
×
[
n
−
(
a
+
b
)
+
1
]
[
m
−
(
a
+
b
)
+
1
]
\large ans(n,m)=\sum_{a=0}^{n-1}\sum_{b=1}^{\min(n,m)-a}(a^2+b^2-2cnt(a,b))\times [n-(a+b)+1][m-(a+b)+1]
ans(n,m)=a=0∑n−1b=1∑min(n,m)−a(a2+b2−2cnt(a,b))×[n−(a+b)+1][m−(a+b)+1]
然后就是推式子环节。
首先可以把含gcd的因式分开考虑（part1）：
设
L
=
a
+
b
L=a+b
L=a+b，记
N
=
min
(
n
,
m
)
N=\min(n,m)
N=min(n,m)
a
n
s
1
=
∑
L
=
1
N
(
n
−
L
+
1
)
(
m
−
L
+
1
)
∑
a
=
1
L
gcd
(
a
,
L
−
a
)
=
∑
L
=
1
N
(
n
−
L
+
1
)
(
m
−
L
+
1
)
∑
a
=
1
L
∑
d
∣
a
,
d
∣
L
ϕ
(
d
)
=
∑
d
=
1
N
ϕ
(
d
)
∑
L
′
=
1
⌊
N
d
⌋
(
n
−
L
′
d
+
1
)
(
m
−
L
′
d
+
1
)
L
′
=
∑
d
=
1
N
ϕ
(
d
)
F
0
(
⌊
N
d
⌋
)
+
∑
d
=
1
N
ϕ
(
d
)
d
×
F
1
(
⌊
N
d
⌋
)
+
∑
d
=
1
N
ϕ
(
d
)
d
2
×
F
2
(
⌊
N
d
⌋
)
\begin{aligned} \large ans1&\large =\sum_{L=1}^N(n-L+1)(m-L+1)\sum_{a=1}^{L}\gcd(a,L-a)\\ &\large =\sum_{L=1}^N(n-L+1)(m-L+1)\sum_{a=1}^{L}\sum_{d|a,d|L}\phi(d) \\ &\large =\sum_{d=1}^N\phi(d)\sum_{L'=1}^{\lfloor \frac{N}{d} \rfloor} (n-L'd+1)(m-L'd+1)L' \\ &\large =\sum_{d=1}^N\phi(d)F_0(\lfloor \frac{N}{d} \rfloor)+ \sum_{d=1}^N\phi(d)d\times F_1(\lfloor \frac{N}{d} \rfloor)+ \sum_{d=1}^N\phi(d)d^2\times F_2(\lfloor \frac{N}{d} \rfloor)\\ \end{aligned}
ans1=L=1∑N(n−L+1)(m−L+1)a=1∑Lgcd(a,L−a)=L=1∑N(n−L+1)(m−L+1)a=1∑Ld∣a,d∣L∑ϕ(d)=d=1∑Nϕ(d)L′=1∑⌊dN⌋(n−L′d+1)(m−L′d+1)L′=d=1∑Nϕ(d)F0(⌊dN⌋)+d=1∑Nϕ(d)d×F1(⌊dN⌋)+d=1∑Nϕ(d)d2×F2(⌊dN⌋)
此处直接上 数论分块+杜教筛。
part2：
记
C
1
=
−
(
n
+
1
)
−
(
m
+
1
)
,
C
2
=
(
n
+
1
)
(
m
+
1
)
,
G
c
(
n
)
=
∑
i
=
1
n
i
c
\large C_1=-(n+1)-(m+1),C_2=(n+1)(m+1),G_c(n)=\sum_{i=1}^{n}i^c
C1=−(n+1)−(m+1),C2=(n+1)(m+1),Gc(n)=∑i=1nic
a
n
s
2
=
∑
L
=
1
N
(
L
2
+
C
1
L
+
C
2
)
∑
a
=
1
L
(
L
2
−
2
L
+
2
a
2
−
2
a
L
)
=
∑
L
=
1
N
(
L
2
+
C
1
L
+
C
2
)
(
L
3
−
2
L
2
+
2
G
2
(
L
)
−
2
L
G
1
(
L
)
)
\begin{aligned} \large ans2&\large =\sum_{L=1}^N(L^2+C_1L+C_2)\sum_{a=1}^{L}(L^2-2L+2a^2-2aL) \\ &\large =\sum_{L=1}^N(L^2+C_1L+C_2)(L^3-2L^2+2G_2(L)-2LG_1(L)) \\ \end{aligned}
ans2=L=1∑N(L2+C1L+C2)a=1∑L(L2−2L+2a2−2aL)=L=1∑N(L2+C1L+C2)(L3−2L2+2G2(L)−2LG1(L))
part1,part2都用到了n以内自然数的幂的和。补补公式推导：
( n + 1 ) c + 1 − n c + 1 = ∑ i = 0 c a i n i , ( 1 ) \large (n+1)^{c+1}-n^{c+1}=\sum_{i=0}^{c} a_in^i,(1) (n+1)c+1−nc+1=∑i=0caini,(1)
可以推出
( n + 1 ) c + 1 − 1 = ∑ i = 1 c a i G i ( n ) + n , ( 2 ) \large (n+1)^{c+1}-1=\sum_{i=1}^{c}a_iG_i(n)+n,(2) (n+1)c+1−1=∑i=1caiGi(n)+n,(2)
通过（2）式，如果求得 G i ( x ) , i ∈ [ 1 , c ) G_i(x),i\in[1,c) Gi(x),i∈[1,c) 则可以得到 G c ( x ) G_c(x) Gc(x)的表达式。
可以看（归纳）出 G c ( x ) G_c(x) Gc(x) 为一个 c + 1 c+1 c+1次的多项式。
公式：
G 1 ( x ) = x ( x + 1 ) 2 = x 2 + x 2 = x ( x + 1 ) 2 G 2 ( x ) = x ( x + 1 ) ( 2 x + 1 ) 6 = 2 x 3 + 3 x 2 + x 6 = x ( x + 1 ) ( 2 x + 1 ) 6 G 3 ( x ) = x 4 + 2 x 3 + x 2 4 = x 2 ( x + 1 ) 2 4 G 4 ( x ) = 6 x 5 + 15 x 4 + 10 x 3 − x 30 = x ( x + 1 ) ( 2 x + 1 ) ( 3 x 2 + 3 x − 1 ) 30 G 5 ( x ) = x 2 ( x + 1 ) 2 ( 2 x 2 + 2 x − 1 ) 12 \large G_1(x)=\frac{x(x+1)}{2}=\frac{x^2+x}{2}=\frac{x(x+1)}{2}\\ \large G_2(x)=\frac{x(x+1)(2x+1)}{6}=\frac{2x^3+3x^2+x}{6}=\frac{x(x+1)(2x+1)}{6}\\ \large G_3(x)=\frac{x^4+2x^3+x^2}{4}=\frac{x^2(x+1)^2}{4}\\ \large G_4(x)=\frac{6x^5+15x^4+10x^3-x}{30}=\frac{x(x+1)(2x+1)(3x^2+3x-1)}{30}\\ \large G_5(x)=\frac{x^2(x+1)^2(2x^2+2x-1)}{12} G1(x)=2x(x+1)=2x2+x=2x(x+1)G2(x)=6x(x+1)(2x+1)=62x3+3x2+x=6x(x+1)(2x+1)G3(x)=4x4+2x3+x2=4x2(x+1)2G4(x)=306x5+15x4+10x3−x=30x(x+1)(2x+1)(3x2+3x−1)G5(x)=12x2(x+1)2(2x2+2x−1)
这样part2就可以解决了。
Tip
part1中，
∑
d
=
1
N
ϕ
(
d
)
d
2
\large\sum_{d=1}^N\phi(d)d^2
∑d=1Nϕ(d)d2 这样的式子一般可以直接卷上
I
d
2
Id^2
Id2 把phi函数系数变为常数：
(
(
I
d
c
ϕ
)
∗
I
d
c
)
(
n
)
=
∑
d
∣
n
ϕ
(
d
)
d
c
×
(
n
d
)
c
=
n
c
+
1
\large ((Id^c\phi)*Id^c)(n)=\sum_{d|n}\phi(d)d^c\times (\frac{n}{d})^c=n^{c+1}
((Idcϕ)∗Idc)(n)=d∣n∑ϕ(d)dc×(dn)c=nc+1
然后卷出来的函数也是很好用于杜教筛的。
code
挑了好久。。。不过跑得飞快，十分满意。
#include<cstdio>
#include<cmath>
#define ll long long
#define mo 998244353
using namespace std;
inline ll ksm(ll x,int y)
{
ll z=1;
for (;y;y>>=1,x=x*x%mo)
if (y&1)z=z*x%mo;
return z;
}
const ll inv3=332748118,inv5=598946612,inv6=166374059;
const int lim=1e6;
int n,m,N,sq;ll C1,C2;
int pri[lim/10];
int phi[3][lim+5],s[3][lim+5];
inline int min(int a,int b){return a<b?a:b;}
inline void init()
{
N=min(n,m);
sq=sqrt(N);
C1=(-(n+1)-(m+1))%mo,
C2=(ll)(n+1)*(m+1)%mo;
s[0][1]=phi[0][1]=
s[1][1]=phi[1][1]=
s[2][1]=phi[2][1]=1;
for (int i=2;i<=lim;++i)
{
if (!phi[0][i])phi[0][pri[++pri[0]]=i]=i-1;
phi[2][i]=(ll)(phi[1][i]=(ll)phi[0][i]*i%mo)*i%mo;
if ((s[0][i]=s[0][i-1]+phi[0][i])>=mo)s[0][i]-=mo;
if ((s[1][i]=s[1][i-1]+phi[1][i])>=mo)s[1][i]-=mo;
if ((s[2][i]=s[2][i-1]+phi[2][i])>=mo)s[2][i]-=mo;
for (int j=1;j<=pri[0]&&i*pri[j]<=lim;++j)
{
if (i%pri[j]==0)
{
phi[0][i*pri[j]]=phi[0][i]*pri[j];
break;
}
phi[0][i*pri[j]]=phi[0][i]*(pri[j]-1);
}
}
}
inline ll G1(int x)
{
return (x&1)?(ll)(x+1>>1)*x%mo:(ll)(x>>1)*(x+1)%mo;
}
inline ll G2(int x)
{
return (ll)x*(x+1)%mo*(2*x+1)%mo*inv6%mo;
}
inline ll G3(int x){return G1(x)*G1(x)%mo;}
inline ll G4(int x)
{
return G2(x)*(3ll*(x+1)*x%mo-1+mo)%mo*inv5%mo;
}
inline ll G5(int x)
{
return G3(x)*(2ll*x*(x+1)%mo-1)%mo*inv3%mo;
}
inline ll F0(int x){return C2*G1(x)%mo;}
inline ll F1(int x){return C1*G2(x)%mo;}
inline ll F2(int x){return G3(x);}
struct box{
ll m1[1005],m2[1005];
inline ll get(int x){return x<=sq?m1[x]:m2[N/x];}
inline void set(int x,ll v){x<=sq?m1[x]=v:m2[N/x]=v;}
}b[3];
inline ll S0(int x)
{
if (x<=lim)return s[0][x];
if (b[0].get(x))return b[0].get(x);
ll re=G1(x);
for (int i=1;i<=x/i;++i)
{
if (i>1)
re=(re-S0(x/i))%mo;
if (i>=x/i)break;
re=(re-S0(i)*((x/i)-(x/(i+1))))%mo;
}
b[0].set(x,re=(re+mo)%mo);
return re;
}
inline ll S1(int x)
{
if (x<=lim)return s[1][x];
if (b[1].get(x))return b[1].get(x);
ll re=G2(x);
for (int i=1;i<=x/i;++i)
{
if (i>1)
re=(re-(ll)i*S1(x/i))%mo;
if (i>=x/i)break;
re=(re-S1(i)*(G1(x/i)-G1(x/(i+1))))%mo;
}
b[1].set(x,re=(re+mo)%mo);
return re;
}
inline ll S2(int x)
{
if (x<=lim)return s[2][x];
if (b[2].get(x))return b[2].get(x);
ll re=G3(x);
for (int i=1;i<=x/i;++i)
{
if (i>1)
re=(re-(ll)i*i%mo*S2(x/i))%mo;
if (i>=x/i)break;
re=(re-S2(i)*(G2(x/i)-G2(x/(i+1))))%mo;
}
b[2].set(x,re=(re+mo)%mo);
return re;
}
inline ll ans1()
{
ll re=0;
for (int i=1;i<=N/i;++i)
{
// d=i
re=(re+phi[0][i]*F0(N/i))%mo;
re=(re+phi[1][i]*F1(N/i))%mo;
re=(re+phi[2][i]*F2(N/i))%mo;
if (i>=N/i)break;
// n/(i+1)<d<=n/i
re=(re+F0(i)*(S0(N/i)-S0(N/(i+1))))%mo;
re=(re+F1(i)*(S1(N/i)-S1(N/(i+1))))%mo;
re=(re+F2(i)*(S2(N/i)-S2(N/(i+1))))%mo;
}
return (re+mo)%mo;
}
inline ll ans2()
{
return (G5(N)*2ll*inv3%mo+
G4(N)*(C1*2ll*inv3%mo-2+mo)%mo+
G3(N)*(C2*2ll*inv3%mo+C1*(mo-2)%mo+inv3)%mo+
G2(N)*(C1*inv3%mo+C2*(mo-2)%mo)%mo+
G1(N)*(C2*inv3%mo)%mo)%mo;
}
int main()
{
freopen("square.in","r",stdin);
freopen("square.out","w",stdout);
scanf("%d%d",&n,&m);
init();
printf("%lld\n",(ans1()*2+ans2())%mo);
return 0;
}