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−1​b=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∑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∑⌊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=1n​ic
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+C1​L+C2​)a=1∑L​(L2−2L+2a2−2aL)=L=1∑N​(L2+C1​L+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=0c​ai​ni,(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=1c​ai​Gi​(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)2​G4​(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;
}

标签：mo,11.23,T4,好题,large,re,int,sum,ll
来源： https://blog.csdn.net/weixin_45010817/article/details/110304508