4815: [Cqoi2017]小Q的表格

4815: [Cqoi2017]小Q的表格

Time Limit: 20 Sec Memory Limit: 512 MB
Submit: 496 Solved: 209
[Submit][Status][Discuss]
Description

小Q是个程序员。
作为一个年轻的程序员，小Q总是被老C欺负，老C经常把一些麻烦的任务交给小Q来处理。每当小Q不知道如何解决
时，就只好向你求助。为了完成任务，小Q需要列一个表格，表格有无穷多行，无穷多列，行和列都从1开始标号。
为了完成任务，表格里面每个格子都填了一个整数，为了方便描述，小Q把第a行第b列的整数记为f(a,b),为了完成
任务，这个表格要满足一些条件：(1)对任意的正整数a,b,都要满足f(a,b)=f(b,a)；(2)对任意的正整数a,b,都要
满足b×f(a,a+b)=(a+b)*f(a,b)。为了完成任务，一开始表格里面的数很有规律，第a行第b列的数恰好等于a*b，
显然一开始是满足上述两个条件的。为了完成任务，小Q需要不断的修改表格里面的数，每当修改了一个格子的数
之后，为了让表格继续满足上述两个条件，小Q还需要把这次修改能够波及到的全部格子里都改为恰当的数。由于
某种神奇的力量驱使，已经确保了每一轮修改之后所有格子里的数仍然都是整数。为了完成任务，小Q还需要随时
获取前k行前k列这个有限区域内所有数的和是多少，答案可能比较大，只需要算出答案mod1,000,000,007之后的结
果。
Input

第1行包含两个整数m,n,表示共有m次操作，所有操作和查询涉及到的行编号和列编号都不超过n。
接下来m行，每行4个整数a,b,x,k表示把第a行b列的数改成x，然后把它能够波及到的所有格子全部修改，保证修改
之后所有格子的数仍然都是整数，修改完成后计算前k行前k列里所有数的和。
1<=m<=10000,1<=a,b,k<=N<=4*10^6,0<=x<=10^18
Output

输出共m行，每次操作之后输出1行，表示答案mod 1,000,000,007之后的结果。

Sample Input

3 3

1 1 1 2

2 2 4 3

1 2 4 2
Sample Output

9

36

14

一开始,表格的前3行前3列如图中左边所示,前2次操

作后表格没有变化,第3次操作之后的表格如右边所示。

1 2 3 2 4 6

2 4 6 4 4 12

3 6 9 6 12 9

HINT

Source

对于$f(a,b) = f(b,a)$，显然$gcd(a,b) = gcd(b,a)$
对于$b\ast f(a,a+b) = (a + b)\ast f(a,b)$由更相减损法，得到$gcd(a,a+b) = gcd(a,b)$
因此，每次修改等价于对所有坐标的最大公约数为一个定值的点修改
对于任意两个坐标的最大公约数相同的点
由约束条件得到，它们的比例肯定是一个定值
因此，任意的$f(a,b)$在修改过后总是能写作$a\ast b \ast k$的形式
重新定义$f(d)$为最大公约数为$d$的点的这个系数
那么$Ans = \sum_{i=1}^k\sum_{j=1}^k i\ast j\ast f(gcd(i,j))$
大力莫比乌斯反演，得到
$Ans = \sum_{d=1}^n f(d)\ast d^2 \sum_{t=1}^{ \lfloor \frac nd \rfloor}\mu(t)\ast t^2 \left(\sum_{i=1}^{\lfloor \frac n{dt} \rfloor} i\right)^2$
定义$g(n) = \sum_{i=1}^n \mu(i)\ast i^2 \left(\sum_{j=1}^{\lfloor \frac ni \rfloor}j\right)^2$
那么$Ans = \sum_{d=1}^n f(d)\ast d^2 \ast g\left(\left \lfloor \frac nd \right \rfloor \right)$
假如能预处理好所有$g(n)$，那么上述公式能在$O(\sqrt n)$计算出来
注意到$\left \lfloor \frac ni \right \rfloor - \left \lfloor \frac {n-1}i \right \rfloor = 1$成立当且仅当$i|n$
那么$g(n) - g(n - 1) = \sum_{d|n} \mu(d) \ast d^2 \ast \left(\frac nd\right)^3 = n^3\sum_{d|n}\frac {\mu(d)} d$
这个式子右边那个东西显然是个积性函数，线性筛解决即可
因此能在$O(n)$递推出所有的$g(n)$
对于修改$f(n)$时，因为还需要时刻查询某段的和
若使用树状数组维护前缀和然后差分，复杂度就多个$log$无法接受
可以牺牲一下修改时间，分块维护$f$的前缀和，这样查询就能$O(1)$了
综上，总复杂度为$O(m\sqrt n)$

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
//#include<ctime>
using namespace std;

