HDU 6607 Easy Math Problem（杜教筛 + min

Easy Math Problem

推式子

∑i=1n∑j=1ngcd(i,j)Klcm(i,j)[gcd(i,j)∈prime]∑i=1n∑j=1ngcd(i,j)K−1ij[gcd(i,j)∈prime]∑d∈primendK+1∑i=1nd∑j=1ndij[gcd(i,j)==1]对∑i=1n∑j=1nij[gcd(i,j)==1]化简2(∑i=1niiϕ(i)+[i==1]2)−1=∑i=1ni2ϕ(i)∑d∈primendK+1∑i=1ndi2ϕ(i)\sum_{i = 1} ^{n} \sum_{j = 1} ^{n} gcd(i, j) ^ K lcm(i, j)[gcd(i, j) \in prime]\\ \sum_{i = 1} ^{n} \sum_{j = 1} ^{n} gcd(i, j) ^{K - 1} ij[gcd(i, j) \in prime]\\ \sum_{d \in prime} ^{n} d ^ {K + 1} \sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{\frac{n}{d}}ij[gcd(i, j) == 1]\\ 对\sum_{i = 1} ^{n} \sum_{j = 1} ^{n}ij[gcd(i, j) == 1]化简\\ 2(\sum_{i = 1} ^{n} i \frac{i\phi(i) + [i == 1]}{2}) - 1 = \sum_{i = 1} ^{n} i ^ 2 \phi(i) \\ \sum_{d \in prime} ^{n} d ^ {K + 1} \sum_{i = 1} ^{\frac{n}{d}} i ^ 2 \phi(i)\\ i=1∑nj=1∑ngcd(i,j)Klcm(i,j)[gcd(i,j)∈prime]i=1∑nj=1∑ngcd(i,j)K−1ij[gcd(i,j)∈prime]d∈prime∑ndK+1i=1∑dnj=1∑dnij[gcd(i,j)==1]对i=1∑nj=1∑nij[gcd(i,j)==1]化简2(i=1∑ni2iϕ(i)+[i==1])−1=i=1∑ni2ϕ(i)d∈prime∑ndK+1i=1∑dni2ϕ(i)

接下来就是快了的数论分块了，首先我们用min_25得到区间[l,r][l, r][l,r]部分的质数的贡献，然后再通过杜教筛得到后面，两者相乘然后累加即可，对于∑i=1nik+1\sum\limits_{i = 1} ^{n} i ^{k + 1}i=1∑nik+1这一部分求和，显然我们可以通过拉格朗日插值来求解。

