/**  Antonio Bucciol    [ webmaster@freshbuc.com ]
/*        per la SCP    [ http://www.scpvolterra.it ]
/*  presenta DICOT, programma che effettua una stima dell'area tra la funzione 
/*  sin(x)/x e l'asse X in un range compreso tra LOW_LIMIT e HIGH_LIMIT
/*
/*  Compilazione:
/*     mpicxx senx_su_x.cpp -o senx_su_x --fast-math -O3
/*  Esecuzione: mpiexec -n <numero_processi> ./senx_su_x <numero_suddivisioni>
/*  All'aumentare del numero di suddivisioni, aumenta la precisione del calcolo
/****************************************************************************/
// inclusione librerie:
#include <mpi.h>	// libreria MPI
#include <limits.h>
#include <math.h>	// libreria di funzioni matematiche
#include <stdio.h>	// libreria standard I/O
#include <stdlib.h>	// libreria standard del C++

/* Modificare le seguenti costanti per cambiare l'intervallo
   dell'asse X da considerare per il calcolo dell'area  */
#define LOW_LIMIT -50
#define HIGH_LIMIT 50

int main(int argc, char *argv[])
{
		// variabili per le lib MPI
	double elapsed_time;	// tempo effettivo di esecuzione
	int p;		// numero di processi allocati
	int id;		// id del processo (da 0 a p-1)

		// variabili specifiche dell'algoritmo
	int i;		// contatore
	double x;	// conterrà i valori di x per cui andrò a calcolare la funzione sin(x)/x
	double h;	// altezza della funzione sin(x)/x per i vari valori di x
	double area_locale = 0;	 // area calcolata dal singolo processo
	double area_totale = 0;  // area totale ricavata dalla somma di tutte le `area_locale'
	
	MPI_Init(&argc, &argv);  // inizializzazione delle lib MPI
	
	MPI_Barrier(MPI_COMM_WORLD);   // aspetta che tutti i processi arrivino a questo punto
	elapsed_time = -MPI_Wtime();   // inizializzazione variabile tempo
	
	MPI_Comm_rank(MPI_COMM_WORLD, &id);	// vedere dispense MPI di Sartorello
	MPI_Comm_size(MPI_COMM_WORLD, &p);	// per spiegazione dettagliata

	if(argc != 2)
		{
		if (!id) printf("Command line: %s <n>\n", argv[0]);
		MPI_Finalize();
		exit(1);
		}

	unsigned long int n = strtoul(argv[1], NULL, 10);	// conversione in long int del numero `n'
	
	if(!id)	{
		printf("Dividing total range by %lu\n", n);
		printf("Max value usable on this box: %lu\n", ULONG_MAX);
		}

	
	int width = HIGH_LIMIT - LOW_LIMIT;	// larghezza intervallo totale
	double local_width = width / p;		// larghezza intervallo locale


	if(p>width)	// il numero di processi non deve superare la larghezza del range
		{
		if (!id) printf("Too many processes!\n");
		MPI_Finalize();
		exit(1);		
		}
	if(local_width == 0)
	    {
	    if (!id) printf("Machine 0 reached! Please decrease partitions number\n");
		MPI_Finalize();
		exit(1);		
		}
	
	double low_value = local_width * id + 1 + LOW_LIMIT;	// valore minimo del range locale
	double high_value = local_width * (id+1) + LOW_LIMIT;	// valore massimo del range locale

	if(id == 0)			// se è il primo processo
		low_value = LOW_LIMIT;		// allora low_value è uguale a LOW_LIMIT
	
	if(id + 1 == p)			// se è l'ultimo processo
		high_value = HIGH_LIMIT;	// allora high_value è uguale a HIGH_LIMIT
	
	
	double margin = (double) width / n;	// distanza tra ogni altezza == base di ogni rettangolo
	
	if(!id) printf("Margin: %.20e\n", margin);
	
	if(id != 0)
		low_value -= 1 - margin;	/* low_value = low_value - 1, ovvero high_value del
					processo precedente, sommato `margin', ovvero il primo
					punto di quel subrange in cui vado a calcolare la funzione */
			
	x = low_value;		// x parte da low_value
	
	//printf("   [%d] Low: %f    High: %f   (%f)\n", id, low_value, high_value, (high_value-low_value));

	if(!id) printf("Domain suddivision completed\n\n");
	
	for(i=1; i<=(local_width/margin); i++)	// per il numero di `margin' in cui è suddivisibile `local_width':
		{
		x += margin;	// aumento `x' di `margin'
		if(x==0)	// se `x' si trova sull'origine, punto per cui la funzione sin(x) non esiste,
		     continue;	// allora salto al prossimo ciclo senza effettuare i calcoli.
		
		h = sin(x) / x;			// salvo in `h' il risultato dell'operazione sin(x)/x
		area_locale += margin * h;	// sommo all'area locale quella del rettangolo con l'altezza attuale
		}
	
		// Sommo `area_locale' di ogni programma, memorizzando il dato in `area_totale'
	MPI_Reduce(&area_locale, &area_totale, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
	
	elapsed_time += MPI_Wtime();	// stoppo il tempo di esecuzione

	//printf("[%d] Area locale: %.20e\n", id, area_locale);

	if(!id) {
		printf("\tArea totale: %.30e\n", area_totale);
	 	printf("\tTotal elapsed time: %10.6f\n", elapsed_time);
		}
	MPI_Finalize();		// terminazione delle funzioni della lib MPI
	return 0;
}