11.23A T4方（思维好题----杜教筛（总复习））

11.23A T4方（思维好题----杜教筛）

Description

给一个$n\times m$的网格（一共$(n+1)\times (m+1)$个格点），定义格点正方形的权值时它完全包含的格子数，求网格中所有正方形（可以是斜的）的权值和。

Solution

容易发现一个正方形的权值就是其面积减去两条斜边穿过的格子。

面积易求。穿过的格子怎么办？让人没有头绪。
考虑穿过的格子坐标满足的条件？那么线段应该与格子有两条边相交，限制比较复杂。

但在思考的时候我们可以发现 格子是否被穿过 似乎可以与 边的相交 联系起来。
于是有这样的思路：当线段穿过一条网格边时，它就进入了一个新格子，且它不会再次回到它。而经过一条横边则意味着线段在纵轴上升高了1，纵边同理。

不过上述方式有一个前提—经过的是网格边；那么经过格点时呢，也是进入了一个新格子，为了方便，我们在上一种情况中一起计算，再减去重复的就好。

如果设斜边在坐标轴上的投影分别为 $a,b$ ，那么有：
$$cnt(a,b)=a+b-\gcd (a,b)$$
则可以通过枚举 $a,b$ 得到答案：
$$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]$$
然后就是推式子环节。

首先可以把含gcd的因式分开考虑（part1）：

设 $L=a+b$，记 $N=\min (n,m)$
$$\begin{array}{rl}ans1& =\sum _{L=1}^N(n-L+1)(m-L+1)\sum _{a=1}^L\gcd (a,L-a)\\ & =\sum _{L=1}^N(n-L+1)(m-L+1)\sum _{a=1}^L\sum _{d|a,d|L}\varphi (d)\\ & =\sum _{d=1}^N\varphi (d)\sum _{L^{\prime }=1}^{\lfloor \frac{N}{d}\rfloor }(n-L^{\prime }d+1)(m-L^{\prime }d+1)L^{\prime }\\ & =\sum _{d=1}^N\varphi (d)F_0(\lfloor \frac{N}{d}\rfloor )+\sum _{d=1}^N\varphi (d)d\times F_1(\lfloor \frac{N}{d}\rfloor )+\sum _{d=1}^N\varphi (d)d^2\times F_2(\lfloor \frac{N}{d}\rfloor )\end{array}$$

此处直接上 数论分块+杜教筛。

part2：

记$C_1=-(n+1)-(m+1),C_2=(n+1)(m+1),G_c(n)={\sum }_{i=1}^n{i}^c$
$$\begin{array}{rl}ans2& =\sum _{L=1}^N(L^2+C_{1}L+C_2)\sum _{a=1}^L(L^2-2L+2a^2-2aL)\\ & =\sum _{L=1}^N(L^2+C_{1}L+C_2)(L^3-2L^2+2G_2(L)-2L{G}_1(L))\end{array}$$
part1,part2都用到了n以内自然数的幂的和。补补公式推导：

我们用 $G_i(n)$ 表示 $n$ 以内自然数的 $i$ 次幂的和。

$(n+1{)}^{c+1}-n^{c+1}={\sum }_{i=0}^c{a}_i{n}^i,(1)$

$n^{c+1}-(n-1{)}^{c+1}={\sum }_{i=0}^c{a}_i(n-1{)}^i,(2)$

$...$

$2^{c+1}-1^{c+1}={\sum }_{i=0}^c{a}_i\times 1^i,(n)$

将左右式分别求和，可以推出：

$(n+1{)}^{c+1}-1={\sum }_{i=0}^c{a}_i{G}_i(n)$

通过上式，如果求得 $G_i(x),i\in [1,c)$ 则可以得到$G_c(x)$的表达式。

可以看（归纳）出 $G_c(x)$ 为一个$c+1$次的多项式。

公式：
$$\begin{gathered} G_1(x)=\frac{x^2+x}{2}=\frac{x(x+1)}{2} \\ {G}_2(x)=\frac{2x^3+3x^2+x}{6}=\frac{x(x+1)(2x+1)}{6} \\ {G}_3(x)=\frac{x^4+2x^3+x^2}{4}=\frac{x^2(x+1{)}^2}{4} \\ {G}_4(x)=\frac{6x^5+15x^4+10x^3-x}{30}=\frac{x(x+1)(2x+1)(3x^2+3x-1)}{30} \\ {G}_5(x)=\frac{x^2(x+1{)}^2(2x^2+2x-1)}{12} \end{gathered}$$
这样 $part2$ 就可以解决了。

Tip

$part1$ 中，${\sum }_{d=1}^N\varphi (d)d^2$ 这样的式子一般可以直接卷上 $I{d}^2$ 把 $phi$ 函数系数变为常数：
$$((I{d}^c\varphi )\ast I{d}^c)(n)=\sum _{d|n}\varphi (d)d^c\times (\frac{n}{d}{)}^c=n^{c+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;
}