Thread: greatest common divisor is a prime number (pairs of numbers problem in C)

Not only you can boost the performance of primalty test using the square root criteria, but you can boost the GCD algorithm as well using binary tricks.
Here's a GCC implementation.

Code:

#include <math.h>

int isPrime( unsigned int n )
{
  unsigned int sqrt_, i;
  
  // 1, 2 and 3 are always primes (special cases).
  // if you want to be mathematically correct, ZERO isn't a prime. Add:
  //  if ( ! n ) return 0;
  if ( n <= 3 )
    return 1;
  
  // if n is even, isn't prime for sure!
  if ( !( n & 1 ) )
    return 0;
  
  // Calculates the upper limit and must be odd.
  // Yep, I'm using sqrt() because double has 53 bits long precision
  // and unsigned int has 32. Don't use sqrtf()!
  sqrt_ = sqrt( n );
  if ( !( sqrt_ & 1 ) )
    sqrt_--;
  
  // We need to check only odd divisors starting from 3.
  for ( i = 3; i <= sqrt_; i += 2 )
  {
    // if n is divisible by any of them, not a prime!
    if ( !( n % i ) )
      return 0;
  }
  
  // if we get here, it is a prime.
  return 1;
}

// Some people prefer to use the XOR method, but my measures showed this
// tends to be faster.
#define swapu(a,b) { unsigned int t = (a); (a) = (b); (b) = t; }

unsigned int gcd( unsigned int u, unsigned int v )
{
  int shift;

  if ( u == 0 )
    return v;

  if ( v == 0 )
    return u;

  // Count the minimum bits we need to shift back,
  // since we must obtain the GREATEST common divisor.
  shift = __builtin_ctz( u | v );

  // Drop the zeroed trailing bits to improve performance.
  u >>= __builtin_ctz( u );

  do
  {
    // Drop the zeroed trailing bits to improve performance.
    v >>= __builtin_ctz( v );

    // v must be greater than u...
    if ( v < u )
      swapu( u, v );

    v -= u;
  } while ( v );

  // At this point u holds the GCD, but we must
  // right shift back...
  return u << shift;
}

Don't forget to use optimization options '-O2 -ffast-math -march=native'.

Originally Posted by flp1969

Not only you can boost the performance of primalty test using the square root criteria, but you can boost the GCD algorithm as well using binary tricks.
Here's a GCC implementation.

Code:

int isPrime( unsigned int n )
{
  unsigned int sqrt_, i;
  
  // 1, 2 and 3 are always primes (special cases).
  // if you want to be mathematically correct, ZERO isn't a prime. Add:
  //  if ( ! n ) return 0;
  if ( n <= 3 )
    return 1;
  
  // if n is even, isn't prime for sure!
  if ( !( n & 1 ) )
    return 0;
  
  // Calculates the upper limit and must be odd.
  // Yep, I'm using sqrt() because double has 53 bits long precision
  // and unsigned int has 32. Don't use sqrtf()!
  sqrt_ = sqrt( n );
  if ( !( sqrt_ & 1 ) )
    sqrt_--;
  
  // We need to check only odd divisors starting from 3.
  for ( i = 3; i <= sqrt_; i += 2 )
  {
    // if n is divisible by any of them, not a prime!
    if ( !( n % i ) )
      return 0;
  }
  
  // if we get here, it is a prime.
  return 1;
}

Don't forget to use optimization options '-O2 -ffast-math -march=native'.

Minor point: one is not prime!

Another optimization is to use multiplication in place of the sqrt() call. It just tends to run faster. You can also take advantage of the fact that any prime p > 3 is always in the form 6x +- 1. So if n % 6 doesn't equal 1 or 5 then it can't possibly be prime.

Originally Posted by flp1969

Don't forget to use optimization options '-O2 -ffast-math -march=native'.

