fastest method to determine if two numbers are coprime

Question

I am working on a mathematical problem that involves coprime integers. I wrote a computer program that allows me to search for the numbers I am looking for. However I am looking at a large set of integers and I have to compare many pairs of numbers and determine if they are coprime. My program has to do this so many times that any reduction in the calculation time for each pair of numbers would significantly reduce the overall run time of the program.

I am currently using the Euclidean algorithm.

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

Even though the Euclidean algorithm is a very efficient method, I don't know if it's the fastest. I don't need the gcd of the two numbers I just need to determine if they are coprime or not.

My pseudocode:

while (B ≠ 0)

{T = B

B = A mod T

A = T}

where A and B are the pair of integers in question and T is a dummy variable used to hold a value to be used in a later calculation. For those who haven't written a computer program before, while means that the process in the brackets will repeat until the condition in the parentheses is false. At the end of the loop the variable A becomes the gcd of the original two integers. if A=1 the two numbers are coprime if A>1 then the numbers are not coprime.

Even though the program above is simple, it is an iterative process and I'm looking for a method that only needs one or two steps.

Thanks in advance!

Answer 1

Since you care about the speed of IsCoprime(n, m), I'm assuming that you're calling that function many times. So let's assume that we're minimizing the time the following code takes to run (K is some constant):

for (int n = 1; n <= K; ++n)
  for (int m = 1; m <= K; ++m)
    IsCoprime(n, m);

You can precompute prime factors all ints up to K, which takes $O(K \log K)$, much smaller than the number of invocations of IsCoprime above, which is $O(K^2)$.

Then for each int up to K, compute its 32-bit "prime factor signature" sig[n].

• Bit 0 is set if n is divisible by 2
• Bit 1 is set if n is divisible by 3
• Bit 2 is set if n is divisible by 5
• ...
• Bit 30 is set if n is divisible by 113
• Bit 31 is set if n is divisible by some prime > 113

Now, the first step you do in IsCoprime(n, m) is to compute binary and of the signatures of its arguments: X = sig[n] & sig[m]. If X is 0, then n and m are coprime. Otherwise, if some bit other than bit 31 of X is set, then n and m are not coprime. If neither of these conditions hold, you need to run some GCD algorithm to get the answer. However, you'll need to call GCD in a small number of cases.

One improvement you can make is having each of bits 16..31 represent multiple prime factors. For example, bit 16 would be set if n is divisible by any of the primes between 47 and 97, or something like that.

And of course you can always use more bits, like 64.

When I tried these techniques for K roughly 30,000, I got an over 10x speed improvement over (binary) GCD.

EDIT (2022-07-11)

Here's a concrete implementation of this idea: https://gist.github.com/zielaj/cfdeca92fc0b8371164e4ca521d83587

For K = 0x4000 (16K), it runs in 0.337s on my laptop, vs 5.798s for binary GCD.

This code uses 64-bit signatures. The lower 32 bits are dedicated to the first 32 primes. Subsequent primes get 2 "random" upper bits each, an approach inspired by Bloom filters: https://en.wikipedia.org/wiki/Bloom_filter

Here's the code snippet for the actual coprimality testing:

  inline bool AreCoprime(uxint n, uxint m) {
    if (n == 0 || m == 0) return n + m == 1;
    uxint sig = signatures[n] & signatures[m];

    // If the signatures are disjoint, then there are no prime factors
    // in common, so the numbers are coprime.
    if (sig == 0) return true;

    // If the signatures overlap somewhere in the first 32 bits, then
    // there's a common prime factor because each of these bits is
    // dedicated to a single prime.
    if ((sig & uniquebits) != 0) return false;

    // If there's at most one bit overlap (in the second half of the
    // signature), then there are no prime factors in common because
    // each of the primes in the second half gets 2 bits.
    if ((sig & (sig - 1)) == 0) return true;

    // If we're still undecided, fall back to the GCD computation.
    return bgcd(n, m) == 1;
  }

Note that if you just want to generate pairs of small coprime integers, you can also use Stern-Brocot tree, which is O(1) per pair. However, when I tried that, using an admittedly unoptimized straightforward recursive implementation, it was still slower than the code above. https://en.wikipedia.org/wiki/Stern%E2%80%93Brocot_tree