题目大意
给定一个N*M的网格。对一个顶点为格点的正方形R(不一定与格线平行),计算出其中有多少个单位格被R完全包含(记作F(R))。求所有正方形的F(R)之和。
题解
首先画一个“勾股图”:
粉色是我们的正方形(不和格线平行),设其外接正方形的边长为L,四周直角三角形的短边为a,则长边为L-a。设其中完整包含了F(L,a)个单位格。
我们这么计算第二部分:一红一白两个直角三角形拼成的矩形内有(L-a)*a个单位格,而对角线穿过了L-gcd(L,a)格,再除以2,就得到了一个三角形内的单位格数。将其乘以4,总共加起来得到:
没错,“斜线下方的格点数”没有简单公式,但“单位格数”有!
可以发现,F(L,a)即使当a>L-a时也成立,所以我们可以直接用a枚举到底,整道题目的答案就是:
这分成两部分:第一部分是一个关于L的多项式,第二部分是gcd(L,a)的前缀和。
至于怎么得到这个多项式呢,如果你数学好可以手推……作为一名数学技能荒废一年的咸鱼选手,我就直接用vector当低配版多项式类,让电脑推了……
前者很简单:L不超过10^6,暴力递推即可。后者也差不多,只需要先递推出来phi函数,然后枚举gcd(a,L)的值,就可以把Σgcd的表都打出来,详见代码。
代码
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<vector>
using namespace std;
typedef long long LL;
typedef vector<LL> Poly;
const LL MOD=1000000007;
const int H=6;
LL gcd(LL a,LL b){
return !b?a:gcd(b,a%b);
}
LL lcm(LL a,LL b){
return a*b/gcd(a,b);
}
LL quickpow(LL a,LL n){
LL ans=1;
while(n){
if(n&1) ans=(ans*a)%MOD;
a=(a*a)%MOD;
n>>=1;
}
return ans;
}
LL inverse(LL a){
return quickpow(a,MOD-2);
}
void print(Poly a){
for(int i=0;i<a.size();i++){cout<<a[i]<<" ";}cout<<endl;
}
Poly operator * (Poly a,LL b){
for(int i=0;i<a.size();i++){
a[i]*=b;
a[i]%=MOD;
}
return a;
}
Poly operator / (Poly a,LL b){
return a*inverse(b);
}
Poly operator + (Poly a,Poly b){
Poly c;
int n=max(a.size(),b.size());
c.resize(n);
for(int i=0;i<a.size();i++){
c[i]+=a[i];
c[i]%=MOD;
}
for(int i=0;i<b.size();i++){
c[i]+=b[i];
c[i]%=MOD;
}
return c;
}
Poly operator * (Poly a,Poly b){
Poly c;
c.resize(a.size()+b.size()-1);
for(int i=0;i<a.size();i++){
for(int j=0;j<b.size();j++){
c[i+j]+=a[i]*b[j]%MOD;
c[i+j]%=MOD;
}
}
return c;
}
const int SIZEN=1000010;
LL calc_presum_with(Poly a,LL ps[]){
LL ans=0;
for(int i=0;i<a.size();i++){
ans+=(a[i]*ps[i])%MOD;
ans%=MOD;
}
return ans;
}
Poly give_Fsum_L(void){//the component without gcd
static LL A[10];
A[0]=1;Poly ans(A,A+1);
A[0]=0,A[1]=1;ans=ans*Poly(A,A+2);
A[0]=1,A[1]=1;ans=ans*Poly(A,A+2);
A[0]=1,A[1]=2;ans=ans*Poly(A,A+2);
ans=ans/3;
A[0]=0,A[1]=0,A[2]=-3;ans=ans+Poly(A,A+3);
return ans;
}
LL G[SIZEN]={0};
LL phi[SIZEN];
LL P[SIZEN][H],PG[SIZEN][H];
void work(LL N,LL M){
static LL A[10];
if(N>M) swap(N,M);
//N是较小者
A[0]=M+1,A[1]=-1;Poly D(A,A+2);
A[0]=N+1,A[1]=-1;D=D*Poly(A,A+2);
Poly F=give_Fsum_L();
Poly F1=D*F;
LL ans=0;
ans+=calc_presum_with(F1,P[N]);
ans%=MOD;
A[0]=2;
Poly F2=D*Poly(A,A+1);
ans+=calc_presum_with(F2,PG[N]);
ans%=MOD;
ans=(ans+MOD)%MOD;
printf("%I64d\n",ans);
}
void prepare(void){
for(int i=1;i<SIZEN;i++){//精妙的求phi方法
phi[i]+=i;
for(int j=i+i;j<SIZEN;j+=i){
phi[j]-=phi[i];
}
}
for(int g=1;g<SIZEN;g++){//求gcd(x,i)的前缀和
G[g]+=g;
G[g]%=MOD;
for(int x=2*g,d=2;x<SIZEN;x+=g,d++){
G[x]+=(g*phi[d]%MOD);
G[x]%=MOD;
}
}
for(int i=0;i<H;i++){P[0][i]=PG[0][i]=0;}
for(LL x=1;x<SIZEN;x++){
LL p=1;
for(int i=0;i<H;i++){
P[x][i]=(P[x-1][i]+p)%MOD;
PG[x][i]=(PG[x-1][i]+(p*G[x])%MOD)%MOD;
p=(p*x)%MOD;
}
}
}
int main(){
prepare();
LL T,N,M;
scanf("%I64d",&T);
while(T--){
scanf("%I64d%I64d",&N,&M);
work(N,M);
}
return 0;
}