首页 > 其他分享 >min25筛学习笔记

min25筛学习笔记

时间:2023-08-03 13:34:00浏览次数:38  
标签:prime right limits int min25 笔记 学习 sum left

min25筛

min25筛用于求一类数论函数的前缀和,适用于函数在素数处的取值可以用一个关于此素数的多项式来表示的数论函数。

处理质数部分

这部分我们需要解决\(\sum\limits_{p \subseteq prime}f(p)\),这里简单起见,假设\(f(p)=p^t\)

用\(s_i\)表示前\(i\)个质数之和,用\(LPF_i\)表示\(i\)的最小质因子,\(p_i\)表示第\(i\)个质数,设$$g(n,i)=\sum\limits_{k=1}^n\left[ LPF_k>p_i \vee k \subseteq prime \right]k^t$$

因此有:$$g(n,0)=\sum\limits_{k=1}nkt$$

假设\(\sqrt{n}\)内有\(cnt\)个质数,那么

\[g(n,cnt)=\sum\limits_{k=1}^n[k \subseteq prime]k^t \]

因此我们的目标就是求解\(g(n,cnt)\),接下来考虑转移,即\(g(n,i)\)和\(g(n,i-1)\)的关系。注意到:

\[g(n,i-1)=\sum\limits_{k=1}^n\left[ LPF_k>p_i \vee k \subseteq prime \right]k^t+\sum\limits_{k=1}^n\left[ LPF_k=p_i \wedge k \not\subseteq prime \right]k^t \]

因此有转移:

  • \(p_i*p_i>n\)

    \[g(n,i)=g(n,i-1) \]

  • \(p_i*p_i\leq n\)

\[\begin{aligned} g(n,i)&=g(n,i-1)-\sum\limits_{k=1}^n\left[ LPF_k=p_i \vee k \not\subseteq prime \right]k^t\\ &=g(n,i-1)-p_i^t\cdot \sum\limits_{k=1}^{\left\lfloor \frac{n}{p_i} \right\rfloor}[LPF_k\ge p_i]k^t\\ &=g(n,i-1)-p_i^t\cdot \left[ g\left(\left\lfloor \frac{n}{p_i} \right\rfloor,i-1\right)- s_{i-1}\right] \end{aligned} \]

处理完整问题

设\(S(n,i)=\sum\limits_{k=1}^n[LPF_k>p_i]f(k)\),则我们的目标是求解\(S(n,0)\)

考虑如何转移:

  • \(p_i*p_i > n\)

\[S(n,i-1)=g(n,i-1)-s_{i-1} \]

  • \(p_i*p_i\leq n\)

\[\begin{aligned} S(n,i-1)&= \sum\limits_{p_j^e\leq n\wedge j\ge i}f(p_j^e)\cdot\left([e>1]+\sum\limits_{k=1}^{\left\lfloor \frac{n}{p_j^e}\right\rfloor}[LPF_k>p_j]f(k) \right)+g(n,cnt)-s_{i-1}\\ &=g(n,cnt)-s_{i-1}+\sum\limits_{p_j^e\leq n\wedge j\ge i}f(p_j^e)\cdot \left( S\left(\left\lfloor \frac{n}{p_j^e} \right\rfloor,j\right)+[e>1] \right){\kern 5pt}\left( 1 \right) \\ &=S(n,i)+\sum\limits_{p_i^e\leq n}f(p_i^e)\cdot \left( S\left(\left\lfloor \frac{n}{p_i^e} \right\rfloor,i\right)+1 \right){\kern 86pt}\left( 2 \right) \end{aligned} \]

可以直接按照\(\left( 1 \right)\)式暴力求解,也可以按照\(\left( 2 \right)\)式仿照质数部分DP递推求解。

时间复杂度分析

对于质数部分的求解,我们只关注\(t \subseteq \left\{ \left\lfloor \frac{n}{k} \right\rfloor|k \subseteq \mathbb{N}^*\right\}\)处\(g\left( t,i \right)\)的取值,而只有当\(p_i*p_i \leq t\)时才会产生转移。因此复杂度可以估计为:

