Author Message
parabolis

Joined: 02 Aug 2004
Posts: 4 Posted: Mon Aug 02, 2004 12:18 pm    Post subject: bool isPrime(uint n) 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 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 = p + 1;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         q = ((SMALL_PRIMES & m) != 0);         q = p + 7;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         q = ((SMALL_PRIMES[i] & m) != 0);         q = p + 11;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         q = ((SMALL_PRIMES[i] & m) != 0);         q = p + 13;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         q = ((SMALL_PRIMES[i] & m) != 0);         q = p + 17;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         q = ((SMALL_PRIMES[i] & m) != 0);         q = p + 19;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         q = ((SMALL_PRIMES[i] & m) != 0);         q = p + 23;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         q = ((SMALL_PRIMES[i] & m) != 0);         q = p + 29;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         q = ((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]   Nick

Joined: 03 Aug 2004
Posts: 2 Posted: Tue Aug 03, 2004 1:38 pm    Post subject: 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 )   Nick

Joined: 03 Aug 2004
Posts: 2 Posted: Tue Aug 03, 2004 2:41 pm    Post subject: 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. 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    parabolis

Joined: 02 Aug 2004
Posts: 4 Posted: Tue Aug 03, 2004 3:44 pm    Post subject: Thanks Nick. I have fixed most of the bugs you pointed out and a couple others that became obvious after your comments.

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 q; // check for primality     bit 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 = p + 1;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         b = ((SMALL_PRIMES[p+1] & m) != 0);         q = p + 7;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         b = ((SMALL_PRIMES[p+7] & m) != 0);         q = p + 11;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         b = ((SMALL_PRIMES[p+11] & m) != 0);         q = p + 13;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         b = ((SMALL_PRIMES[p+13] & m) != 0);         q = p + 17;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         b = ((SMALL_PRIMES[p+17] & m) != 0);         q = p + 19;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         b = ((SMALL_PRIMES[p+19] & m) != 0);         q = p + 23;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         b = ((SMALL_PRIMES[p+23] & m) != 0);         q = p + 29;            i = q >>> 4;            m = 1 << ((q & 0xE) >>> 1);         b = ((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; }   jcc7

Joined: 22 Feb 2004
Posts: 657
Location: Muskogee, OK, USA Posted: Tue Aug 03, 2004 4:41 pm    Post subject: 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.     parabolis

Joined: 02 Aug 2004
Posts: 4 Posted: Tue Aug 03, 2004 9:40 pm    Post subject: Further optimization 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.   Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First
 All times are GMT - 6 Hours Page 1 of 1