Bzoj2154: Crash的数字表格
2018-01-09 17:49
169 查看
题意
\(求ans=\sum_{i=1}^{n}\sum_{j=1}^{n}lcm(i, j)\)n,m<=10^7
Sol
\(原式=\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{i*j}{gcd(i, j)}\)假设n < m,
\(则ans=\sum_{d=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor} d*i*j*[gcd(i,j)==1]\)
枚举d,看里面的
\(令x=\lfloor\frac{n}{d}\rfloor, y=\lfloor\frac{m}{d}\rfloor\)
\(设f(k)=\sum_{i=1}^{x}\sum_{i=1}^{y}i*j*[gcd(i, j)==1]\)
\(设g(i)=\sum_{i|d} f(d)\)即表示gcd是i及i的倍数的数对的乘积和
\(就是i^2*\frac{(\lfloor\frac{x}{i}\rfloor + 1)*\lfloor\frac{x}{i}\rfloor}{2}*\frac{(\lfloor\frac{y}{i}\rfloor + 1)*\lfloor\frac{y}{i}\rfloor}{2}\)
然后就可以莫比乌斯反演求出f数组,从而得到答案
但是暴力求显然跑不过10^7的点
所以可以用两个数论分块+前缀和优化,记得取模,会爆
代码
不用数论分块,暴力# include <bits/stdc++.h> # define RG register # define IL inline # define Fill(a, b) memset(a, b, sizeof(a)) using namespace std; typedef long long ll; const int _(1e7), MOD(20101009); IL ll Read(){ char c = '%'; ll x = 0, z = 1; for(; c > '9' || c < '0'; c = getchar()) if(c == '-') z = -1; for(; c >= '0' && c <= '9'; c = getchar()) x = x * 10 + c - '0'; return x * z; } int prime[_], num, mu[_]; bool isprime[_]; IL void Prepare(){ mu[1] = 1; for(RG int i = 2; i <= _; ++i){ if(!isprime[i]){ prime[++num] = i; mu[i] = -1; } for(RG int j = 1; j <= num && i * prime[j] <= _; ++j){ isprime[i * prime[j]] = 1; if(i % prime[j]) mu[i * prime[j]] = -mu[i]; else{ mu[i * prime[j]] = 0; break; } } } } IL ll Calc(RG ll n, RG ll m){ RG ll f = 0, g; for(RG ll i = 1; i <= n; ++i){ RG ll x = n / i, y = m / i; g = i * i % MOD * (x * (x + 1) >> 1) % MOD * (y * (y + 1) >> 1) % MOD; (f += 1LL * mu[i] * g % MOD) %= MOD; } return (f + MOD) % MOD; } int main(RG int argc, RG char *argv[]){ Prepare(); RG ll n = Read(), m = Read(), ans = 0; if(n > m) swap(n, m); for(RG ll d = 1; d <= n; ++d) (ans += d * Calc(n / d, m / d) % MOD) %= MOD; printf("%lld\n", ans); return 0; }
用数论分块
# include <bits/stdc++.h> # define RG register # define IL inline # define Fill(a, b) memset(a, b, sizeof(a)) using namespace std; typedef long long ll; const int _(1e7 + 10), MOD(20101009); IL ll Read(){ char c = '%'; ll x = 0, z = 1; for(; c > '9' || c < '0'; c = getchar()) if(c == '-') z = -1; for(; c >= '0' && c <= '9'; c = getchar()) x = x * 10 + c - '0'; return x * z; } int prime[_], num, mu[_], squ[_], s[_]; bool isprime[_]; IL void Prepare(){ mu[1] = 1; RG int maxn = 1e7; for(RG int i = 2; i <= 1e7; ++i){ if(!isprime[i]){ prime[++num] = i; mu[i] = -1; } for(RG int j = 1; j <= num && i * prime[j] <= maxn; ++j){ isprime[i * prime[j]] = 1; if(i % prime[j]) mu[i * prime[j]] = -mu[i]; else{ mu[i * prime[j]] = 0; break; } } } } IL ll Calc(RG ll n, RG ll m){ RG ll f = 0, g, j; for(RG ll i = 1; i <= n; i = j + 1){ RG ll x = n / i, y = m / i; j = min(n / (n / i), m / (m / i)); g = (x * (x + 1) >> 1) % MOD * ((y * (y + 1) >> 1) % MOD) % MOD; (f += 1LL * (squ[j] - squ[i - 1]) % MOD * g % MOD) %= MOD; } return (f + MOD) % MOD; } int main(RG int argc, RG char *argv[]){ Prepare(); RG ll n = Read(), m = Read(), ans = 0, j; if(n > m) swap(n, m); for(RG int i = 1; i <= n; ++i) squ[i] = ((squ[i - 1] + 1LL * mu[i] * i * i % MOD) % MOD + MOD) % MOD, s[i] = (s[i - 1] + i) % MOD; for(RG ll d = 1; d <= n; d = j + 1){ j = min(n / (n / d), m / (m / d)); (ans += 1LL * ((s[j] - s[d - 1]) % MOD + MOD) % MOD * Calc(n / d, m / d) % MOD) %= MOD; } printf("%lld\n", ans); return 0; }
相关文章推荐
- bzoj2154 Crash的数字表格(莫比乌斯反演)
- bzoj2154 Crash的数字表格
- 莫比乌斯反演套路三、四--BZOJ2154: Crash的数字表格 && BZOJ2693: jzptab
- bzoj2154 Crash的数字表格
- bzoj2154 Crash的数字表格(反演)
- BZOJ2154 Crash的数字表格
- bzoj2154: Crash的数字表格
- BZOJ2154: Crash的数字表格
- bzoj2154 Crash的数字表格
- bzoj2154: Crash的数字表格
- bzoj2154: Crash的数字表格
- Bzoj2154: Crash的数字表格
- BZOJ2154 Crash的数字表格
- bzoj2154 Crash的数字表格
- bzoj2154: Crash的数字表格
- BZOJ2154: Crash的数字表格
- [BZOJ2154] Crash的数字表格
- 【BZOJ2154】Crash的数字表格 莫比乌斯反演
- 2693: jzptab/2154: Crash的数字表格 莫比乌斯反演
- Bzoj2154 Crash的数字表格 乘法逆元+莫比乌斯反演(TLE)