代码

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>#define mp make_pair
#define pb push_back
#define endl '\n'
#define mid (l + r >> 1)
#define lson rt << 1, l, mid
#define rson rt << 1 | 1, mid + 1, r
#define ls rt << 1
#define rs rt << 1 | 1using namespace std;typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pii;const double pi = acos(-1.0);
const double eps = 1e-7;
const int inf = 0x3f3f3f3f;inline ll read() {ll f = 1, x = 0;char c = getchar();while(c < '0' || c > '9') {if(c == '-')    f = -1;c = getchar();}while(c >= '0' && c <= '9') {x = (x << 1) + (x << 3) + (c ^ 48);c = getchar();}return f * x;
}const int N = 1e6 + 10, mod = 1e9 + 7, inv2 = (mod + 1) >> 1, inv6 = (mod + 1) / 6;namespace Djs {int prime[N], cnt;ll phi[N];bool st[N];void init() {phi[1] = 1;for(int i = 2; i < N; i++) {if(!st[i]) {prime[++cnt] = i;phi[i] = i - 1;}for(int j = 1; j <= cnt && 1ll * i * prime[j] < N; j++) {st[i * prime[j]] = 1;if(i % prime[j] == 0) {phi[i * prime[j]] = phi[i] * prime[j];break;}phi[i * prime[j]] = phi[i] * (prime[j] - 1);}}for(int i = 1; i < N; i++) {phi[i] = (phi[i - 1] + 1ll * i * i % mod * phi[i] % mod) % mod;}}ll calc1(ll n) {n %= mod;return (n * (n + 1) % mod * inv2 % mod) * (n * (n + 1) % mod * inv2 % mod) % mod;}ll calc2(ll n) {n %= mod;return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;}unordered_map<ll, ll> ans_s;ll S(ll n) {if(n < N) return phi[n];if(ans_s.count(n)) return ans_s[n];ll ans = calc1(n);for(ll l = 2, r; l <= n; l = r + 1) {r = n / (n / l);ans = ((ans - (calc2(r) - calc2(l - 1)) * S(n / l) % mod) % mod + mod) % mod;}return ans_s[n] = ans;}
}namespace Lagrange {const int N = 110;ll fac[N], pre[N], suc[N], inv[N], prime[N], sum[N], mu[N], n, k, cnt;bool st[N];ll quick_pow(ll a, int n) {ll ans = 1;while(n) {if(n & 1) ans = ans * a % mod;a = a * a % mod;n >>= 1;}return ans;}void init1() {fac[0] = inv[0] = 1;for(int i = 1; i < N; i++) {fac[i] = 1ll * fac[i - 1] * i % mod;}inv[N - 1] = quick_pow(fac[N - 1], mod - 2);for(int i = N - 2; i >= 1; i--) {inv[i] = 1ll * inv[i + 1] * (i + 1) % mod;}}void init() {for(int i = 1; i < N; i++) {st[i] = 0;}cnt = 0;sum[1] = 1;for(int i = 2; i < N; i++) {if(!st[i]) {prime[cnt++] = i;sum[i] = quick_pow(i, k);}for(int j = 0; j < cnt && i * prime[j] < N; j++) {st[i * prime[j]] = 1;sum[i * prime[j]] = 1ll * sum[i] * sum[prime[j]] % mod;if(i % prime[j] == 0) break;}}for(int i = 1; i < N; i++) {sum[i] = (sum[i] + sum[i - 1]) % mod;}}ll solve(ll n) {n %= mod;ll ans = 0;pre[0] = suc[k + 3] = 1;for(int i = 1; i <= k + 2; i++) pre[i] = 1ll * pre[i - 1] * (n - i) % mod;for(int i = k + 2; i >= 1; i--) suc[i] = 1ll * suc[i + 1] * (n - i) % mod;for(int i = 1; i <= k + 2; i++) {ll a = 1ll * pre[i - 1] * suc[i + 1] % mod, b = 1ll * inv[i - 1] * inv[k + 2 - i] % mod;if((k + 2 - i) & 1) b *= -1;ans = ((ans + 1ll * sum[i] * a % mod * b % mod) % mod + mod) % mod;}return (ans - 1 + mod) % mod;}
}namespace Min_25 {int prime[N], id1[N], id2[N], cnt, m, k, T;ll g[N], sum[N], a[N], n;bool st[N];int ID(ll x) {return x <= T ? id1[x] : id2[n / x];}void init(ll x, int y) {n = x, k = y;cnt = 0, m = 0;T = sqrt(n + 0.5);for(int i = 2; i <= T; i++) {if(!st[i]) {prime[++cnt] = i;sum[cnt] = (sum[cnt - 1] + Lagrange::quick_pow(i, k)) % mod;}for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {st[i * prime[j]] = 1;if(i % prime[j] == 0) {break;}}}for(ll l = 1, r; l <= n; l = r + 1) {r = n / (n / l);a[++m] = n / l;if(a[m] <= T) id1[a[m]] = m;else id2[n / a[m]] = m;g[m] = Lagrange::solve(a[m]);}for(int j = 1; j <= cnt; j++) {for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {g[i] = ((g[i] - Lagrange::quick_pow(prime[j], k) * (g[ID(a[i] / prime[j])] - sum[j - 1]) % mod) % mod + mod) % mod;}}for(int i = 1; i <= T; i++) {st[i] = 0;}}ll solve(ll x) {if(x <= 1) return 0;return g[ID(x)];}
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);Djs::init();Lagrange::init1();int T = read();while(T--) {ll n = read(), k = read() + 1;Lagrange::k = k;Lagrange::init();Min_25::init(n, k);ll ans = 0;for(ll l = 1, r; l <= n; l = r + 1) {r = n / (n / l);ans = (ans + (Min_25::solve(r) - Min_25::solve(l - 1)) * Djs::S(n / l) % mod) % mod;}printf("%lld\n", (ans % mod + mod) % mod);}return 0;
}

