[CQOI2017] 小Q的表格

[题目链接]

[思路要点]

[egin{align} ans &= sum_{i=1}^{k}sum_{j=1}^kf(i,j)\ &= sum_{d=1}^kf(d)cdot sum_{i=1}^{left[frac k d ight]} sum_{j=1}^{left[frac k d ight]}[(i,j)==1]cdot ij \ end{align} ]

令 $$ g(lim, x)=sum_{i=1}^{lim} sum_{j=1}^{lim}[(i,j)==x]cdot ij$$，则 (ans=sum_{d=1}^kf(d)cdot gleft(left[frac k d ight],1 ight))

不妨设 (lim) 已确定，并且此时 (g[i]=g(lim,i))，考虑 (g[i]) 怎么求

运用莫比乌斯反演，设

[G[i] = sum_{i|d} g[d] ]

则

[g[i] = sum _{i|d} muleft(left[frac d i ight] ight) G(d) ]

而

[egin{align} G[x] &= sum _{x|d} g[d]\ &= sum_{x|d} sum_{i=1}^{lim}sum_{j=1}^{lim}left[left(i,j ight)==d ight]cdot d^2 \ &= sum_{i=1}^{lim}sum_{j=1}^{lim}left[x mid left(i,j ight) ight]cdot x^2 \ &= x^2cdotsum_{i=1}^{left[frac {lim} x ight]} sum_{j=1}^{left[frac {lim} x ight]} ij end{align} ]

因此我们设 (S[n]=sum_{i=1}^{n} sum_{j=1}^{n} ij=(frac {n(n+1)} 2)^2)

则有 (G[x]=x^2 cdot Sleft[frac {lim} x ight])

那么 (g[i]= sum_{imid d}muleft(left[frac d i ight] ight)cdot d^2 cdot Sleft[frac {lim} d ight])

[g[1] = sum_{d = 1} ^ {lim} mu (d) cdot d^2 cdot Sleft[frac {lim} d ight] ]

此时令 (g(lim)=g[1])，则

[egin{align} ans &= sum_{d=1}^k f(d) cdot gleft(left[frac k d ight],1 ight)\ &=sum_{d=1}^k f(d) cdot gleft(left[frac k d ight] ight) end{align} ]

此时问题基本得到解决，我们现在需要在 (Theta (sqrt k)) 的时间内解决上述式子的求值

运用数论分块，只要维护 (f) 的前缀和以及 (g) 数组的值即可

对于每次操作一共询问 (O(sqrt k)) 次，修改 (1) 次，大致需要 (O(sqrt k)) 修改 (O(1)) 查询的数据结构，可以使用分块维护，分别维护块内前缀和以及块尾维护从 (1) 到这里的总和。

对于 (g) 的求值，可以使用差分

[egin{align} g(n) &= sum_{i=1} ^n mu(i) cdot i^2 cdot Sleft[frac n i ight]\ g(n - 1) &=sum_{i=1} ^{n-1} mu(i) cdot i^2 cdot Sleft[frac {n-1} i ight]\ g(n)-g(n-1) &= sum_{i=1}^{n-1} mu(i) cdot i^2 cdot left(Sleft[frac n i ight] - Sleft[frac{n-1} i ight] ight)+mu(n)cdot n^2 cdot S[1] end{align} ]

对于后面的 (mu(n)cdot n^2 cdot S[1]) 可以直接加上，前面的式子考虑 (left(Sleft[frac n i ight] - Sleft[frac{n-1} i ight] ight)) 什么时候不为 (0)

发现 (left(Sleft[frac n i ight] - Sleft[frac{n-1} i ight] ight) e 0) 当且仅当 (i mid n)，并且此时 (left[frac n i ight] - left[frac {n - 1} i ight]=1)

而由于 (S[i]=left(frac {i(i+1)} 2 ight)^2)

故

[egin{align} S[x] - S[x - 1] &= frac {x^2cdot (x+1)^2}4 - frac {(x-1)^2cdot x^2} 4\ &= x^3 end{align} ]

因此

[egin{align} g(n)-g(n-1) &= sum_{imid n} mu(i) cdot i^2 cdot left(Sleft[frac n i ight] - Sleft[frac{n-1} i ight] ight)+mu(n)cdot n^2 cdot S[1]\ &= sum _{i mid n} mu(i) cdot i^2 cdot left(left[frac n i ight] ight)^3 + mu(n) cdot n^2 cdot S[1]\ &= n^2 sum_{imid n} mu(i) cdot left[frac n i ight]\ &= n^2 varphi(n) + mu(n) cdot n^2 cdot S[1] end{align} ]

所以 (g) 数组可以线性递推解出

因此总复杂度 (Theta (n+msqrt n))

[代码]

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
const int mod=1e9+7;
const int N=5e6;
int a[N],s[N],ss[N],prime[N],p[N],phi[N],cnt,f[N],bl[N],blo,n,m;
void make(int n)
{
    phi[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!p[i])prime[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++)
        {
            int x=i*prime[j];p[x]=1;
            if(i%prime[j])phi[x]=phi[i]*(prime[j]-1);
            else{phi[x]=phi[i]*prime[j];break;}
        }
    }
    for(int i=1;i<=n;i++)f[i]=(f[i-1]+1ll*i*i%mod*phi[i])%mod;

    blo=sqrt(n);
    for(int i=1;i<=n;i++)bl[i]=(i-1)/blo+1,a[i]=1ll*i*i%mod;
    for(int i=1;i<=bl[n];i++)
    {
        int l=(i-1)*blo+1,r=min(i*blo,n);s[l]=a[l];
        for(int j=l+1;j<=r;j++)s[j]=(s[j-1]+a[j])%mod;
        ss[i]=(ss[i-1]+s[r])%mod;
    }
}
int query(int n){return (s[n]+ss[bl[n]-1])%mod;}
void update(int x,int w)
{
    int blx=bl[x],r=min(blx*blo,n);
    int tw=(w-a[x]+mod)%mod;a[x]=w;
    for(int i=x;i<=r;i++)s[i]=(s[i]+tw)%mod;
    for(int i=blx;i<=bl[n];i++)ss[i]=(ss[i]+tw)%mod;
}
int solve(int n)
{
    int i=1,lt=sqrt(n);
    int ans=0;

    for(;i<=lt;i++)ans=(ans+1ll*a[i]*f[n/i])%mod;
    for(;i<=n;i=lt+1)
    {
        lt=n/(n/i);
        ans=(ans+1ll*(query(lt)-query(i-1))*f[n/i])%mod;
    }
    return (ans+mod)%mod;
}
int gcd(int x,int y){return y?gcd(y,x%y):x;}
int main()
{
    scanf("%d%d",&m,&n);
    make(n);
    for(int i=1;i<=m;i++)
    {
        int u,v,k;long long x;
        scanf("%d%d%lld%d",&u,&v,&x,&k);
        int d=gcd(u,v);
        update(gcd(u,v),(x/(u/d)/(v/d))%mod);
        printf("%d
",solve(k));
    }
}