題目描述

小Q是個程式設計師。

作為一個年輕的程式設計師，小Q總是被老C欺負，老C經常把一些麻煩的任務交給小Q來處理。每當小Q不知道如何解決時，就只好向你求助。

為了完成任務，小Q需要列一個表格，表格有無窮多行，無窮多列，行和列都從1開始標號。為了完成任務，表格裡面每個格子都填了一個整數，為了方便描述，小Q把第a行第b列的整數記為f(a,b)。為了完成任務，這個表格要滿足一些條件：

(1)對任意的正整數a,b，都要滿足f(a,b)=f(b,a)；

(2)對任意的正整數a,b，都要滿足b×f(a,a+b)=(a+b)×f(a,b)。

為了完成任務，一開始表格裡面的數很有規律，第a行第b列的數恰好等於a×b，顯然一開始是滿足上述兩個條件的。為了完成任務，小Q需要不斷的修改表格裡面的數，每當修改了一個格子的數之後，為了讓表格繼續滿足上述兩個條件，小Q還需要把這次修改能夠波及到的全部格子裡都改為恰當的數。由於某種神奇的力量驅使，已經確保了每一輪修改之後所有格子裡的數仍然都是整數。為了完成任務，小Q還需要隨時獲取前k行前k列這個有限區域內所有數的和是多少，答案可能比較大，只需要算出答案mod1,000,000,007之後的結果。

題解

這題思路太神了。

乍一看需要維護一個二維的東西看起來很麻煩，但實際上是可以轉化為一維上的問題的。

先來看第一個限制，關於對角線對稱的數字相等，這樣我們的棋盤大小變為原來的一半。

然後看第二個限制，挪一下位置，變成了

f(a,a+b)/(a+b)=f(a,b)/b

f(a,a+b)/((a+b)*a)=f(a,b)/(a*b)

我們發現了規律，f(a,b)是和a*b正相關的，而且這種關係非常像輾轉相減，那麼一直減到最後就會變成f(gcd,gcd)。

於是我們得出了結論，每個位置(a,b)的值之和gcd(a,b)相關。

於是它變成了一個序列上的問題，每次修改只需要修改一個數。

然後考慮答案。

ans=∑f(g,g)∑∑i*j*(gcd(i,j)==g)  i,j,g<=k

後面的那個東西∑∑i*j*(gcd(i,j)==g)經過推導後可以得到∑i*i*Φ(i) i<=k/g

然後就可以除法分塊了。

但我們還需要支援單點修改，區間求和，可以用樹狀陣列，但是這題的資料範圍比較特別，所以用分塊搞就可以了，根號修改，O(1)查詢。

程式碼

#include<iostream>
#include<cstdio>
#include<cmath>
#define N 4000009
using namespace std;
typedef  
 
long long ll;
const int mod=1e9+7;
ll n,m,sum[N],pre[N],phi[N],prime[N],f[N],val[N],n1,be[N];
bool vis[N];
inline ll rd(){
    ll x=0;char c=getchar();bool f=0;
    while(!isdigit(c)){if(c=='-')f=1;c=getchar();}
    while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
    return f?-x:x;
} 
inline void work(ll x,ll y){
    pre[x]=(pre[x]+y)%mod;
    for(int i=x+1;be[i]==be[i-1];++i)pre[i]=(pre[i]+y)%mod;
    for(int i=be[x];i<=be[n];++i)sum[i]=(sum[i]+y)%mod;
}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
inline ll query(int x){return (pre[x]+sum[be[x]-1])%mod;}
inline void prework(){
    phi[1]=1;
    for(int i=2;i<=n;++i){
        if(!vis[i])prime[++prime[0]]=i,phi[i]=i-1;
        for(int j=1;j<=prime[0]&&i*prime[j]<=n;++j){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;}
            phi[i*prime[j]]=phi[i]*phi[prime[j]];
        }
    }
    for(int i=1;i<=n;++i)f[i]=(f[i-1]+phi[i]*val[i])%mod;
} 
int main(){
    m=rd();n=rd();n1=sqrt(n);
    for(int i=1;i<=n;++i){
        be[i]=(i-1)/n1+1;val[i]=1ll*i*i%mod;
        (sum[be[i]]+=val[i])%=mod;
        if(be[i]==be[i-1])pre[i]=(pre[i-1]+val[i])%mod;
        else pre[i]=val[i];    
    }
    for(int i=2;i<=be[n];++i)(sum[i]+=sum[i-1])%=mod;
    prework();ll a,b,x,k;
    while(m--){
        a=rd();b=rd();x=rd();k=rd();ll g=gcd(a,b);ll ans=0,r;
        x=x/((a/g)*(b/g));x%=mod;work(g,(x-val[g]+mod)%mod);val[g]=x;
        for(ll l=1;l<=k;l=r+1){
            r=k/(k/l);
            ans+=(query(r)-query(l-1))*f[k/l]%mod;
            ans=(ans%mod+mod)%mod;
        }
        printf("%lld\n",ans);
    }
    return 0;
}