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

  1. #16
    Registered User
    Join Date
    Feb 2019
    Posts
    773
    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'.
    Last edited by flp1969; 10-28-2019 at 04:39 PM.

  2. #17
    Registered User Sir Galahad's Avatar
    Join Date
    Nov 2016
    Location
    The Round Table
    Posts
    199
    Quote Originally Posted by flp1969 View Post
    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.

  3. #18
    misoturbutc Hodor's Avatar
    Join Date
    Nov 2013
    Posts
    1,795
    Quote Originally Posted by flp1969 View Post
    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.
    Last edited by Hodor; 10-29-2019 at 04:19 AM.

  4. #19
    Registered User
    Join Date
    Feb 2019
    Posts
    773
    Quote Originally Posted by Sir Galahad View Post
    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 6x1 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...
    Last edited by flp1969; 10-29-2019 at 04:59 AM.

  5. #20
    C++ Witch laserlight's Avatar
    Join Date
    Oct 2003
    Location
    Singapore
    Posts
    28,162
    Quote 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.

    Quote 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.

    Quote 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.
    Last edited by laserlight; 10-29-2019 at 05:34 AM.
    Quote Originally Posted by Bjarne Stroustrup (2000-10-14)
    I get maybe two dozen requests for help with some sort of programming or design problem every day. Most have more sense than to send me hundreds of lines of code. If they do, I ask them to find the smallest example that exhibits the problem and send me that. Mostly, they then find the error themselves. "Finding the smallest program that demonstrates the error" is a powerful debugging tool.
    Look up a C++ Reference and learn How To Ask Questions The Smart Way

  6. #21
    Registered User
    Join Date
    Feb 2019
    Posts
    773
    Quote Originally Posted by laserlight View Post
    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.
    Neither I said that... All I'm saying is there are a lot of x values, using this criteria, where they are NOT primes... 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...

    There are another attempts, of course (here), but I never seen this n mod 6...

  7. #22
    Registered User
    Join Date
    Aug 2019
    Location
    inside a singularity
    Posts
    309
    > So if n % 6 doesn't equal 1 or 5 then it can't possibly be prime.

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

    No, this check isn't going to give you correct results. You need to check till near the square root to get the correct result. However, you may add this check to eliminate the numbers that give you remainders other than 1 or 5 instead of checking further in a loop to make your code more efficient and I guess that's what you intended to imply. For the ones that do give you 1 or 5 as remainder, you'll have to check through a loop anyway...

  8. #23
    Registered User
    Join Date
    Feb 2019
    Posts
    773
    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?

  9. #24
    Registered User Sir Galahad's Avatar
    Join Date
    Nov 2016
    Location
    The Round Table
    Posts
    199
    Quote Originally Posted by flp1969 View Post
    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;
    }

  10. #25
    C++ Witch laserlight's Avatar
    Join Date
    Oct 2003
    Location
    Singapore
    Posts
    28,162
    Quote 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.

    Quote 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".

    Quote 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.
    Last edited by laserlight; 10-29-2019 at 12:59 PM.
    Quote Originally Posted by Bjarne Stroustrup (2000-10-14)
    I get maybe two dozen requests for help with some sort of programming or design problem every day. Most have more sense than to send me hundreds of lines of code. If they do, I ask them to find the smallest example that exhibits the problem and send me that. Mostly, they then find the error themselves. "Finding the smallest program that demonstrates the error" is a powerful debugging tool.
    Look up a C++ Reference and learn How To Ask Questions The Smart Way

  11. #26
    Registered User
    Join Date
    Feb 2019
    Posts
    773
    Quote Originally Posted by Sir Galahad View Post
    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)...
    Last edited by flp1969; 10-29-2019 at 01:44 PM.

  12. #27
    Registered User
    Join Date
    Aug 2019
    Location
    inside a singularity
    Posts
    309
    Quote Originally Posted by laserlight View Post
    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".
    I did mention something else ahead of "...err, 25 is prime cuz remainder is 1, 35 is prime cuz remainder is 1, 49 is prime, 55 is prime...." that completely agrees with what you said...

  13. #28
    Registered User Sir Galahad's Avatar
    Join Date
    Nov 2016
    Location
    The Round Table
    Posts
    199
    Quote Originally Posted by flp1969 View Post
    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.

  14. #29
    C++ Witch laserlight's Avatar
    Join Date
    Oct 2003
    Location
    Singapore
    Posts
    28,162
    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.
    Quote Originally Posted by Bjarne Stroustrup (2000-10-14)
    I get maybe two dozen requests for help with some sort of programming or design problem every day. Most have more sense than to send me hundreds of lines of code. If they do, I ask them to find the smallest example that exhibits the problem and send me that. Mostly, they then find the error themselves. "Finding the smallest program that demonstrates the error" is a powerful debugging tool.
    Look up a C++ Reference and learn How To Ask Questions The Smart Way

  15. #30
    Registered User
    Join Date
    Feb 2019
    Posts
    773
    Quote Originally Posted by laserlight View Post
    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)...
    Last edited by flp1969; 10-29-2019 at 03:00 PM.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. code for getting the GCD (greatest common divisor)
    By Noobpoint in forum C Programming
    Replies: 1
    Last Post: 03-07-2012, 07:13 AM
  2. The greatest common divisor (GCD) help..
    By damha in forum C Programming
    Replies: 4
    Last Post: 04-09-2011, 05:18 AM
  3. Greatest Common Divisor.....
    By muran_pling in forum C++ Programming
    Replies: 10
    Last Post: 12-18-2006, 05:02 AM
  4. Greatest Common Divisor problem
    By fenixataris182 in forum C++ Programming
    Replies: 8
    Last Post: 07-12-2005, 07:55 PM
  5. Greatest common divisor
    By wiz23 in forum C++ Programming
    Replies: 5
    Last Post: 04-13-2005, 04:50 PM

Tags for this Thread