/*********************************************************
*   CRIVELLO DI ERATOSTENE
*   input: L - limite entro il quale ricercare i primi.
*   output: vettore contenente tutti i primi tra 2 e L.
*     
*   Languasco & Zaccagnini, Introduzione alla Crittografia, Hoepli.
*   Capitolo 6.3 pag. 165-166  
*********************************************************/  
 
{Erat(L) = local(a, B, B1, i , j , k , l, primi, primiout);

/* inizializziamo un vettore di L-1 componenti
 		 * (1 non e' primo) con gli interi tra 2 e L
 		 * e definiamo un altro vettore della stessa
 		 * dimensione in cui memorizzeremo i primi
 		 */
 		 
 		B=L-1;
 		a=vector(B);
 		primi=vector(B);
		for(i=1, B, a[i] = i+1);
	 
		/* intendiamo come marcato un intero la cui
 		 * corrispondente componente di a e' zero.
 		 * La componente viene marcata a zero se l'intero
 		 * contenuto e' un multiplo di uno degli interi 
 		 * precedenti
 		 */

        i = 1; l=1; B1=floor(sqrt(L));
        while(i <= B1,
        
                /* "saltiamo" gli interi gia' marcati */
                 
                while( a[i] == 0, i = i + 1);
                
                
                /* alla fine del while, l'intero contenuto
                * nella posizione a[i] deve essere un primo
                * perche' non e' diviso da nessun intero piu'
                * piccolo. Allora lo memorizzo nella prima
                * posizione disponibile del vettore primi
                */
                
                primi[l]=a[i]; 
                l=l+1;
              
                /* marchiamo adesso i multipli dell'intero
                * contenuto in a[i]=i+1, ossia azzeriamo le 
                * posizioni corrispondenti agli interi 
                * j(i+1), per j che varia fino a L/i+1.
                * queste posizioni hanno indice j(i+1)-1
                */
				              
                for(j=2, floor(L/(i+1)), a[(i+1)*j-1] = 0);
                
                /* passiamo ora a studiare la primalita'
                * dell'intero successivo incrementando la 
                * variabile che governa il while piu' esterno
                */
                
                i=i+1; 
        ); 
         
        /* inseriamo ora nel vettore dei primi i 
        * primi piu' grandi di radice di L
        */  
         
        for (j=i+1,B, 
                     if ((a[j]<>0), primi[l]=a[j];l=l+1));
        
        /* Sappiamo anche  quanti sono i primi minori o uguali
        * a L (sono l-1) e prepariamo l'output del vettore dei primi
        */
        
        print("Il numero di primi minori o uguali a ", L, " e': ", l-1);
        
        primiout=vector(l-1);
        for(j=1,l-1,primiout[j]=primi[j]);
       
        /* restituiamo in output il vettore dei primi
        */
       
        print("I numeri primi minori o uguali a ", L, " sono dati dal vettore: "); 
        return(primiout);
}