第五节 min25 筛
|
176
|
0
|
第二章 数论函数

1605 字
|
20 分钟

min25 筛分为两部分。第一部分可以求质数点值的前缀和（例如：n 以内的质数之和），第二部分可以处理积性函数 f(x) 前缀和，要求 f(p) 是一个低阶多项式（如：f(p)=p^2+p），且 f(p^c) 可快速求值。

以下介绍均以洛谷模版题为例：洛谷P5325 ：f(p^k)=p^k(p^k-1)

Contents

第一部分-求质数点的前缀和

设 min(p) 为 i 的最小质因子，p_j 为从小到大的第 j 个质数，设：

g(n,j)=\sum_{i=1}^{n}[i\in prime\ ||\ min(p)>p_j]f(i)

设 F(j) 为质数处点值等于 f(j) 的完全积性函数，则有转移方程：

g(n,j)=
\begin{cases}
g(n,j-1)& {p_j^2>n}\cr g(n,j-1)-F(j)\left(g(\frac{n}{p_j},j-1)-g(p_j-1,j-1)\right)&{p_j^2\leq n}
\end{cases}

因此 g(n,0) 即为除去 1 后所有点值视为质数的前缀和，也即转移的开始；g(n,tot) （ tot 为转移的轮次，也即 n 以内质数数目）即为所求质数点值之和。采用滚动数组省略第二维。

筛质数个数

例题：loj6235 区间素数个数

求 1-n 之间素数个数，2\leq n\leq 10^{11}

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 1e6 + 50;
ll n;

int num[N], prime[N], tot;
void init() // 省略线性筛素数模版

struct Min25 {
    ll g[N << 1], id1[N << 1], id2[N << 1], val[N << 1];
    ll n, lim, sp[N];

    void init() {
        for (int i = 1; i <= tot; i++) {
            sp[i] = sp[i - 1] + 1;
        }
    }
    void build(ll _n) {
        n = _n; lim = sqrt(n);
        int m = 0;
        for(ll l = 1, r; l <= n; l = r + 1) {
            r = n / (n / l);
            val[++m] = n / l;
            if (val[m] <= lim) id1[val[m]] = m;
            else id2[n/val[m]] = m;
            g[m] = val[m] - 1;
        }
        for (int j = 1; j <= tot; j++) {
            for (int i = 1; 1LL * prime[j] * prime[j] <= val[i]; i++) {
                ll tmp = val[i] / prime[j];
                if (tmp <= lim) tmp = id1[tmp];
                else tmp = id2[n / tmp];
                g[i] -= g[tmp] - sp[j - 1];
            }
        }
    }
}Sn;

int main() {
    scanf("%lld", &n);
    init();
    Sn.init(); Sn.build(n);
    printf("%lld\n", Sn.g[1]);
    return 0;
}

g 的意义

有下面两道例题：

例题1：P6021 质数前缀统计

题目大意：

设 S(n) 为 n 以内质数的 k 次方和。求：

\sum_{i=1}^{\lfloor\sqrt n\rfloor}i^2S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \pmod p

0 \le k \le 10，1 \le N \le 4\times {10}^{10}，{10}^9

解题思路：

首先我们需要能够快速求出 \sum_{i=1}^{x}i^k。这个式子的通项可以通过拉格朗日差值得到。然后我们需要处理出 i^2 前缀和的通项。

题目的式子显然应当通过数论分块来做。在 min25 筛的过程中，我们注意到，以下几个变量的意义：

val[m]=n/l;

\begin{cases}
id1[n/l]=m,&n/l\le lim\cr id2[n/(n/l)]=m,&n/l>lim
\end{cases}

g[m]=(n/l\ 以内质数点值之和)

我们就可以很容易得到 n/l 以内质数点值之和。（这里有个卡常技巧，预处理出所有质数的 k 次方值 powp[i]！省略掉快速幂的一个log）

struct Min25 {
    ll g[N << 1], id1[N << 1], id2[N << 1], val[N << 1];
    ll n, lim, sp[N];

