再探欧式筛——一种泛用性更强的欧拉筛法/线性筛法实现
一、引言
欧式筛/欧拉筛法/线性筛法（Euler Sieve）是一种能够在 \(O(n)\) 时间复杂度内，处理 \([1,n]\) 内质数的方法。
其相比埃氏筛/埃拉托斯特尼筛法（Eratosthenes Sieve）的 \(O(n\log\log n)\) 时间复杂度，主要的优化在于欧式筛保证了所有正整数 \(n\) 均只被其最小质因数 \({minp}_n\) 筛到，从而保证了线性复杂度。
此外，欧式筛同时还具有线性时间内预处理欧拉函数 \(\boldsymbol \varphi(n)\)、莫比乌斯函数 \(\boldsymbol \mu(n)\) 等简单积性函数的能力。但对于除数个数函数 \(\boldsymbol d(n)\) （或记为 \(\boldsymbol {\sigma_0}(n)\)）、除数和函数 \(\boldsymbol \sigma(n)\) 等比较复杂的积性函数，则实现起来比较复杂。
本文提出了一个能处理上述复杂函数，泛用性更强的一种欧拉筛实现；并通过实验证明，这个实现的效率也非常高。
二、相关基础内容
2.1 筛法
我们称，一个只有 \(1\) 和其本身是自己因数的正整数 \(n\) 为质数。一个比较简单的求解 \([1,n]\) 范围内质数的方式为：
对已知的质数，划去所有其倍数。未被划去的数即为质数。参考代码如下：
bool nprime[MAXN];
for(int i=2; i<=n; ++i) if(!nprime[i])
for(int j=i+i; j<=n; j+=i)
nprime[j]=1;
此即为埃氏筛法，时间复杂度可证明为 \(O(n\log \log n)\)。
为保证时间复杂度为线性，我们需要保证每个数 \(n\) 仅被自身的最小质因数 \({minp}_n\) 筛到。当然，这个描述的实现方法就和动态规划中的填表法一样，不容易理解；我们也可以参考动态规划的思想，转化为刷表法实现。参考代码如下：
vector<int> pr;
bool nprime[MAXN];
for(int i=2; i<=n; ++i) {
if(!nprime[i])
pr.push_back(i);
for(int p : pr)
if(p*i>n)
break;
else {
nprime[p*i]=1;
if(i%p==0) //再大的就不是枚举最小质因数了
break;
}
}
于是，我们就实现了 \(O(n)\) 的欧式筛法。
2.2 积性函数和狄利克雷卷积
我们称，对任意满足 \(\gcd(a,b)=1\) 的 \(a,b\)，都满足 \(\boldsymbol f(ab)=\boldsymbol f(a)\boldsymbol f(b)\) 的函数为积性函数，用粗体表示。
很显然，根据积性函数的性质可以有：
其中 \(p_i\) 是互不相同的质数。因此，一个比较简单的描述积性函数的方法就是，直接描述它在质数幂次的结果
通过上式可以很简单地推出，积性函数 \(\boldsymbol f(n)\) 和积性函数 \(\boldsymbol g(n)\) 的点积 \(\boldsymbol f(n)\cdot \boldsymbol g(n)\) 也是积性函数。为此，我们在后文将直接使用 \((\boldsymbol f\cdot \boldsymbol g)(n)\) 来表示这种点积后的积性函数。
此外，我们还能定义狄利克雷卷积运算，称积性函数 \(\boldsymbol f(n)\) 和积性函数 \(\boldsymbol g(n)\) 的狄利克雷卷积为：
通过定义式也可以证明，狄利克雷卷积也是积性函数。
2.3 简单积性函数的线性筛法
我们定义积性函数：欧拉函数 \(\boldsymbol \varphi(n)\) 表示 \([1,n]\) 中与 \(n\) 互质的正整数个数。
此篇文章 中证明了，欧拉函数在质数幂次 \(p^e(e>0)\) 处的结果为：\(\boldsymbol \varphi(p^e)=p^{e-1}(p-1)\)。
此外，另一个比较常见的积性函数是莫比乌斯函数 \(\boldsymbol \mu(n)\)，其在质数幂次 \(p^e(e>0)\) 处的结果为：
这两个函数的一个共性是：其在质数幂次处的结果可以很容易地整体乘法递推。我们假设质数 \(p<{minp}_n,e>0\)，则从已知 \(\boldsymbol f(p^e\cdot n)\) 递推出 \(\boldsymbol f(p^{e+1}\cdot n)\) 是很容易的：
因此，我们可以拓展一下欧式筛，在刷表的时候，确认一下当前枚举的质数是不是自己的最小质因数，即可根据积性性质或上述递推性质分别转移了。
vector<int> pr;
int phi[MAXN], mu[MAXN];
bool nprime[MAXN];
for(int i=2; i<=n; ++i) {
phi[1]=1;
mu[1]=1;
if(!nprime[i]) {
pr.push_back(i);
phi[i]=i-1;
mu[i]=-1;
}
for(int p : pr)
if(p*i>n)
break;
else {
nprime[p*i]=1;
if(i%p==0) {//递推转移
phi[p*i]=p*phi[i];
mu[p*i]=0;
break;
}
else {//积性转移
phi[p*i]=phi[p]*phi[i];
mu[p*i]=mu[p]*mu[i];
}
}
}
这个方法的使用很广泛，然而当积性函数不满足可乘法递推性质时，这个实现会使得递推转移部分的代码写得非常复杂；甚至若维护不好，可能还会在维护过程中失去线性时间复杂度。
例如，对于除数个数函数，其在质数幂次 \(p^e(e>0)\) 处的结果为 \(\boldsymbol {\sigma_0}(p^e)=e+1\)。若强行维护，则 \(\boldsymbol {\sigma_0}(n\cdot p^{e+1})={e+2\over e+1}\cdot \boldsymbol {\sigma_0}(n\cdot p^e)\)。若想到了维护最小质因数的幂次数组，还尚且能 \(O(n)\) 转移，否则还要在每次花费 \(O(e)\) 的复杂度处理 \({e+2\over e+1}\)，时间复杂度为：
其中 \(\text{Prime}\) 表示 \([1,n]\) 中的质因数集合。
很显然，这个实现已经不够理想了，若是更复杂的积性函数 \(\displaystyle \boldsymbol {\sigma_k}(p^e)=\sum_{i=0}^e p^{ki}\)，维护起来只会更加复杂。
三、再探欧拉筛
3.1 优化转移
上述欧拉筛造成积性函数只方便递推转移的关键在于，我们在递推过程中只能方便地表示每个数的最小质因数。而在大部分的积性函数中，我们要知道 \(\boldsymbol f(p^e\cdot n)(p<{minp}_n,e>0)\)，是需要知道最小质因数的幂次 \(p^e\)。
这就启示我们能不能额外维护最小质因数的幂次 \({pe}_n\)，这样我们即可通过查表快速计算积性函数 \(\boldsymbol f(n)=\boldsymbol f({pe}_n)\cdot \boldsymbol f({n\over {pe_n}})\)。
事实上，\({pe}_n\) 是很好处理的。我们在枚举不超过 \({minp}_n\) 的质数 \(p\) 时，若 \(p<{minp}_n\)，则 \({pe}_{n\cdot p}=p\)，否则 \({pe}_{n\cdot p}={pe}_n\cdot p\)；这在欧式筛中是很容易实现的。
并且，由于维护了 \({pe}_n\)，我们可以不必再维护 \({nprime}_n\) 来判断一个数是否不是质数。参考代码如下：
vector<int> pr;
//...要筛的的积性函数
int pe[MAXN];
//f[1]=1;
for(int i=2, x; i<=n; ++i) {
if(!pe[i]) {
pr.push_back(i);
pe[i]=i;
//...要筛的积性函数
}
for(int p : pr)
if((x=p*i)>n)//记忆 x=p*i，避免后续重复乘法
break;
else {
if(i%p==0) {
pe[x]=pe[i]*p;
//f[x]=f[x]*f[x/pe[x]]
break;
}
else {//积性转移
pe[x]=p;
//f[x]=f[i]*f[p];
}
}
}
3.2 质数部分的处理
若仔细思考上述代码，很快就会发现一个问题：
虽然 \(\boldsymbol f(n)=\boldsymbol f({pe}_n)\cdot \boldsymbol f({n\over {pe_n}})\) 在大部分情况下都是对的，但按照该写法，\(n=p^e\) 就会有：
但是，这个转移式并不能求出在质数幂次时的积性函数值！
因此，关于质数幂次处的积性函数值，必须在筛到质数的时候同步处理：当我们枚举到质数 \(p\) 时，我们枚举 \([1,n]\) 范围内所有 \(p\) 的幂次 \(p^e\)，根据定义计算其结果。
一个比较朴素的实现方法是：通过不停让一个初始为 \(p\) 的变量乘以 \(p\)，若超过 \(n\) 了则直接退出。
但是这个实现方法是非常危险的，当我们筛到的质数 \(p\) 接近 \(n\) 时，我们进行一次乘法则会直接达到 \(n^2\) 的规模，很容易发生溢出。
除了开 long long/unsigned long long 以外，还有一个比较巧妙但常数较小的方法如下：
我们只要在维护 \(p^e\) 的同时再维护 \(p^{e-1}\) ，通过 \(p^{e-1}\) 是否大于 \(\lfloor{n\over p}\rfloor\) 即可；而由于我们又需要在后续判断 \(i\cdot p>n\) ，这也可以同样等价于 \(i>\lfloor{n\over p}\rfloor\)，这一次除法很好地解决了溢出问题和后续多次乘法的问题。
3.3 处理质数幂次的时间复杂度
很显然，若积性函数满足原来的可乘法递推性时，它在 \(\boldsymbol f(n)\) 转移到 \(\boldsymbol f(n\cdot p)\) 的时间复杂度是 \(O(1)\) 的，因此它从 \(\boldsymbol f(p^e)\) 转移到 \(\boldsymbol f(p^{e+1})\) 显然也是 \(O(1)\) 的。因此这一部分的时间复杂度和所有质因数的幂次数显然是等价的，而该数显然不超过 \(n\)，因此时间复杂度仍然为 \(O(n)\)。
我们不妨考虑不满足可乘法递推性时，均摊到质数幂次的时间复杂度容许为多高？不妨假设所有质数幂次的均摊复杂度为 \(T_2(n)\)，其中筛到各质数时需要先花费均摊为 \(T_1(n)\) 的时间复杂度预处理，则：
根据质数定理，\(\pi(n)\sim O({n\over \ln n})\)；为保证总复杂度为 \(O(n)\)，故 \(T_1(n)=O(\log n)\)；而 \(T_2(n)\) 则比较繁琐：
因此，同理可得，\(T_2(n)=O(\log n)\) 即可保证线性复杂度
3.4 一些细节优化
当然，我们知道，维护 \(n\) 中最小质因数 \({pe}_n\) 和维护 \(n\) 中除了最小质因数的部分 \({fr}_n\) 是等价的。那两个到底谁的效率更高呢？
我们分别考虑枚举的 \(p<{minp}_n\) 和 \(p={minp}_n\) 的情况：
当 \(p<{minp}_n\) 时，
\(
\begin{cases}
{pe}_{n\cdot p}&=p\\
{fr}_{n\cdot p}&=n\\
\end{cases}
\)，两者效率近似；
当 \(p={minp}_n\) 时，
\(
\begin{cases}
{pe}_{n\cdot p}&={pe}_n\cdot p\\
{fr}_{n\cdot p}&={fr}_n\\
\end{cases}
\)，我们发现两者均需要查表，但 \({fr}_n\) 减少了一次的乘法操作。
因此，维护 \({fr}_n\) 在最小质因数的幂次大于 \(1\) 时，均减少了一次的乘法操作。通过打表，我们可以看到 \(n\) 在不同范围时的，乘法优化次数：
| \(\log_{10} n\) | 最小质因数幂次大于 \(1\) 的数 |
|---|---|
| \(5\) | \(33009\) |
| \(6\) | \(330078\) |
| \(7\) | \(3300900\) |
| \(8\) | \(33009533\) |
可以观察到，这个方法基本上优化了 \({n\over 3}\) 次的乘法。
3.5 积性函数解耦
如果待筛的积性函数非常多，这容易造成这种欧式筛在质数的地方处理得过于冗长。这不易于维护。
然而实际上，由于我们已知了 \({fr}_n\)，我们已经具备了筛任意积性函数的能力。我们只需要先预处理质数幂次的积性函数值，再利用积性性质分别统计不同质数幂次的贡献。这样，我们就实现了欧式筛和积性函数线性求解的解耦。
额外需要注意的是，由于我们在后续筛积性函数的时候还是需要知道，各个质数有多少个幂在 \([1,n]\) 中；因此我们可以在预处理质数的时候同时保存这个值。
vector<pair<int, int> > pr;
int fr[MAXN];
for(int i=2, x=n/2, e, p1, p2; i<=n; x=n/(++i)) {
if(!pe[i]) {
for(e=0, p1=1, p2=i; p1<=x; ++e, p1=p2, p2*=i)
fr[p2]=1;
pr.emplace_back(i, e);
}
for(auto [p, e] : pr)
if(p>x)
break;
else if(i%p)
fr[p*i]=p;
else {
fr[p*i]=fr[i];
break;
}
}
int f[MAXN];
f[1]=1;
for(auto [p, e] : pr)
for(int i=1, pe=p; i<=e; ++i, pe*=p)
f[pe]=calc(p, i);//算积性函数结果
for(int i=2; i<=n; ++i)
f[i]=f[fr[i]]*f[i/fr[i]];
代码中，我们用 calc(p, i) 来计算 \(p^i\) 的积性函数结果；这里我们假设了积性函数的运算永远不会溢出，不然后续的代码需要加上取模等操作。
3.6 更广义的积性函数
事实上，我们这种的筛法不一定需要限制积性函数 \(\boldsymbol f\) 的值域为整数。实际上，我们只需要保证：
- \(\boldsymbol f(1)\) 是运算的单位元；
- \(p_1<p_2<\cdots<p_m\) 则 \(\displaystyle \boldsymbol f(\prod_{i=1}^m p_i^{e_i})=\boldsymbol f(p_1^{e_1})\cdot \boldsymbol f(p_2^{e_2}) \cdots \boldsymbol f(p_m^{e_m})\) 或 \(\displaystyle \boldsymbol f(\prod_{i=1}^m p_i^{e_i})=\boldsymbol f(p_m^{e_m})\cdot \boldsymbol f(p_{m-1}^{e_{m-1}}) \cdots \boldsymbol f(p_1^{e_1})\)；其运算符合结合律。
需要格外注意的是，若这个运算不符合交换律，
f[fr[i]]和f[i/fr[i]]的乘法必须严格按照顺序。
一个广义的积性函数是：积性函数 \(\boldsymbol f\) 是好多个朴素积性函数连接组成的向量，而运算是 Hadamard 积 \(\odot\)，即各个元素分别相乘。例如 \(\boldsymbol f(n)={(\boldsymbol \varphi(n), \boldsymbol \mu(n))}^T\)，则：
这启示我们可以将多个积性函数封装为一个结构体，统一筛出（某种意义上说，这个方法也对访存比较友好）。
当然，积性函数的值域显然也可以是矩阵、置换之类的东西；但鉴于没有例题，因此本人仅构造一个（可能不具有实际意义的）例子：
设某函数 \(\displaystyle g(\prod_{i=1}^m p_i^{e_i})=\sum_{i=1}^m e_i\cdot \boldsymbol \varphi (\prod_{j=1\wedge j\neq i}^m p_j^{e_j})\)，要求线性求解。
做法是考虑到 \(\displaystyle g(p^e \cdot q^m)=e\cdot \boldsymbol \varphi(q^m)+m\cdot \boldsymbol \varphi(p^e)\)，这和多项式乘法有点类似；于是我们定义形式多项式 \(\boldsymbol f(n,\chi)=\boldsymbol \varphi(n)+g(n)\cdot \chi\)，运算是模 \(\chi^2\) 意义下的多项式乘法，则：
因此，我们只要定义广义积性函数 \(\boldsymbol f(n)=\left<\boldsymbol \varphi(n), g(n)\right>\)，运算为模 \(\chi^2\) 意义下的多项式乘法，即可通过欧式筛法求解。
四、例题与实验结果
接下来，本文会根据题目所求积性函数从简单到复杂，分别进行一些实现。通过例题的方式讲解如何在 \(O(\log n)\) 时间内对不同的积性函数，在质数幂次处的结果进行维护转移。
此外，本段将对本文实现方法和除本文实现方法外最快的实现方法进行对比，以验证该方法对效率的影响。
4.1 朴素的积性函数
4.1.1 欧拉函数
简单观察即可知道要求的是：
\(\displaystyle \sum_{i=1}^n\sum_{j=1}^m [2\gcd(i,j)-1]=2\sum_{d=1}^n\boldsymbol \varphi(d)(n/d)(m/d)-nm\)
其中 \((n/d)\) 表示 \(\lfloor{n\over d}\rfloor\)
按照乘法递推的方法即可筛欧拉函数。一个实现方法是：在递推函数中，维护静态变量 lst，表示上一次返回的结果；通过传入参数 p,e 可以快速判断结果的变化，例如：
static uint lst;//unsigned
if(e==1) return lst=p-1;
else return lst=lst*p;
该实现方法的评测 记录 为 31ms；而除本文外，目前最快的 记录 为 24ms。本文耗时仅为其约 1.375 倍。
4.1.2 除数和函数
考虑积性函数在质数幂次位置的结果：
因此，可以同上维护静态变量 lst 表示上次的输出结果，也是容易转移的。
static uint lst;
if(e==1) return lst=fpow(p, k)-1;
else return lst=lst*fpow(p, k);
但实际上，由于每次反复计算 fpow(p, k) 会造成额外的 \(O(\log k)\) 开销，这里可以再额外引入静态变量 pk 进行维护。
static uint pk = -1, lst = -1;
if (e == 1)
pk = fpow(p, k), lst = 1 + pk;
else
lst = ((ull)lst * pk + 1) % P;
return lst;
该实现方法的评测 记录 为 1819ms；而除本文外，目前最快的 记录 为 1646ms。本文耗时约为其 1.105 倍。
4.2 较为复杂的积性函数
4.2.1 \(\boldsymbol \varphi*\boldsymbol \mu\) 的欧拉筛
考虑积性函数在质数幂次位置的结果：
故 \(e=1\) 时 \((\boldsymbol\varphi*\boldsymbol \mu)(p^e)=p-1-1=p-2\)；
\(e=2\) 时 \((\boldsymbol\varphi*\boldsymbol \mu)(p^e)=p^2-p-p+1=p^2-2p+1\)
\(e>2\) 时 \((\boldsymbol\varphi*\boldsymbol \mu)(p^e)=p\boldsymbol \varphi(p^{e-1})-p\boldsymbol \varphi(p^{e-2})=p\cdot (\boldsymbol\varphi*\boldsymbol \mu)(p^{e-1})\)。
通过上面的分类讨论，很容易想到，同样是静态变量 lst 维护上次的结果，然后根据 \(e\) 进行分类讨论即可。
static uint lst=-1;
if(e==1) return lst=p-2;
else if(e==2) return lst=(p-1)*(p-1);
else return lst*=p;
该实现方法的评测 记录 为 876ms；而同等方法下，除本文外，目前最快的 记录 为 1295ms。本文耗时约为其 0.676 倍。
4.2.2 \(\boldsymbol \varphi*\boldsymbol \varphi\) 的欧拉筛
考虑积性函数在质数幂次位置的结果：
故 \(e=1\) 时 \((\boldsymbol\varphi*\boldsymbol \varphi)(p^e)=p-1+p-1=2(p-1)\)；
\(e=2\) 时 \((\boldsymbol\varphi*\boldsymbol \varphi)(p^e)=p(p-1)+p(p-1)+(p-1)(p-1)=(3p-1)(p-1)\)
\(e>2\) 时 \((\boldsymbol\varphi*\boldsymbol \varphi)(p^e)=p\cdot (\boldsymbol\varphi*\boldsymbol \varphi)(p^{e-1})+(p-1)\cdot \boldsymbol\varphi(p^{e-1})\)。
欧拉函数的筛法我们之前已经进行了讨论，而 \(\boldsymbol \varphi*\boldsymbol \varphi\) 的筛法通过分类讨论，可以想到法一的做法。此外，其转移过程中还额外需要知道 \(\boldsymbol \varphi\) 的结果，这恰好是我们同时求的结果（当然，若不需要求，可以增加这个积性函数，或作为静态变量）。
static puu lst;//pair<ull, ull> where ull=unsigned long long
if(e==1) return lst=puu((p-1)*2, p-1);
else if(e==2) return lst=puu((p-1)*(3*p-1), (p-1)*p);
ull phi=lst.se*p;
ull phi2=p*lst.fi+(p-1)*lst.se;
return lst=puu(phi2, phi);
该实现方法的评测 记录 为 533ms；而除本文外，目前最快的 记录 为 654ms。本文耗时约为其 0.815 倍。
此处因为找不到同等方法，而本题做法较多。故仅与全场最快方法进行了对比。
4.2.3 \(\boldsymbol {id^k}*\boldsymbol \mu\) 的欧拉筛
考虑积性函数在质数幂次位置的结果：
故 \(e=1\) 时 \((\boldsymbol{{id}^k}*\boldsymbol \mu)(p^e)=p^k-1\)；
\(e>1\) 时 \((\boldsymbol{{id}^k}*\boldsymbol \mu)(p^e)=p^k\cdot (\boldsymbol{{id}^k}*\boldsymbol \mu)(p^{e-1})\)。
此处维护 lst 的考量和欧拉函数类似，而维护 pk 可以避免反复求解 fpow(p, k) 和 \(\boldsymbol {\sigma_k}\) 类似。
static uint pk, lst;
if(e==1) {
pk=fpow(p, k);
return lst=(pk+P-1)%P;
}
else
return lst=(ull)lst*pk%P;
该实现方法的评测 记录 为 3396ms；而除本文外，目前最快的 记录 为 1870ms。本文耗时约为其 1.816 倍。
4.3 复杂积性函数的组合
4.3.1 两个复杂积性函数的组合
LOJ #6680. 「hyOI2019」henry_y 的函数
考虑到： \(\lfloor{n\over i}\rfloor>\lfloor{n-1\over i}\rfloor\) 时，\(\exists d\in \mathrm{Z}\wedge {n-1\over i}<d\leq {n\over i}\)；
此时有 \(n-1<id\leq n\)，故当且仅当 \(n=id\)；
也即 \(i\nmid n\) 时 \({\lfloor{n\over i}\rfloor}-{\lfloor{n-1\over i}\rfloor}=0\)。
考虑积性函数在质数幂次位置的结果：
此题前半部分的按位异或可以 \(O(1)\) 求出，后半部分的递推，第一个和欧拉函数类似，第二个则无需额外递推。
static puu lst;//pair<uint, uint>
if (e == 1)
return lst = puu(p + (p ^ 1), 1 + (p ^ 1));
else
return lst = puu(
((ull)lst.fi * p + (p ^ e)) % P,
(lst.se + (p ^ e)) % P
);
该实现方法的评测 记录 为 1459ms；而除本文外，目前最快的 记录 为 1670ms。本文耗时约为其 0.874 倍。
4.3.2 多个复杂积性函数的组合
“山大地纬杯”第十二届山东省ICPC大学生程序设计竞赛(正式赛） I-gcds
此题不确定是否存在更优解法。目前按照本人的解法，过程较为复杂，可以参考本人的 博客园链接 或 知乎链接。
最终的问题落到了求解 \(\boldsymbol {f_{t, i}}=\boldsymbol {id^{t+1}} * (\boldsymbol \mu\cdot \boldsymbol {id^t})*\boldsymbol {id^i}\) 上，其中 \(\left<t,i\right>\in \{\left<1, 1\right>, \left<1, 2\right>, \left<2, 1\right>, \left<2, 2\right>, \left<2, 3\right>, \left<3, 2\right>, \left<3, 3\right>, \left<3, 4\right>\}\)。
考虑积性函数在质数幂次位置的结果：
此题由于递推式中 \(p^{ie}, \boldsymbol {\sigma_{t+1-i}}(p^e), p^t\) 均不是代求项，我们可以额外拓展 \(4+3+3=10\) 项积性函数，但本题本身已经有 \(8\) 个积性函数代求，未免太过夸张。
实际上，我们只需要额外维护静态变量来表示这些递推项即可：我们额外维护静态变量 sigma[0-4] 表示 \(\boldsymbol {\sigma_0}(p^e)\) 到 \(\boldsymbol {\sigma_4}(p^e)\)，维护 pe[0-4] 表示 \(p^{e\cdot 0}\) 到 \(p^{e\cdot 4}\)，维护 pk[0-4] 表示 \(p^0\) 到 \(p^4\)。
由于积性函数 \(\boldsymbol f_{t,i}(p^e)\) 的递推式中有用到这些函数在 \(p^{e-1}\) 处的结果，但是它们的使用和 \(p^e\) 的结果不交叉；
一个方法是我们每个静态变量再翻倍，分别表示 \(p^e\) 的结果和 \(p^{e-1}\) 的结果；或者是我们先统一统计了 \(p^{e-1}\) 的结果，再更新上述 15 个静态变量，最后再统计 \(p^e\) 的贡献。
const uint farr_t[] = {1, 1, 2, 2, 2, 3, 3, 3};
const uint farr_i[] = {1, 2, 1, 2, 3, 2, 3, 4};
const uint fnum = sizeof(farr_t) / sizeof(farr_t[0]);
typedef array<uint, fnum> farr;
static uint sigma[5], pk[5], pe[5];
if(e==1) {
pe[0]=pk[0]=1; sigma[0]=2;
for(int i=1; i<=4; ++i) {
pe[i]=pk[i]=(ull)pe[i-1]*p%P;
sigma[i]=(pe[i]+1)%P;
}
farr f;
for(int j=0; j<fnum; ++j) {
uint t=farr_t[j], i=farr_i[j];
f[j]=((ull)pe[i]*sigma[t+1-i]-pk[t]+P)%P;
}
return f;
}
else {
farr f;
for(int j=0; j<fnum; ++j) {//统计 p^{e-1} 贡献
uint t=farr_t[j], i=farr_i[j];
f[j]=(P-(ull)pe[i]*sigma[t+1-i]%P*pk[t]%P);
}
for(int i=0; i<=4; ++i) {//更新静态变量
pe[i]=(ull)pe[i]*pk[i]%P;
sigma[i]=((ull)sigma[i]*pk[i]+1)%P;
}
for(int j=0; j<fnum; ++j) {//统计 p^e 贡献
uint t=farr_t[j], i=farr_i[j];
f[j]=(f[j]+(ull)pe[i]*sigma[t+1-i])%P;
}
return f;
}
该实现方法的评测 记录 为 2171ms；而除本文外，目前最快的 记录 为 2823ms。本文耗时约为其 0.769 倍。
4.4 实验小结
| 序号 | 题目 | 所求项 | 本文记录(ms) | 除本文方法的最快记录(ms) | 比率 |
|---|---|---|---|---|---|
| 1 | Luogu P1447 [NOI2010] 能量采集 | \(\boldsymbol \varphi\) | 31 | 24 | 1.375 |
| 2 | LOJ #124. 除数函数求和 1 | \(\boldsymbol {\sigma_k}\) | 1819 | 1646 | 1.105 |
| 3 | BZOJ #4804. 欧拉心算 法一 | \(\boldsymbol \varphi*\boldsymbol \mu\) | 876 | 1295 | 0.676 |
| 4 | BZOJ #4804. 欧拉心算 法二 | \(\boldsymbol \varphi, \boldsymbol \varphi*\boldsymbol \varphi\) | 533 | 654 | 0.815 |
| 5 | BZOJ #4407. 于神之怒加强版 | \(\boldsymbol {id^k}*\boldsymbol \mu\) | 3396 | 1870 | 1.816 |
| 6 | LOJ #6680. 「hyOI2019」henry_y 的函数 | \(\boldsymbol g*\boldsymbol I, \boldsymbol g*\boldsymbol {id}\) | 1459 | 1670 | 0.874 |
| 7 | “山大地纬杯”第十二届山东省ICPC大学生程序设计竞赛(正式赛） I-gcds | \(\boldsymbol {id^{t+1}} * (\boldsymbol \mu\cdot \boldsymbol {id^t})*\boldsymbol {id^i}\) | 2171 | 2823 | 0.769 |
本部分先是通过了例题，讲解了这种新式的欧式筛写法该如何维护；再通过实验，表明了该算法的效率和这些题最好实现的差异。
就实验数据而言，新式的欧式筛实现方法在大部分情况下，并不会过度地劣于最优算法。并且，就结果而言，代筛的积性函数越复杂，新式实现方法的效率比就越高。这也充分体现了新式写法的擅长点——即复杂积性函数，或多个积性函数的组合。
五、结论
本文对欧式筛求积性函数进行了讨论。通过对传统欧式筛的实现方法进行分析，提出了传统欧式筛解决起来较为复杂的积性函数情形。
随后，本文根据此情形提出了解决方法，通过维护 pe 数组或 fr 数组，与预处理质数幂次的方法，解决部分积性函数难以线性筛出的问题。
通过理论分析，本文证明了只需要保证均摊 \(O(\log n)\) 的时间处理质数幂次，即可保证该算法的线性时间复杂度。此外，本文提出了通过维护 fr 数组而非 pe 数组的方法可以提高运行效率。
接下来，本文提出了积性函数求解解耦的概念，随后简要探讨了广义的积性函数可以如何用新式实现方法求解。
最后，本文讨论了新式实现方法在具体题目里的求解实现细节；并通过实验表明，新式实现方法在大部分情形下并不劣于传统实现方法，尤其是在筛复杂积性函数或多个积性函数组合的情形。
附录
实验用的新式欧式筛板子：
typedef unsigned uint;
#define fi first
#define se second
#define eb emplace_back
#define vt vector<T>
#define tmpt template<class T>
struct euler_sieve{
uint n;
vector<pair<uint, uint> > pr;
uint fr[MAXN];
euler_sieve(uint n_):n(n_) {
for(uint i=2, x=n/2, p1, p2, e; i<=n; x=n/(++i)) {
if(!fr[i]) {
for(p1=1, p2=i, e=0; p1<=x; p1=p2, p2*=i, ++e)
fr[p2]=1;
pr.eb(i, e);
}
for(int t=0, T=pr.size(); t<T; ++t) {
uint p=pr[t].fi;
if(p>x)
break;
else if(i%p)
fr[p*i]=i;
else {
fr[p*i]=fr[i];
break;
}
}
}
}
tmpt vt f(function<T(uint, uint)> fp, function<T(T, T)> mul=[](T a, T b) { return a*b; }, T I=1) {
vt v(n+1); v[1]=I;
for(uint i=0, I=pr.size(); i<I; ++i)
for(uint p=pr[i].fi, e=pr[i].se, pe=p, ex=1; ex<=e; ++ex, pe*=p)
v[pe]=fp(p, ex);
for(uint i=2; i<=n; ++i) v[i]=mul(v[i/fr[i]], v[fr[i]]);
return v;
}
}*es;
