/* sCrivellante - v0.5 		Ottimizzazione:  Eliminati messaggi di broadcast (Kpe)	*/
/*						 Eliminati numeri pari	(Kpe&Buc)	*/
/*						 Ottimizzato l'uso in cache (Kpe&Buc)	*/
/*											*/
/*											*/
/*		Per lanciare il programma:						*/
/*		$ mpdboot -n 5								*/
/*		$ mpicxx scrivellante.c -o scrivellante					*/
/*		$ mpiexec -n 5 ./scrivellante <n>					*/
/*											*/
/*	Il programma in questione calcola tutti i numeri primi da 0 ad n, parametro	*/
/*	che viene acquisito da linea di comando, attraverso l'applicazione del metodo	*/
/*	noto come "Crivello di Eratostene".						*/
/*	Per evitare lo scambio di messaggi di broadcast atti a comunicare di volta in 	*/
/*	volta il numero primo di cui marcare i multipli, vengono inizialmente su 	*/
/*	ogni processo calcolati i numeri primi da 3 a radice di n (utilizzando l'array 	*/
/*	priv8) che verranno poi utilizzati nel marcamento in marked.			*/
/*	Dall'esecuzione sono stati eliminati i numeri pari, ottenendo così un 		*/
/*	dimezzamento del dominio su cui effettuare le computazioni.			*/
/*	L'array viene scandito a blocchi di dimensione DIV, in modo tale da ottimizzare	*/
/*	l'uso della cache evitando frequenti accessi in RAM che lo scorrimento di un 	*/
/*	array molto vasto comporterebbe.						*/
	
#include <mpi.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

// Macro
#define BLOCK_LOW(id,p,n) ((id)*(n)/(p))
#define BLOCK_HIGH(id,p,n) (BLOCK_LOW((id)+1,p,n)-1)
#define BLOCK_SIZE(id,p,n) (BLOCK_LOW((id)+1,p,n-2)-BLOCK_LOW(id,p,n-2))

// Costanti
#define DIV 30000			// blocco di numeri da esaminare di volta in volta: dipende dalla cache del sistema