    void init() {
        for (int i = 1; i <= tot; i++) {
            sp[i] = (sp[i - 1] + powp[i]) % P;
        }
    }
    void build(ll _n) {
        n = _n; lim = sqrt(n);
        int m = 0;
        for(ll l = 1, r; l <= n; l = r + 1) {
            r = n / (n / l);
            val[++m] = n / l;
            if (val[m] <= lim) id1[val[m]] = m;
            else id2[n/val[m]] = m;
            g[m] = (qm(k, val[m]) - 1 + P) % P;
        }
        for (int j = 1; j <= tot; j++) {
            for (int i = 1; 1LL * prime[j] * prime[j] <= val[i]; i++) {
                ll tmp = val[i] / prime[j];
                if (tmp <= lim) tmp = id1[tmp];
                else tmp = id2[n/tmp];
                g[i] -= 1LL * powp[j] * (g[tmp] - sp[j-1] + P) % P;
                if(g[i] < 0) g[i] += P;
            }
        }
    }
    ll PrimeSum(ll l) { /* n / l 以内质数点值之和 */
        return n / l <= lim ? g[id1[n / l]] : g[id2[n / (n / l)]];
    }
}Sn;

int solve() {
    ll up = sqrt(n), ret = 0;
    for(ll l = 1, r; l <= up; l = r + 1) {
        r = n / (n / l);
        ll sz = sum(min(up, r)) - sum(l - 1) + P;
        ret += 1LL * sz * Sn.PrimeSum(l) % P;
    }
    return ret % P;
}

例题2：hdu6834. Yukikaze and Smooth numbers

第二部分-求积性函数前缀和

设sp[x]表示小于等于第x个质数处点值的前缀和。设S(n,x)表示求1到n中所有最小质因子大于p_x的函数值之和。答案就是S(n,0)。考虑转移，将满足条件的数分成两部分，第一部分是大于p_x的质数，也就是g(n)-sp[x]，另一部分是最小质因子大于p_x的合数，我们可以枚举最小质因子。转移方程：

S(n,x)=g(n)−sp_x+\sum_{p_k^e\le n\ \And\ k>x}f({p_k}^e)\left(S\left(\frac{n}{{p_k}^e},k\right)+[e\not = 1]\right)

最后不要忘记加上 1 处的点值。

下面分情况来讨论不同的多项式：

无系数多项式

洛谷P5325

f(p^k)=p^k(p^k-1)

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int P = 1e9 + 7;
const int N = 1e5 + 50;
const int inv6 = 166666668;
ll n;