\[\sum\limits_{i=1}^{\sqrt{n}}\frac{\sqrt{i}}{\ln{\sqrt{i}}}+\sum\limits_{i=1}^{\sqrt{n}}\frac{\sqrt{\frac{n}{i}}}{\ln{\sqrt{\frac{n}{i}}}}=O\left( \int_1^{\sqrt{n}}\frac{\sqrt{\frac{n}{x}}}{\log_2\sqrt{\frac{n}{x}}}dx \right)=O\left( \frac{n^{\frac{3}{4}}}{\log_2{n}} \right) \]

对于第二部分的复杂度分析类似,值得一提的是,按照\(\left( 1 \right)\)式暴力求解的时间复杂度也是对的,而且这种写法更加简洁。

一些例子

  • 求\(\sum\limits_{i=1}^n\sum\limits_{i=1}^nf^k\left( \operatorname{gcd}(i,j) \right)\)其中\(f(x)\)表示\(x\)的次大质因子,重复的质因子计算多次,例如\(f(4)=f(6)=2\)规定\(f(1)=0,f(p)=1\)其中\(p\)为质数。\(n,k\leq 2\times 10^9\)

\[\begin{aligned} &{\kern 10pt}\sum\limits_{i=1}^n\sum\limits_{i=1}^nf^k\left( \operatorname{gcd}(i,j) \right)\\ &=\sum\limits_{d=1}^n\sum\limits_{i=1}^{\left\lfloor \frac{n}{d} \right\rfloor}\sum\limits_{j=1}^{\left\lfloor \frac{n}{d} \right\rfloor}\left[ \operatorname{gcd}(i,j)=1 \right] f^k(d)\\ &=\sum\limits_{d=1}^nf^k(d)\sum\limits_{i=1}^{\left\lfloor \frac{n}{d} \right\rfloor}\mu(i)\cdot \left\lfloor \frac{n}{d\cdot i} \right\rfloor^2\\ &令g(n)=\sum\limits_{i=1}^n\mu(i)\cdot \left\lfloor \frac{n}{i} \right\rfloor^2,则原式=\sum\limits_{d=1}^nf^k(d)\cdot g\left(\left\lfloor\frac{n}{d} \right\rfloor\right)\\ &于是接下来要做的事情是求g的前缀和以及求f^k的前缀和,g的前缀和可以用数论分块+杜教筛做。\\ &然后考虑求f的前缀和,设S(n,i)=\sum\limits_{j=1}^n\left[ LPF_j>p_i \vee j \subseteq prime \right]f^k(j)\\ &若p_i\cdot p_i>n,S(n,i-1)=\sum\limits_{j=1}^n\left[ j \subseteq prime \right]\\ &若p_i\cdot p_i\leq n,S(n,i-1)=S(n,i)+S\left(\left\lfloor \frac{n}{p_i} \right\rfloor,i-1\right)-\sum\limits_{j=1}^{\left\lfloor \frac{n}{p_i} \right\rfloor}\left[ j \subseteq prime \right]+p_i\cdot \left( \left( \sum\limits_{j=1}^{\left\lfloor \frac{n}{p_i} \right\rfloor}[j \subseteq prime] \right)-i+1 \right)\\ &对于n以内质数个数可以直接套用前面min25筛的质数部分 \end{aligned} \]

code:

#include <cstdio>
#include <iostream>
#define int unsigned int

using namespace std;

const int M = 3e5 + 5;

int prime[M], cnt = 0, a[M], g[M << 1], f[M << 1], primeMi[M];
int n, k, v[M << 1], ss = 0, mu[M], mus[M];

inline int getIdx(int o) { return (o > M) ? (M + n / o) : o; }

int mi(int o, int p) {
    int yu = 1;
    while (p) {
        if (p & 1) yu *= o;
        o *= o;
        p >>= 1;
    }
    return yu;
}

int cal(int o) {
    if (o < M) return mu[o];
    if (mus[n / o]) return mus[n / o];
    int yu = 1;
    for (int l = 2, r; l <= o; l = r + 1) {
        r = o / (o / l);
        yu -= (r - l + 1) * cal(o / l);
    }
    return mus[n / o] = yu;
}

