题目描述
小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 imes (a + b)}=frac{f(a,b)}{a imes b})
然后发现这玩意儿跟辗转相减法很类似啊
那(frac{f(a,b)}{a imes b})就应该(=frac{f(a,a \% b)}{a imes (a \% b)})
然后由于(f(a,b)=f(b,a)),所以再把(a)和(a\%b)换过来继续辗转相除
最后就是(frac{f(a,b)}{a imes 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)的上界相同
所以可以化成(2sum_{i=1}^{k/d}isum_{j=1}^{i}j[gcd(i,j)=1])
有一个东西就是(sum_{i=1}^{n}{i[gcd(i,n)=1]}=frac{phi(n) imes n}{2})
证明的话也比较显然,就是对于一个与(n)互质的数(x),那么一定有(n-x)与(x)互质
然后把他们两两配对,正好有(frac{phi(n)}{2})对
那么再看原式就化成了(sum_{i=1}^{k/d}phi(i) imes i^2)
这个东西可以(O(n))预处理出来
所以(ans=sum_{d=1}^{k}f(d)sum_{i=1}^{k/d}phi(i) imes 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
",(ans % mod + mod) % mod) ;
}
return 0 ;
}