FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

bool isPrime(uint n)

 
Post new topic   Reply to topic     Forum Index -> Deimos
View previous topic :: View next topic  
Author Message
parabolis



Joined: 02 Aug 2004
Posts: 4

PostPosted: Mon Aug 02, 2004 12:18 pm    Post subject: bool isPrime(uint n) Reply with quote

Here you are. Let me know if you find bugs.
Code:

bool isPrime(uint n)
{
    if ((n & 1) == 0) return (n == 2);

    if (n < 0x10000)
    {
        uint i = n >>> 4;
        uint mask = 1 << ((n & 0xE) >>> 1);
        return ((SMALL_PRIMES[i] & mask) != 0);
    }

    uint i, j, m, r;
    uint[8] q;

    r = n / 3;
    if (n == r * 3) return false;
    r = n / 5;
    if (n == r * 5) return false;

    for (uint p=0; p<0x10000; p+=30)
    {
        q[0] = p + 1;
           i = q[0] >>> 4;
           m = 1 << ((q[0] & 0xE) >>> 1);
        q[0] = ((SMALL_PRIMES[0] & m) != 0);

        q[1] = p + 7;
           i = q[1] >>> 4;
           m = 1 << ((q[1] & 0xE) >>> 1);
        q[1] = ((SMALL_PRIMES[i] & m) != 0);

        q[2] = p + 11;
           i = q[2] >>> 4;
           m = 1 << ((q[2] & 0xE) >>> 1);
        q[2] = ((SMALL_PRIMES[i] & m) != 0);

        q[3] = p + 13;
           i = q[3] >>> 4;
           m = 1 << ((q[3] & 0xE) >>> 1);
        q[3] = ((SMALL_PRIMES[i] & m) != 0);

        q[4] = p + 17;
           i = q[4] >>> 4;
           m = 1 << ((q[4] & 0xE) >>> 1);
        q[4] = ((SMALL_PRIMES[i] & m) != 0);

        q[5] = p + 19;
           i = q[5] >>> 4;
           m = 1 << ((q[5] & 0xE) >>> 1);
        q[5] = ((SMALL_PRIMES[i] & m) != 0);

        q[6] = p + 23;
           i = q[6] >>> 4;
           m = 1 << ((q[6] & 0xE) >>> 1);
        q[6] = ((SMALL_PRIMES[i] & m) != 0);

        q[7] = p + 29;
           i = q[7] >>> 4;
           m = 1 << ((q[7] & 0xE) >>> 1);
        q[7] = ((SMALL_PRIMES[i] & m) != 0);

        for( j = 0; j <= 8; j++ )
        {
            if (q[j])
            {
                r = n / q[j];
                if (n == r * q[j]) return false;
                if (r < q[j]) return true;
            }
        }
    }
    return true;
}
[/code]
Back to top
View user's profile Send private message
Nick



Joined: 03 Aug 2004
Posts: 2

PostPosted: Tue Aug 03, 2004 1:38 pm    Post subject: Reply with quote

Looks very good! Personally think it could need some comenting, though...

Also, you should replace
Code:

r = n / q;
if ( n == r * q ) return false;

with
Code:

if ( n ? q == 0 ) return false;

This gives a whooping 10? speed increase! (Not of the entire function, only the part you replace Smile )
Back to top
View user's profile Send private message
Nick



Joined: 03 Aug 2004
Posts: 2

PostPosted: Tue Aug 03, 2004 2:41 pm    Post subject: Reply with quote

Oh and I found a few trivial bugs. First, you store the current divider in q[i] for use later, but then overwrite it with the bool indicating if it is a prime. Second, in the first block in the loop, it should read SMALL_PRIMES[i], not SMALL_PRIMES[0]. Lastly, the inner loop should stop at 7 (the loop condition should be j < 8, not j <= 8.) Other than that it works great.

Also, if you do the change I mentioned in the last post, it will be more efficient to compare q[j] with a precalculated square root of n in the very last if statement instead of calculating r = n / q.

Hope this helps Wink
Back to top
View user's profile Send private message
parabolis



Joined: 02 Aug 2004
Posts: 4

PostPosted: Tue Aug 03, 2004 3:44 pm    Post subject: Reply with quote

Thanks Nick. I have fixed most of the bugs you pointed out and a couple others that became obvious after your comments.

SMALL_PRIMES[i] was bad to start with because i is not used.

I added a comment that explains why the loop works the way it does but otherwise it is not obvious to me where any other comments belong. Any suggestions or bug sightings?

Code:

