ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

min25筛 学习笔记

2021-08-06 02:01:23  阅读:240  来源: 互联网

标签:prime le limits min25 笔记 学习 sum rm ll


之前做题要用到min25,就断断续续地学了几下,用后即忘,简直就是浪费时间。不如现在好好记下来,巩固一下记忆。

在找博客学习过程中发现了一个写得非常好的博客:Min-25筛学习笔记 | LNRBHAW,配合Min_25 筛 - OI Wiki (oi-wiki.org)食用,效果很好。

作用和适用范围

min25可以以低于线性的复杂度计算出一类积性函数\(f\)​的前缀和(不一定是局限于积性函数,见题目Hasse Diagram)。

要求对于素数\(p\)​,\(f(p)\)​是一个关于\(p\)​的低次多项式,并且\(f(p^c)\)​可以快速计算出来。例如\(f(p)=p^2-p\)​,\(f(p^k)=p^k(p^k-1)\)​这样的函数。

实现

定义:

  • \({\rm lpf}(n)\)​为\(n\)​的最小质因数。例如\({\rm lpf}(6)=2\)​,特别地,\({\rm lpf}(1)=1\)​。

  • \(p_k\)代表第\(k\)大的质数。

  • \(S(n,k)=\sum\limits_{i=2}^{n}{[p_k<{\rm lpf}(i)]f(i)}\)

  • \(F_{\rm prime}(n)=\sum\limits_{p_i\le n}{f(p_i)}\)

这样,\(S(n,1)+f(1)\)就是我们要求的前缀和。

\(S(n,k)\)中可以分为两个部分分别计算,一个是素数的部分,一个合数部分。

  • 素数:即为\(S(n,k)=F_{\rm prime}(n)-F_{\rm prime}(p_{k-1})\)

  • 合数:通过枚举每个\(i\)​的最小质因子及其次数就可以求得递推式(枚举最小质因子保证了不重不漏)。

素数部分

因为\(f(p)=\sum{a_kp^k}\)​​,有

\(F_{\rm prime}(n)=\sum\limits_{p_i\le n}{f(p_i)}=\sum\limits_{p_i\le n}{\sum{a_kp_i^k}}=\sum{a_k\sum\limits_{p_i\le n}{p_i^k}}=\sum{a_k\sum\limits_{p_i\le n}{g(p_i)}}\)

令\(g(p)=p^k\)​​,可以发现\(g(p)\)​​是完全积性函数(注意\(g(p)\)​​系数是1,否则就不是完全积性函数了)。所以现在要求的是\(F_{\rm prime}(n)=\sum\limits_{p_i\le n}{g(p_i)}\)​​

设函数\(G(n,i)\)为用前\(i\)个质数做埃式筛剩下的\(g\)值之和,这样\(F_{\rm prime}(n)=G(n,\lfloor\sqrt{n}\rfloor)\)

它由有两个部分组成:

  • \(p_i^2>n\),筛完了,\(G(n,i)=G(n,i-1)\)

  • \(p_i^2\le n\)​,\(p_i\)的​在区间\([p_i,\lfloor \frac{n}{p_i}\rfloor]\)​且没被前\(i-1\)​个素数筛掉的倍数会被筛掉,要减去。而这些倍数位于\(G(\lfloor \frac{n}{p_i}\rfloor,i-1)\)​中,除此之外,\(G(\lfloor \frac{n}{p_i}\rfloor,i-1)\)​还包含了\(p_1,p_2,...,p_{i-1}\)​这些值没被筛去,要从中减去\(F_{\rm prime}(p_{i-1})\)​

综上有

\(G(n,i)=G(n,i-1)-[p_i^2\le n]g(p_i)(G(\lfloor \frac{n}{p_i}\rfloor,i-1)-F_{\rm prime}(p_{i-1}))\)​​

注意\(g(p)=p^k,p{\rm 是质数}\),是完全积性函数,故解析式为\(g(x)=x^k,x\in Z^+\),可以O(1)求出其\(g\)前缀和\(G(n,0)\),\(F_{\rm prime}(p_{i-1})\)也可以在O(\(\sqrt{n}\))预处理出来。

可以发现都是形如\(\lfloor \frac{n}{p_i}\rfloor\),可以使用数论分块+递推预处理出\(2\sqrt{n}\)项\(F_{\rm prime}\)值,分别为\(F_{\rm prime}(i),i=1,...,\sqrt{n}\)和\(F_{\rm prime}(\frac{n}{i}),i=1,...,\sqrt{n}\)。这\(2\sqrt{n}\)个就是下面要用的值。

合数部分

枚举最小质数可得递推式如下:

