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

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

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

    return 0;
}