const int maxm = 2005;
const int maxn = 4E6 + 4;
typedef long long LL;
const LL mo = 1000000007;

int n,m,Sqrt,S,tot,pri[maxn],bel[maxn],pos[maxn],g[maxn],h[maxn];
bool not_pri[maxn];
LL s1[maxm],s2[maxm][maxm],f[maxn];

inline int gcd(int x,int y) {return !y ? x : gcd(y,x % y);}
inline int Mul(const LL &x,const LL &y) {return x * y % mo;}
inline int Add(const int &x,const int &y) {return x + y < mo ? x + y : x + y - mo;}
inline int Dec(const int &x,const int &y) {return x - y >= 0 ? x - y : x - y + mo;}

inline int ksm(int x,int y)
{
    int ret = 1;
    for (; y; y >>= 1)
    {
        if (y & 1) ret = Mul(ret,x);
        x = Mul(x,x);
    }
    return ret;
}

inline int getint()
{
    char ch = getchar(); int ret = 0;
    while (ch < '0' || '9' < ch) ch = getchar();
    while ('0' <= ch && ch <= '9')
        ret = ret * 10 + ch - '0',ch = getchar();
    return ret;
}

inline LL getLL()
{
    char ch = getchar(); LL ret = 0;
    while (ch < '0' || '9' < ch) ch = getchar();
    while ('0' <= ch && ch <= '9')
        ret = ret * 10LL + 1LL * (ch - '0'),ch = getchar();
    return ret;
}

void Pre_Work()
{
    g[1] = h[1] = f[1] = 1;
    Sqrt = sqrt(n); S = n / Sqrt + 1;
    bel[1] = 1 / Sqrt; pos[1] = 1 % Sqrt;
    s1[bel[1]] += 1LL; s2[bel[1]][pos[1]] = 1LL;

    for (int i = 2; i <= n; i++)
    {
        if (!not_pri[i])
            pri[++tot] = i,h[i] = Dec(1,ksm(i,mo - 2));
        for (int j = 1; j <= tot; j++)
        {
            int Nex = i * pri[j];
            if (Nex > n) break;
            not_pri[Nex] = 1;
            if (i % pri[j] == 0)
            {
                h[Nex] = h[i]; break;
            }
            h[Nex] = Mul(h[i],h[pri[j]]);
        }
        int tmp = Mul(i,i);
        g[i] = Add(g[i - 1],Mul(Mul(tmp,i),h[i]));
        if (pos[i - 1] == Sqrt - 1)
            pos[i] = 0,bel[i] = bel[i - 1] + 1;
        else pos[i] = pos[i - 1] + 1,bel[i] = bel[i - 1];
        f[i] = 1; s1[bel[i]] += 1LL * tmp; s2[bel[i]][pos[i]] = tmp;
    }

    for (int i = 0; i <= S; i++)
        for (int j = 1; j < Sqrt; j++)
            s2[i][j] += s2[i][j - 1];
    for (int i = 1; i <= S; i++) s1[i] += s1[i - 1];
}

inline void Modify(int k,LL d)
{
    d = Mul(d,Mul(k,k));
    for (int i = bel[k]; i <= S; i++) s1[i] += d;
    for (int i = pos[k]; i < Sqrt; i++) s2[bel[k]][i] += d;
}

inline LL Query(const int &k)
{
    LL ret = bel[k] ? s1[bel[k] - 1] : 0;
    return ret + s2[bel[k]][pos[k]];
}

inline int Mod(LL x) {x %= mo; return x >= 0 ? x : x + mo;}

inline void Solve()
{
    int a,b,x,N,G,Ans = 0;
    a = getint(); b = getint();
    x = getLL() % mo; N = getint();
    G = gcd(a,b); x = Mul(x,ksm(Mul(a,b),mo - 2));
    Modify(G,x - f[G]); f[G] = x;

    for (int i = 1,last; i <= N; i = last + 1)
    {
        int tmp = N / i; last = N / tmp;
        int A = Mod(Query(last) - Query(i - 1));
        Ans = Add(Ans,Mul(A,g[tmp]));
    }
    printf("%d\n",Ans);
}

int main()
{
    #ifdef DMC
        freopen("DMC.txt","r",stdin);
        freopen("test.txt","w",stdout);
    #endif

    m = getint(); n = getint();
    Pre_Work();
    //cerr << (double)(clock()) / CLOCKS_PER_SEC << endl;
    while (m--) Solve();
    //cerr << (double)(clock()) / CLOCKS_PER_SEC << endl;
    return 0;
}