3
int prime(unsigned long long n){
    unsigned val=1, divisor=7;
    if(n==2 || n==3) return 1; //n=2, n=3 (special cases).
    if(n<2 || !(n%2 && n%3)) return 0; //if(n<2 || n%2==0 || n%3==0) return 0;
    for(; divisor<=n/divisor; val++, divisor=6*val+1) //all primes take the form 6*k(+ or -)1, k[1, n).
        if(!(n%divisor && n%(divisor-2))) return 0; //if(n%divisor==0 || n%(divisor-2)==0) return 0;
    return 1;
}

The code above is something a friend wrote up for getting a prime number. It seems to be using some sort of sieving, but I'm not sure how it exactly works. The code below is my less awesome version. I would use sqrt for my loop, but I saw him doing something else (probably sieving related) and so I didn't bother.

int prime( unsigned long long n ){
    unsigned i=5;
    if(n < 4 && n > 0)
        return 1;
    if(n<=0 || !(n%2 || n%3))
        return 0;
    for(;i<n; i+=2)
        if(!(n%i)) return 0;
    return 1;
}

My question is: what exactly is he doing?

Jonathan Leffler
  • 730,956
  • 141
  • 904
  • 1,278
Painguy
  • 565
  • 6
  • 13
  • 25
  • 1
    i've asked my friend already, but I didn't quite understand his explanation. – Painguy Apr 23 '12 at 04:38
  • 3
    @dan04: **"My question is what exactly is he doing?"** Tip: Whenever a statement ends with a Question mark(`?`) it is a Question. – Alok Save Apr 23 '12 at 04:41
  • Not really sieving. Just a slightly more sophisticated version of trial division that avoids testing with numbers that couldn't possibly be divisors (an extended version of your code not trying to test odd numbers against even divisors, since you know they won't divide evenly). – Jerry Coffin Apr 23 '12 at 04:48
  • 1
    There seems to be a problem with the `prime()` function as posted - it identifies some multiples of 5 as prime. – Michael Burr Apr 23 '12 at 04:53

5 Answers5

5

Your friend's code is making use of the fact that for N > 3, all prime numbers take the form (6×M±1) for M = 1, 2, ... (so for M = 1, the prime candidates are N = 5 and N = 7, and both those are primes). Also, all prime pairs are like 5 and 7. This only checks 2 out of every 3 odd numbers, whereas your solution checks 3 out of 3 odd numbers.

Your friend's code is using division to achieve something akin to the square root. That is, the condition divisor <= n / divisor is more or less equivalent to, but slower and safer from overflow than, divisor * divisor <= n. It might be better to use unsigned long long max = sqrt(n); outside the loop. This reduces the amount of checking considerably compared with your proposed solution which searches through many more possible values. The square root check relies on the fact that if N is composite, then for a given pair of factors F and G (such that F×G = N), one of them will be less than or equal to the square root of N and the other will be greater than or equal to the square root of N.


As Michael Burr points out, the friend's prime function identifies 25 (5×5) and 35 (5×7) as prime, and generates 177 numbers under 1000 as prime whereas, I believe, there are just 168 primes in that range. Other misidentified composites are 121 (11×11), 143 (13×11), 289 (17×17), 323 (17×19), 841 (29×29), 899 (29×31).

Test code:

#include <stdio.h>

int main(void)
{
    unsigned long long c;

    if (prime(2ULL))
        printf("2\n");
    if (prime(3ULL))
        printf("3\n");
    for (c = 5; c < 1000; c += 2)
        if (prime(c))
            printf("%llu\n", c);
    return 0;
}

Fixed code.

The trouble with the original code is that it stops checking too soon because divisor is set to the larger, rather than the smaller, of the two numbers to be checked.

static int prime(unsigned long long n)
{
    unsigned long long val = 1;
    unsigned long long divisor = 5;

    if (n == 2 || n == 3)
        return 1;
    if (n < 2 || n%2 == 0 || n%3 == 0)
        return 0;
    for ( ; divisor<=n/divisor; val++, divisor=6*val-1)
    {
        if (n%divisor == 0 || n%(divisor+2) == 0)
            return 0;
    }
    return 1;
}

Note that the revision is simpler to understand because it doesn't need to explain the shorthand negated conditions in tail comments. Note also the +2 instead of -2 in the body of the loop.

Community
  • 1
  • 1
