【BZOJ4815】小Q的表格（CQOI2017）-数论+分块

测试地址：小Q的表格
做法：本题需要用到数论+分块。
这道题可真是个神题……之所以标题笼统地写了数论，是因为它用到的知识实在是有点多……这里给大家，也算给自己理一遍思路，以免忘了。
首先分析题目所给二元函数的条件：
$f(a,b)=f(b,a)$
$b\cdot f(a,a+b)=(a+b)\cdot f(a,b)$
将第二个式子中的$b$换成$b-a$，得到：$(b-a)\cdot f(a,b)=b\cdot f(a,b-a)$
变换一下形式，得到：$f(a,b)=\frac{b}{b-a}\cdot f(a,b-a)$
这个东西有什么用呢？我们继续得到：$f(a,b-a)=\frac{b-a}{b-2a}\cdot f(a,b-2a)$
于是有：$f(a,b)=\frac{b}{b-2a}\cdot f(a,b-2a)$
观察这个式子，容易用数学归纳法得出$f(a,b)=\frac{b}{b \mod a}\cdot f(a,b\mod a)$的结论。注意由于$a=0$或$b=0$时$f(a,b)$没有意义，这里若$b\mod a=0$，我们就把它的值当成$a$，下面都用这样的方式进行表示。
又由于$f(a,b)=f(b,a)$，我们可以把$a$和$b\mod a$换过来再进行一遍这样的过程。我们不断这样换来换去，发现两个数的变动和辗转相除法中两个数的变化非常相似，那么令$d=\gcd(a,b)$，最后我们可以得到$f(a,b)=x\cdot f(d,d)$，因此我们知道，只要确定了$f(d,d)$，整个表格就随之确定了。于是现在的问题是，系数$x$是个什么东西。
我们知道$f(a,b)=\frac{b}{b \mod a}\cdot f(a,b\mod a)$，把右边$f$中的两项换过来再求一遍，得到：
$f(a,b)=\frac{b}{b\mod a}\cdot \frac{a}{a\mod(b\mod a)}\cdot f(b\mod a,a\mod(b\mod a))$
如果再往下算一步，系数会乘上一个分子为$b\mod a$的分数，和$\frac{b}{b\mod a}$的分母相互抵消，那么我们就得到了一个新的分母，太长就不写出来了，但一定和$f$中的一项相同。于是我们发现，随着等式右边两项变化成$x,y$，系数会同时变化成$\frac{a}{x}\cdot \frac{b}{y}$，因此当最后$x=y=d$时，有$f(a,b)=\frac{a}{d}\cdot \frac{b}{d}\cdot f(d,d)$，这样我们就求出了系数。
推导了上面一大堆东西之后，修改操作就很显然了，把所有修改都变成针对$f(d,d)$的修改即可。现在主要是询问操作，我们要求：
$ans=\sum_{i=1}^n\sum_{j=1}^nf(i,j)$
显然把上面算出的式子带进去，把$d$提出来得：
$ans=\sum_{d=1}^nf(d,d)\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor}[\gcd(i,j)=1]ij$
令$g(n)=\sum_{i=1}^{n}\sum_{j=1}^{n}[\gcd(i,j)=1]ij$，显然式子就变成了一个数论分块的形式，那么现在有两个问题要处理：
首先我们要解决快速求$f(d,d)$前缀和的问题。因为询问一共有$m\sqrt n$次，我们不能承受$O(\log n)$的询问复杂度，怎么办呢？注意到修改只有$m$次，所以我们分块，对每个点维护块内前缀和，再对每个块维护块的前缀和，即可做到修改$O(\sqrt n)$，询问$O(1)$的复杂度。
然后我们要解决预处理$g(n)$的问题。我们发现直接求这个东西非常不好求，怎么想都是$O(\sqrt n)$的复杂度，于是我们考虑递推，即考虑$g(n)$与$g(n-1)$相比多了什么。
显然，它多出的贡献是$2n\sum_{i=1}^n[\gcd(i,n)=1]i$。注意当$n=1$时不能这样求，因为对于$n>1$，$\gcd(n,n)=n\ne 1$，不会有重复的问题，但当$n=1$时就会有算重，特判一下即可。
对于$\sum_{i=1}^n[\gcd(i,n)=1]i$这个式子，它是在算在$n$之内与$n$互质的所有数之和。注意到这样的数总是成对出现的，如果$x$与$n$互质，那么$n-x$也与$n$互质，所以我们有：
$\sum_{i=1}^n[\gcd(i,n)=1]i=\frac{\varphi(n)}{2}n$
虽然当$n=1$时显然不成立，但它带进原式中时恰好把原来的特判给省掉了，于是我们找到了递推式：
$g(n)=g(n-1)+n^2\varphi(n)$
注意到$n^2\varphi(n)$是积性函数，线性筛预处理即可，时间复杂度为$O(n)$。
经过了上面一系列漫长的讨论，我们终于完全地解决了这个问题，时间复杂度为$O(m\sqrt n)$。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1000000007;
int m;
ll n,blocksiz,prime[2000010],g[4000010];
ll val[4000010],sum[4000010],blocksum[10010];
bool vis[4000010]={0};

void calc_g()
{
    prime[0]=0;
    g[0]=0,g[1]=1;
    for(ll i=2;i<=n;i++)
    {
        if (!vis[i])
        {
            prime[++prime[0]]=i;
            g[i]=i*i%mod*(i-1)%mod;
        }
        for(int j=1;j<=prime[0]&&i*prime[j]<=n;j++)
        {
            vis[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                g[i*prime[j]]=g[i]*prime[j]%mod*prime[j]%mod*prime[j]%mod;
                break;
            }
            g[i*prime[j]]=g[i]*prime[j]%mod*prime[j]%mod*(prime[j]-1)%mod;
        }
    }
    for(ll i=1;i<=n;i++)
        g[i]=(g[i]+g[i-1])%mod;
}

void init()
{
    scanf("%d%lld",&m,&n);

    blocksiz=sqrt(n);
    blocksum[0]=0;
    for(ll i=1,nowblock=0;i<=n;i++)
    {
        val[i]=i*i%mod;
        if (i%blocksiz) sum[i]=(sum[i-1]+val[i])%mod;
        else sum[i]=val[i];
        if (i%blocksiz==blocksiz-1)
            blocksum[(i+1)/blocksiz]=(blocksum[i/blocksiz]+sum[i])%mod;
    }
    calc_g();
}

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

void modify(ll pos,ll x)
{
    ll add=x-val[pos];
    val[pos]=x;
    for(int i=pos;;i++)
    {
        sum[i]+=add;
        if (i%blocksiz==blocksiz-1) break;
    }
    for(int i=pos/blocksiz+1;i<=n/blocksiz+1;i++)
        blocksum[i]+=add;
}

ll query(ll l,ll r)
{
    l--;
    ll ans=0;
    if (l>0) ans=ans-blocksum[l/blocksiz]-sum[l];
    ans=ans+blocksum[r/blocksiz]+sum[r];
    ans=(ans%mod+mod)%mod;
    return ans;
}

void work()
{
    ll a,b,k,x;
    for(int i=1;i<=m;i++)
    {
        scanf("%lld%lld%lld%lld",&a,&b,&x,&k);

        ll d=gcd(a,b);
        x/=a/d*b/d;
        modify(d,x);

        ll ans=0;
        for(ll i=k;i>=1;i=k/(k/i+1))
        {
            ll l=k/(k/i+1)+1,r=i;
            ans=(ans+query(l,r)*g[k/i])%mod;
        }
        printf("%lld\n",ans);
    }
}

int main()
{
    init();
    work();

    return 0;
}