I'd be wary of using -ffast-math as it enables "unsafe" optimisations that can (and does) break IEEE compliance. It also disables the proper setting of errno (edit: actually it stops setting errno completely I think). It rearranges stuff (expressions) that programmers might deliberately put in a specific order for the most accurate result. It does a heap of other potentially unsafe things as well and it's probably not needed in 99.9% of cases. If you're programming a game or a demo then maybe (only maybe) it'll make a noticeable difference although I've never seen it. In short, and in general, the disadvantages far outweigh the potential advantages and I would not use it.

Originally Posted by Sir Galahad

Minor point: one is not prime!

Another optimization is to use multiplication in place of the sqrt() call. It just tends to run faster. You can also take advantage of the fact that any prime p > 3 is always in the form 6x +- 1. So if n % 6 doesn't equal 1 or 5 then it can't possibly be prime.

Good point... 1 isn't prime as well...
How to use multiplication instead of square root?
I didn't get it the '6x ± 1' tip. If this is true, than criptographic algorithms like RSA would be very easy to break... Do you have any proof?

119 or 121, for example, aren't primes (6*20 ± 1): 17 * 7 = 119 and 11 * 11 = 121.

Even testing for the remainder for 1 and 5 you get a lot of non primes in UINT_MAX/6 range:

Code:

/* prime.c */
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <errno.h>
#include <limits.h>
#include <math.h>

static int isPrime( unsigned int );

int main( void )
{
  unsigned int n, m, o;

  for (n = 1; n < UINT_MAX/6; n++)
  {
    int p1, p2;
    unsigned int r;

    m = 6*n-1; 
    o = 6*n+1; 
    r = n % 6;

    p1 = isPrime(m);
    p2 = isPrime(o);

    // Shows all values in form 6x±1 which aren't primes
    if ( ! p1 && ! p2 && ( r == 1 || r == 5 ) )
      printf( "%u or %u, x=%u (remainder=%u)\n", m, o, n, r );
  }
}

int isPrime( unsigned int n )
{
  unsigned int sqrt_, i;

  // To be mathematically correct, zero AND one aren't primes.
  if ( n < 2 )
    return 0;

  // 2 and 3 are always primes (special cases).
  if ( n <= 3 )
    return 1;

  // if n is even, isn't prime for sure!
  if ( !( n & 1 ) )
    return 0;

  // Calculates the upper limit and must be odd.
  sqrt_ = sqrt( n );
  if ( !( sqrt_ & 1 ) )
    sqrt_--;

  // We need to check only odd divisors starting from 3.
  for ( i = 3; i <= sqrt_; i += 2 )
  {
    // if n is divisible by any of them, not a prime!
    if ( !( n % i ) )
      return 0;
  }

  // if we get here, it is a prime.
  return 1;
}

Here the 20 first such values:

Code:

185 or 187, x=31 (remainder=1)
245 or 247, x=41 (remainder=5)
425 or 427, x=71 (remainder=5)
473 or 475, x=79 (remainder=1)
533 or 535, x=89 (remainder=5)
581 or 583, x=97 (remainder=1)
713 or 715, x=119 (remainder=5)
833 or 835, x=139 (remainder=1)
869 or 871, x=145 (remainder=1)
893 or 895, x=149 (remainder=5)
1001 or 1003, x=167 (remainder=5)
1073 or 1075, x=179 (remainder=5)
1145 or 1147, x=191 (remainder=5)
1157 or 1159, x=193 (remainder=1)
1253 or 1255, x=209 (remainder=5)
1265 or 1267, x=211 (remainder=1)
1337 or 1339, x=223 (remainder=1)
1505 or 1507, x=251 (remainder=5)
1517 or 1519, x=253 (remainder=1)
1589 or 1591, x=265 (remainder=1)

This code will take a long time to process them all...

Originally Posted by flp1969

If this is true, than criptographic algorithms like RSA would be very easy to break

