[CQOI2017]小Q的表格

题目描述

小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之后的结果。

输入输出格式

输入格式：

输入文件第1行包含两个整数m,n，表示共有m次操作，所有操作和查询涉及到的行编号和列编号都不超过n。

接下来m行，每行4个整数a,b,x,k，表示把第a行b列的数改成x，然后把它能够波及到的所有格子全部修改，保证修改之后所有格子的数仍然都是整数，修改完成后计算前k行前k列里所有数的和。

输出格式：

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

输入输出样例

输入样例#1：

3 3
1 1 1 2
2 2 4 3
1 2 4 2

输出样例#1：

9
36
14

输入样例#2:

4 125
1 2 4 8
1 3 9 27
1 4 16 64
1 5 25 125

输出样例#2：

2073
316642
12157159
213336861

---

题解

首先这个表格是对称的
然后可以发现ta给的一个奇怪的限制:
\(b×f(a,a+b)=(a+b)×f(a,b)\)
那么根据ta给的式子可以化成\(\frac{f(a , a + b)}{a + b}=\frac{f(a , b)}{b}\)
发现似乎两边同时乘上\(\frac{1}{a}\)很舒服
\(\frac{f(a , a + b)}{a \times (a + b)}=\frac{f(a,b)}{a\times b}\)
然后发现这玩意儿跟辗转相减法很类似啊
那\(\frac{f(a,b)}{a\times b}\)就应该\(=\frac{f(a,a \% b)}{a \times (a \% b)}\)
然后由于\(f(a,b)=f(b,a)\),所以再把\(a\)和\(a\%b\)换过来继续辗转相除
最后就是\(\frac{f(a,b)}{a\times b}=\frac{f(d , d)}{d^2}(d=gcd(a,b))\)
我们设\(f(d)=f(d,d)\)
那么题目就是求\(ans=\sum_{i=1}^{k}\sum_{j=1}^{k}f(gcd(i,j))\)
先简单化简一下式子\(ans=\sum_{d=1}^{k}\sum_{i=1}^{k}\sum_{j=1}^{k}f(d)[gcd(i,j)=d]\)
\(ans=\sum_{d=1}^{k}f(d)\sum_{i=1}^{k/d}\sum_{j=1}^{k/d}{ij[gcd(i,j)=1]}\)
然后化简一波后面的那个东西
\(\sum_{i=1}^{k/d}\sum_{j=1}^{k/d}ij[gcd(i,j)=1]\)
可以发现\(i,j\)的上界相同
所以可以化成\(2\sum_{i=1}^{k/d}i\sum_{j=1}^{i}j[gcd(i,j)=1]\)
有一个东西就是\(\sum_{i=1}^{n}{i[gcd(i,n)=1]}=\frac{\phi(n) \times n}{2}\)
证明的话也比较显然,就是对于一个与\(n\)互质的数\(x\),那么一定有\(n-x\)与\(x\)互质
然后把他们两两配对,正好有\(\frac{\phi(n)}{2}\)对
那么再看原式就化成了\(\sum_{i=1}^{k/d}\phi(i) \times i^2\)
这个东西可以\(O(n)\)预处理出来
所以\(ans=\sum_{d=1}^{k}f(d)\sum_{i=1}^{k/d}\phi(i) \times i^2\)
对于每次修改,为了平衡操作和修改时间,用分块来做到\(O(\sqrt n)\)修改,\(O(1)\)查询

代码

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
# define LL long long
const int M = 4000005 ;
const int N = 3005 ;
const int mod = 1e9 + 7 ;
using namespace std ;

inline LL read() {
    char c = getchar() ; LL x = 0 , w = 1 ;
    while(c>'9'||c<'0') { if(c=='-') w = -1 ; c = getchar() ; }
    while(c>='0'&&c<='9') { x = x*10+c-'0' ; c = getchar() ; }
    return x*w ;
} 

