洛咕 P3700 [CQOI2017]小Q的表格

神仙题orz

首先推一下给的两个式子中的第二个

(bcdot F(a,a+b)=(a+b)cdot F(a,b))

先简单的想，(F(a,a+b))和(F(a,b))会相互影响

可以换一种角度想，(F(a,b-a))和(F(a,b))会相互影响((b>a))

那么可以从(F(x,y))一路推下去

(F(x,y)=F(x,y-x)=F(x,y-2x)=cdots=F(x,ymod x))

（注意这里的( ext{mod})结果是0的话就没有办法再减了）

这时横坐标比纵坐标大了，利用题目给的式子1，swap横纵坐标

(F(x,y)=F(x,ymod x)=F(ymod x,x)=F(ymod x,xmod(ymod x))=cdots)

总结一下，如果继续这样推下去，当横、纵坐标相等时会就不能再减了

刚才是怎么推的呢，就是当x,y不等的时候每次都把x对y取模然后交换x,y

是不是很熟悉，就是求gcd的过程

那么可以推出来，(gcd(x,y))相等的会相互影响

具体影响多少是显然的，(F(x,y)=F(gcd(x,y),gcd(x,y)) imes frac{xy}{gcd(x,y)^2})

所以只要维护(F(i,i))的值就行了，下面设(F(i,i)=F(i))

现在要算(sum_{i=1}^nsum_{j=1}^nF(i,j))

(ans=sum_{i=1}^nsum_{j=1}^nF(gcd(i,j)))

考虑枚举(gcd(i,j))，

(ans=sum_{k=1}^nF(k)sum_{i=1}^nsum_{j=1}^nfrac{ij}{k^2}[gcd(i,j)=k])

(ans=sum_{k=1}^nF(k)sum_{i=1}^{n/k}sum_{j=1}^{n/k}ij[gcd(i,j)=1])

设(g(n)=sum_{i=1}^{n}sum_{j=1}^{n}ij[gcd(i,j)=1])，(ans=sum_{k=1}^nF(k)g(n/k))

可以只求一半，(g(n)=2sum_{i=1}^{n}sum_{j=1}^{i}ij[gcd(i,j)=1]-sum_{i=1}^{n}i^2[gcd(i,i)=1])

右边那一块是显然的

(g(n)=2sum_{i=1}^{n}sum_{j=1}^{i}ij[gcd(i,j)=1]-1)

设(h(n)=sum_{i=1}^{n}in[gcd(i,n)=1])，(g(n)=2sum_{i=1}^{n}h(i)-1)

怎么求(h)呢，显然可行的(i)有(varphi(n))种，如果(gcd(i,n)=1)，那么(gcd(n-i,n)=1)

所以(i)和(n-i)可以配对，加起来是(n)，一共(varphi(n)/2)对，(h(n)=frac{varphi(n) imes n^2}{2})，注意特判(h(n)=1)，因为不能配对

然后求出了(h)就可以求出(g)，(F)每次只会修改一个。要求的(ans=sum_{k=1}^nF(k)g(n/k))显然(n/k)只有根号种取值，数论分块即可，(F)树状数组前缀和

#include<bits/stdc++.h>
#define il inline
#define vd void
#define int ll
#define mod 1000000007
typedef long long ll;
il int gi(){
    int x=0,f=1;
    char ch=getchar();
    while(!isdigit(ch)){
        if(ch=='-')f=-1;
        ch=getchar();
    }
    while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
    return x*f;
}
int g[4000010],h[4000010];
int pri[4000010],pr,phi[4000010];
bool yes[4000010];
int n,m,f[4000010];
il vd update(int x,int p){while(x<=n)f[x]=(f[x]+p)%mod,x+=x&-x;}
il int query(int x){int ret=0;while(x)ret=(ret+f[x])%mod,x-=x&-x;return ret;}
signed main(){
    m=gi(),n=gi();
    int x,y,k,t,G;
    for(int i=1;i<=n;++i)update(i,1ll*i*i%mod);
    phi[1]=1;
    for(int i=2;i<=n;++i){
        if(!yes[i])phi[i]=i-1,pri[++pr]=i;
        for(int j=1;j<=pr&&i*pri[j]<=n;++j){
            yes[i*pri[j]]=1;
            if(i%pri[j]==0){
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
    }
    for(int i=1;i<=n;++i)h[i]=1ll*phi[i]*i%mod*i%mod;
    g[1]=1;
    for(int i=2;i<=n;++i)g[i]=(g[i-1]+h[i])%mod;
    while(m--){
        x=gi(),y=gi();G=std::__gcd(x,y);
        t=(gi()/(x/G)/(y/G))%mod;
        update(G,(0ll+t-query(G)+query(G-1)+mod)%mod);
        k=gi();
        int ans=0;
        for(int i=1;i<=k;++i){
            int j=k/(k/i);
            ans=(ans+1ll*(query(j)-query(i-1)+mod)*g[k/i])%mod;
            i=j;
        }
        printf("%lld
",ans);
    }
    return 0;
}