/*
Integracao pelo metodo de Montecarlo

Funcao:  f(x)=-(x^4/2)+2*x^3-x^2+x+5

Objetivo: Fazer a integração definida entre 0 e 3 que é a própria area entre a curva e o eixo x da curva. O resultado deve ser 26,7.

Metodo: Inserir randomicamente pontos entre os pontos (0,0) e (3,14), sendo (x,y). A area da curva é proporcionar a quantidade de pontos sob a curva. Portanto AreaCurva=42*(PtIN/PtTot).
*/


#include <stdio.h>
#include <time.h>
#include <sys/timeb.h>
#include <stdlib.h>
#include <math.h>
#include <mpi/mpi.h>

#define MAX 10000

#define XMIN 0
#define XMAX 3
#define YMIN 0
#define YMAX 14

#define PART 1000

main(int argc, char *argv[])
{
	/*x e y: pontos randômicos / area:area final / areap:area parcial / result: resultado da função / 
	parcial: parte dada a cada processador / parte: subdivião da parte dada a cada procesador
	xini e xfim:valores iniciais e finais de cada parte / xmin e xmax: valores iniciais e finais da parte dada a cada processador*/
		double  x, y, area, areap, result, parcial, parte, *p_area, *recebe, xini, xfim, xmin, xmax;
	//i:conta num pontos por parte / j:conta partes / pttin:pontos dentro / pttot:pontos totais
		long int i, j, p, ptin=0, pttot=0;
	//nodolocal: nodo local / nnodos:numero de nodos / NPROC:numero de procesadores / nodo: nodos utilizados para calcular partes da curva
		int nodolocal, nnodos, NPROC, nodo[]={0,1,2,3,4};
	//usado para iniciar a variável randômica
	struct timeb tempoatual;
				
		NPROC=sizeof(nodo)/sizeof(int);
		
		//inicializacao da variavel randomica
		ftime(&tempoatual);
		srand(tempoatual.millitm);
		
		//inicializacao da biblioteca mpi
		MPI_Init(&argc,&argv);
		MPI_Comm_size(MPI_COMM_WORLD, &nnodos);
		MPI_Comm_rank(MPI_COMM_WORLD, &nodolocal);
		
		parte=(double)(XMAX-XMIN)/(double)NPROC;              //parte da curva dada a cada processador(range do x)
		parcial=parte/(double)PART;                           //subdivisao da parte que eh dada a cada processador
		xmin=XMIN;                                            //x minimo e maximo de cada trexo dado a cada processador
		xmax=parte+xmin;
		for (p=0;p<NPROC;p++)                                 //varre os nodos distribuindo as tarefas(partes da curva)
		{
				xini=xmin;
				xfim=xmin+parcial;
				
				if(nodolocal==nodo[p])                              //bloco paralelizado e executado por cada processador
				{
						areap=0;
						p_area=(double *)malloc(sizeof(double)*PART);
						if (p_area==NULL)
						{
								printf("ERRO 1!");
								exit(1);
						}
						
						for (j=0;j<PART;j++)                              //particionamento da parcela da curva dada a cada processador
						{
								for(i=0;i<MAX;i++)                              //numero maximo de pontos em cada parte
								{
										//calculo do x e y randomico
										x=((double)random())/RAND_MAX;
										x=(x*parcial)+xini;
										y=((double)rand())/RAND_MAX;
										y=y*(YMAX-YMIN)+YMIN;
										
										//calculo do valor da funcao que está sendo integrada
										result=-(pow(x,4.0)/2)+2.0*pow(x,3.0)-pow(x,2.0)+x+5.0;
										
										if (result>YMAX)                              //verifica se o valor y da curva esta fora do valor do YMAX estipulado
										{
												printf("Curva fora da area estipulada por YMAX (x=%f/y=%f)! O valor YMAX deve ser aumentado.\n",x,result);
												exit(1);
										}
										
										if (y<=result)                                //verifica se o ponto randomico esta abaixo ou acima da curva, se esta abaixo contabiliza como ponto 'in'
												ptin++;
								}
								pttot=i;
								//calcula a area parcial em funcao da relacao entre os pontos que ficaram abaixo da curva e os pontos totais
								p_area[j]=(YMAX-YMIN)*parcial*(double)ptin/(double)pttot;     
								
								//pega a proxima parte da curva para dar ao outro nodo processar
								xini=xfim;
								xfim=xini+parcial;
								pttot=0;
								ptin=0;
						}
						
						for (j=0;j<PART;j++)
								areap+=p_area[j];
						
		//printf("Nodo %d: Integral Definida entre %.2f e %.2f = %f (resultado parcial)\n",nodolocal, (double)xmin, (double)xmax, areap);
						
						free(p_area);
				}
				
				xmin=xmax;
				xmax=xmin+parte;
		}

	//alocação da variável que vai receber os dados dos nodos escravos
		recebe=(double *)malloc(sizeof(double)*NPROC);
		if (recebe==NULL)
		{
				printf("ERRO 1!");
				exit(1);
		}

		/*receber os dados calculados nos nodos para serem somados no nodo master       
		areap: dado que cada nodo escravo envia para o nodo master 
		recebe: variável que vai receber os dados, cada posicao(indice) da variável recebe é o número do nodo que gerou o dado*/
		int erro=MPI_Gather( &areap, 1, MPI_DOUBLE, recebe, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
		
		if (erro!=0)
				printf("Nodo %d: erro=%d\n", nodolocal, erro);
		
		if(nodolocal == 0)//bloco executado pelo nodo master
		{
				area=0;
				for (p=0;p<NPROC;p++)
						area+=recebe[p];   //soma todos os valores parciais
				
				printf("\n\tNodo %d: Integral Definida entre %.2f e %.2f = %f\n\n",nodolocal, (double)XMIN, (double)XMAX, area);
		}
		
		free(recebe);
		
		MPI_Finalize();
}