线性时间Euler的Totient函数计算

Question

在浏览了一下之后,我发现这个代码用于使用Sieve of Eratostenes在线性时间内计算Euler的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

编码既不在这里也不在那里,但算法本身就很棒.但它与Eratosthenes的筛子几乎没有关系.这种方法有点让人想起Sundaram的Sieve,因为它系统地产生多个质数以标记复合材料; 它比Sundaram的更好,因为每个复合材料只被划掉一次(没有透支).同样地,计算每个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)=披(M*P)=披(M)*披(P)=披(M)*(P - 1)

(3)k = m*p,m = m'*p ^ n和m'和p互质和p prime:phi(k)= phi(m)*p

这两个相关的数学事实在Euler的产品公式中列在维基百科关于欧拉函数的文章中.情况1在外环中处理,情况2是then内环中条件的分支,情况3是else也终止内环的分支.情况2和3为较小数字的预先存在的phi值构建复合材料的phi值,所有这些都最终来自素数的phi值(在外环中设置).

算法的亮度在于它安排要完成的工作的方式,使得每个值只计算一次并且在需要时已经计算.它为复合实现的递归基于通过分解其最小的素因子来有效地分割每个复合:m = m'*p.这有利于根据上述情况2和3计算复合材料的phi,并且它导致用于生产没有重复的复合材料的简单规则.除此之外,外部循环向内部循环呈现所有可能的m'(尽管对于i> N/2,内部循环从不接受任何而外部循环仅用于收集剩余的质数).然后内环产生所有复合物,它们是m'的乘积和不超过其自身主要因子最小值的素因子.

该代码已根据从phi函数的OEIS页面检索的前100,000个phi值的列表进行验证.它是明确地写出来展示算法的工作原理; 在用于程序之前,必须对其进行审查,调整和强化(特别是防止溢出).

附录 -如果任何人无意中跳过DanaJ的评论:所述is_composite[]阵列是多余的,因为在给定的非compositeness(素性)i可以通过测试来确定phi[i]在外环为零.原因是复合材料的phi值向上传播 - 它们是在早期迭代期间计算的,当i它们是其中一个因素时.另一种推理方法是,is_composite[m]只有true在计算和存储相应的phi值时才设置,并且这些值永远不能为零.因此,外环中的测试变为if (phi[i] == 0).而实施者(而不是"纯粹的鉴赏家")可能想更仔细地看待DanaJ的评论......