It's also true that any prime > 2 is in the form 2k+1 for some non-negative integer k, but you can easily see why that doesn't mean primality testing is easy: there are also plenty of composite numbers in that form (i.e., all primes > 2 are odd, but plenty of odd numbers are composite!)

So, "if n % 6 doesn't equal 1 or 5 then it can't possibly be prime" is similiar to saying "if n is not odd then it can't possibly be prime": you can rule out composite numbers that way, but by itself it is insufficient to prove primality.

Originally Posted by flp1969

119 or 121, for example, aren't primes (6*20 ± 1): 17 * 7 = 119 and 11 * 11 = 121.

The statement is that "any prime p > 3 is always in the form 6x +- 1", not "any integer in the form 6x +- 1 is prime". Therefore, examples of composite numbers in this form are not counterexamples to the statement. Rather, you have to find a prime > 3 that is not of that form.

Originally Posted by flp1969

Do you have any proof?

I've seen it used by people other than Sir Galahad, but I suppose that's not proof

We could observe that any integer >= 6 could be expressed in the form 6k+0, 6k+1, 6k+2, 6k+3, 6k+4, or 6k+5, for some integer k > 0. Clearly, 6k+0, 6k+2, and 6k+4 are divisible by 2, and hence cannot be prime. 6k+3 is divisible by 3, and hence cannot be prime. Therefore, only those of the form 6k+1 and 6k+5 could be prime, i.e., if an integer >= 6 is prime, it must be of one of these two forms. But 6k+5 is congruent to 6k-1 under modulo 6, hence the statement expresses this as 6k+-1. Furthermore, 5 is also prime yet is in the form 6k-1, hence the range was specified as being > 3 rather than >= 6.

Originally Posted by flp1969

And, correct me if I'm wrong... testing n, getting the remainder of the division by 6 (if it is 1 or 5) only checks if n isn't divisible by 3 (since my routine already ruled out even numbers, except 2). Why 6 then? And if somebody is testing by divisors of 3, how about 5, 7, 11, 13, 17, ... ? Do you want a huge sequence of ifs?

Valid point. But still a little cheaper to use the 6x +- 1 identity. With the sqrt() call removed it might look something like this:

Code:

int prime(unsigned n)
{
 if(n <= 3)
  return n >= 2;
 unsigned d = n % 6;
 if(d != 1 && d != 5)
  return 0; 
 for(unsigned f = 5; f * f <= n; f += 2)
  if(n % f == 0)
   return 0;
 return 1;
}

Originally Posted by flp1969

It is not a reliable criteria because of this first reason and because you will need to test for primalty AGAIN using another ad hoc criteria...

Think about how the fact that all primes > 2 are of the form 2k+1 is commonly used in primality testing: we simply check for 2 as a possible divisor, then generate odd numbers as potential factors to test, because we know that all possible prime factors of the number we're testing is odd. It's not optimal compared to checking with primes only, but it's more efficient compared to checking with all integers less <= the square root of the number we're testing.

Likewise, making use of the 6k+-1 form, we can test for 2 and 3 as potential factors, then generate numbers in the form 6k+-1 to test.

If we're generating a list of primes through brute force testing rather than a sieve, likewise we can make use of the fact that they will be of this form by only generating candidates to test that are of this form.

The way I see it, the larger the multiplier we choose, the more this becomes like incorporating a sieve into the brute force test.

Originally Posted by Zeus_

...err, 25 is prime cuz remainder is 1, 35 is prime cuz remainder is 1, 49 is prime, 55 is prime....

Sigh. As per my previous post, "all primes greater than n are of this form" is not the same as "all integers greater than n of this form are prime".

Originally Posted by flp1969

Why 6 then? And if somebody is testing by divisors of 3, how about 5, 7, 11, 13, 17, ... ? Do you want a huge sequence of ifs?

I think the reason why 6 is used is because its the smallest example of this such that the multiplier is not a prime: obviously we can use forms with a larger multiplier instead, but depending on the range we have in mind, at some point it's just easier to (directly use a sieve to?) generate a list of primes for optimal brute force primality testing.

