7/17 D 题题解
x义x
·
·
个人记录
题目大意
给出一个 N\times M 的网格，可以从 (1,1) 向网格中的某个整点发出一条光线，光线撞到四壁会反弹，当光线射入四角之一时结束，求对于每一个终点，可能的光线轨迹的数量。
## 题解
下面为了方便认为网格左上角是 $(0,0)$，右下角是 $(n,m)$。
枚举互质的 $i,j$，我们考虑它的终点。有一个著名 trick，把反射光按镜面对称后可以与入射光拼成一条直线，所以你可以想象一排无限对称复制的网格：
$$\begin{matrix}A&B&A&B&...\\C&D&C&D&...\\A&B&A&B&...\\C&D&C&D&...\\\vdots&\vdots&\vdots&\vdots&\ddots\end{matrix}$$
假设光线射到 $(an,bm)$（需要有 $\gcd(a,b)=1$，否则必然射不到）。那么它射到的是 $A,B,C,D$ 中的哪一个由 $(a\bmod 2,b\bmod 2)$ 决定。
发现有 $a=\dfrac{im}{\gcd(im,jn)},b=\dfrac{jn}{\gcd(im,jn)}$。考虑他们 $\bmod 2$ 的结果，我们发现如果设 $\text{cnt}(x)$ 是 $x$ 的质因子中 2 的数量，则有
$$\begin{cases}\texttt{射到}B&(\text{cnt}(i)+\text{cnt}(m)-\text{cnt}(j)-\text{cnt}(n)>0)\\\texttt{射到}D&(\text{cnt}(i)+\text{cnt}(m)-\text{cnt}(j)-\text{cnt}(n)=0)\\\texttt{射到}C&(\text{cnt}(i)+\text{cnt}(m)-\text{cnt}(j)-\text{cnt}(n)<0)\end{cases}$$
注意 $A$ 是必然射不到的。
设这个函数为 $\text{calc}(i,j)$。射到 $X$ 的方案数即为
$$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=1][\text{calc}(i,j)=X]$$
（这里忽略了射向 $(0,1),(1,0)$ 的两条直线，记得特判。）
反演，同时因为 $\text{cnt}(id)=\text{cnt}(i)+\text{cnt}(d)$，所以有
$$\sum_{d=1}^{\min(n,m)}\mu(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}[\text{calc}(i,j)=X]$$
杜教筛出 $\mu$ 的前缀和即可数论分块。
考虑统计后面那个东西，把 $\text{cnt}(x)=i$ 的 $x$ 塞进桶里再处理一下前后缀和即可 $\log n$ 统计。
参考代码：
```cpp
#include<bits/stdc++.h>
#define ll long long
using namespace std;
ll N,sqN,M,sqM;
const int len=5000000,LEN=5000000;
bool flg[len+5];int pri[len+5],cnt;
ll mu[len+5],smu[len+5];
void init(){
mu[1]=1;
for(int i=2;i<=len;i++){
if(!flg[i]) pri[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*pri[j]<=len;j++){
flg[i*pri[j]]=1;
if(i%pri[j]) mu[i*pri[j]]=-mu[i];
else {mu[i*pri[j]]=0;break;}
}
}
for(int i=1;i<=len;i++) smu[i]=smu[i-1]+mu[i];
}
bool visN[LEN];ll ans_muN[LEN];
ll inline ret_muN(ll x){
if(x<=sqN) return smu[x];
return ans_muN[N/x];
}
void SolveN(ll x){
if(x<=sqN) return;
if(visN[N/x]) return;
visN[N/x]=1;
ans_muN[N/x]=1;
for(ll l=2,r;l<=x;l=r+1){
r=(x/(x/l));
SolveN(x/l);
ans_muN[N/x]-=(r-l+1)*ret_muN(x/l);
}
}
bool visM[LEN];ll ans_muM[LEN];
ll inline ret_muM(ll x){
if(x<=sqM) return smu[x];
return ans_muM[M/x];
}
void SolveM(ll x){
if(x<=sqM) return;
if(visM[M/x]) return;
visM[M/x]=1;
ans_muM[M/x]=1;
for(ll l=2,r;l<=x;l=r+1){
r=(x/(x/l));
SolveM(x/l);
ans_muM[M/x]-=(r-l+1)*ret_muM(x/l);
}
}
int cntn,cntm;
int ans[4];
const int p=998244353;
ll valN[50],valM[50];
ll SvalN[50],SvalM[50];
ll TvalN[50],TvalM[50];
int main(){
init();
scanf("%lld%lld",&N,&M);
N--;M--;sqN=pow(N,2.0/3.0);sqM=pow(M,2.0/3.0);
memset(visN,0,sizeof(visN));
SolveN(N);
memset(visM,0,sizeof(visM));
SolveM(M);
ll N1=N;while(N1%2==0) cntn++,N1/=2;
ll M1=M;while(M1%2==0) cntm++,M1/=2;
for(ll l=1,r,lstans=0,nowans;l<=N&&l<=M;l=r+1){
if(N/(N/l)<M/(M/l)) r=N/(N/l),nowans=ret_muN(r);
else r=M/(M/l),nowans=ret_muM(r);
for(ll i=48;i>=0;i--)
valN[i]=((N/l)>>i)-((N/l)>>(i+1)),
valM[i]=((M/l)>>i)-((M/l)>>(i+1));
SvalM[0]=valM[0];
for(int i=1;i<=48;i++) SvalM[i]=SvalM[i-1]+valM[i];
TvalM[48]=valM[48];
for(int i=47;i>=0;i--) TvalM[i]=TvalM[i+1]+valM[i];
for(int i=0;i<=48;i++){
if(i+cntm-cntn>=0&&i+cntm-cntn<=48) ans[3]=(ans[3]+1LL*((nowans-lstans)%p+p)%p*(valN[i]%p)%p*(valM[i+cntm-cntn]%p)%p)%p;
if(i+cntm-cntn>0)
if(i+cntm-cntn>48) ans[2]=(ans[2]+1LL*((nowans-lstans)%p+p)%p*(valN[i]%p)%p*(SvalM[48]%p)%p)%p;
else ans[2]=(ans[2]+1LL*((nowans-lstans)%p+p)%p*(valN[i]%p)%p*(SvalM[i+cntm-cntn-1]%p)%p)%p;
if(i+cntm-cntn<48)
if(i+cntm-cntn<0) ans[1]=(ans[1]+1LL*((nowans-lstans)%p+p)%p*(valN[i]%p)%p*(TvalM[0]%p)%p)%p;
else ans[1]=(ans[1]+1LL*((nowans-lstans)%p+p)%p*(valN[i]%p)%p*(TvalM[i+cntm-cntn+1]%p)%p)%p;
}
lstans=nowans;
}
printf("0\n%d\n%d\n%d\n",ans[1]+1,ans[2]+1,ans[3]);
return 0;
}
```