int find(int o) {
    int l = 1, r = ss, res = 0;
    while (l <= r) {
        int mid = (l + r) >> 1;
        if (o <= v[mid] / o) {
            res = mid;
            l = mid + 1;
        } else
            r = mid - 1;
    }
    return res;
}

signed main() {
    scanf("%u%u", &n, &k);
    mu[1] = 1;
    for (int i = 2; i < M; i++) {
        if (!a[i]) {
            prime[++cnt] = i;
            primeMi[cnt] = mi(prime[cnt], k);
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && i * prime[j] < M; j++) {
            a[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            } else {
                mu[i * prime[j]] = -mu[i];
            }
        }
    }
    mu[0] = 0;
    for (int i = 1; i < M; i++) mu[i] += mu[i - 1];
    for (int l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        v[++ss] = n / l;
        g[getIdx(n / l)] = n / l - 1;
    }
    for (int i = 1; i <= cnt; i++) {
        for (int j = 1; j <= ss; j++) {
            int t = v[j];
            if (prime[i] > t / prime[i]) break;
            g[getIdx(t)] -= g[getIdx(t / prime[i])] - i + 1;
        }
    }
    for (int j = 1; j <= ss; j++) {
        f[getIdx(v[j])] = g[getIdx(v[j])];
    }
    for (int i = cnt; i >= 1; i--) {
        int idx = find(prime[i]);
        for (int j = idx; j >= 1; j--) {
            int t = v[j];
            if (prime[i] > t / prime[i]) continue;
            int p1 = getIdx(t);
            int p2 = getIdx(t / prime[i]);
            f[p1] += f[p2] - g[p2] + primeMi[i] * (g[p2] - i + 1);
        }
    }
    int ans = 0;
    f[0] = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        int t = n / l, yu = 0;
        r = n / t;
        for (int l1 = 1, r1; l1 <= t; l1 = r1 + 1) {
            r1 = t / (t / l1);
            yu += (cal(r1) - cal(l1 - 1)) * (t / l1) * (t / l1);
        }
        ans += (f[getIdx(r)] - f[getIdx(l - 1)]) * yu;
    }
    cout << ans << endl;
    return 0;
}
  • 对于正整数\(n\),定义如下操作每次选择\(p_i,p_j\)满足\(p_i|n \wedge p_j|n\wedge p_i\neq p_j\)将\(n\)变为\(\frac{n}{p_i\cdot p_j}\)。求有多少n以内的正整数满足可以通过若干次操作变为\(1\)。\(n\leq 10^{11}\)

\[\begin{aligned} &首先可以通过若干次操作变为1的正整数一定满足质因数的次数之和为偶数且质因数的最大幂次不超过质因数总幂次的一半。\\ &设f(n,j,m,s)=\sum\limits_{i=1}^n\left[ LPF_i>p_j\wedge 当前最大质因子次数为m,质因子此数之和为s 且加入i后(m,s)能到达合法状态\right]\\ &则我们要求f(n,0,0,0)\\ &则有转移f(n,j,m,s)=[在插入一个新质因子后,(m,s)能到达合法状态]\cdot \left( \left( \sum\limits_{i=1}^n[i \subseteq prime] \right) -j\right)+\\ &\sum\limits_{i=j+1}\sum\limits_{a=1}\left( f\left( \left\lfloor \frac{n}{p_i^a} \right\rfloor,i,\max(m,a),s+a \right) +[a>1\wedge 插入p_i^a后(m,s)能到达合法状态] \right) \end{aligned} \]

code:

#pragma GCC optimize(2)
#include <chrono>
#include <cstdio>
#include <ctime>
#include <iostream>

using namespace std;

const int M = 317005;

int prime[M], cnt = 0, a[M], ss = 0;
long long g[M << 1], n, v[M << 1];

int getIdx(long long o) { return (o > M) ? (n / o + M) : o; }

int check(int sum, int mx) { return ((mx << 1) <= sum) && ((sum & 1) == 0); }