HDU 6607 Easy Math Problem（杜教筛 + min_25 + 拉格朗日插值）相关推荐

[复习]莫比乌斯反演,杜教筛,min_25筛
[复习]莫比乌斯反演,杜教筛,min_25筛莫比乌斯反演做题的时候的常用形式: \[\begin{aligned}g(n)&=\sum_{n|d}f(d)\\f(n)&=\sum_ ...
「LibreOJ Round #11」Misaka Network 与求和（杜教筛 + Min_25）
#572. 「LibreOJ Round #11」Misaka Network 与求和推式子 ∑i=1n∑j=1nf(gcd(i,j))k∑d=1nf(d)k∑i=1nd∑j=1nd[gcd(i,j ...
【LOJ#572】Misaka Network 与求和（莫比乌斯反演/杜教筛/min_25筛）
[LOJ#572]Misaka Network 与求和 https://www.cnblogs.com/cjyyb/p/10170630.html 看到次大质因子就可以想到是min_25筛了,然后只需 ...
[51NOD1847]奇怪的数学题（杜教筛+min_25筛+第二类斯特林数）
f(x)f(x)f(x)表示xxx的次大约数,有f(x)=xx的最小质因数f(x)=\frac{x}{x的最小质因数}f(x)=x的最小质因数x,那么 ∑i=1n∑j=1nsgcd(i,j)k=∑i ...
Loj #572. 「LibreOJ Round #11」Misaka Network 与求和（莫比乌斯反演 + 杜教筛 + min_25筛（递推版））
直接反演一下:∑i=1n∑i=1nf(gcd(i,j))k\sum_{i = 1}^n\sum_{i = 1}^nf(gcd(i,j))^ki=1∑ni=1∑nf(gcd(i,j))k=∑d=1n ...
1847 奇怪的数学题（杜教筛 + Min_25 + 第二类斯特林数）
1847 奇怪的数学题推式子 ∑i=1n∑j=1nsgcd(i,j)k∑d=1nsgcd(d)k∑i=1nd∑j=1nd[gcd(i,j)=1]∑d=1nsgcd(d)k(2∑i=1ndϕ(i)−1 ...
Lady Layton with Math（杜教筛）
Lady Layton with Math ∑i=1n∑j=1nϕ(gcd(i,j))∑d=1nϕ(d)∑i=1n∑j=1n[gcd(i,j)=d]∑d=1nϕ(d)∑i=1nd∑j=1nd[gcd( ...
Easy Math（ACM-ICPC 2018 徐州赛区网络预赛）（递归 + 杜教筛）
Easy Math 推式子 ∑i=1mμ(in)∑i=1mμ(i×nd×d),d是n的一个质因子i,d互质项有(−∑i=1mμ(i×nd)),由于减去了多余的非互质项,所以加上,−∑i=1mμ(i×n ...
2018 ICPC 徐州网络赛 D. Easy Math（思维，反演，杜教筛）
整理的算法模板合集: ACM模板点我看算法全家桶系列!!! 实际上是一个全新的精炼模板整合计划题目链接 https://nanti.jisuanke.com/t/A2003 Problem 计算 ...

HDU 6607 Easy Math Problem（杜教筛 + min_25 + 拉格朗日插值）

Easy Math Problem

推式子

代码

HDU 6607 Easy Math Problem（杜教筛 + min_25 + 拉格朗日插值）相关推荐

最新文章

热门文章