algorithm - 快速迭代 GCD

Question

我有 GCD(n, i) 其中 i=1 在循环中增加 1 到 n。是否有任何算法可以比天真的增加和使用欧几里得算法计算 GCD 更快地计算所有 GCD？

PS我注意到如果n是素数，我可以假设从1到n-1的数字会给出1，因为素数对它们来说是互质的。除了素数之外的其他数字有什么想法吗？

Answer 1

概括

gcd 的可能答案由 n 的因数组成。

您可以如下有效地计算这些。

算法

首先将n分解为素数的乘积，即n=p1^n1*p2^n2*..*pk^nk。

然后，您可以遍历 n 的所有因子，并为 n 的每个因子将 GCD 数组在该位置的内容设置为因子。

如果您确保这些因素以合理的顺序（例如排序）完成，您应该会发现多次写入的数组条目最终将被写入最高值（这将是 gcd）。

代码

下面是一些 Python 代码来为数字 1400=2^3*5^2*7 执行此操作：

prime_factors=[2,5,7]
prime_counts=[3,2,1]
N=1
for prime,count in zip(prime_factors,prime_counts):
    N *= prime**count

GCD = [0]*(N+1)
GCD[0] = N
def go(i,n):
    """Try all counts for prime[i]"""
    if i==len(prime_factors):
        for x in xrange(n,N+1,n):
            GCD[x]=n
        return
    n2=n
    for c in xrange(prime_counts[i]+1):
        go(i+1,n2)
        n2*=prime_factors[i]
go(0,1)    
print N,GCD

Answer 2

C++ 实现，在 O(n * log log n) 中工作（假设整数的大小为 O(1)）：

#include <cstdio>
#include <cstring>
using namespace std;

void find_gcd(int n, int *gcd) {
  // divisor[x] - any prime divisor of x
  //              or 0 if x == 1 or x is prime
  int *divisor = new int[n + 1];
  memset(divisor, 0, (n + 1) * sizeof(int));

  // This is almost copypaste of sieve of Eratosthenes, but instead of
  // just marking number as 'non-prime' we remeber its divisor.
  // O(n * log log n)
  for (int x = 2; x * x <= n; ++x) {
    if (divisor[x] == 0) {
      for (int y = x * x; y <= n; y += x) {
        divisor[y] = x;
      }
    }
  }

  for (int x = 1; x <= n; ++x) {
    if (n % x == 0) gcd[x] = x;
    else if (divisor[x] == 0) gcd[x] = 1; // x is prime, and does not divide n (previous line)
    else {
      int a = x / divisor[x], p = divisor[x]; // x == a * p
      // gcd(a * p, n) = gcd(a, n) * gcd(p, n / gcd(a, n))
      // gcd(p, n / gcd(a, n)) == 1 or p
      gcd[x] = gcd[a];
      if ((n / gcd[a]) % p == 0) gcd[x] *= p;
    }
  }
}

int main() {
  int n;
  scanf("%d", &n);
  int *gcd = new int[n + 1];
  find_gcd(n, gcd);
  for (int x = 1; x <= n; ++x) {
    printf("%d:\t%d\n", x, gcd[x]);
  }
  return 0;
}

Answer 3

二进制 GCD 算法：
https ://en.wikipedia.org/wiki/Binary_GCD_algorithm

比欧几里得算法快：
https ://en.wikipedia.org/wiki/Euclidean_algorithm

我在 C 中为类型“__uint128_t”（在 Intel i7 Ubuntu 上使用 gcc）实现了“gcd()”，基于迭代 Rust 版本：
https ://en.wikipedia.org/wiki/Binary_GCD_algorithm#Iterative_version_in_Rust

使用“__builtin_ctzll()”可以有效地确定尾随 0 的数量。我针对 gmplib "mpz_gcd()" 对两个最大的 128 位斐波那契数的 100 万个循环（它们导致最大迭代次数）进行了基准测试，并看到了 10% 的减速。利用 u/v 值只会减小的事实，我在“<=UINT64_max”时切换到 64 位特殊情况“_gcd()”，现在看到 gmplib 的加速比为 1.31，详情请参阅：
https ://www.raspberrypi.org /forums/viewtopic.php?f=33&t=311893&p=1873552#p1873552

inline int ctz(__uint128_t u)
{
  unsigned long long h = u;
  return (h!=0) ?      __builtin_ctzll( h )
                : 64 + __builtin_ctzll( u>>64 );
}

unsigned long long _gcd(unsigned long long u, unsigned long long v)
{
  for(;;) {
    if (u > v) { unsigned long long a=u; u=v; v=a; }

    v -= u;

    if (v == 0)  return u;

    v >>= __builtin_ctzll(v);
  }
}

__uint128_t gcd(__uint128_t u, __uint128_t v)
{
       if (u == 0) { return v; }
  else if (v == 0) { return u; }

  int i = ctz(u);  u >>= i;
  int j = ctz(v);  v >>= j;
  int k = (i < j) ? i : j;

  for(;;) {
    if (u > v) { __uint128_t a=u; u=v; v=a; }

    if (v <= UINT64_MAX)  return _gcd(u, v) << k;

    v -= u;

    if (v == 0)  return u << k;

    v >>= ctz(v);
  }
}