BZOJ4815: [Cqoi2017]小Q的表格

BZOJ4815

$b×f(a,a+b)=(a+b)*f(a,b)$很像辗转相减法。。
那么每次修改点$(a,b)$的值，会修改所有满足$gcd(i,j)==gcd(a,b)$的点$(i,j)$的值。
记$d=gcd(a,b)$，那么$f_{i,j}=x*i/d*j/d$$(gcd(i,j)==d$

那么可以转化为，每次修改对角线上的值。记$Num_i$为$f{i,i}$

$$ans=\sum_{d=1}^nNum_d\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)==d]i/d*j/d$$
$$ans=\sum_{d=1}^nNum_d\sum_{I=1}^{n/d}\sum_{j=1}^{n/d}[gcd(i,j)==1]i*j$$
令 $G(n)=\sum_{I=1}^{n}\sum_{j=1}^{n}[gcd(i,j)==1]i*j$
根据 BZOJ2154: Crash的数字表格 的做法
$$G(n)=\sum_{i=1}^n μ(i)*i^2*S^2(n/i)$$
$$S(i)=\frac {i(i+1)}{2}$$
$n\sqrt n$求 $G(n)$无法接受，复杂度还是高了一些。
考虑递推。当且仅当 $n|d$时， $\lfloor \frac{n}{d}\rfloor$比 $\lfloor\frac{n-1}{d}\rfloor$大1
那么 $G(n)=G(n-1)+\sum_{i|n}μ(i)i^2(S^2(\lfloor\frac{n}{i}\rfloor)-S^2(\lfloor\frac{n-1}{i}\rfloor)$
令 $H(n)=\sum_{i|n}μ(i)i^2(S^2(\lfloor\frac{n}{i}\rfloor)-S^2(\lfloor\frac{n-1}{i}\rfloor)=n^3\sum_{i|n}μ(i)\frac{1}{i}$
令 $h(n)=\sum_{i|n}μ(i)\frac{1}{i}$
发现 $h(n)$是一个积性函数。
证明一下：
对于两个互质的数 $x,y$
记 $x$的因子集为$A$，大小为 $s1$
记 $y$的因子集为$B$，大小为 $s2$
那么 $x*y$的因子集 $C$为$A*B$，大小为 $s1*s2$
举个例子 $x=10,y=21$
$A=1,2,5,10\space$ $s1=4$
$B=1,3,7,21\space$ $s2=4$
那么 $$C= \left( \begin{matrix} 1&2&5&10\\ 3&6&15&30\\ 7&14&35&70\\ 21&42&105&210 \end{matrix} \right)$$

显然$μ(a)\frac{1}{a}*μ(b)\frac{1}{b}=μ(a*b)\frac{1}{a*b}$
那么相当于$h(xy)$为$A$中每个$μ(a)*\frac {1}{a}$乘以$B$中每个$μ(b)*\frac{1}{b}$
然后就可以线性筛出$h(n)$

$$ans=\sum_{d=1}^n Num_d G(n/d)$$

然后就可以下底函数分块来求了。但是有修改操作,如果用树状数组，修改复杂度为$O(log_n)$，但是一次查询的复杂度为$O(log_n\sqrt n)$，会$T$掉。
那么可以考虑牺牲一下修改的复杂度，用分块来维护前缀和，修改变成$O(\sqrt n)$，查询变成$O(1)$了。

然后这个题终！于！做！完！了！
$19s$算是给卡过去了。。好像$S^2(n)=\sum_{i=1}^{n}φ(i)*i^2$但是不理解为什么啊$QAQ$，如有神犇愿意解答感激不尽！

【代码】

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <queue>
#define N 4000005
#define M 600005
#define Mod 1000000007
#define INF 0x7fffffff
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const ull base=31;

ll read()
{
    ll x=0,f=1;char ch=getchar();
    while(!isdigit(ch)){if(ch=='-') f=-1;ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
    return x*f;
}

int m,n,block,Mx;
int Prime[N],bl[N],R[N];bool Not_Prime[N]={1,1};
ll f[N],h[N],g[N],Inv[N],sum[2005][2005],Sum[2005];

ll SSqr(int x){
    return 1LL*x*x%Mod*x%Mod;
}

ll gcd(ll a,ll b){
    return b==0?a:gcd(b,a%b);
}

void Pre()
{
    Inv[1]=h[1]=g[1]=1;block=sqrt(m);
    for(register int i=2;i<=m;i++)
    {
        Inv[i]=(Mod-Mod/i)*Inv[Mod%i]%Mod;
        if(!Not_Prime[i]) Prime[++Prime[0]]=i,h[i]=(1-Inv[i]+Mod)%Mod;
        for(register int j=1;j<=Prime[0]&&i*Prime[j]<=m;j++)
        {
            Not_Prime[i*Prime[j]]=1;
            if(i%Prime[j]==0) {
                h[i*Prime[j]]=h[i];
                break;
            }
            h[i*Prime[j]]=h[i]*h[Prime[j]]%Mod;
        }
        g[i]=(g[i-1]+SSqr(i)*h[i]%Mod)%Mod;
    }
    for(register int i=1;i<=m;i++) {
        f[i]=1LL*i*i;
        bl[i]=(i-1)/block+1;
        Sum[bl[i]]=(Sum[bl[i]]+f[i])%Mod;
        sum[bl[i]][i-R[bl[i]-1]]=(sum[bl[i]][i-R[bl[i]-1]-1]+f[i])%Mod;
        if(i==n||(i%block==0)) R[bl[i]]=i;
    } 
    Mx=bl[m];
    for(int i=1;i<=Mx;i++)
        Sum[i]=(Sum[i]+Sum[i-1])%Mod;
}

void Update(ll d)
{
    int Bl=bl[d];
    for(int i=d;i<=R[Bl];i++)
        sum[Bl][i-R[Bl-1]]=(sum[Bl][i-1-R[Bl-1]]+f[i])%Mod;
    Sum[Bl]=(Sum[Bl-1]+sum[Bl][R[Bl]-R[Bl-1]])%Mod;
    for(int i=Bl+1;i<=Mx;i++)
        Sum[i]=Sum[i-1]+sum[i][R[i]-R[i-1]];
}

ll Get_Sum(int x,int y)
{
    int sx=bl[x],sy=bl[y];
    if(sx==sy) return (sum[sx][y-R[sx-1]]-sum[sx][x-1-R[sx-1]]+Mod)%Mod;
    return ((sum[sx][R[sx]-R[sx-1]]-sum[sx][x-1-R[sx-1]]+sum[sy][y-R[sy-1]]+Sum[sy-1]-Sum[sx])%Mod+Mod)%Mod;
}

void Solve(int x)
{
    int pos;
    ll ans=0;
    for(int i=1;i<=x;i=pos+1)
    {
        pos=(x/(x/i));
        ans=(ans+Get_Sum(i,pos)*g[x/i]%Mod)%Mod;
    }   
    printf("%lld\n",ans);
}

int main()
{
    n=read(),m=read();
    Pre();
    for(int i=1;i<=n;i++)
    {
        static int a,b,k;ll x;
        a=read(),b=read(),x=read(),k=read();x%=Mod;
        int GCD=gcd(a,b);
        f[GCD]=x*Inv[a/GCD]%Mod*Inv[b/GCD]%Mod;
        Update(GCD);
        Solve(k);
    }
    return 0;
}