Jonathan Leffler
  • 730,956
  • 141
  • 904
  • 1,278
  • divisor * divisor <= n might overflow because the calculated number would become to large correct? – Painguy Apr 23 '12 at 05:18
  • Yes, overflow is why the division is safer but slower. – Jonathan Leffler Apr 23 '12 at 05:20
  • Thanks for posting the fix. I wish I could give another +1 for the revisions to the negated conditions, which gave me a headache trying to follow, even though the comments were right there. I'm guessing that the less readable versions were used in a misguided micro optimization attempt. – Michael Burr Apr 23 '12 at 06:04
  • @MichaelBurr: yes, I think that the negated stuff was an attempt at micro-optimization, and the use of one line instead of two was probably similarly motivated (readability was not a concern, that's for sure). – Jonathan Leffler Apr 23 '12 at 06:11
  • It could honestly be a literal translation of how he was thinking about it in his head. "Return true if it's 2 or 3. Return false if it's 1, or if a test that n%2 and n%3 are nonzero fails. ..." –  Apr 23 '12 at 06:27
  • @Hurkyl: You're right; imputing reasons from the code is hard and maybe I should give them the benefit of the doubt. (The jibe about 'one line instead of two' was intended to be humorous.) However, given the tail comments, someone (friend or OP) felt the need to translate the negative conditions into the corresponding positive conditions. And that indicates that someone (other than just Michael and me) finds the negatives too hard to understand for comfort. – Jonathan Leffler Apr 23 '12 at 06:31
  • 1
    Nice answer. I would do away with `val` in the loop altogether and just `divisor += 6`, seems more obvious to me. – Daniel Fischer Apr 23 '12 at 09:06
2

He's checking for the basis 6k+1/6k-1 as all primes can be expressed in that form (and all integers can be expressed in the form of 6k+n where -1 <= n <= 4). So yes it is a form of sieving.. but not in the strict sense.

For more: http://en.wikipedia.org/wiki/Primality_test

hkf
  • 4,440
  • 1
  • 30
  • 44
1

In case the 6k+-1 portion is confusing, note that you can perform some factorization of most forms of 6k+n and some are obviously composite and some need to be tested.

Consider numbers:
6k + 0 -> composite
6k + 1 -> not obviously composite
6k + 2 -> 2(3k+1) --> composite
6k + 3 -> 3(2k+1) --> composite
6k + 4 -> 2(3k+2) --> composite
6k + 5 -> not obviously composite

I've not seen this little trick before, so it's neat, but of limited utility since a sieve of Eratosthenese is more efficient for finding many small prime numbers, and larger prime numbers benefit from faster, more intelligent, tests.

sarnold
  • 102,305
  • 22
  • 181
  • 238
  • Thanks for posting an explanation that even someone like me, whose math skills have regressed to middle school algebra levels, can completely comprehend. – Michael Burr Apr 23 '12 at 05:25
  • @sarnold: You greatly underestimate its utility. This trick lets you implement the sieve of Eratosthenes nearly **three times faster** than the straight implementation, because it means you only have to sieve on 1/3 of the numbers. Carrying out this idea a few more steps (e.g. mod 30, or mod 210, or more) is **central** to the fastest algorithms for generating all small prime numbers. Search for "wheel sieve". –  Apr 23 '12 at 06:32
  • @Hurkyl: _three times faster_ is compared against checking every number, including even numbers. No one will do that... Compared against checking every odd number, this only checks 2/3 odd numbers, which is certainly a useful improvement -- 50% speedup. I had wondered whether or not the benefits of the approach would be _diminishing returns_ or _incredibly useful_ with larger modulo terms -- thanks for the report that this is in fact as useful as it sounds. :) – sarnold Apr 24 '12 at 21:45
0
#include<stdio.h>
int main()
{
    int i,j;
    printf("enter the value :");
    scanf("%d",&i);

    for (j=2;j<i;j++)
    {
        if (i%2==0 || i%j==0)
        {
            printf("%d is not a prime number",i);
            return 0;
        }
        else
        {
            if (j==i-1)
            {
                printf("%d is a prime number",i);
            }
            else
            {
                continue;
            }
        }
    }
}
-1
#include<stdio.h>

int main()
{
   int n, i = 3, count, c;

   printf("Enter the number of prime numbers required\n");
   scanf("%d",&n);

   if ( n >= 1 )
   {
      printf("First %d prime numbers are :\n",n);
      printf("2\n");
    }

   for ( count = 2 ; count <= n ;  )
   {
      for ( c = 2 ; c <= i - 1 ; c++ )
      {
     if ( i%c == 0 )
        break;
  }
      if ( c == i )
      {
         printf("%d\n",i);
         count++;
      }
      i++;
   }

   return 0;
}
Joya
  • 1
  • 1