ANSI C: Test di numeri primi con la Wheel Factorization

Crivello di eratostene: animazioneCosa sono i numeri primi e perché tanto interesse intorno a loro? Un numero n > 0 si definisce primo se è divisibile solo per 1 e per sé stesso con resto 0. L’interesse è dovuto a due motivi, principalmente:

1) I numeri primi sono le basi di tutta l’aritmetica come si può evincere dal Teorema Fondamentale della matematica che recita:“Ogni intero si fattorizza in modo unico come prodotto di numeri primi” e questo li rende indubbiamente affascinanti.
2) I numeri primi molto grandi (ordine di 21024), sono alla base del più utilizzato sistema di crittografia: l’RSA questo invece, attira molto gli esperti in sicurezza informatica…!

Come si trovano i numeri primi?

Fu Eratostene di Cirene (276 – 194 a.c.), il primo a proporre un metodo per “setacciare” i numeri primi, detto appunto: crivello di Eratostene. Esistono diversi algoritmi che risolvono il crivello di Eratostene, questo che propongo mi sembra il più efficiente. Supponiamo dunque di voler trovare tutti i numeri primi da 1 a n, l’algoritmo in pseudocodice è il seguente:


for i := 2 to n-1 do
 a[i] := 1
for i := 2 to n-1 do
 if a[i]
  for j := i to (j*i)<n do
   a[i*j] := 0
for i := 2 to n-1 do
 if a[i]
  print i

In pratica si procede da i=2 a n, eliminando tutti i multipli di i
(nella immagine in alto a sinistra è possibile vedere un’animazione del processo). Va detto che ancora oggi, pur non essendo l’algoritmo più efficiente in assoluto, è molto utilizzato quando si lavora su numeri relativamente piccoli, per la sua estrema semplicità.

Il metodo della ruota di fattorizzazione (Wheel Factorization)

Già il crivello di Eratostene permette un notevole incremento di prestazioni rispetto ad un processo di “forza bruta” specialmente se l’algoritmo è opportunamente ottimizzato. Il metodo della ruota risulta ancora più efficiente, e vediamo il perché:
Si scelgono inizialmente come numeri primi di partenza: 2, 3, 5 per ottenere una ruota di lunghezza 30. La lunghezza della ruota è determinta dal prodotto dei numeri primi che decidiamo di utilizzare. Nel nostro caso 2 * 3 * 5 = 30. A questo punto troviamo tutti gli interi che non sono multipli di 2, 3 e 5, il vettore risultate sarà: s[] = {1, 7, 11, 13, 17, 19, 23, 29 }. Abbiamo ottenuto un insieme di numeri che, non sono multipli di 2, 3, 5 e anche se sommati alla lunghezza della ruota o ai suoi multipli (30, 60, 90, ecc), non daranno mai multipli di 2, 3, 5.
Il concetto è più chiaro osservando la figura sotto, le colonne gialle della tabella contengono i soli numeri che andremo a testare, i rimanenenti non sono sicuramente primi, esclusi i numeri di partenza (2, 3, 5). Se provassimo a visualizzare ciascuna riga della tabella come una circonferenza, potremmo immaginare la tabella come una ruota in cui le colonne gialle rappresentano i raggi dove cercare i numeri primi.

Tabella di fattorizzazione
Tabella di fattorizzazione

Il codice è tratto dall’articolo: “Prime Number Determination Using Wheel Factorization” di rickoshay originariamente scritto in C#, e da me convertito in ANSI C


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#include <stdio .h>
#include <stdlib .h>
 
unsigned long FirstDivisor(unsigned long candidatePrime);
int isPrime(unsigned long primeSuspect);
 