bool notp[M] ;
int n , pnum , p[M / 2] , phi[M] , psum[M] ;
int Tag[N] , bel[M] , f[M] , sz , cnt , lpos[N] , rpos[N] , sum[M] ;
int gcd(int a , int b) {
    if(b == 0) return a ;
    return gcd(b , a % b) ;
}
inline int Fpw(int Base , int k) {
    int temp = 1 ;
    while(k) {
        if(k & 1) temp = 1LL * temp * Base % mod ;
        Base = 1LL * Base * Base % mod ; k >>= 1 ;
    }
    return temp ;
}
inline int Inv(int x) {
    return Fpw(x , mod - 2) ;
}
inline void PreSolve(int n) {
    phi[1] = 1 ;
    for(int i = 2 ; i <= n ; i ++) {
        if(!notp[i]) {
            p[++pnum] = i ;
            phi[i] = i - 1 ;
        }
        for(int j = 1 ; 1LL * i * p[j] <= n ; j ++) {
            notp[i * p[j]] = true ;
            if(i % p[j] == 0) {
                phi[i * p[j]] = phi[i] * p[j] ;
                break ;
            }
            phi[i * p[j]] = phi[i] * phi[p[j]] ;
        }
    }
    for(int i = 1 ; i <= n ; i ++) {
        psum[i] = ((psum[i - 1] + 1LL * i * i % mod * phi[i] % mod) % mod + mod) % mod ;
        f[i] = 1LL * i * i % mod ;
        sum[i] = (sum[i - 1] + f[i]) % mod ;
    }
}

inline void Get_Unit() {
    sz = sqrt(n) ;
    int l = 1 , r = sz ;
    while(l <= n) {
        r = min(l + sz - 1 , n) ; 
        ++ cnt ; lpos[cnt] = l ; rpos[cnt] = r ;
        for(int i = l ; i <= r ; i ++) bel[i] = cnt ;
        l = r + 1 ;
    }
}
void Change(int l , int r , int k) {
    for(int i = lpos[bel[l]] ; i <= rpos[bel[l]] ; i ++) {
        sum[i] = (sum[i] + Tag[bel[l]]) % mod ;
        if(i >= l && i <= r) sum[i] = (sum[i] + k) % mod ;
    }
    Tag[bel[l]] = 0 ;
    for(int i = bel[l] + 1 ; i < bel[r] ; i ++) Tag[i] = (Tag[i] + k) % mod ;
    if(bel[l] == bel[r]) return ;
    for(int i = lpos[bel[r]] ; i <= rpos[bel[r]] ; i ++) {
        sum[i] = (sum[i] + Tag[bel[l]]) % mod ;
        if(i >= l && i <= r) sum[i] = (sum[i] + k) % mod ;
    }
    Tag[bel[r]] = 0 ;
}
inline int query(int l , int r) {
    return (((sum[r] + Tag[bel[r]]) % mod - (sum[l - 1] + Tag[bel[l - 1]]) % mod) % mod + mod) % mod ;
}
int main() {
    int Case = read() ; n = read() ;
    PreSolve(n) ; Get_Unit() ;
    int a , b , k , d , dlt ; LL x , ans ;
    while(Case --) {
        ans = 0 ;
        a = read() ; b = read() ; x = read() % mod ; k = read() ;
        d = gcd(a , b) ;
        dlt = ((1LL * x * d % mod * d % mod * Inv(1LL * a * b % mod) % mod - f[d]) % mod + mod) % mod ;
        f[d] = 1LL * x * d % mod * d % mod * Inv(1LL * a * b % mod) % mod ;
        Change(d , n , dlt) ;
        for(int l = 1 , r ; l <= k ; l = r + 1) {
            r = k / (k / l) ;
            ans = (ans + 1LL * query(l , r) * psum[k / l] % mod) % mod ;
        }
        printf("%lld\n",(ans % mod + mod) % mod) ;
    }
    return 0 ;
}