\(S(n, k)=\sum\limits_{i\ge k}\sum\limits_{c>1}^{p_i^{c+1}\le n}{(f(p_i^c)S(\frac{n}{p_i^c},i+1)+f(p_i^{c+1})})\)​​

由于\(f(x)\)是积性函数,\(f(p_i^c)S(\frac{n}{p_i^c},i+1)\)相当于最小质因数为\(p_i^c\)的合数对应的值的和,最后还要补上一个\(f(p_i^{c+1})\)。

合并

最终结果即为:

\(S(n, k)=\sum\limits_{i\ge k}\sum\limits_{c>1}^{p_i^{c+1}\le n}{(f(p_i^c)S(\frac{n}{p_i^c},i+1)+f(p_i^{c+1})})+F_{\rm prime}(n)-F_{\rm prime}(p_{k-1})\)

直接递归计算即可。

细节

  • 实现中将\(F_{\rm prime}(i),i=1,...,\sqrt{n}\)这小于\(\sqrt{n}\)的映射到\(\rm id1[i]\),\(F_{\rm prime}(\frac{n}{i}),i=1,...,\sqrt{n}\)​这些大于\(n\)的映射到\(\rm id2[i]\),相当于hashmap。

  • 求素数部分用的是非递归写法,直接递推上去即可。

例题

P5325 【模板】Min_25筛

#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 10;
const int M = 1e9 + 7;
const int INF = 0x3f3f3f3f;
typedef long long ll;


ll num[N], ans1[N], ans2[N];
ll id1[N], id2[N];
ll pri[N];
ll m, n;
int cnt;
bool isnp[N];

void init() {
    isnp[1] = 1;
    pri[cnt++] = 1;
    for(int i = 2; i < N; i++) {
        if(!isnp[i]) {
            pri[cnt++] = i;
            for(int j = 2 * i; j < N; j += i)
                isnp[j] = 1;
        }
    }
}
ll r2 = 500000004;
ll r6 = 166666668;

ll sum1(ll x) { // 求和,用于求初始值S(n,0)
    x %= M;
    return x * (x + 1) % M * r2 % M;
}

ll sum2(ll x) {
    x %= M;
    return x * (x + 1) % M * (2 * x + 1) % M * r6 % M;
}

ll f1(ll x) { // 一次幂
    x %= M;
    return x;
}

ll f2(ll x) { // 二次幂
    x %= M;
    return x * x % M;
}

ll f(ll x) { // 真正的
    x %= M;
    return x * (x - 1) % M;
}

int ID(ll x) { // hash映射
    if(x <= m) return id1[x];
    return id2[n / x];
}

void solve(ll n, ll (*f)(ll), ll (*sum)(ll), ll ans[]) {
    static ll sump[N]; // 预处理素数和
    int cur = 0;
    int mx = 0;
    m = 0;
    while(m * m <= n) {
        if(pri[mx] <= m) {
            mx++;
            sump[mx] = (sump[mx - 1] + f(pri[mx])) % M;
        }
        m++;
    }
    ll i = 1;
    while(i <= n) { // 数论分快求出映射关系和初始值
        num[cur] = n / i;
        ans[cur] = (sum(num[cur]) - f(1)) % M;
        if(num[cur] <= m) id1[num[cur]] = cur;
        else id2[n / num[cur]] = cur;
        cur++;
        i = n / (n / i) + 1;
    }
    for(int i = 1; i < mx; i++) { // 直接递推
        for(int j = 0; j < cur && 1ll * pri[i] * pri[i] <= num[j]; j++) {
            ans[j] -= f(pri[i]) * (ans[ID(num[j] / pri[i])] - sump[i - 1]) % M;
            ans[j] %= M;
        }
    }
}

ll F(ll n, int k) {
    ll res = 0;
    for(int i = k; 1ll * pri[i] * pri[i] <= n; i++) {
        ll p = pri[i];
        while(p * pri[i] <= n) {
            res += (f(p) * F(n / p, i + 1) % M + f(p * pri[i])) % M;
            res %= M;
            p = p * pri[i];
        }
        
    }
    res += (ans2[ID(n)] - ans1[ID(n)]) - ((k - 1) ? (ans2[ID(pri[k - 1])] - ans1[ID(pri[k - 1])]) : 0) % M;
    return res % M;
}



int main() {
    init();
    while(cin >> n) {
        solve(n, f1, sum1, ans1); // 素数部分
        solve(n, f2, sum2, ans2);
        ll ans = (F(n, 1) + 1 + M) % M; // 合数+素数统计
        cout << (ans % M + M) % M << endl;
    }
}

标签:prime,le,limits,min25,笔记,学习,sum,rm,ll
来源: https://www.cnblogs.com/limil/p/15106764.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有