martedì 23 agosto 2011

i numeri primi di project euler? OpenMP è un buon alleato

Oggi è una giornata molto calda, non avevo niente da fare e ho ripreso in mano project euler. Fino ad oggi avevo evitato tutti gli esercizi con i numeri primi, e quindi ho deciso di cominciare: esercizio 7, tocca a te!

L'esercizio sette di project euler recita:

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.
What is the 10 001st prime number?


Come sempre nei problemi di project euler le richieste possono sembrare esagerate per essere completate nel giro di un minuto, ma ce la si può fare :D

Avevo bisogno di risultati pronti per la libreria dinamica che stavo sviluppando e... Wolfram|Alpha ha sempre una risposta!: ebbravo il motore computazionale di conoscenza di Stephen, potendo sfruttare Mathematica per lui è fin troppo semplice :D

Comunque, necessitavo di una libreria di verifica dei numeri primi... l'ho scritta, sfruttando OpenMP.
OpenMP è un'interessante libreria, disponibile per C/C++ e Fortran, che permette di sfruttare il parallelismo senza modificare pesantemente o dover riscrivere l'applicazione, troverete (probabilmente) un mio articolo su questa libreria nel numero di settembre di UnderAttHack, per maggiori dettagli leggetelo :D


//codice della DLL, test.cpp
//il file.h contiene solamente i prototipi delle due funzioni
//compilare con g++ -shared -fPIC -fopenmp test.cpp -o test.so
#define DLLEXPORT extern "C"
#include
DLLEXPORT int isParallelPrime(long long int n)
{
     int isprime=1;
     if((n&1)==0) isprime=0;
     long long int limit=(long long int)pow(n, 0.5)+1;
     #pragma omp parallel for
     for(long long int i=3; i        if (!(n%i))
        {
           #pragma omp atomic
           isprime*=0;
        }

     return isprime;
}

DLLEXPORT int isSerialPrime(long long int n)
{
     int isprime=1;
     if ((n&1)==0) isprime=0;
     for(long long int i=3; i<(long long int)pow(n, 0.5)+1; i+=2)
        if (!(n%i))
            isprime=0;
     return isprime;
}


il codice Python che risolve effettivamente il problema 7

#!/usr/bin/python
import sys
from ctypes import cdll

#usage: python ex7.py 10001
if __name__=="__main__":
    n=int(sys.argv[1])
    mydll=cdll.LoadLibrary("./test.so")
    nprimes=0
    candidate_prime=1
    while nprimes:

        if mydll.isParallelPrime(candidate_prime)
            nprimes+=1
        if nprimes2:
            candidate_prime+=2
        elif nprimes            candidate+=1
    #prints result.
    print nprimes, n, i


voglio mostrarvi i tempi (lies, damn lies and benchmarks) di esecuzione di ex7.py (AMD Phenom II X4 955 su Fedora 15 con kernel 2.6.40)

[alfateam123@alfateam123 testDll]$ time python ex7.py 10001 #usando isParallelPrime
10001 10001 104743

real    0m0.227s
user    0m0.707s
sys     0m0.038s
[alfateam123@alfateam123 testDll]$ time python ex7.py 10001 #usando isSerialPrime
10001 10001 104743

real    0m0.713s
user    0m0.695s
sys     0m0.012s


In conclusione: non avevo niente di meglio da fare che mostrarvi il frutto delle mie fatiche mentali, non mi interessa che le apprezziate o meno. Divertitevi. Ah, e non perdetevi il mio (probabile) articolo su UnderAttHack! :D