int main(int argc, char *argv[])
{
	int count = 0;		// Contatore dei numeri primi in marked
	double elapsed_time;
	int first;		// Indice del primo multiplo del numero primo
	int global_count;	// Contatore di tutti i numeri primi calcolati nei vari processi (servirà solo al processo 0)
	int high_value;		// Valore massimo processato dal processo corrente
	int i;
	int id;			// Rank del processo
	int index = 0;		// Indice che scorrerà marked
	int indice = 0;		// Indice che scorrerà priv8
	int low_value;		// Valore minimo processato nel processo corrente
	char *marked;		// Array di char che rappresenterà l'insieme dei numeri da cui ricavare i primi
	char *priv8;		// Array di char che servirà a calcolare tutti i numeri primi tra 3 e radice di n
	int n;			// Estremo superiore
	int p;			// Numero di processi impiegati nell'esecuzione del programma
	int proc0_size;		// Dimensione dell'array assegnato al processo 0
	int prime;		// Numero primo di volta in volta preso da priv8 per marcarne i multipli su marked
	int size;		// Dimensione dell'array assegnato al processo
	int primo;		// Numero primo di volta in volta preso da priv8 per marcarne i multipli

	
	
	MPI_Init(&argc, &argv);
	
	MPI_Barrier(MPI_COMM_WORLD);
	elapsed_time = -MPI_Wtime();
	
	MPI_Comm_rank(MPI_COMM_WORLD, &id);
	MPI_Comm_size(MPI_COMM_WORLD, &p);
	
	if(argc != 2)
		{
		if (!id) printf("Command line: %s <m>\n", argv[0]);
		MPI_Finalize();
		exit(1);
		}
	n = atoi(argv[1]);
	int size2 = (int) ((sqrt((double) n))/2)-2;	// Dimensione di priv8
	
	low_value = 3 + BLOCK_LOW(id,p,n-2);
	high_value = 3 + BLOCK_HIGH(id,p,n-2);
	if(low_value % 2 == 0)			// non esistono valori pari nei
		low_value +=1;			// nostri array, quindi neppure
	if(high_value % 2 == 0)			// i valori estremi potranno essere
		high_value -=1;			// tali, ergo li ritocchiamo
	
	size = (high_value-low_value)/2;

	proc0_size = (n - 1) / p;
	int low_index = ((low_value-1)/2)-1;	// Indice del primo elemento che dovrà elaborare il processo
	
	int low_subvalue = low_value;		// Indice del primo elemento del sottoinsieme di array assegnato al processo
	int high_subvalue = low_value + DIV;	// Indice dell'ultimo elemento del sottoinsieme di array assegnato al processo
//	int low_subindex = ((low_subvalue-1)/2)-1;
	if (high_subvalue>high_value) high_subvalue=high_value;

/* parte di controllo dei valori trovati
	printf("[%d] Low value: %d\n", id, low_value);
	printf("[%d] Low subvalue: %d High sub_value: %d\n",id,low_subvalue,high_subvalue);
	printf("[%d] High value: %d\n", id, high_value);
	printf("[%d] Size: %d\tSize2: %d\n", id, size, size2);
	printf("[%d] Proc0 size: %d\n", id, proc0_size);
	printf("[%d] Low value = %d , low index = %d\n", id, low_value, low_index);
*/
	
	if( (3 + proc0_size) < (int) sqrt((double) n) )
		{
		if(!id) printf("Too many processes\n");
		MPI_Finalize();
		exit(1);
		}
	
	marked = (char *) malloc(size+1);
	priv8 = (char *) malloc(size2+1);
	
	if(marked == NULL || priv8 == NULL)
		{
		printf("Cannot allocate enough memory\n");
		MPI_Finalize();
		marked[i] = 0;
		}
		
	for(i=0; i<=size; i++)	
		marked[i] = 0;
		
/*	for(i=0; i<size; i++) 
		printf("[%d] %d\n",id,(i+low_index)*2+3); */
		
	for(i=0; i<=size2; i++)
		priv8[i] = 0;
	
	// Algoritmo per la ricerca sequenziale di tutti i numeri primi compresi tra 3 e radice di n
	do {
	   primo=((indice)*2)+3;
	   for (i=((primo*primo)/2)-1;i<=size2; i+=primo)
	      priv8[i]=1;
	   while(priv8[++indice]);
	} while(indice <= size2);
	
	int high_subindex = ((high_subvalue-1)/2)-1-low_index;	// Indice dell'ultimo elemento del sottoinsieme di array assegnato al processo
	
while(low_subvalue < high_value)
	{
	index = 0;
	prime = 3;
//	printf("[%d] Low subvalue: %d High sub_value: %d\n",id,low_subvalue,high_subvalue);
	while(prime * prime <= high_subvalue)
	    {
	    if(prime * prime > low_subvalue)
	    	 first = ((prime * prime) / 2 - 1) - low_index;		// nota: uguale a ((prime*prime-3) / 2) - low_index
	    else {
	    	  if((low_subvalue % prime) == 0)
		  	first = ((low_subvalue/2)-1) - low_index;
		  else  {
		  	first = (low_subvalue + prime - (low_subvalue % prime));
			if(first % 2 == 0)
				first+=prime;
			first = ((first/2)-1) - low_index;
			}
		 }
//	    printf("[%d] Prime: %d First: %d High_subvalue/2: %d\n",id,prime,(first+low_index)*2+3,high_subvalue/2);
//	    printf("[%d] Prime: %d First: %d High_subindex: %d\n",id,prime,first,((high_subvalue-1)/2)-1);
	    for(i=first; i<=high_subindex; i+=prime)	// forse high_balue/2-1
	    	{
	    	 marked[i] = 1;
//		 printf("[%d] Markato: %d\n",id,(i+low_index)*2+3);
		} 
    	    while(priv8[++index]);
	    prime = (index+1)*2+1;
	    }
	  low_subvalue += DIV;
	  high_subvalue += DIV;
	 if (high_subvalue > high_value)
	 	high_subvalue = high_value;
//	  low_subindex = ((low_subvalue-1)/2)-1;
	high_subindex = ((high_subvalue-1)/2)-1-low_index;
	 }
	 
	 
	   
	if(id==0) count = 1;		// Poiché anche 2, non incluso nell'array,è un numero primo
	
	// Conteggio dei numeri primi trovati localmente
	for(i=0; i<=size; i++)
		if(! marked[i])
			{
			count++;
//			printf("[%d] %d° num primo: Indice %d = %d\n", id,count, i, (i+low_index)*2+3);
			}	// decommentare la riga precedente per visualizzare i numeri primi trovati
//	printf("[%d] Numeri primi: %d\n",id,count);
			
	// Conteggio totale dei numeri primi trovati da ogni processo eseguito dal processo 0
	MPI_Reduce(&count, &global_count, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
		
	elapsed_time += MPI_Wtime();
	
	if(!id) {
		printf("%d primes are less than or equal to %d\n", global_count, n);
	 	printf("Total elapsed time: %10.6f\n", elapsed_time);
		}
	MPI_Finalize();
	return 0;
}