【模板】min25 筛
使用条件:
- 求积性函数单点前缀和。
- 若 \(p\) 为质数,则 \(f(p)\) 是低次多项式(或者可以拆成常数个完全积性函数,同时它们的前缀和都可以计算)
- 若 \(p\) 为质数,则要求 \(f(p^c)\) 能在比较低的时间计算。
- \(n\leq 10^{10}\)(?)
算法流程建议直接看 OI-wiki 而不是我在这里复读一遍。你就找个时间认真地逐行读一下,注意一定是认真地读,搞清楚每个等号的由来。
代码中除了 \(idl, idg\) 记录了 \(x, n/x\) 对应的 id 之外还记了一个 \(idp_k=id(pri_k)\)。筛质数点值前缀和对我的代码负优化,没加,或者可能我写的有问题。\(id\) 的写法也很牛,你需要提前处理根号,不能做太多无用的除法。
点击查看代码
本题 \(f(p^k)=p^k(p^k-1)\)
#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define endl "\n"
#define debug(...) void(0)
#endif
using LL = long long;
template <class T>
using must_int = enable_if_t<is_integral<T>::value, int>;
template <unsigned umod>
struct modint {/*{{{*/
static constexpr int mod = umod;
unsigned v;
modint() : v(0) {}
template <class T, must_int<T> = 0>
modint(T x) {
x %= mod;
v = x < 0 ? x + mod : x;
}
modint operator+() const { return *this; }
modint operator-() const { return modint() - *this; }
friend int raw(const modint &self) { return self.v; }
friend ostream& operator<<(ostream& os, const modint &self) {
return os << raw(self);
}
modint &operator+=(const modint &rhs) {
v += rhs.v;
if (v >= umod) v -= umod;
return *this;
}
modint &operator-=(const modint &rhs) {
v -= rhs.v;
if (v >= umod) v += umod;
return *this;
}
modint &operator*=(const modint &rhs) {
v = 1ull * v * rhs.v % umod;
return *this;
}
modint &operator/=(const modint &rhs) {
assert(rhs.v);
return *this *= qpow(rhs, mod - 2);
}
template <class T, must_int<T> = 0>
friend modint qpow(modint a, T b) {
modint r = 1;
for (; b; b >>= 1, a *= a)
if (b & 1) r *= a;
return r;
}
friend modint operator+(modint lhs, const modint &rhs) { return lhs += rhs; }
friend modint operator-(modint lhs, const modint &rhs) { return lhs -= rhs; }
friend modint operator*(modint lhs, const modint &rhs) { return lhs *= rhs; }
friend modint operator/(modint lhs, const modint &rhs) { return lhs /= rhs; }
bool operator==(const modint &rhs) const { return v == rhs.v; }
bool operator!=(const modint &rhs) const { return v != rhs.v; }
};/*}}}*/
using mint = modint<int(1e9 + 7)>;
LL sqr(LL x) { return x * x; }
template <int N>
struct siever {
int pri[N + 10], cnt;
siever() : cnt(0) {
fill(pri + 1, pri + N + 1, true);
pri[0] = 1, pri[1] = false;
for (int i = 1; i <= N; i++) {
if (pri[i]) pri[++cnt] = i;
for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
pri[i * pri[j]] = false;
if (i % pri[j] == 0) break;
}
}
}
};
mint ID1(mint x) { return x * (x + 1) * 500000004; }
mint ID2(mint x) { return x * (x + 1) * (2 * x + 1) * 166666668; }
LL n, v[200010], lim;
mint f(mint pc) { return pc * (pc - 1); }
int m, idl[100010], idg[100010], idp[100010];
int id(LL x) { return x <= lim ? idl[x] : idg[n / x]; }
mint g1[200010], g2[200010], *g = g2;
siever<100010> sie;
mint dfs(LL n, int k) {
if (n < sie.pri[k]) return 0;
mint ans = g[id(n)] - g[idp[k - 1]];
for (int i = k; sqr(sie.pri[i]) <= n; i++) {
int pi = sie.pri[i];
LL now = pi;
for (int c = 1; now * pi <= n; c++, now *= pi) {
ans += f(now) * dfs(n / now, i + 1) + f(now * pi);
}
}
return ans;
}
int main() {
#ifndef LOCAL
cin.tie(nullptr)->sync_with_stdio(false);
#endif
cin >> n;
lim = sqrtl(n);
for (LL l = 1, r; l <= n; l = r + 1) {
v[++m] = n / l, r = n / v[m];
(v[m] <= r ? idl[v[m]] : idg[r]) = m;
g1[m] = ID1(v[m]) - 1, g2[m] = ID2(v[m]) - 1;
}
idp[0] = id(1);
for (int k = 1; sqr(sie.pri[k]) <= n; k++) {
idp[k] = id(sie.pri[k]);
int p = sie.pri[k];
for (int i = 1; i <= m && v[i] >= sqr(p); i++) {
g1[i] -= p * (g1[id(v[i] / p)] - g1[idp[k - 1]]);
g2[i] -= sqr(p) * (g2[id(v[i] / p)] - g2[idp[k - 1]]);
}
}
for (int i = 1; i <= m; i++) g[i] -= g1[i];
cout << dfs(n, 1) + 1 << endl;
return 0;
}
LOJ6053【简单的函数】:这道题就是说 \(p\) 为质数时 \(f(p)=p\oplus 1=p+(-1)^{[p\neq 2]}\),就是说 \(f(2)=3\) 而其它 \(f(p)=p-1\),\(f(2)\) 要特判。这个特判要判的东西其实很少,首先 \(g\) 的部分根本不需要管,然后最后用 \(g\) 的时候 mint ans = g[id(n)] - g[idp[k - 1]];
再把 \(2\) 的多出来一点点的答案加进去,实际上根本没啥要改的。
可以认为就是:\(g\) 部分就是算 lambda x: x ** k
这个函数的质数点处前缀和,然后你把它组合一下让它摇身一变变成原函数的质数点处前缀和。然后再暴力 dfs,去把每个数搜出来,按照质因子从小到大的顺序,对于质数用刚才的 \(g\) 直接前缀和求出,剩下的合数让它递归去。
本文来自博客园,作者:caijianhong,转载请注明原文链接:https://www.cnblogs.com/caijianhong/p/18407205/template-min25