快速迭代GCD

Question

我有GCD( n，i)，其中i=1在循环中递增1到n。有没有什么算法可以计算所有的GCD比朴素递增更快，并使用欧几里德算法计算GCD？

我注意到，如果n是素数，我可以假设从1到n-1的数等于1，因为素数对它们来说是共同质数。除了质数之外，还有其他数字的想法吗？

Answer 1

C++实现，工作在O(n * log )(假设整数的大小为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 2

摘要

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

您可以按如下方式高效地计算这些值。

算法

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

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

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

代码

以下是对数字1400=2^3*5^2*7执行此操作的一些Python代码：

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 3

二进制GCD算法：

https://en.wikipedia.org/wiki/Binary_GCD_algorithm

比欧几里德算法更快：

https://en.wikipedia.org/wiki/Euclidean_algorithm

我用C语言实现了"__uint128_t“类型的"gcd()”(在英特尔i7 Ubuntu上和gcc一起)，基于iterative版本：

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);
  }
}