闲话 22.12.16

闲话

字符串能不能滚出 OI 啊?😅

所以今天写点小清新的杜教筛。
凭借已经做过 6/9 的历史,我成功在一晚上切穿了杜教筛专题,快来膜我!

昨天闲话里发生了很奇妙的事情
我不知道为什么
这里 @ 一下 jjdw 是有意义的吧(

杜教筛

发现自己还没写过杜教筛的内容,那就权当帮同学预习吧(

首先这玩意的前置知识是 Dirichlet 卷积,记这个复合为 \(*\),正常的点值乘法为 \(\times\)。其实推导式子是不用莫反的。

假设我们有一个积性函数 \(f\),定义它的前缀和函数为 \(S(f, n) = \sum_{i=1}^n f(i)\)。现在需要解决的是快速算出 \(S(f, n)\)。由 duls 总结的杜教筛对一部分具有特殊性质的 \(f\) 提供了一种 \(O(n^{2/3})\) 时间复杂度的解决方案。

具体地,我们构造 \(g,h\text{ s.t. } h = f * g\)。假设我们容易对给定的多个 \(k \le n \text{ s.t. } \exists x, k = \left\lfloor \frac nx \right\rfloor\) 求得 \(S(g, k)\)\(S(h, k)\),我们就可以应用杜教筛。杜教筛的推导如下:

\[\begin{aligned} S(h, n) & = \sum_{x=1}^n h(x) \\ & = \sum_{x = 1}^n \sum_{d|x} g(d) f(\frac xd) \\ & = \sum_{d = 1}^n g(d) \sum_{x = 1}^{\left\lfloor \frac n d \right\rfloor} f(x) \\ & = \sum_{d = 1}^n g(d) S(f, \left\lfloor \frac n d \right\rfloor) \\ & = g(1) S(f, n) + \sum_{d = 2}^n g(d) S(f, \left\lfloor \frac n d \right\rfloor) \\ g(1) S(f, n) & = S(h, n) - \sum_{d = 2}^n g(d) S(f, \left\lfloor \frac n d \right\rfloor) \\ S(f, n) & = \frac{1}{g(1)} \left( S(h, n) - \sum_{d = 2}^n g(d) S(f, \left\lfloor \frac n d \right\rfloor) \right) \end{aligned}\]

所以就可以按照上面的方式递归了!(一般 \(g(1)\)\(1\) 所以可以忽略)

那复杂度是多少呢?应用分析 \(Min25\) 筛时用到的渐进方式可以得到

\[\begin{aligned} & \sum_{i=1}^{\sqrt n} O(\sqrt i) + \sum_{i=1}^{\sqrt n} O(\sqrt \frac ni) \\ = \ & O\left(\int_{1}^{\sqrt n} x^{\frac 12} \text dx\right) + O\left(\int_{1}^{\sqrt n} (\frac n x)^{\frac 12} \text dx\right) \\ = \ & O\left(n^{3/4} \right) + O\left(n^{3/4} \right) \end{aligned}\]

因此不预处理的复杂度是 \(O\left(n^{3/4} \right)\) 的。这里默认记忆化了。

欸你刚才不是说 \(2/3\) 吗?咋变慢了?
别着急,现在考虑预处理。我们预处理前 \(B\) 个位置的值,且有 \(B \ge \sqrt n\),复杂度就变成了

\[\begin{aligned} \sum_{i=B}^{\sqrt n} O(\sqrt \frac ni) = O\left(\frac{n}{\sqrt B}\right) \end{aligned}\]

\(B = n^{2/3}\) 得到总时间复杂度 \(O(n^{2/3})\)

总结一下,杜教筛是一类特殊的递归方案。我们首先选取两个易求前缀和的函数,随后将所求函数 \(f\) 的前 \(n^{2/3}\) 个点值求出来,做前缀和,最后记忆化递归,在过程中模拟如上的式子就能得到最终答案了。
还发现,由于杜教筛的记忆化性质,其可以在一次求解中得到 \(n\) 对应 \(O(\sqrt n)\) 个数论分块点值处的前缀和。因此其时间复杂度应当均摊分析,用于数论分块中时总复杂度没有退化。
重复式子:构造 \(g,h\text{ s.t. } h = f * g\)

\[S(f, n) = \frac{1}{g(1)} \left( S(h, n) - \sum_{d = 2}^n g(d) S(f, \left\lfloor \frac n d \right\rfloor) \right) \]

给出几个常见的数论函数:

  • 欧拉函数 \(\varphi(n) = \sum_{(d,n) =1} 1\)
  • 莫比乌斯函数 \(\mu(n)\)
  • 恒等函数 \(\text I(n) = 1\)
  • 元函数 \(e(n) = [n = 1]\)
  • 单位函数 \(\text{id} (n) = n\)
  • 约数个数函数 \(\text{d} (n) = \sum_{d | n} 1\)
  • 题面里给的函数 \(f(n) = \small{我怎么知道}\)

这里记 \(k\) 个数论函数 \(f\) 卷起来的函数为 \(f_k\)

给出几个常见的构造:

  • $\mu * \text I = e \ \to \ f = (f * \text I) * \mu \ $ 该变换即莫比乌斯反演。
  • \(\varphi * \text I = \text {id}\ \to \ \text{id} * \mu = \varphi\)
  • \(\text{id} * \text{id} = \text{id} \times \text d\)

结合例题来理解吧。

例题

首先是板子。

【模板】杜教筛

给定正整数 \(n\),求 \(S(\varphi, n)\)\(S(\mu, n)\)

\(n \le 10^{11}\)

照上式模拟即可。具体看代码。

code
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define rep(i,a,b) for ( register int (i) = (a); (i) <= (b); ++(i) )
#define pre(i,a,b) for ( register int (i) = (a); (i) >= (b); --(i) )
using namespace std;
const int N = 5e6 + 10, L = 5e6;
int T, n;

int prime[N], cnt, phi[N], mu[N];
inline void sieve() {
    phi[1] = mu[1] = 1;
    rep(i,2,L) {
        if (phi[i] == 0) prime[++cnt] = i, phi[i] = i - 1, mu[i] = - 1;
        rep(j,1,cnt) {
            int tmp = i * prime[j];
            if (tmp > L) break;
            if (i % prime[j] == 0) {
                phi[tmp] = prime[j] * phi[i]; mu[tmp] = 0;
                break;
            } else 
                phi[tmp] = phi[prime[j]] * phi[i]; mu[tmp] = -mu[i];
        }
    } 
	rep(i,2,L) phi[i] = phi[i-1] + phi[i], mu[i] = mu[i-1] + mu[i];
}

int _phi[N], _mu[N];
inline int getPhi(int x) {
	if (x <= L) return phi[x];
	if (~_phi[n / x]) return _phi[n / x];
	int ret = 1ll * x * (x + 1) >> 1;
	for (int l = 2, r; l <= x; l = r+1) {
		r = x / ( x / l );
		ret -= (r - l + 1) * getPhi(x / l);
	} 
	return _phi[n / x] = ret;
}

inline int getMu(int x) {
	if (x <= L) return mu[x];
	if (~_mu[n / x]) return _mu[n / x];
	int ret = 1;
	for (int l = 2, r; l <= x; l = r + 1) {
		r = x / ( x / l );
		ret -= (r - l + 1) * getMu(x / l);
	} 
	return _mu[n / x] = ret;
}

signed main() {
	sieve(); cin >> T;
	while (T--) {
		cin >> n;
		memset(_phi, -1, sizeof _phi);
		memset(_mu, -1, sizeof _mu);
		cout << getPhi(n) << ' ' << getMu(n) << '\n';
	} 
}



简单的最大公约数

给定 \(n,k\),求

\[\sum^k_{a_1=1}\sum^k_{a_2=1}\sum^k_{a_3=1}\dots\sum^k_{a_n=1}\gcd(a_1,a_2,a_3,\dots,a_n)\pmod{ 1000000007} \]

\(n,k \le 10^{11}\)

这里是早就写好了的题解



DZY Loves Math IV

给定 \(n,m\),求

\[\sum_{i=1}^n \sum_{j=1}^m \varphi(ij) \pmod{1000000007} \]

\(n \le 10^5, m \le 10^9\)

考虑 \(n\) 比较小,我们可以对每个 \(n\) 计算答案后求和。

现在的问题是求解

\[S(n, m) = \sum_{i=1}^m \varphi(in) \]

假设 \(n = pq\),其中 \(p\) 是将 \(n\) 的所有质因子拿出一个来求积的得数。我们有

\[\begin{aligned} S(n, m) & \ = \sum_{i=1}^m \varphi(in) \\ & \ = q\sum_{i=1}^m \varphi(ip) \\ & \ = q\sum_{i=1}^m \varphi(\frac{p}{\gcd(i, p)}) \varphi(i) \gcd(i, p) \\ & \ = q\sum_{i=1}^m \varphi(\frac{p}{\gcd(i, p)}) \varphi(i) \sum_{d|(i, p)} \varphi(d) \qquad & (\varphi * \text{I} = \text{id} ) \\ & \ = q\sum_{i=1}^m \varphi(i) \sum_{d|(i, p)} \varphi\left(\frac{dp}{\gcd(i, p)} \right)\qquad & (\varphi 是积性的) \\ & \ = q\sum_{i=1}^m \varphi(i) \sum_{d|(i, p)} \varphi\left(\frac{p}{d} \right) \\ & \ = q\sum_{d|p} \varphi\left(\frac{p}{d} \right)\sum_{i=1}^{m / d} \varphi(id) \\ & \ = q\sum_{d|p} \varphi\left(\frac{p}{d} \right) S(d, \left\lfloor\frac md\right\rfloor) \end{aligned}\]

递归计算即可。边界不再赘述。

由于最终是 \(\left\lfloor\frac md\right\rfloor\) 的形式,所以第二维的取值数是 \(O(n\sqrt m)\) 的。每一层要枚举的总状态数是\(O(\sqrt n)\) 的,因此这一部分的复杂度是 \(O(n(\sqrt n + \sqrt m))\)。考虑杜教筛的记忆化形式,边界本质上只需要进行一次杜教筛即可全部求出来。因此总时间复杂度为 \(O(n(\sqrt n + \sqrt m) + m^{2/3})\)
理论复杂度能跑 \(10^{11}\),本机 O2 加持下跑了 1.61s。

Submission



神犇和蒟蒻

给定 \(n,m\),求

\[\sum_{i=1}^n \mu(i^2) \qquad \sum_{i=1}^n \varphi(i^2) \]

\(n \le 10^9\)

没啥,放这题就是整点水清凉一下(

\(\mu\) 那个显然答案就是 \(1\)
\(\varphi\)……哦想起来了,放这题是因为想告诉你们 \(\varphi(n^k) = n^{k-1}\varphi(n)\)。挺好用的。

\(f(n) = \varphi(n) \times n\)\(f * \text{id} = \text{id}_2\)。剩下的看代码。

Submission



Lucas的数论

给定 \(n\),求

\[\sum_{i=1}^n \sum_{j=1}^n \text d(ij) \pmod{1000000007} \]

\(n \le 10^9\)

还记得 约数个数和 那题吗?
直接套那个式子

\[\begin{aligned} &\sum_{i=1}^n \sum_{j=1}^n \text d(ij) \\ = & \ \sum_{i=1}^n \sum_{j=1}^n \sum_{x|i} \sum_{y|j} [gcd(x, y) = 1] \\ = & \ \sum_{i=1}^n \sum_{j=1}^n \sum_{x|i} \sum_{y|j} \sum_{d|(x, y)}\mu(d) \\ = & \ \sum_{x = 1} \sum_{y = 1} \sum_{d|(x, y)}\mu(d) \left\lfloor \frac n{x} \right\rfloor \left\lfloor \frac n{y} \right\rfloor \\ = & \ \sum_{d = 1}\mu(d) \sum_{x = 1} \left\lfloor \frac n{xd} \right\rfloor \sum_{y = 1} \left\lfloor \frac n{yd} \right\rfloor \end{aligned}\]

然后对 \(d\) 进行数论分块,对 \(\mu\) 采用杜教筛,对 \(x,y\) 采用二次数论分块。

根据杨卓凡第一小定理(?)可知总时间复杂度为 \(O(n^{3/4})\)

Submission


留着别的题明天写吧。


upd: 对杜教筛的巨大常数优化:非递归

模板题实现

Submission.

#include <bits/stdc++.h>
using namespace std;
#define inline __attribute__((__gnu_inline__, __always_inline__, __artificial__)) inline
using pii = pair<int,int>; using vi = vector<int>; using vp = vector<pii>; using ll = long long; 
using ull = unsigned long long; using db = double; using ld = long double; using lll = __int128_t;
template<typename T1, typename T2> T1 max(T1 a, T2 b) { return a > b ? a : b; }
template<typename T1, typename T2> T1 min(T1 a, T2 b) { return a < b ? a : b; }
#define multi int _T_; cin >> _T_; for (int TestNo = 1; TestNo <= _T_; ++ TestNo)
#define timer cerr << 1. * clock() / CLOCKS_PER_SEC << '\n';
#define iot ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define rep(i,s,t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i,s,t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
#define eb emplace_back
#define pb pop_back
const int N = 4e6 + 10;
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;
ll n, th, mp[N], lp[N], ilp[N], lv[N], cntl;
ll phimp[N], mump[N];
int T[N], cT;
db iv[N];

ll phi[N], mu[N], pri[N], cntp;
void sieve(int N) {
    phi[1] = mu[1] = 1;
    rep(i,2,N) {
        if (!phi[i]) pri[++ cntp] = i, phi[i] = i - 1, mu[i] = -1;
        rep(j,1,cntp) {
            if (i * pri[j] > N) break;
            if (i % pri[j] == 0) {
                phi[i * pri[j]] = phi[i] * pri[j], mu[i * pri[j]] = 0;
                break;
            } else {
                phi[i * pri[j]] = phi[i] * (pri[j] - 1), mu[i * pri[j]] = - mu[i];
            }
        }
    }
    rep(i,1,N) phi[i] += phi[i - 1], mu[i] += mu[i - 1];
}

#define Sid(n) (1ll * (n) * ((n) + 1) / 2)
#define SI(n) (n)
#define Se(n) (1)

signed main() {
    cin >> cT;
    rep(i,1,cT) cin >> T[i], n = max(n, T[i]);
    th = min(N - 5, (int)pow(n, 2./3));
    sieve(th);
    rep(i,1,(int)sqrt(n)) iv[i] = 1. / i;
    // cerr << (ll)(1e10) / 10000 << ' ' << iv[10000] << ' ' << (int)((ll)(1e10) * iv[10000]) << endl;
    rep(i,1,cT) {
        n = T[i];
        if (n <= th) {
            cout << phi[n] << ' ' << mu[n] << '\n';
            continue;
        }
        cntl = 0;
        for (ll l = 1, r; n / l > th; l = r + 1) {
            r = n / (n / l);
            mp[lp[++ cntl] = l] = n / l;
            rep(i,l,r) lv[i] = l;
        }
        pre(i,cntl,1) {
            ll nc = mp[lp[i]], sqnc = (int)sqrt(nc);
            phimp[i] = Sid(nc) + 1ll * SI(sqnc) * phi[sqnc];
            mump[i] = Se(nc) + 1ll * SI(sqnc) * mu[sqnc];
            rep(j,1,sqnc) {
                int it = int(nc * iv[j] + 1e-8); // 警钟敲烂:如果 nc mod j = 0 则需要 + eps 防丢精
                phimp[i] -= 1ll * (phi[j] - phi[j - 1]) * SI(it);
                mump[i] -= 1ll * (mu[j] - mu[j - 1]) * SI(it);
                if (j > 1) if (it <= th) {
                    phimp[i] -= phi[it];
                    mump[i] -= mu[it];
                } else {
                    phimp[i] -= phimp[lv[lp[i] * j]];
                    mump[i] -= mump[lv[lp[i] * j]];
                }
            }
            // cerr << "Partial answer : " << nc << ' ' << sqnc << ' ' << phimp[i] << ' ' << mump[i] << endl;
        } 

        cout << phimp[1] << ' ' << mump[1] << '\n';
    }
}
posted @ 2022-12-16 21:09  joke3579  阅读(165)  评论(1编辑  收藏  举报