bool isPrime(uint n)
{
    if ((n & 1) == 0) return (n == 2);

    if (n < 0x10000)
    {
        uint i = n >>> 4;
        uint mask = 1 << ((n & 0xE) >>> 1);
        return ((SMALL_PRIMES[i] & mask) != 0);
    }


    uint i, j, m. sqrtN = (n / n) + 1;
    uint[8] q; // check for primality
    bit[8] b;  // remembers if q[i] is prime

    if (n ? 3 == 0) return false;
    if (n ? 5 == 0) return false;

    // Make use of the fact that since 30 = 2*3*5
    // we only need to check the following every
    // iteration: 1,7,11,13,17,19,23,29
    for (uint p=0; p<0x10000; p+=30)
    {
        q[0] = p + 1;
           i = q[0] >>> 4;
           m = 1 << ((q[0] & 0xE) >>> 1);
        b[0] = ((SMALL_PRIMES[p+1] & m) != 0);


        q[1] = p + 7;
           i = q[1] >>> 4;
           m = 1 << ((q[1] & 0xE) >>> 1);
        b[1] = ((SMALL_PRIMES[p+7] & m) != 0);


        q[2] = p + 11;
           i = q[2] >>> 4;
           m = 1 << ((q[2] & 0xE) >>> 1);
        b[2] = ((SMALL_PRIMES[p+11] & m) != 0);


        q[3] = p + 13;
           i = q[3] >>> 4;
           m = 1 << ((q[3] & 0xE) >>> 1);
        b[3] = ((SMALL_PRIMES[p+13] & m) != 0);


        q[4] = p + 17;
           i = q[4] >>> 4;
           m = 1 << ((q[4] & 0xE) >>> 1);
        b[4] = ((SMALL_PRIMES[p+17] & m) != 0);


        q[5] = p + 19;
           i = q[5] >>> 4;
           m = 1 << ((q[5] & 0xE) >>> 1);
        b[5] = ((SMALL_PRIMES[p+19] & m) != 0);


        q[6] = p + 23;
           i = q[6] >>> 4;
           m = 1 << ((q[6] & 0xE) >>> 1);
        b[6] = ((SMALL_PRIMES[p+23] & m) != 0);


        q[7] = p + 29;
           i = q[7] >>> 4;
           m = 1 << ((q[7] & 0xE) >>> 1);
        b[7] = ((SMALL_PRIMES[p+29] & m) != 0);


        for( j = 0; j < 8; j++ )
        {
            if (b[j])
            {
                if (n ? q[j] == 0) return false;
                if (sqrtN < q[j]) return true;
            }
        }
    }
    return true;
}
Back to top
View user's profile Send private message
jcc7



Joined: 22 Feb 2004
Posts: 657
Location: Muskogee, OK, USA

PostPosted: Tue Aug 03, 2004 4:41 pm    Post subject: Reply with quote

Sounds cool. In your latest posting, there's a subtle typographical error (a period should be a comma):
Code:
uint i, j, m. sqrtN = (n / n) + 1;
should be
Code:
uint i, j, m, sqrtN = (n / n) + 1;


As to the commenting issue, you might add a comment around the beginning mentioning that values under 0x10000 are hard-coded in the table. The rest of it seems well-commented to me.
Back to top
View user's profile Send private message AIM Address
parabolis



Joined: 02 Aug 2004
Posts: 4

PostPosted: Tue Aug 03, 2004 9:40 pm    Post subject: Further optimization Reply with quote

First I want to thank jcc7 for the typo and comment suggestion. I fixed them in my code. I am not sure they are worth reposting the whole function again however. It is rather long...

I posted this at digitalmars but thought it might be more appropriate here:

I was just looking at mango's primes.d and figured out that the longest stretch of non-primes in a ushort is 52. I also gave some thought to a way to speed up the primeGreaterEqual and related functions and realized that a division and clever switch statement would allow the same speed up for these functions.

(The problem being that starting from 0 it is simple to maintain a window of 30, however it is trickier to fall into a window of 30 in the general case)

So given no two ushort primes have a distance greater than 60 you can guarantee that they will be found within 16 attempts. In mango's implementation there is a table of 1000 shorts and the next prime is found by using a binary search which has a worst case around 10 but the distribution of cases is the opposite of the distribution of prime distance - the majority of cases in bsearch are worst cases.

================================
Note: mango is apparently removing the primes.d file they had in 9.2 so the example that suggests the optimization no longer exists. However it would again roughly double the speed of the code for the getPrimeXXX functions.
Back to top
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic     Forum Index -> Deimos All times are GMT - 6 Hours
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group