long long S(long long o, int p, int sum, int mx) {
    if (o < prime[p + 1]) return 0;
    long long yu = (check(sum + 1, max(mx, 1))) ? (g[getIdx(o)] - p) : 0;
    for (int i = p + 1; i <= cnt; i++) {
        if (prime[i] > o / prime[i]) break;
        int t = 1;
        for (long long j = prime[i]; j <= o; j *= prime[i], t++) {
            yu += S(o / j, i, sum + t, max(mx, t));
            yu += check(sum + t, max(mx, t)) * (t > 1);
        }
    }
    return yu;
}

signed main() {
    for (int i = 2; i < M; i++) {
        if (!a[i]) prime[++cnt] = i;
        for (int j = 1; j <= cnt && i * prime[j] < M; j++) {
            a[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
        }
    }
    scanf("%lld", &n);
    for (long long l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        g[getIdx(n / l)] = n / l - 1;
        v[++ss] = n / l;
    }
    for (int i = 1; i <= cnt; i++) {
        for (int j = 1; j <= ss; j++) {
            long long t = v[j];
            if (prime[i] > t / prime[i]) break;
            g[getIdx(t)] -= g[getIdx(t / prime[i])] - i + 1;
        }
    }
    cout << S(n, 0, 0, 0) << endl;
    return 0;
}
  • 求\(\sum\limits_{i=1}^n\sum\limits_{j=1}^mf(i\cdot j)\),令积性函数\(f(n)\)满足\(f(p^c)=p^{\operatorname{gcd}(c,k)}\),其中\(p\)为给定常数,\(p\)为素数,\(c\)为正整数。\(n\leq 10^5,m\leq 10^{10},k\leq 10^9\)

\[\begin{aligned} &令F(n,x)=\sum\limits_{i=1}^nf(i\cdot x),则ans=\sum\limits_{i=1}^nF(m,i)\\ &我们仿照min25筛的思想,每次提取出一个质数的贡献,假设提取质因子p的贡献,p在x中次数为c\\ &F(n,x)=\sum\limits_{p^e\leq n}f(p^e)\cdot \left( F\left( \left\lfloor \frac{n}{p^e} \right\rfloor,\frac{x}{p^c} \right) -F\left( \left\lfloor \frac{n}{p^{e+1}} \right\rfloor,\frac{x}{p^{c-1}} \right) \right)\\ &这样看起来状态数有O(n\cdot \sqrt{m}),但如果每次都取x的最大质因子,状态数和转移数都在可接受的范围内\\ &极限数据下状态数仅有10^6级别,转移数在10^7次级别。\\ &接下来需要解决的问题是F(n,1)\\ &F(n,1)=\sum\limits_{i=1}^nf(i),注意到f和id在质数处的取值都一样,f的前缀和可以很方便的用Powerful Number解决。\\ &设f=id*h,则S_f(n)=\sum\limits_{d \subseteq PN}h(d)\cdot \frac{\left\lfloor \frac{n}{d} \right\rfloor\cdot (\left\lfloor \frac{n}{d}+1 \right\rfloor)}{2},h(p^c)=p^{\operatorname{gcd}(c,k)}-\sum\limits_{t=1}^cp^t\cdot h(p^{c-t}) \end{aligned} \]

#include <assert.h>

#include <algorithm>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <map>
#include <vector>
#define int long long

using namespace std;

const int mo = 1e9 + 7;
const int inv2 = (mo + 1) >> 1;
const int M = 1e5 + 5;
const int N = 2e6 + 5;

int n, m, prime[200005], cnt = 0, vis[N + 5], k, tmp = 0, pns[1000005],
                         cs[N + 5];
vector<int> h[10005];
vector<pair<int, int>> pn;
vector<int> pp[100005];
map<int, int> mp[100005];

inline int gcd(int o, int p) { return (!p) ? o : gcd(p, o % p); }

inline int getIdx(int o) { return (o <= M) ? o : (M + m / o); }

inline int Mod(int o) { return (o >= mo) ? (o - mo) : o; }

int mi(int o, int p) {
    int yu = 1;
    while (p) {
        if (p & 1) yu = yu * o % mo;
        o = o * o % mo;
        p >>= 1;
    }
    return yu;
}

vector<array<int, 3>> factor(int o) {
    vector<array<int, 3>> ans;
    for (int i = 1; i <= cnt; i++) {
        if (o % prime[i] == 0) {
            int yu = 0;
            int tt = 1;
            while (o % prime[i] == 0) {
                yu++;
                tt *= prime[i];
                o /= prime[i];
            }
            ans.push_back(array<int, 3>{yu, prime[i], tt});
        }
        if (prime[i] * prime[i] >= o) break;
    }
    if (o > 1) {
        ans.push_back(array<int, 3>{1, o, o});
    }
    return ans;
}

void init() {
    memset(pns, -1, sizeof(pns));
    memset(cs, -1, sizeof(cs));
    for (int i = 2; i <= N; i++) {
        if (!vis[i]) prime[++cnt] = i;
        for (int j = 1; j <= cnt && i * prime[j] <= N; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                break;
            }
        }
    }
    for (int i = 1; i <= cnt; i++) {
        if (prime[i] > 100000) break;
        h[i].push_back(1);
        for (int j = prime[i], t = 1; j <= m; j *= prime[i], t++) {
            int sum = 0;
            for (int l = 0; l < t; l++) {
                sum = Mod(sum + mi(prime[i], t - l) * h[i][l] % mo);
            }
            h[i].push_back(Mod(mi(prime[i], gcd(t, k)) - sum + mo));
        }
    }
}

