重点 (1)
[frac{f(a,a+b)}{a(a+b)}=frac{f(a,b)}{ab}=frac{f(a,b-a)}{a(b-a)}=frac{f(a,b~mod~a)}{a(b~mod~a)}=frac{f(d,d)}{d^2}
]
其中 (d=gcd(a,b+b))
那么
[ans=sum_{i=1}^{k}sum_{j=1}^{k}f(i,j)=sum_{i=1}^{k}sum_{j=1}^{k}sum_{d=1}^{k}[gcd(i,j)==d]frac{ijf(d,d)}{d^2}
]
最终可以得到
[ans=sum_{d=1}^{k}f(d,d)sum_{i=1}^{lfloorfrac{k}{d}
floor}sum_{i=1}^{lfloorfrac{k}{d}
floor}[iperp j]ij
]
重点 (2)
考虑求
[sum_{i=1}^{x}sum_{j=1}^{x}[iperp j]ij=2sum_{i=1}^{x}isum_{j=1}^{i}[iperp j]j-1
]
除了 (i=1) 以外，小于等于 (i) 的与 (i) 互质的数成对存在，那么这些数字的和就是 (frac{varphi(i)i}{2})
所以上面变成
[sum_{i=1}^{x}varphi(i)i^2
]
即
[ans=sum_{d=1}^{k}f(d,d)sum_{i=1}^{lfloorfrac{k}{d}
floor}varphi(i)i^2
]
筛出 (sum_{i=1}^{x}varphi(i)i^2)
维护 (f(d,d)) 的前缀和即可
分块即可暴力树状数组就过了
# include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod(1e9 + 7);
const int maxn(4e6 + 5);
inline int Pow(ll x, int y) {
register ll ret = 1;
for (; y; y >>= 1, x = x * x % mod)
if (y & 1) ret = ret * x % mod;
return ret;
}
inline void Inc(int &x, int y) {
x = x + y >= mod ? x + y - mod : x + y;
}
int n, m, f[maxn], pr[maxn / 10], tot, sum[maxn], ans, lst[maxn];
bitset <maxn> ispr;
inline int Gcd(int a, int b) {
return !b ? a : Gcd(b, a % b);
}
inline void Add(int x, int v) {
for (; x <= n; x += x & -x) Inc(sum[x], v);
}
inline int Query(int x) {
register int ret = 0;
for (; x; x ^= x & -x) Inc(ret, sum[x]);
return ret;
}
int main() {
register int i, j, a, b, k, cur, p, d;
register ll x;
scanf("%d%d", &m, &n), f[1] = 1, ispr[1] = 1;
for (i = 2; i <= n; ++i) {
if (!ispr[i]) pr[++tot] = i, f[i] = i - 1;
for (j = 1; j <= tot && i * pr[j] <= n; ++j) {
ispr[i * pr[j]] = 1;
if (i % pr[j]) f[i * pr[j]] = f[i] * (pr[j] - 1);
else {
f[i * pr[j]] = f[i] * pr[j];
break;
}
}
}
for (i = 1; i <= n; ++i) {
lst[i] = (ll)i * i % mod, Add(i, lst[i]);
f[i] = (ll)lst[i] * f[i] % mod, Inc(f[i], f[i - 1]);
}
while (m) {
ans = 0, --m, scanf("%d%d%lld%d", &a, &b, &x, &k), x %= mod;
d = Gcd(a, b), Add(d, mod - lst[d]);
lst[d] = (ll)d * d % mod * x % mod * Pow((ll)a * b % mod, mod - 2) % mod;
Add(d, lst[d]);
for (p = 0, i = 1, j; i <= k; i = j + 1) {
j = k / (k / i), cur = Query(j);
Inc(ans, (ll)(cur - p + mod) * f[k / i] % mod), p = cur;
}
printf("%d
", ans);
}
return 0;
}