zoukankan      html  css  js  c++  java
  • Min_25 筛学习笔记

    模板

    在本文中用 \(p\) 表示质数,\(P_i\) 表示第 \(i\) 个质数。

    题解

    处理整个函数的前缀和

    题目要我们求的是一个 积性函数 的前缀和。

    我们可以对于 \(n\) 内的每一个质因子进行单独考虑。但是这样子的复杂度显然过大,而这很大一部分原因是因为 \(n\) 以内的质数。

    因此如果我们把所有数分成质数和合数考虑,那么在合数中的不同最小质因子个数就是 \(\frac{\sqrt n}{\log n}\) 级别的了。

    接下来我们设函数 \(S\)\(S(n, j) = \sum\limits_{i = 2} ^ n f(i) (\min\limits_{p | i} p > P_j)\)

    对这个函数进行分类讨论。

    质数

    如果我们已经得到质数前缀和函数 \(q(x) = \sum\limits_{i = 1}^{n} f(i) (i \in Prime)\)

    那么 \(S(n,j)\) 在质数处的和就很好算,为 \(q(n) - \sum\limits_{i = 1}^{j} f(P_i)\) (简单容斥)。

    合数

    在合数处答案考虑枚举每一个数的最小质因数及他的次数,答案为 \(\sum\limits_{k > j} \sum\limits_{e=1}^{P_k^e \le n} f(P_k^e) (S(\frac{n}{P_k ^ e}, k) + \left[e>1\right])\)

    后面的 \(e>1\) 表示如果 \(k\) 次数不为 \(1\),那么就还要再加 \(f(P_k^e)\),因为 \(S\) 是不算为 \(1\) 的数。质数的 \(1\) 次方已经在前面枚举质数的地方算过了因此在 \(e = 1\) 时不再计算。

    合起来

    \[S(n, j) = q(n) - \sum\limits_{i = 1}^{j} f(P_i) + \sum\limits_{k>j} \sum\limits_{e=1}^{P_k^e \le n} f(P_k^e)[S(\frac{n}{P_k^e}, k) + (e>1)] \]

    最终求的答案是 \(S(n, 0) + f(1)\), 还要加上 \(f(1)\) 因为 \(S\) 函数不会计算到 \(1\)

    处理质数处的前缀和

    我们现在的目标是计算 \(q(n)\)

    完全积性函数 \(f'\) 在质数处取值与 \(f\) 相同。

    考虑一步一步把合数的贡献给消除。设一个函数 \(g\),满足 \(g(n,j)=\sum\limits_{i=2}^n f'(i)(i \in prime | \min\limits_{p|i} p > P_j)\)

    考虑计算 \(g(n, j-1) - g(n, j)\)

    \(P_j^2 > n\) 时, \(g(n, j - 1)\) 显然没有多余贡献, 即 \(g(n, j) = g(n, j-1)\)

    否则如果数 \(x\)\(g(n, j - 1)\) 中有多余贡献的数,\(x\) 一定有质因子 \(P_j\)。可以先把有质因子 \(P_j\) 的数都除以 \(P_j\)

    这些数剩下的所有质因子都大于等于 \(P_j\),因此是 \(g(\frac{n}{P_j}, j - 1) - \sum\limits_{i=1}^{j-1}f'(P_i)\) (注意 \(g(\frac{n}{P_j}, j-1)\) 会多算一些质数,因此后面需要减 \(\sum\limits_{i=1}^{j-1}f'(P_i)\))。

    因此得出状态转移方程:

    \[g(n, j)=\begin{cases}g(n, j-1) (P_j^2 > n)\\g(n, j-1) - f'(P_j)(g(\frac{n}{P_j}, j-1) - \sum\limits_{i=1}^{j-1}f'(P_i))(P_j^2 \le n)\end{cases} \]

    这样可以方便地计算 \(q(n) = g(n, |P|)\)

    \(g\) 函数需要的地方都是 \(g(\left\lfloor\frac{n}{k}\right\rfloor, ...)\) 的形式。

    众所周知 \(\left\lfloor\frac{n}{k}\right\rfloor\) 的取值只有 \(2\sqrt{n}\) 种情况。让 \(id1_k\) 记录 \(k\le\sqrt n\)\(id2_k\) 记录 \(\frac{n}{k} \le \sqrt n\),可以只枚举这 \(2\sqrt n\) 种情况。

    在处理质数处的完全积性函数时, \(p^k(p^k-1)\) 可以用 \(p^{2k}\)\(p^k\) 表示(尽管在合数的情况答案变了,但是在质数处没变,因此答案不变)。

    于是就可以为 \(g\) 函数赋初值(\(g(n, 0) = \sum\limits_{i = 2}^{n} f'(i)\))。


    总结一下,Min_25 筛能用的前提:质数处的 \(f(p)\) 值是关于 \(p\) 的多项式,质数次方处的 \(f(p^e)\) 值可以快速计算。

    代码

    #include<bits/stdc++.h>
    #define L(i, j, k) for(int i = j, i##E = k; i <= i##E; i++)
    #define R(i, j, k) for(int i = j, i##E = k; i >= i##E; i--)
    #define ll long long
    #define ull unsigned long long
    #define db double
    #define pii pair<int, int>
    #define mkp make_pair
    #define ll long long
    using namespace std;
    const int N = 2e5 + 7;
    const int mod = 1e9 + 7;
    bool Prime[N];
    int tot;
    ll p[N], n, sq, w[N], m, sum1[N], sum2[N], g1[N], g2[N], id1[N], id2[N];
    // id1 id2 如文中所示 w[k] 为第 k 个需要处理的数 
    // 文中的 g 为 g2 - g1, g1 : f'(k) = k, g2: f'(k) = k ^ 2
    // sum1 : 质数前缀和 ; sum2: 质数平方前缀和
    void xxs() { // 线性筛
    	L(i, 2, sq)  {
    		if(!Prime[i]) ++tot, p[tot] = i;
    		for(int j = 1; p[j] * i <= sq && j <= tot; j++) {
    			Prime[p[j] * i] = 1;
    			if(i % p[j] == 0) break;
    		}
    	}
    	L(i, 1, tot) {
    		sum1[i] = (sum1[i - 1] + p[i]) % mod;
    		sum2[i] = (sum2[i - 1] + p[i] * p[i] % mod) % mod;
    	}
    }
    ll getid(ll x) {
    	if(x <= sq) return id1[x];
    	else return id2[n / x];
    }
    ll f1(ll x) { // 1 ~ x 前缀和
    	x %= mod;
    	return x * (x + 1) / 2 % mod;
    }
    ll f2(ll x) { // 1 ~ x 前缀平方和
    	x %= mod;
    	return x * (x + 1) % mod * (2 * x % mod + 1) % mod * 166666668 % mod;
    }
    ll S(ll x, int j) {
    	if(p[j] > x) return 0;
    	ll Ans = ((g2[getid(x)] - g1[getid(x)] + mod) % mod - (sum2[j] - sum1[j] + mod) % mod + mod) % mod;
    	for(int i = j + 1; i <= tot && p[i] * p[i] <= x; i ++) 
    		for(ll e = 1, sp = p[i]; sp <= x; sp *= p[i], e ++) 
    			Ans = (Ans + sp % mod * (sp % mod - 1) % mod * (S(x / sp, i) + (e > 1)) % mod) % mod;
    	return Ans;
    }
    int main() {
    	scanf("%lld", &n), sq = sqrt(n), xxs();
    	for(ll l = 1, r; l <= n; l = r + 1) {
    		r = (n / (n / l)), w[++m] = n / l;
    		g1[m] = f1(w[m]) - 1, g2[m] = f2(w[m]) - 1; // 处理 g(?, 0)
    		if(w[m] <= sq) id1[w[m]] = m;
    		else id2[n / w[m]] = m;
    	}
    	L(i, 1, tot) { // dp 处理 g 函数
    		for(int j = 1; j <= m && p[i] * p[i] <= w[j]; j++) {
    			g1[j] = (g1[j] - p[i] * (g1[getid(w[j] / p[i])] - sum1[i - 1]) % mod + mod) % mod;
    			g2[j] = (g2[j] - p[i] * p[i] % mod * (g2[getid(w[j] / p[i])] - sum2[i - 1]) % mod + mod) % mod;
    		}
    	}
    	printf("%lld\n", (S(n, 0) + 1) % mod);
    	return 0;
    }
    

    祝大家学习愉快!

  • 相关阅读:
    webpy使用笔记(一)
    如何衡量离散程度
    Hash哈希(二)一致性Hash(C++实现)
    Hash哈希(一)
    sqlmap使用笔记
    Windows7 IIS7.5 HTTP Error 503 The service is unavailable 另类解决方案
    [转]IP动态切换脚本
    全国各地电信DNS服务器地址
    比较好的汉字拼音化类
    c#读取INI文件类
  • 原文地址:https://www.cnblogs.com/zkyJuruo/p/13445596.html
Copyright © 2011-2022 走看看