void pre(int o, int p, int q) {
    cs[o] = p % mo;
    for (int j = q; j <= cnt; j++) {
        if (o * prime[j] > N) break;
        for (int t = prime[j] * o, s = 1; t <= N; t *= prime[j], s++) {
            pre(t, p * mi(prime[j], gcd(s, k)) % mo, j + 1);
        }
    }
}

void getPN(int o, int p, int q) {
    pn.push_back(make_pair(o, p));
    for (int j = q; j <= cnt; j++) {
        if (prime[j] > 100000) break;
        if (prime[j] * prime[j] > m / o) break;
        for (int t = prime[j] * prime[j] * o, s = 2; t <= m;
             t *= prime[j], s++) {
            getPN(t, p * h[j][s] % mo, j + 1);
        }
    }
}

inline int s2(int o) { return o * (o + 1) % mo * inv2 % mo; }

int calH(int o) {
    int t = getIdx(o);
    if (pns[t] != -1) return pns[t];
    int ans = 0;
    for (auto [l, r] : pn) {
        if (l > o) break;
        ans = Mod(ans + r * s2(o / l % mo) % mo);
    }
    return pns[t] = ans;
}

int cal(int o, vector<array<int, 3>> p, int q) {
    if (o == 0) return 0;
    if (o == 1) return cs[q];
    if (p.empty()) return calH(o);
    if (o * q <= N) return pp[q][o];
    if (mp[q].find(o) != mp[q].end()) return mp[q][o];
    tmp++;
    int ans = 0;
    int q1 = q / p.back()[2];
    int q2 = q1 * p.back()[1];
    for (int i = 1, t = 0; i <= o; i *= p.back()[1], t++) {
        vector<array<int, 3>> p1 = p;
        vector<array<int, 3>> p2 = p;
        p1.pop_back();
        p2.pop_back();
        p2.push_back(array<int, 3>{1, p.back()[1], p.back()[1]});
        ans =
            (ans + mi(p.back()[1], gcd(k, t + p.back()[0])) *
                       (cal(o / i, p1, q1) - cal(o / i / p.back()[1], p2, q2)) %
                       mo) %
            mo;
        ans = Mod(ans + mo);
    }
    return mp[q][o] = ans;
}

signed main() {
    scanf("%lld%lld%lld", &n, &m, &k);
    init();
    pre(1, 1, 1);
    getPN(1, 1, 1);
    sort(pn.begin(), pn.end());
    for (int i = 1; i <= n; i++) {
        pp[i].push_back(0);
        for (int j = 1; i * j <= N; j++) {
            int yu = Mod(pp[i][pp[i].size() - 1] + cs[i * j]);
            pp[i].push_back(yu);
        }
    }
    int ans = 0;
    for (int i = 1; i <= n; i++) {
        ans = (ans + cal(m, factor(i), i)) % mo;
    }
    cout << ans << endl;
    return 0;
}

