Description

小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
之後的結果。

Input

第1行包含兩個整數m,n,表示共有m次操作，所有操作和查詢涉及到的行編號和列編號都不超過n。
接下來m行，每行4個整數a,b,x,k表示把第a行b列的數改成x，然後把它能夠波及到的所有格子全
部修改，保證修改之後所有格子的數仍然都是整數，修改完成後計算前k行前k列裏所有數的和。
1<=m<=10000,1<=a,b,k<=N<=4*10^6,0<=x<=10^18

Output

輸出共m行，每次操作之後輸出1行，表示答案mod 1,000,000,007之後的結果。

Sample Input

3 3

1 1 1 2

2 2 4 3

1 2 4 2

Sample Output

9

36

14

一開始,表格的前3行前3列如圖中左邊所示,前2次操

作後表格沒有變化,第3次操作之後的表格如右邊所示。

1 2 3 2 4 6

2 4 6 4 4 12

3 6 9 6 12 9

簡要題解

觀察表格，可發現
\[ \begin{equation*} f(a,b)=f(gcd(a,b),gcd(a,b)) \cdot \frac{a}{gcd(a,b)} \cdot \frac{b}{gcd(a,b)} \end{equation*} \]
於是我們只需記下\(f(i,i)\)的值即可
利用歐拉函數計算答案，分塊維護

想法

手動模擬表格，至少花了30分鐘我才發現
\[ \begin{equation*} f(a,b)=f(gcd(a,b),gcd(a,b)) \cdot \frac{a}{gcd(a,b)} \cdot \frac{b}{gcd(a,b)} \end{equation*} \]

由上面的式子可知，我們只需記下所有\(f(i,i)\)，其他所有\(f\)值就都可以求出來了
查詢時
\[ \begin{equation*} ans=\sum\limits_{i=1}^k f(i,i) \cdot [\sum\limits_{x=1}^{k/i} \sum\limits_{y=1}^{k/i} xy \cdot (gcd(x,y)==1)] \end{equation*} \]
我們記
\[ \begin{equation*} \begin{aligned} g(n)&=\sum\limits_{x=1}^n \sum\limits_{y=1}^n xy \cdot (gcd(x,y)==1)\ &=2\sum\limits_{x=1}^n \sum\limits_{y=1}^x xy \cdot (gcd(x,y)==1) - \sum\limits_{x=1}^n x^2 \cdot (gcd(x,x)==1) \ &=2\sum\limits_{x=1}^n \sum\limits_{y=1}^x xy \cdot (gcd(x,y)==1) - 1 \end{aligned} \end{equation*} \]
再記
\[ \begin{equation*} \begin{aligned} h(n) &=\sum\limits_{y=1}^n ny \cdot (gcd(n,y)==1) \ &=n \sum\limits_{y=1}^n y \cdot (gcd(n,y)==1) \\end{aligned} \end{equation*} \]
這個式子可以理解為，1~n中所有與n互質的數的和 × n
有歐拉函數 \(\varphi(n)\) 表示n以內與n互質的數的個數
我們知道，若x與n互質，則n-x也與n互質。
也就是說與n互質的數可兩兩配對，和為n。
那麽
\[ \begin{equation*} h(n)=\frac{\varphi(n) \cdot n^2}{2} \end{equation*} \]
故
\[ \begin{equation*} g(n)=2 \sum\limits_{x=1}^n h(x) - 1 \end{equation*} \]
那麽
\[ \begin{equation*} ans=\sum\limits_{i=1}^k f(i,i) \cdot g(k/i) \end{equation*} \]
由於我們知道 \(k/i\) 只有 \(\sqrt{n}\) 種取值
於是就可以愉快地分塊 \(f(i,i)\) 啦！

設計一個 \(O(\sqrt{n})\) 修改， \(O(1)\) 查詢的分塊結構就好了……

代碼

註：調了若幹次才終於對…心好累

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>

#define P 1000000007

using namespace std;

typedef long long ll;
const int N = 4000005;
const int SN = 4005;

ll phi[N];
int prime[N],pnum;
int n,sn,block;
void getphi(){
    phi[1]=1;
    for(int i=2;i<=n;i++) phi[i]=i-1;
    for(int i=2;i<=n;i++){
        if(phi[i]==i-1) prime[pnum++]=i;
        for(int j=0;j<pnum && (ll)prime[j]*i<=n;j++){
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
        phi[i]=(phi[i-1]+(ll)i*i%P*phi[i])%P; //g(i)
    }
}

ll num[N],add[SN];
ll query(int x){
    if(!x) return 0;
    int bl=(x-1)/sn+1;
    return (num[x]+add[bl])%P;
}
void change(int x,ll c){
    int bl=(x-1)/sn+1;
    for(int i=x;i<=min(bl*sn,n);i++) num[i]=(num[i]+c)%P;
    for(int i=bl+1;i<=block;i++) add[i]=(add[i]+c)%P;  
}
int gcd(int a,int b) { return b ? gcd(b,a%b) : a; }

int main()
{
    int m;
    scanf("%d%d",&m,&n);
    sn=(int)sqrt(n); block=(n-1)/sn+1;
    for(int i=1;i<=n;i++) num[i]=(num[i-1]+(ll)i*i%P)%P;
    getphi();
    
    int a,b,k,g;
    ll ans=0,x,X;
    while(m--){
        scanf("%d%d%lld%d",&a,&b,&x,&k);
        if(a<b) swap(a,b);
        g=gcd(a,b);
        X=(x/(1ll((ll)a/g)*(b/g)))%P;  //罪魁禍首！！就是這裏！！！
        change(g,(X-(query(g)-query(g-1)+P)%P+P)%P);
        ans=0;
        for(int l=1,r;l<=k;l=r+1){
            r=k/(k/l);
            ans=(ans+(query(r)-query(l-1)+P)%P*phi[k/l]%P)%P;    
        }
        printf("%lld\n",ans);
    }
    
    return 0;
}

[bzoj4815] [洛谷P3700] [Cqoi2017] 小Q的表格