[CF 449E] Jzzhu and Squares

作法：

考慮一個斜的正方形，設他下面頂點指向右邊頂點的向量為$(x,y)$（也就是如果下端點是$(0,0)$，那麼剩下三個點的座標為$(x,y),(x-y,x+y),(-y,x)$），我們想要知道這個正方形裡面包含了幾個單位正方形，而一個著名的結論是：給定一個$n\times m$的方格表，那麼他的對角線會經過$n+m-gcd(n,m)$個格點。因此就可以用這個結論算出這個正方形包住了幾個單位正方形了：把這個正方形補成邊平形座標軸的正方形，並把他切成四個$x\times y$的長方形，還有中間邊長$|x-y|$的正方形。那麼就可以知道在四個長方形中分別包含了$\frac{xy-x-y+gcd(x,y)}{2}$個所求格子了，再加上中間的正方形$(x-y)^2$個，就共有$(x-y)^2+2(xy-x-y+gcd(x,y))=x^2+y^2-2(x+y-gcd(x,y))$個。並且注意到這種正方形因為要用長$x+y$的正方形包住他，所以會在原本給定的範圍裡出現$(n-x-y+1)(m-x-y+1)$次，其中$x+y\leq min(n,m)$。並且注意到我們也要算$x=0$的情況（上面通式也會對），但是不能也算$y=0$的，否則會重複算到。因此就可以把答案寫成：
$\displaystyle \sum_{x\geq 0,y>0,x+y\leq min(n,m)} (n-x-y+1)(m-x-y+1)(x^2+y^2-2(x+y-gcd(x,y)))$
考慮把$x+y$值相同的一起算，令$t=x+y$，那麼就可以改寫成：
$\displaystyle \sum_{t=1}^{min(n,m)}\sum_{x=0}^{t-1} (n-t+1)(m-t+1)(x^2+(t-x)^2-2t+2gcd(t,x))$
其中用到了簡單的$gcd$的性質：$gcd(x,t-x)=gcd(x,t)$。我們可以把所求的$gcd$的部份拆出來，也就是改寫成
$\displaystyle \sum_{t=1}^{min(n,m)}\sum_{x=0}^{t-1} (n-t+1)(m-t+1)(x^2+(t-x)^2-2t)+$$\displaystyle 2\sum_{t=1}^{min(n,m)}\sum_{x=0}^{t-1} (n-t+1)(m-t+1)gcd(t,x)$
首先第一項其實就是一個$t$的多項式，只是係數會有$n,m$之類的，因此我們可以先預處理好$t^i$的前綴和陣列，其中$i=1,...,6$，就可以在$O(1)$算出第一項了（或是借住電腦的力量直接把通式爆出來XD）。再來則是第二項，把$(n-t+1)(m-t+1)$乘開後會得到其實我們主要需要的東西就是$\displaystyle \sum_{x=0}^{t-1} gcd(t,x)$，還有$\displaystyle \sum_{x=0}^{t-1} t\cdot gcd(t,x)$和$\displaystyle \sum_{x=0}^{t-1} t^2\cdot gcd(t,x)$的前綴和陣列，而我們只要有辦法對每個$t$都算出$\displaystyle \sum_{x=0}^{t-1} gcd(t,x)$的值就可以了。記我們要求的這個陣列為$f$，我們看某個數$i$在哪些$gcd(t,x)$的項出現了，首先$i$要整除$t$，因此只有index 被$i$整除的$f$值要加上一些$i$。至於要加上多少個$i$，因為$gcd(t,x)=i$，所以我們也可以設$x=y\cdot i$，其中$y<\frac{t}{i}$。這樣就可以兩邊同除以$i$，變成$gcd(\frac{t}{i},y)=1$，因此我們需要知道$1,...,\frac{t}{i}-1$中有幾個數和$\frac{t}{i}$互質，這正是歐拉函數 phi ，因此只要先預處理好每個數的 phi 值，再用他算出$f$陣列的值就可以了。

code :

#include<bits/stdc++.h>
#define LL long long
#define MOD 1000000007
using namespace std;
const int maxn=1000000+10 ;
 
LL pw(LL x,int n)
{
    if(n<=1) return n ? x : 1 ;
    LL t=pw(x,n/2) ;
    if(n%2) return (t*t%MOD)*x%MOD ;
    else return t*t%MOD ;
}
LL inv(LL x)
{
    return pw(x,MOD-2) ;
}
 
int phi[maxn] ;
LL f[maxn],s1[maxn],s2[maxn],s3[maxn] ;
main()
{
    for(int i=1;i<maxn;i++)
    {
        phi[i]+=i ;
        for(int j=2*i;j<maxn;j+=i) phi[j]-=phi[i] ;
    }
    for(int i=1;i<maxn;i++)
    {
        f[i]=(f[i]+i)%MOD ;
        s1[i]=(s1[i-1]+f[i])%MOD ;
        s2[i]=(s2[i-1]+f[i]*i)%MOD ;
        s3[i]=(s3[i-1]+f[i]*i%MOD*i)%MOD ;
        for(int j=2*i,k=2;j<maxn;j+=i,k++) f[j]+=i*phi[k] ;
    }
 
    int T ; scanf("%d",&T) ;
    while(T--)
    {
        LL n,m,r ; scanf("%lld%lld",&n,&m) ;
        r=min(m,n) ;
        LL X=r*(r+1)%MOD*inv(180)%MOD ;
        LL Y=6*m*(5*n*((r*r-3*r-1)%MOD)%MOD-4*r*r%MOD*r%MOD+14*r*r%MOD-4*r-6)%MOD
                 +(r-1)*(n*((-24*r*r%MOD+60*r+36)%MOD)%MOD+
                         (20*r*r%MOD)*r%MOD-60*r*r%MOD-5*r+30)%MOD ;
        LL ans=X*Y%MOD ;
        ans+=((m+1)*(n+1)%MOD*s1[r]%MOD-(m+n+2)*s2[r]%MOD+s3[r])*2%MOD + MOD ;
        printf("%lld\n",((ans%MOD)+MOD)%MOD) ;
    }
}