【LOJ】#3020. 「CQOI2017」小 Q 的表格

#3020. 「CQOI2017」小 Q 的表格

这个的话求出来(g = gcd(a,b))

会修改所有gcd为g的位置

我们要求((g,g))这个位置的数一定是(g^{2})的倍数

之后的(gcd(a,b) == g)的位置也会满足

$frac{f(g,g)}{g^{2}} = frac{f(a,b)}{ab} $

注意(frac{f(g,g)}{g^{2}})，(frac{f(a,b)}{ab})都不一定是整数，但是(f(g,g))和(f(a,b))是

这样的话我们求出(f(g,g))的值，很容易得到

(ans = sum_{d = 1}^{k} f(d,d) sum_{i = 1}^{lfloor frac{k}{d} floor} sum_{j = 1}^{lfloor frac{k}{d} floor} [gcd(i,j) == 1] i cdot j)

这个时候不要用(mu) ！不要用(mu)！不要用(mu)！

我们这么考虑，就是枚举一个(i)，然后乘上(i)以内和(i)互质的数的和，由于顺序可以交换，所以要乘上1

(i)以内和(i)互质的数的公式是(frac{i imes varphi(i)}{2})

(1)的话除外，不过由于外面乘了一个2，所以这个式子在这里对1成立

所以

(ans = sum_{d = 1}^{k} f(d,d) sum_{i = 1}^{lfloor frac{k}{d} floor} i^{2}varphi(i))

后面的可以预处理，前面的话要更新前缀和，用分块实现(O(sqrt{N}))维护和(O(1))查询即可

#include <bits/stdc++.h>
#define fi first
#define se second
#define pii pair<int,int>
#define mp make_pair
#define pb push_back
#define space putchar(' ')
#define enter putchar('
')
#define eps 1e-10
#define ba 47
#define MAXN 4000005
//#define ivorysi
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef double db;
template<class T>
void read(T &res) {
    res = 0;T f = 1;char c = getchar();
    while(c < '0' || c > '9') {
	if(c == '-') f = -1;
	c = getchar();
    }
    while(c >= '0' && c <= '9') {
	res = res * 10 +c - '0';
	c = getchar();
    }
    res *= f;
}
template<class T>
void out(T x) {
    if(x < 0) {x = -x;putchar('-');}
    if(x >= 10) {
	out(x / 10);
    }
    putchar('0' + x % 10);
}
const int MOD = 1000000007;
int M,N;
int f[MAXN],prime[MAXN],tot,phi[MAXN];
int h[MAXN],bl[2005],br[2005],id[MAXN],cnt,lz[2005];
int sum[MAXN];
bool nonprime[MAXN];
int inc(int a,int b) {
    return a + b >= MOD ? a + b - MOD : a + b;
}
int mul(int a,int b) {
    return 1LL * a * b % MOD;
}
void update(int &x,int y) {
    x = inc(x,y);
}
int gcd(int a,int b) {
    return b == 0 ? a : gcd(b,a % b);
}
int main(){
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    read(M);read(N);
    phi[1] = 1;
    h[1] = 1;
    for(int i = 2 ; i <= N ; ++i) {
	if(!nonprime[i]) {
	    prime[++tot] = i;
	    phi[i] = i - 1;
	}
	for(int j = 1 ; j <= tot ; ++j) {
	    if(prime[j] > N / i) break;
	    nonprime[i * prime[j]] = 1;
	    if(i % prime[j] == 0) {
		phi[i * prime[j]] = phi[i] * prime[j];
		break;
	    }
	    else {
		phi[i * prime[j]] = phi[i] * phi[prime[j]];
	    }
	}
	h[i] = mul(phi[i],mul(i,i));
	update(h[i],h[i - 1]);
    }
    for(int i = 1 ; i <= N ; ++i) {
	sum[i] = mul(i,i);
	f[i] = mul(i,i);
	update(sum[i],sum[i - 1]);
    }
    int S = sqrt(N);
    for(int i = 1 ; i <= N ; ++i) {
	int r = min(N,i + S - 1);
	++cnt;
	bl[cnt] = i;br[cnt] = r;i = r;
	for(int j = bl[cnt] ; j <= br[cnt] ; ++j) id[j] = cnt;
    }
    int a,b,k;
    int64 x;
    for(int i = 1 ; i <= M ; ++i) {
	read(a);read(b);read(x);read(k);
	int g = gcd(a,b);
	int d = (x / (1LL * a / g * b / g)) % MOD;
	int c = inc(d,MOD - f[g]);
	f[g] = d;
	for(int j = g ; j <= br[id[g]] ; ++j) {
	    update(sum[j],c);
	    
	}
	for(int j = id[g] + 1 ; j <= cnt ; ++j) update(lz[j],c);
	int res = 0;
	for(int j = 1 ; j <= k ; ++j) {
	    int r = k / (k / j);
	    int s = inc(sum[r],lz[id[r]]);
	    s = inc(s,MOD - inc(sum[j - 1],lz[id[j - 1]]));
	    update(res,mul(s,h[k / j]));
	    j = r;
	}
	out(res);enter;
    }
    return 0;
}