/** * Computes Highly Wilsonian primes using a * method that is half as dumb as the dumbest * method. * * Optimized for a dual-core CPU (uses 2 threads). * * Compile using * g++ -O3 -pthread bf.cpp * * Author: Igor Naverniouk (abednego@gmail.com), * Date: May 9, 2007. **/ #include <stdio.h> #include <string.h> #include <pthread.h> #include <algorithm> using namespace std; // We will only check primes smaller than N. #define N 100000000 // This should be larger than N/log(N) #define NP 6000000 // The number of primes and the primes themselves. long long np, pr[NP]; // The two threads we'll be running. pthread_t a, b; // Erastosthenes to the rescue! void sieve() { bool *p = new bool[N]; memset( p, -1, N ); np = 0; p[0] = p[1] = false; for( long long i = 2; i < N; i++ ) if( p[i] ) { pr[np++] = i; for( long long j = i * i; j < N; j += i ) p[j] = false; } delete [] p; } // The values of w(n) that we have found so far (or -1 if unknown). long long w[128]; // This computes the w(p) using the dumbest possible method. inline int _( long long p ) { if( p < 5 ) return p; int ans = 4; long long x = 1, m1 = p - 1, h = (p-1)/2; for( long long i = 2; i < h; i++ ) { x = ( x * i ) % p; if( x == 1 || x == m1 ) ans += 2; } if( h + h == p - 1 ) { x = ( x * h ) % p; if( x == 1 || x == m1 ) ans++; } return ans; } // The next prime we should check is pr[nexti]. long long nexti; pthread_mutex_t nexti_mutex = PTHREAD_MUTEX_INITIALIZER; // A thread that checks primes in order. We'll have two of these. void *eater( void *x ) { while( nexti < np ) { // Grab a prime p and compute myw = w(p). pthread_mutex_lock( &nexti_mutex ); long long p = pr[nexti++]; // Print a message every 100,000 numbers, so we know where we are. if( nexti > 10 && pr[nexti - 1] / 100000 > pr[nexti - 2] / 100000 ) { printf( "checking p = %Ld\n", pr[nexti - 1] ); fflush( stdout ); } pthread_mutex_unlock( &nexti_mutex ); int myw = _( p ); // Remember it if it's something we didn't have before. if( w[myw] < 0 ) { w[myw] = p; printf( "%d found w( %Ld ) = %d\n", ( int )x, p, myw ); fflush( stdout ); } } pthread_exit( NULL ); } int main( int argc, char **argv ) { // First, we find some primes. printf( "Sieving...\n" ); fflush( stdout ); sieve(); printf( "Testing the first %Ld primes: 2 -- %Ld\n", np, pr[np - 1] ); fflush( stdout ); memset( w, -1, sizeof( w ) ); nexti = 0; if( argc == 2 ) nexti = lower_bound( pr, pr + np, atoi( argv[1] ) ) - pr; // Then we make two threads to test them in parallel. pthread_create( &a, NULL, eater, &a ); pthread_create( &b, NULL, eater, &b ); pthread_exit( NULL ); }