Eratostene

Eratostene e la misura della circonferenza della terra

Eratostene (276 a.C. – 195 a.C., circa) valutò in modo approssimato la lunghezza L della circonferenza terrestre conoscendo la distanza tra due località A (Siene) e B (Alessandria d'Egitto) che stanno più o meno sullo stesso meridiano (che è lungo complessivamente L) e misurando l'inclinazione dei raggi del sole in B nello stesso giorno in cui in A erano verticali.
Dalle misure dell'ampiezza dell'angolo α e dell'arco s (vedi la figura a lato) ha ottenuto una stima di L: L = s·360°/α.
Per curiosità, pare che Eratostene ottenne la lunghezza di circa 39·10³ km, mentre in realtà essa è di circa 40.1·10³ km (questa è una stima del valore trovato da Eratostene, in quanto l'unità di misura delle lunghezze era lo "stadio", che non si sa precisamente a quale valore in metri corrispondesse).

Il crivello di Eratostene

Ad Eratostene è dovuto anche l'algortimo per calcolare i numeri primi noto come crivello di Eratostene ("crivello" è un sininomo di "setaccio"), che trovi descritto dettagliatamente in WolframAlpha: sieve of Eratosthenes (vedi nota). Con WolframAlpha puoi anche calcolarli: ad esempio con primes <= 100 trovi i numeri primi minori o eguali a 100.
Vediamo come realizzare il "crivello" con R.
L'idea è semplice: cancellare i multipli di 2, poi i multipli di 3, poi quelli di 5 (quelli di 4 li ho già cancellati con i multipli di 2), …

n <- 100
x <- NULL; for (i in 1:n) x[i] <- i
k<- 2; h<- k+k; while(h<= n) {x[h]<- 0; h<- h+k}; z<- NULL; for(i in 1:n) z<- c(z,x[i]); z
  [1]  1  2  3  0  5  0  7  0  9  0 11  0 13  0 15  0 17  0 19  0 21  0 23  0 25  0
 [27] 27  0 29  0 31  0 33  0 35  0 37  0 39  0 41  0 43  0 45  0 47  0 49  0 51  0
 [53] 53  0 55  0 57  0 59  0 61  0 63  0 65  0 67  0 69  0 71  0 73  0 75  0 77  0
 [79] 79  0 81  0 83  0 85  0 87  0 89  0 91  0 93  0 95  0 97  0 99  0
k<- 3; h<- k+k; while(h<= n) {x[h]<- 0; h<- h+k}; z<- NULL; for(i in 1:n) z<- c(z,x[i]); z
  [1]  1  2  3  0  5  0  7  0  0  0 11  0 13  0  0  0 17  0 19  0  0  0 23  0 25  0
 [27]  0  0 29  0 31  0  0  0 35  0 37  0  0  0 41  0 43  0  0  0 47  0 49  0  0  0
 [53] 53  0 55  0  0  0 59  0 61  0  0  0 65  0 67  0  0  0 71  0 73  0  0  0 77  0
 [79] 79  0  0  0 83  0 85  0  0  0 89  0 91  0  0  0 95  0 97  0  0  0
k<- 5; h<- k+k; while(h<= n) {x[h]<- 0; h<- h+k}; z<- NULL; for(i in 1:n) z<- c(z,x[i]); z
  [1]  1  2  3  0  5  0  7  0  0  0 11  0 13  0  0  0 17  0 19  0  0  0 23  0  0  0
 [27]  0  0 29  0 31  0  0  0  0  0 37  0  0  0 41  0 43  0  0  0 47  0 49  0  0  0
 [53] 53  0  0  0  0  0 59  0 61  0  0  0  0  0 67  0  0  0 71  0 73  0  0  0 77  0
 [79] 79  0  0  0 83  0  0  0  0  0 89  0 91  0  0  0  0  0 97  0  0  0
 ...

Ecco come automatizzare meglio il procedimento:

n <- 1000
x <- NULL; for (i in 1:n) x[i] <- i
k <- 2
while(k <=n) { h <- k+k ; while(h <= n) {x[h]<- 0; h<- h+k}; k <- k+1}
# stampo solo i valori diversi da 0:
z <- NULL; for(i in 1:n) if(x[i] > 0) z <- c(z,x[i]); z
  [1]   1   2   3   5   7  11  13  17  19  23  29  31  37  41  43  47  53  59  61
 [20]  67  71  73  79  83  89  97 101 103 107 109 113 127 131 137 139 149 151 157
 [39] 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263
 [58] 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379
 [77] 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491
 [96] 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617
[115] 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743
[134] 751 757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 863 877
[153] 881 883 887 907 911 919 929 937 941 947 953 967 971 977 983 991 997

Nota. Nel sito si osserva che per trovare i numeri primi minori o eguali ad N posso fermarmi a "setacciare" i numeri fino alla parte intera della radice quadrata di N. La cosa può essere dimostrata. Non è difficile farlo, ma è un dettaglio che qui non ci interessa affrontare.