unsigned long FirstDivisor(unsigned long candidatePrime)
{
 const int wheelFactor = 30;
 int aSieve30[] = {7, 11, 13, 17, 19, 23, 29, 31};
 int firstPrimes[] = {2, 3, 5}, i;
 unsigned long sqrtmax, pass;
 
 if (candidatePrime == 1)
 {  
    return 0;
 }
 for (i=0;i<sizeof (firstPrimes)/sizeof(int);i++) {
     if (candidatePrime % firstPrimes[i] == 0) {
        return firstPrimes[i];
     }
 }
 
 sqrtmax  = (unsigned long) sqrt(candidatePrime);
 for (pass=0; pass<sqrtmax; pass+=wheelFactor) {
  for (i=0;i<sizeof(aSieve30)/sizeof(int);i++) {
      if (candidatePrime % (pass + aSieve30[i]) == 0) {
         return pass + aSieve30[i];
      }
  }   
 }
 return candidatePrime; 
}
 
int isPrime(unsigned long primeSuspect)
{
   if (primeSuspect == 0) {
      return 0;
   }
   if (FirstDivisor(primeSuspect) == primeSuspect) {
      return 1;
   }
   else {
      return 0;
   }
}
 
int main(int argc, char *argv[])
{
  unsigned long x, i;
 
  if (argc == 2) {
     x = strtoul(argv[1],NULL,10);
     if (isPrime(x)) {
        printf("%u e' un numero primo\n",x);
     }
     else {
        printf("%u non e' un numero primo\n",x);        
     }     
  }
  else {
   printf ("\nUsage: %s n \nWhere n is an integer number to check for primality\n",argv[0]);
  }
  system("PAUSE");	
  return 0;
}

Download

Il download dell’applicazione (compilata per Win32) è disponibile qui, in formato compresso zip. Sono compresi i sorgenti e i files di progetto per il compilatore DevC++ (ver.: 4.9.9.2). Per l’utilizzo è sufficiente lanciare il programma da una console DOS e passare come parametro il numero da testare.

Conclusioni:

Adottare il metodo della Wheel Factorization per testare la primalità di un numero, consente un incremento di performance pari a (1 – 8 /30) * 100 = 73% con una ruota di dimensioni 30. Le prestazioni possono crescere aumentando le dimensioni della ruota, specialmente quando si ragiona con numeri abbastanza grandi. Il codice proposto utilizza al massimo numeri di 4 byte (unsigned long int = 4,294,967,295) e le performance rimangono accettabili ricompilando per numeri di 8 byte (unsigned long long = 18,446,744,073,709,551,615). In questo caso però sarebbe consigliabile utilizzare una ruota prodotta dai primi 8 numeri primi, come suggerito nell’articolo da cui ho tratto lo spunto per il codice.

Riferimenti ed approfondimenti:

ANSI C: giocare con la numerazione binaria per indovinare l’eta’

tavole per indovinare l'etàDiversi anni fa, ero alla ricerca di spunti divertenti e non troppo banali per insegnare agli allievi di un corso di programmazione in ANSI C
come funziona il sistema di numerazione binaria. Ho trovato su un vecchio libretto di curiosità matematiche un gioco per indovinare l’età attraverso l’utilizzo di tavole “magiche”. Va premesso che, ovviamente le tavole non sono magiche e vedremo il perché, e che per limiti di spazio, il gioco funziona solo per età inferiori ai 32 anni (la mia versione arriva a 64). Il gioco consiste nell’individuare la propria età in ciascuna delle colonne della figura a fianco, e segnare in quale di queste è presente. Mettiamo che individuate la vostra età nelle colonne 1,3,4,5 la vostra età è 29 anni.
Vediamo come sono composte le colonne, partendo da sinistra: la prima colonna è una serie di numeri dispari da 1 a 31, la seconda coppie alterne, la terza quadruple alterne, la quarta ottuple alterne e così via. A questo punto si intuisce che la tabella è un sistema per convertire i numeri decimali in cifre binarie, dove la prima colonna corrisponde a 20 nel sistema posizionale binario,
la seconda a 21, la terza a 22 e così via. Per calcolare l’età, quindi, basta sommare le intestazioni delle colonne dove essa viene individuata.
La stringa individuata, indicando con 1 e con 0 rispettivamente la presenza e l’assenza, e nel caso di 29 anni è 10111 che rovesciata da 11101 ovvero il numero 29 in notazione binaria. E’ necessario capovolgere la stringa perché, nella numerazione posizionale, la prima colonna si deve trovare a destra.

