/**
 * 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 );
}