Luogu P3700「CQOI2017」小Q的表格

为什么我连分块都想不到啊...

---

题意

定义一个矩阵$f$满足

$ f(a,b)=f(b,a)$

$ b·f(a,a+b)=(a+b)·f(a,b)$

初始$ f(a,b)=ab$

有$ m$次修改，每次会将$ f(x,y)$改成$ c$，并修改这个矩阵使得仍然满足以上两条

保证修改后矩阵中每个数均为正整数

你需要在每次修改后求出$ \sum\limits_{i=1}^k\sum\limits_{j=1}^kf(i,j)$

数据范围满足$ m \leq 10000$，涉及到的矩形范围不超过$4·10^6$

---

$Solution $

给定的信息类似一个辗转相减的过程

和求$ gcd(x,y)$的过程非常像

我们不断辗转相减可以得到$ f(x,y)=f(d,d)*\frac{xy}{d^2}$

其中$ d=gcd(x,y)$

考虑求答案的式子的意义

$ ans=\sum\limits_{i=1}^k\sum\limits_{j=1}^kf(i,j)$

$ans=\sum\limits_{k=1}^nf(k,k)\sum\limits_{i=1}^{\frac{n}{k}}\sum\limits_{j=1}^{\frac{n}{k}}ij[gcd(i,j)=1]$

设

$ A(n)=\sum\limits_{i=1}^n\sum\limits_{j=1}^nij[gcd(i,j)=1]$

有

$A(n)=\sum\limits_{i=1}^n i^2\varphi(i)$

原理是若$ x<y,x和y互质则y-x和y互质$

这样可以化简为

$ans=\sum\limits_{k=1}^nf(k,k)A(\frac{n}{k})$

显然$ A(x)$可以预处理，然后数论分块计算

单次$ O(\sqrt{4000000})=O(2000)$

考虑修改

现在问题是：单点修改，区间求和

修改次数只有$ 10000$次，而求和次数达到$ 2*10^7$次

考虑分块

对于每个块我们维护块内前缀和以及前$ i$个块的和

这样查询是$ O(1)$的而修改是$ O(2000)$的

就以一种非常优秀的均摊复杂度解决了这道题

---

$ my \ code$ 

#include<ctime>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#define p 1000000007
#define rt register int
#define ll long long
using namespace std;
inline ll read(){
    ll x = 0; char zf = 1; char ch = getchar();
    while (ch != '-' && !isdigit(ch)) ch = getchar();
    if (ch == '-') zf = -1, ch = getchar();
    while (isdigit(ch)) x = x * 10 + ch - '0', ch = getchar(); return x * zf;
}
void write(ll y){if(y<0)putchar('-'),y=-y;if(y>9)write(y/10);putchar(y%10+48);}
void writeln(const ll y){write(y);putchar('\n');}
int i,j,k,m,n,x,y,z,cnt,blo;
int phi[4000010],ss[4000010],inv[4000010],f[4000010];bool pri[4000010];
int val[4000010],qz1[2010][2010],qz2[2010];
void init(int M){
    //mu[1]=1;
    phi[1]=1;blo=(int)sqrt(M);
    for(rt i=2;i<=M;i++){
        if(!pri[i])phi[i]=i-1,ss[++cnt]=i;
        for(rt j=1;j<=cnt&&i*ss[j]<=M;j++){
            pri[i*ss[j]]=1;
            phi[i*ss[j]]=phi[i]*phi[ss[j]];
            if(i%ss[j]==0){
                phi[i*ss[j]]=phi[i]*ss[j];
                break;
            }
        }
    }
    for(rt i=1;i<=M;i++)val[i]=(val[i-1]+1ll*i*i%p*phi[i]%p)%p;
    for(rt i=1;i<=M;i++){
        f[i]=1ll*i*i%p;
        if((i-1)%blo==0)qz1[(i-1)/blo][0]=f[i];
        else qz1[(i-1)/blo][(i-1)%blo]=(qz1[(i-1)/blo][(i-1)%blo-1]+f[i])%p;
    }
    qz2[0]=qz1[0][blo-1];
    for(rt i=1;i<=blo;i++)qz2[i]=(qz2[i-1]+qz1[i][blo-1])%p;
    inv[0]=inv[1]=1;
    for(rt i=2;i<=M;i++)inv[i]=1ll*inv[p%i]*(p-p/i)%p;
}
int qz(int x){
    if(!x)return 0;int ks=(x-1)/blo;
    if(ks==0)return qz1[0][x-1];
    else return (qz1[ks][(x-1)%blo]+qz2[ks-1])%p;
}
int calc(int n){
    int ans=0;
    for(rt i=1;i<=n;){
        int R=n/(n/i);
        (ans+=1ll*(qz(R)-qz(i-1))*val[n/i]%p)%=p;
        i=R+1;
    }
    return (ans+p)%p;
}
void update(int x,int y){
    int ks=(x-1)/blo;
    for(rt i=x;(i-1)%blo!=0||i==x;i++)(qz1[ks][(i-1)%blo]+=y)%=p;
    for(rt i=ks;i<=blo;i++)(qz2[i]+=y)%=p;
    (f[x]+=y)%=p;
}
int main(){
    n=read();m=read();
    init(m);
    for(rt i=1;i<=n;i++){
        x=read();y=read();int v=read()%p;k=read();
        int gc=__gcd(x,y);
        v=1ll*v*inv[x]%p*inv[y]%p*gc%p*gc%p;
        update(gc,v-f[gc]);
        writeln(calc(k));
    }
    return 0;
}