Questo è il codice in linguaggio C per l’implementazione del gioco. (Per comodità di costruzione dell’algoritmo, le colonne sono state trasposte in righe)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#include 
#include 
 
int main(int argc, char *argv[])
{
  int i,j,a,b,tabella[6][32],anni=0;
  char risposta[6];    
  a=1;  
  tabella[0][0]=1;  
  for (i=0;i&lt;6;i++){      
      printf("\nriga %d:\t",i+1);
      b=1;            
      for (j=0;j&lt;32;j++){ if (j==0){ tabella[i][j]= a;} else {tabella[i][j] = tabella[i][j-1]+1;} if (b&gt;a){
             tabella[i][j] = tabella[i][j] + a;
             b=1;   
          }                           
          printf("%2d ",tabella[i][j]);
          if (!((j+1)%16)) printf("\n\t");
          b++;                   
      }      
      a=a&lt;&lt;1;      
  }    
  printf("\nIn quale di queste righe e' presente la tua eta'?\n");  
  printf("\n(immettere: p.e. 110010 per indicare riga1, riga2, riga5)? ");      
  gets(risposta);      
  for (i=0;i&lt;6;i++){ switch (i){ case 0: if (risposta[i]!='0') anni++; break; case 1: if (risposta[i]!='0') anni+=2; break; case 2: if (risposta[i]!='0') anni+=4; break; case 3: if (risposta[i]!='0') anni+=8; break; case 4: if (risposta[i]!='0') anni+=16; break; case 5: if (risposta[i]!='0') anni+=32; break; default: anni=-1; } } if (anni&gt;0) printf ("\n Tu hai %d anni!\n",anni);
  else printf("\n hai sbagliato ad immettere i dati!\n");
  return 0;
}

output del programma per indovinare l'età

Notate che ho aggiunto la riga 6 (che corrisponde a 25) per avere a disposizione un range di età fra 1 e 63 anni.

Nella figura a fianco, l’output del programma per la mia età…

Più che la reale comprensione del giochino, forse la difficoltà maggiore si incontra nel strutturare l’algoritmo per la costruzione della tabella, io l’ho risolto così: un paio di cicli annidati e una serie di costrutti condizionali. Ho usato anche un operatore bit a bit per l’incremento della variabile a, nei casi di moltiplicazione *2 (shift a sinistra) è estremamente veloce!
Mi aspetto che ci siano altri modi anche più compatti, a voi il divertimento! (fatemi sapere…)

Download (compilato per win32): indovina.exe

Conclusioni:

Vi consiglio di evitare di proporre il giochino a matematici o esperti informatici che scopriranno subito il meccanismo, a parte questo mi sembra un esercizio divertente per rispolverare qualche reminiscenza di C, un linguaggio oggi forse un po’ in disuso ma da sempre un caposaldo della programmazione.

Riferimenti ed approfondimenti:

Calcolo della distanza geodetica tra due punti della superficie terrestre

La TerraDiversi anni fa, quando lavoravo al CNR, insieme a Bruno Fornari, un caro collega ed amico purtroppo scomparso, scrivemmo un programmino di “servizio” che utilizzammo per calcolare matrici di distanze terrestri (geodetiche).
Il programma è veramente semplice, più interessante è capire come si calcola la distanza tra due punti di una superficie sferica. Devo riconoscere che, al momento, ho avuto la sensazione che i miei studi scientifici non fossero serviti a molto, perché non avevo la più pallida idea di come si facesse! (Gli studi dell’OCSE-Pisa dicono il vero…) Per fortuna, i consigli di Bruno, fisico e scienziato di prim’ordine, mi hanno fornito la formula giusta.
Per il nostro scopo, non ci serviva una grande precisione, per cui l’assunto iniziale era che la terra fosse una sfera perfetta (in realtà non lo è, vi rimando ai link sotto per gli approfondimenti).
Dunque, osservando la figura in alto, diciamo che, in base alla trigonometria sferica (teorema di Eulero), tra i lati a, b e p del triangolo sferico ABP vale la relazione:

