c++ - 线性时间欧拉的 Totient 函数计算

Question

经过一番浏览，我找到了这段代码，用于使用 Eratostenes 筛法在线性时间内计算欧拉的 phi 值。但是未能理解这段代码中使用的主要逻辑，特别是内部 for 循环中所做的事情和使用的想法这个循环用于计算 phi 值。如果有人能帮助我理解这段代码，那将会很有帮助。

#define MAXN 3000000
int phi[MAXN + 1], prime[MAXN/10], sz;
bitset <MAXN + 1> mark;

for (int i = 2; i <= MAXN; i++ ){
   if(!mark[i]){
      phi[i] = i-1;
      prime[sz++]= i;
   }
    for (int j=0; j<sz && prime[j]*i <= MAXN; j++ ){
            mark[prime[j]*i]=1;
            if(i%prime[j]==0){
                  phi[i*prime[j]] = phi[i]*prime[j];
                  break;
            }
            else phi[i*prime[j]] = phi[i]*(prime[j]-1 );
   }
}

Answer 1

代码既不在这里也不在那里，但算法本身出色。不过，它与埃拉托色尼筛法关系不大。该方法有点让人联想到 Sieve of Sundaram。因为它系统地产生素数的倍数以标记复合物；它比 Sundaram 的更好，因为每个复合 Material 都恰好划掉一次(no overdraw)。同样，每个 phi 值都只计算并分配一次。

如果先稍微改动一下代码会更容易理解:

enum {  N = 3000000  };
vector<unsigned> primes;
vector<unsigned> phi(N + 1);        // indexed directly with numbers, hence the +1
vector<bool> is_composite(N + 1);   // indexed directly with numbers, hence the +1

phi[1] = 1;   // phi[0] is 0 already; phi[1] needs to be set explicitly

for (unsigned i = 2; i <= N; ++i)
{
   if (not is_composite[i])  // it's a prime
   {
      phi[i] = i - 1;        // a prime is coprime to all numbers before it
      primes.push_back(i);
   }

   for (auto current_prime: primes)
   {
      unsigned ith_multiple = i * current_prime;

      if (ith_multiple > N)
         break;

      is_composite[ith_multiple] = true;

      if (i % current_prime)  // i and current_prime are coprime -> phi is multiplicative in this case
      {
         phi[ith_multiple] = phi[i] * (current_prime - 1);  // (current_prime - 1) == phi[current_prime]
      }
      else
      {
         phi[ith_multiple] = phi[i] * current_prime;  // based on phi(p^(k+1)) / phi(p^k) == p

         break;
      }
   }
}

phi(k) 值的计算必须考虑三种不同的情况:

(1) k 是质数:phi(k) = k - 1，这是平凡的

(2) k = m * p 其中 m 和 p 互质且 p 质数:phi(k) = phi(m * p) = phi(m) * phi(p) = phi(m) * (p - 1)

(3) k = m * p with m = m' * p^n and m' and p coprime and p prime: phi(k) = phi(m) * p

两个相关的数学事实列在 Euler's product formula 下在关于 Euler's totient function 的维基百科文章中.案例1在外循环中处理，案例2是then内循环条件的分支，案例 3 是 else也终止内部循环的分支。案例 2 和案例 3 从先前存在的较小数字的 phi 值构建复合 Material 的 phi 值，所有这些最终都来自素数的 phi 值(在外循环中设置)。

该算法的出色之处在于它安排要完成的工作的方式，这样每个值都只计算一次，并且在需要时已经计算过了。它为复合 Material 实现的递归是基于通过分解出其最小素因数来有效地拆分每个复合 Material :m = m' * p。这有助于根据上述情况 2 和 3 计算复合 Material 的 phi，并且它导致生成没有重复的复合 Material 的简单规则。除其他外，外循环将所有可能的 m' 呈现给内循环(尽管对于 i > N/2 内循环从不接受任何并且外循环旋转仅用于收集剩余的素数)。然后，内循环生成所有复合 Material ，这些复合 Material 是 m' 和不超过其自身素数中最小素数的素数的乘积。

代码已根据 list of the first 100,000 phi values 进行了验证从 OEIS page for the phi function 检索.它是专门为展示算法的工作原理而编写的；在程序中使用之前，必须对其进行审查、调整和强化(特别是防止溢出)。

附录 - 以防有人无意中跳过 DanaJ 的评论:is_composite[]数组是多余的，因为给定的 i 的非复合性(素性)可以通过测试确定phi[i]在外循环中为零。原因是复合 Material 的 phi 值向上传播 - 它们是在 i 时的早期迭代中计算的。是他们的因素之一。另一种推理方式是 is_composite[m]仅设置为 true当计算和存储相应的 phi 值时，这些值永远不会为零。因此，外循环中的测试变为 if (phi[i] == 0) .实现者(而不是“纯粹的”鉴赏家)可能想更仔细地研究 DanaJ 的评论......