标签:prime,right,limits,int,min25,笔记,学习,sum,left
From: https://www.cnblogs.com/clapp/p/17603088.html

相关文章

  • (一)flask学习笔记
    1、flask路由(用了装饰器)@app.route('/login',methods=["GET","POST"])deflogin():pass2、路由参数@app.route('/login',methods=["GET","POST"],endpoint='login')deflogin():pass  ......
  • Markdown学习
    标题(#+空格+标题名字一级标题)(##+空格+标题名字二级标题)(###+空格+标题名字三级标题)字体(加粗)(斜体)(加粗斜体)(划线)引用好好学习天天向上分割线图片超链接百度列表第一点第二点第三点第一点第二点第三点表格姓名班级学号马迪雅602代......
  • Markdown学习
    标题(#+空格+标题名字一级标题(##+空格+标题名字二级标题(###+空格+标题名字三级标题字体加粗斜体加粗斜体划线引用好学学习,天天向上分割线图片超链接百度列表第一点第二点第三点第一点第二点第三点表格姓名年龄性别曾老师25男......
  • 算法笔记(二)—— 认识N(logN)的排序算法
    递归行为的时间复杂度估算整个递归过程是一棵多叉树,递归过程相当于利用栈做了一次后序遍历。对于master公式,T(N)表明母问题的规模为N,T(N/b)表明每次子问题的规模,a为调用次数,加号后面表明,除去调用之外,剩余语句的复杂度是多少,算出d。根据上次三个判断公式进行算法时间复杂度计算......
  • 小柏实战学习FineBI(图文教程一)
    前言:一定要知道百度,必应,谷歌这个三个网站,这三个不知道的话也要会使用ChatGPT,并且要学会看报错信息,学会优雅的提问.  本节课主题:FineBI的下载,安装,配置. 零:官网填写信息,获取试用码:https://www.finebi.com/  一:下载:https://www.finebi.com/product/downlo......
  • 【Linux】Kali Linux 渗透安全学习笔记(2) - OneForAll 简单应用
    OneForAll(以下简称“OFA”)是一个非常好用的子域收集工具,可以通过一级域名找到旗下的所有层级域名,通过递归的方式我们很容易就能够知道此域名下的所有域名层级结构,对于进一步通过域名推测站点功能起到非常重要的作用。声明:本文测试的站点为自家站点仅做学习使用,不存在侵犯网络......
  • C学习(一)基本概念
    《C语言程序设计:现代方法》第2章,2.1C程序转为机器码,需要3个步骤:预处理:预处理器preprocessor,执行#开头的命令/指令,类似于编辑器,可添加修改程序编译:编译器compiler,.c-->.exe/.out机器指令/目标代码【Windows是.exe,Linux是.out】链接:链接器linker,把编译器产生的目标代码和其......
  • 关于菜鸡学习RHEL8的一些小笔记--->stratis和vdo
    #注:stratis和vdo目前都是属于redhat的预览技术,并没有实际投入到生产环境stratis精简卷(适用于海量应用场景,只需关注精简池容量,无需去管文件系统):stratis(redhat8的新功能)会吧磁盘放在一个精简配置的共享池子里面(精简存储池),stratis文件系统也是没有固定的大小,也不会提前分配没有......
  • Qt+GDAL开发笔记(二):在windows系统msvc207x64编译GDAL库、搭建开发环境和基础Demo
    前言  上一篇使用mingw32版本的gdal,过程曲折,为更好的更方便搭建环境,在windows上msvc方式对于库比较友好。<br>大地坐标简介概述  大地坐标(Geodeticcoordinate)是大地测量中以参考椭球面为基准面的坐标,地面点P的位置用大地经度L、大地纬度B和大地高H表示。原理  当点在......
  • JavaScript学习 -- RSA算法应用实例及公钥私钥的生成方法
    正文:RSA算法是一种非对称加密算法,用于加密、解密和数字签名等场景。本文将介绍如何在JavaScript中使用RSA算法,并提供一个实际的案例,同时也会说明如何生成公钥和私钥。首先,确保您已经引入了jsencrypt库。以下是一个使用RSA算法进行加密和解密的示例,同时也包含了公钥和私钥的生成方法......