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

#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;}