int prime[N], num[N], tot;
int sp1[N << 1], sp2[N << 1];
void init() {
    for(int i = 2; i < N; i++) {
        if(!num[i]) {
            prime[++tot] = i;
        }
        for(int j = 1, x; j <= tot && (x = prime[j] * i) < N; j++) {
            num[x] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}

struct Min25 {
    ll val[N << 1], id1[N << 1], id2[N << 1], g1[N << 1], g2[N << 1];
    ll n, lim, sp1[N], sp2[N];
    // g1为i的转移滚动数组，g2为i^2的转移滚动数组
    // sp1为质数处i的前缀和，sp2为质数处i^2的前缀和

    void init() {
        for (int i = 1; i <= tot; i++) {
            sp1[i] = (sp1[i - 1] + prime[i]) % P;
            sp2[i] = (sp2[i - 1] + 1LL * prime[i] * prime[i]) % P;
        }
    }
    void build(ll _n) {
        n = _n; lim = sqrt(n);
        int m = 0;
        for(ll l = 1, r; l <= n; l = r + 1) {
            r = n / (n / l);
            val[++m] = n / l;
            if (val[m] <= lim) id1[val[m]] = m;
            else id2[n / val[m]] = m;
            g1[m] = g2[m] = val[m] % P;
            g1[m] = 1LL * g1[m] * (g1[m] + 1) / 2 % P - 1;
            // i的前缀和，别忘记刨除1
            g2[m] = g2[m] * (g2[m] + 1) % P * (2 * g2[m] + 1) % P * inv6 % P - 1;
            // i^2的前缀和，别忘记刨除1^2=1
        }
        for (int j = 1; j <= tot; j++) {
            for (int i = 1; 1LL * prime[j] * prime[j] <= val[i]; i++) {
                ll tmp = val[i] / prime[j];
                if (tmp <= lim) tmp = id1[tmp];
                else tmp = id2[n / tmp];
                g1[i] -= 1LL * prime[j] * (g1[tmp] - sp1[j - 1] + P) % P;
                g2[i] -= 1LL * prime[j] * prime[j] % P * (g2[tmp] - sp2[j - 1] + P) % P;
                if (g1[i] < 0) g1[i] += P;
                if (g2[i] < 0) g2[i] += P;
            }
        }
    }
    ll S(ll x, int y) {
        if (prime[y] > x || x == 1) return 0;
        ll tmp = (x <= lim ? id1[x] : id2[n/x]);
        ll ret = (g2[tmp] - g1[tmp] + P - (sp2[y] - sp1[y]) + P) % P;
        // 初始值：g - sp
        for (int i = y + 1; i <= tot && 1LL * prime[i] * prime[i] <= x; i++) {
            ll pe = prime[i];
            for(int e = 1; pe <= x; e++, pe *= prime[i]) {
                ll tmp = pe % P;
                (ret += (tmp * (tmp - 1) % P * (S(x / pe, i) + (e != 1)))) %= P;
            }
        }
        return ret;
    }
    int solve(ll n) {
        return (S(n, 0) + 1) % P;
    }
}Sn;

int main() {
    init();
    scanf("%lld", &n);
    Sn.init(); Sn.build(n);
    printf("%d\n", Sn.solve(n));
    return 0;
}

多项式含有系数

SP34096 DIVCNTK

题目大意：

求

S_3(n) = \sum_{i=1}^n \sigma_0(i^3)

其中 \sigma_0(x) 为 x 的约数个数。

题目分析：

设 f(x)=\sigma_0(x^k)，则有 f(p)=k+1，f(p^c)=kc+1。

#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long ll;
const int N = 1e5 + 50;
ll n, k, lim;

int num[N], prime[N], tot;
void init() {
    for(int i = 2; i < N; i++) {
        if(!num[i]) {
            prime[++tot] = i;
        }
        for(int j = 1, x; j <= tot && (x = prime[j] * i) < N; j++) {
            num[x] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}

ll g[N << 1], id1[N << 1], id2[N << 1], val[N << 1];
ll sp[N << 1];

struct Min25 {
    ll g[N << 1], id1[N << 1], id2[N << 1], val[N << 1];
    ll n, lim, sp[N];

    void init() {
        for (int i = 1; i <= tot; i++) {
            sp[i] = sp[i - 1] + (k + 1);
        }
    }
    void build(ll _n) {
        n = _n; lim = sqrt(n);
        int m = 0;
        for(ll l = 1, r; l <= n; l = r + 1) {
            r = n / (n / l);
            val[++m] = n / l;
            if(val[m] <= lim) id1[val[m]] = m;
            else id2[n / val[m]] = m;
            g[m] = (k + 1) * (val[m] - 1);
            // 这里减去1的点值为k+1！
        }
        for(int j = 1; j <= tot; j++) {
            for(int i = 1; 1LL * prime[j] * prime[j] <= val[i]; i++) {
                ll tmp = val[i] / prime[j];
                if (tmp <= lim) tmp = id1[tmp];
                else tmp = id2[n/tmp];
                g[i] -= (g[tmp] - sp[j - 1]);
                /* 特别注意！这里本来是两部分相乘，即(1 * (g-sp))，
                    其中1为函数值。保持答案为无系数结果再乘4，
                    由于前面的g和预处理的sp都已经乘上系数
                    所以前面的1不再乘上系数4！*/
            }
        }
    }
    ll S(ll x, int y) {
        if (x < prime[y] || x == 1) return 0;
        ll tmp = (x <= lim ? id1[x] : id2[n/x]);
        ll ret = g[tmp] - sp[y];
        for(int i = y + 1; i <= tot && 1LL * prime[i] * prime[i] <= x; i++) {
            ll pe = prime[i];
            for(int e = 1; pe <= x; e++, pe *= prime[i]) {
                ret += (k * e + 1) * (S(x / pe, i) + (e != 1));
            }
        }
        return ret;
    }

    ll solve(ll n) {
        return S(n, 0) + 1;
    }
}Sn;

int main() {
    init();
    int T; scanf("%d", &T);
    while(T--)
    {
        scanf("%lld %lld", &n, &k);
        Sn.init(); Sn.build(n);
        printf("%llu\n", Sn.solve(n));
    }
    return 0;
}

存在特殊点的情况

loj6053 简单的函数

题目大意：

求 \sum_{i=1}^nf(i)，其中 f(x) 有以下性质：

  1. $f(1)=1$
  2. $f(p^c)=p\oplus c$
  3. $f(ab)=f(a)f(b)$（a和b互质）

题目分析：

在质数处，有 p\oplus 1=p-1，但是还存在一个特殊点：2\oplus 1=3，使得其并非一个完全的多项式。但我们可以先暂时忽略这个特殊点，第一步的 build() 函数和 init() 中预处理 sp 数组均按照原来的方法做。在求 S(n,x) 的过程中，我们把它分成了两部分：第一部分是所有大于 p_x 的质数，另一部分是最小质因子大于 p_x 的合数。而求合数的值是一个递归的过程，因此我们只需妥善处理好质数的情况即可。

我们观察到，第一部分的值为 g[n]-sp[x]，我们发现当 x\geq 1 时，2 这个特殊点正好被抵消掉；当 x=0 时，sp[x]=0，应该等于 3 的点值被我们处理成了 1，因此ret+=2 即可。

ll S(ll x, int y) {
    if (x < prime[y] || x == 1) return 0;
    ll tmp = (x <= lim ? id1[x] : id2[n / x]);
    ll ret = ((g1[tmp] - g2[tmp] + P) - (sp[y] - y) + P) % P;
    if (y == 0) ret += 2;
    for (int i = y + 1; i <= tot && 1LL * prime[i] * prime[i] <= x; i++) {
        ll pe = prime[i];
        for (int e = 1; pe <= x; e++, pe *= prime[i]) {
            ret += (prime[i] ^ e) * (S(x / pe, i) + (e != 1)) % P;
        }
        ret %= P;
    }
    return ret;
}

筛 phi 和 mu

注意到 \varphi(p)=p-1，且 \varphi(p^c)=(p-1)p^{c-1} 均可快速求，可用 min25 筛；

注意到 \mu(p)=-1，且 \mu(p^c)=-1\cdot [c=1] 均可快速求，可用 min25 筛。

暂无评论

发送评论 编辑评论 


				

			

						

				

					

						

							

								
							

							
						

					

				

				

					

						

							

								
							

							
						

					

				

				

					

						

							

								
							

							
							
													

					

				

			

			
			

				

					
				

			

				

											

							
							
						

																				
					
											
						

	

		
|´・ω・)ノ
ヾ(≧∇≦*)ゝ
(☆ω☆)
（╯‵□′）╯︵┴─┴
￣﹃￣
(/ω＼)
∠( ᐛ 」∠)＿
(๑•̀ㅁ•́ฅ)
→_→
୧(๑•̀⌄•́๑)૭
٩(ˊᗜˋ*)و
(ノ°ο°)ノ
(´இ皿இ｀)
⌇●﹏●⌇
(ฅ´ω`ฅ)
(╯°A°)╯︵○○○
φ(￣∇￣o)
ヾ(´･ ･｀｡)ノ"
( ง ᵒ̌皿ᵒ̌)ง⁼³₌₃
(ó﹏ò｡)
Σ(っ °Д °;)っ
( ,,´･ω･)ﾉ"(´っω･｀｡)
╮(╯▽╰)╭
o(*////▽////*)q
＞﹏＜
( ๑´•ω•) "(ㆆᴗㆆ)
😂
😀
😅
😊
🙂
🙃
😌
😍
😘
😜
😝
😏
😒
🙄
😳
😡
😔
😫
😱
😭
💩
👻
🙌
🖕
👍
👫
👬
👭
🌚
🌝
🙈
💊
😶
🙏
🍦
🍉
😣
Source: github.com/k4yt3x/flowerhd
	

	

		
颜文字
Emoji
小恐龙
花!
	


									

			

			
			
		

	






上一篇
下一篇