Originally Posted by Sir Galahad

Code:

int prime(unsigned n)
{
 if(n <= 3)
  return n >= 2;
 unsigned d = n % 6;
 if(d != 1 && d != 5)
  return 0; 
 for(unsigned f = 5; f * f <= n; f += 2)
  if(n % f == 0)
   return 0;
 return 1;
}

Hummmm.... interesting approach!! I have to study a little bit more to try to catch any flaws... awesome!
For now I can only see three potential problems... doing 'f*f' can result in an overflow if we got an value for 'n' bit enough. Let's say we want to check if 4294967293 is a prime or not, the maximum value of f, without causing an overflow, is 65535, which will result in f*f being 4294836225, but for f=65537, f*f will wrap to 131073. This way, for big 'n', well get, probably, a lot more comparisons (and with wrong values) than using a square root. This is easy to solve (I think), all you have to do is use a larger precision for f (unsigned long long, for instance)...

The second one is that I did a square root calculation ONCE. Here you do multiple multiplications inside the loop. Let's say MUL takes 4 cycles and we have 32000 iterations. This will loose 128000 clock cycles only to calculate f*f. It can be a little worse. Remember MUL instruction multiplies EAX with another operand resulting in EDX:EAX (IMUL is more flexible)...

The third is when you check for even values, you are using the loop... let's say I want to check if 200 is prime, in this case n > 3 and d=2... My routine checks if the value itself is even and stops there if isn't 2...

I believe your routine is way slower than mine in a long run (except if that test with n%6 catches the non primal value)...

Originally Posted by flp1969

For now I can only see three potential problems... doing 'f*f' can result in an overflow if we got an value for 'n' bit enough. Let's say we want to check if 4294967293 is a prime or not, the maximum value of f, without causing an overflow, is 65535, which will result in f*f being 4294836225, but for f=65537, f*f will wrap to 131073. This way, for big 'n', well get, probably, a lot more comparisons (and with wrong values) than using a square root. This is easy to solve (I think), all you have to do is use a larger precision for f (unsigned long long, for instance)...

The second one is that I did a square root calculation ONCE. Here you do multiple multiplications inside the loop. Let's say MUL takes 4 cycles and we have 32000 iterations. This will loose 128000 clock cycles only to calculate f*f. It can be a little worse. Remember MUL instruction multiplies EAX with another operand resulting in EDX:EAX (IMUL is more flexible)...

The third is when you check for even values, you are using the loop

Makes sense. Not sure about that last one though. Even values would be rejected by the test before the loop.

flp1969: following your example from post #16, we're probably looking at something along these lines:

Code:

int isPrime(unsigned int n)
{
    if (n < 2)
        return 0;
    if (n <= 3)
        return 1;
    if (n % 2 == 0 || n % 3 == 0)
        return 0;

    unsigned int sqrt_ = (unsigned int)sqrt(n);

    // We only need to check divisors of the forms 6k-1 and 6k+1, for k > 0
    for (unsigned int i = 5; i <= sqrt_; i += 4) // +4 because +2 in the loop body
    {
        if (n % i == 0)
            return 0;
        i += 2
        if (n % i == 0)
            return 0;
    }
    return 1;
}

I guess it's also possible to construct this more carefully with a check for the square root after each modulo test, but the idea remains the same: we're skipping numbers with an alternating +4 instead of always doing +2 (or +1) because we know that the possible prime factors to test with are of those two forms.

Originally Posted by laserlight

flp1969: following your example from post #16, we're probably looking at something along these lines:

This is essentially the same routine I posted before, except you check for divisors by 3 in the beginning (starting the checks in the loop, for divisors by 5)... A better way, which can consume a lot of memory, is to check for the previously found primes less or equal to the square rooot of n (a modified version of the sieve of Erasthotenes)...