cos p = cos a cos b + sen a sen b cos φ

Ora, dette lat(A), lon(A), lat(B), lon(B), la latitudine e la longitudine dei punti A e B e, considerato che:

  • a = 90° – lat(B)
  • b = 90° – lat(A)
  • φ = lon(A)lon(B)

abbiamo tutti i dati per calcolare la lunghezza del lato p considerando il raggio della Terra approssimabile a R = 6371 km.

Ecco dunque il codice in linguaggio ANSI C:

dis_geod.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include 
/* Questa funzione calcola la distanza tra due punti 
sulla superficie terrestre, date le coordinate in
latitudine e longitudine espresse in
gradi decimali */
double disgeod (double latA, double lonA,
                double latB, double lonB)
{
      /* Definisce le costanti e le variabili */
      const double R = 6371;
      const double pigreco = 3.1415927;
      double lat_alfa, lat_beta;
      double lon_alfa, lon_beta;
      double fi;
      double p, d;
      /* Converte i gradi in radianti */
      lat_alfa = pigreco * latA / 180;
      lat_beta = pigreco * latB / 180;
      lon_alfa = pigreco * lonA / 180;
      lon_beta = pigreco * lonB / 180;
      /* Calcola l'angolo compreso fi */
      fi = fabs(lon_alfa - lon_beta);
      /* Calcola il terzo lato del triangolo sferico */
	  p = acos(sin(lat_beta) * sin(lat_alfa) + 
        cos(lat_beta) * cos(lat_alfa) * cos(fi));
      /* Calcola la distanza sulla superficie 
      terrestre R = ~6371 km */
      d = p * R;
      return(d);
}

Ed ecco il codice per il programma di interfaccia utente:

prodisge.c

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include 
#include "dis_geod.h"
double disgeod(double latA, double lonA, 
              double latB, double lonB);
void main (void)
{
	float latA, lonA, latB, lonB, distanza;
	printf("\n Inserisci la latitudine del punto A :\t");
	scanf ("%f", &amp;latA);
	printf("\n Inserisci la longitudine del punto A :\t");
	scanf ("%f", &amp;lonA);
	printf("\n Inserisci la latitudine del punto B :\t");
	scanf ("%f", &amp;latB);
	printf("\n Inserisci la longitudine del punto B :\t");
	scanf ("%f", &amp;lonB);
	distanza =  disgeod(latA, lonA, latB, lonB);
	printf("\n La distanza fra A e B e' : %f  km\n", distanza);
}

Ed ecco infine il risultato del calcolo della distanza fra Roma (lat: +41.91;, lon: +12.45) e Milano (lat: +45.48 lon: +09.18)

output del programma disgeod

Se provate il calcolo con altri programmi, potrete avere risultati un poco differenti, questo è dovuto alle approssimazioni usate per il raggio della Terra, per il valore di Π e per la conversione dei gradi sessagesimali in decimali.
Ovviamente, al programma originale fu aggiunto il codice necessario per costruire le matrici di distanze, che qui, per brevità, ho tagliato.

Download: (compilato per win32) prodisge.exe

Conclusioni:

Ormai, che i GPS ci parlano amichevolmente in auto, e Google Earth ci da tutte le informazioni che vogliamo con pochi click, programmini del genere sembrano della preistoria. Ma, a furia di avere sempre la “pappa pronta”, non ci annichiliremo un po’ troppo le meningi? Trovo che ogni tanto sarebbe bene fare qualche ripasso, tanto per capire anche come funziona ciò che si utilizza passivamente!

Riferimenti ed approfondimenti: