#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <mpi.h>

/// Math part
int fact( int x );
int CNM( int n, int m );
double** MatrixCreate( int dim );
void MatrixDelete( double** m, int dim);
void MatrixCopy( double** dest, double** src, int dim);
void MatrixMul( double** m1, double** m2, int dim, double** res );
void MatrixPow( double** m1, int dim, int power, double** res );
/// Calculation part
int GetK( double** matr, int dim, double delta );
void Getf( double** matr, int dim, int K, double *f, int numproc, int myid );
double GetM( double *f, int count );
double GetD( double *f, int count );
double* ParallelMfromN( int procCount, int M, double *f, int K );
double* ParallelAND( int procCount, double *f, int K );
double* ParallelOR( int procCount, double *f, int K );
double* Consecutive( int procCount, double **f, int K );
/// Other stuff
void InitMyMatrix( double** matrix, int dim );
/// Entry point
int main( int argc, char **argv, char **envp )
{
	// Всякий разный MPI шит
	MPI_Init( &argc, &argv );
	int numprocs,myid;
	MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
	MPI_Comm_rank(MPI_COMM_WORLD,&myid);
	if( myid == 0 ){
		printf( "proc count = %i\nID = %i\n", numprocs, myid );
		fflush( stdout );
	} else {
		printf( "Start Child proc id = %i\n", myid);
	}

	double start, end;
	start= MPI_Wtime( );



	int dim = 8, I = 5, M = 2, L = 7;
	double delta = 0.000001;
	double **mymatrix = MatrixCreate( dim );
	InitMyMatrix( mymatrix, dim );
	int K = GetK( mymatrix, dim, delta );
	double *f = new double[K];

	// Эта функция работает через библиотеку MPI
	Getf( mymatrix, dim, K, f, numprocs, myid);

	if (myid == 0) {	// Поэтому все последующие вычисления буду проходить в главном потоке
						// он-же поток с id == 0

		double *fMN = ParallelMfromN( I, M, f, K );
		double *fAND = ParallelAND( L, f, K );
		double **f3 = new double*[3];
		f3[0] = new double[K];
		f3[1] = new double[K];
		f3[2] = new double[K];
		for( int i = 0; i < K; i++ )
		{
			f3[2][i] = fMN[i];
			f3[0][i] = f[i];
			f3[1][i] = fAND[i];
		}
		double *fRes = Consecutive( 3, f3, K );
		double getM = GetM( fRes, 3 * K );
		double getD = GetD( fRes, 3 * K );
		// The end

		// Cчитаем время выполнения и финализируем MPI
		end = MPI_Wtime( );
		MPI_Finalize( );
		printf("Time = %f\n\n", end-start );


		double sum = 0;				   //some stuff
		double *fsum = fMN;
		for( int i = 0; i < K; i++ )
		{
			sum += fsum[i];
			printf( "f[%i]\t=\t%e\n", i + 1, fsum[i] );
		}
		printf( "\n=====[%f]=====\n\n", sum );
		sum = 0;
		for( int i = 0; i < 3 * K; i++ )
			sum += fRes[i];
		printf( "=====\nsum = %f\n", sum );
		printf( "=====\nM = %f\nD = %f\n", getM, getD );
		fflush( stdout );
		MatrixDelete( mymatrix, dim );
		delete f;
		delete fMN;
		delete fAND;
		delete f3[0];
		delete f3[1];
		delete f3[2];
		delete[] f3;
		delete fRes;
		return 0;
	}
	// Это выполняют все остальные потоки кроме главного
	MPI_Finalize();
	delete f;
	return 0;
}
/// Math part
int fact( int x )
{
	int res = 1;
	for( int i = 2; i <= x; i++ )
		res *= i;
	return res;
}
int CNM( int n, int m )
{ return fact( n ) / (fact( m ) * fact( n - m )); }
double** MatrixCreate( int dim )
{
	double** matrix = new double*[dim];
	for( int i = 0; i < dim; i++ )
		matrix[i] = new double[dim];
	return matrix;
}
void MatrixDelete( double** m, int dim)
{
	for( int i = 0; i < dim; i++ )delete m[i];
	delete[] m;
}
void MatrixCopy( double** dest, double** src, int dim)
{
	for( int i = 0; i < dim; i++ )
		for( int j = 0; j < dim; j++ )
			dest[i][j] = src[i][j];
}
void MatrixMul( double** m1, double** m2, int dim, double** res )
{
	double** tmpBuf = MatrixCreate( dim );
	double sum;
	for( int i = 0; i < dim; i++ )
		for( int j = 0; j < dim; j++ )
		{
			sum = 0;
			for( int k = 0; k < dim; k++ )
				sum += m1[i][k] * m2[k][j];
			tmpBuf[i][j] = sum;
		}
	MatrixCopy( res, tmpBuf, dim );
	MatrixDelete( tmpBuf, dim );
}

void MatrixPow( double** m1, int dim, int power, double** res )
{
	if( power == 0 )
	{
		for( int i = 0; i < dim; i++ )
			for( int j = 0;  j < dim; j++ )
				res[i][j] = (i == j ? 1 : 0);
		return;
	}
	double** m2 = m1;
	if( m1 == res )
	{
		m2 = MatrixCreate( dim );
		MatrixCopy( m2, m1, dim );
	}
	MatrixCopy( res, m1, dim );
	for( int i = 1; i < power; i++ )
		MatrixMul( res, m2, dim, res );
}
/// Calculation part
int GetK( double** matr, int dim, double delta )
{
	double** tmpMatr = MatrixCreate( dim );
	int k = 0;
	do{
		k++;
		MatrixPow( matr, dim, k, tmpMatr );
	}while( (1 - tmpMatr[0][dim - 1]) >= delta );
	MatrixDelete( tmpMatr, dim );
	return k;
}
/**
##########################################################################################
		РЕАЛИЗОВАНА С ПОМОЩЬЮ БИБЛИОТЕКИ MPI.
##########################################################################################
**/
void Getf( double** matr, int dim, int K, double *f , int numprocs, int myid)
{
	double sum;

	double** m1 = MatrixCreate( dim );
	double** m2 = MatrixCreate( dim );
	sum = 0;
	double* fm1 = new double[K];
	double* fm2 = new double[K];

	/**
	##########################################################################################
	ЦИКЛ ВИДА:

		for (int cnt = 0; cnt < MAX_CNT; cnt++)

	ПЕРЕД ВОЗВЕДЕНИЕМ В СТЕПЕНЬ ДЛЯ ТОГО, ЧТОБЫ БОЛЕЕ ЧЕТКО ВИДЕТЬ
	РАЗЛИЧИЯ В ПРОИЗВОДИТЕЛЬНОСТИ В MPI РЕЖИМЕ
	##########################################################################################
	**/

	int MAX_CNT = 100;

	for( int i = 1; i <= K; i++ )
	{
		if (numprocs == 1) {	// Все выполняется в одном процессе
			for (int cnt = 0; cnt < MAX_CNT; cnt++)MatrixPow( matr, dim, i, m1 );
			for (int cnt = 0; cnt < MAX_CNT; cnt++)MatrixPow( matr, dim, i - 1, m2 );

			fm1[i-1] = m1[0][dim - 1];
			fm2[i-1] = m2[0][dim - 1];

			//if (i == 1) printf("ONE THREAD MODE\n");

		} else {				// Процессов несколько
			if (myid == 0) {	// Первую матрицу возводим в главном процессе
				for (int cnt = 0; cnt < MAX_CNT; cnt++)MatrixPow( matr, dim, i, m1 );
				fm1[i-1] = m1[0][dim - 1];

				//if (i == 1)printf("MPI MODE; ID = %i; START = %f\n", myid, MPI_Wtime( ));
				//if (i == K)printf("MPI MODE; ID = %i; END = %f\n", myid, MPI_Wtime( ));
			}
			if (myid == 1) { // Вторую матрицу в отдельном специальном процессе
				for (int cnt = 0; cnt < MAX_CNT; cnt++)MatrixPow( matr, dim, i - 1, m2 );
				fm2[i-1] = m2[0][dim - 1];

				//if (i == 1)printf("MPI MODE; ID = %i; START = %f\n", myid, MPI_Wtime( ));
				//if (i == K)printf("MPI MODE; ID = %i; END = %f\n", myid, MPI_Wtime( ));
			}
		}
	}

	if (numprocs != 1) {
		double start, end;
		start= MPI_Wtime( );
		// И синхронизируемся с главным процессом
		MPI_Bcast(fm2, K , MPI_DOUBLE, 1, MPI_COMM_WORLD);
		end = MPI_Wtime( );
		printf("ID=%i\tTime = %f\n\n", myid, end-start );
	}

	// Все вспомогательные вычисления только в главном потоке
	if (myid == 0) {

		printf("STUFF TIME = %f\n", MPI_Wtime( ));

		for( int i = 1; i <= K; i++ ) {
			if (myid == 0) {	// Только в главном процессе
				f[i - 1] = fm1[i-1] - fm2[i-1]; //m1[0][dim - 1] - m2[0][dim - 1];
				sum += f[i - 1];
			}
		}
	}

	MatrixDelete( m1, dim );
	MatrixDelete( m2, dim );
}

double GetM( double *f, int count )
{
	double res = 0;
	for( int i = 0; i < count; i++ )
		res += i * f[i];
	return res;
}
double GetD( double *f, int count )
{
	double M = GetM( f, count );
	double res = 0;
	for( int i = 0; i < count; i++ )
		res += pow( i - M, 2 ) * f[i];
	return res;
}
double* ParallelMfromN( int procCount, int M, double *f, int K )
{
	double *fAND = ParallelAND( M, f, K );
	double *fMN = ParallelOR( CNM( procCount, M ), fAND, K );
	return fMN;
}
double* ParallelAND( int procCount, double *f, int K )
{
	double *fAnd = new double[K];
	double S, M, res;
	for( int k = 0; k < K; k++ )
	{
		S = 0;
		for( int n = 1; n <= procCount; n++ )
		{
			M = 1;
			for( int i = 1; i <= procCount; i++ )
			{
				if( n == i )
					res = 1;
				else
				{
					res = 0;
					for( int ki = 0; ki < K; ki++ )
						if( (i < n && ki < k) || (i >= n && ki <= k) )
							res += f[ki];
				}
				M *= res;
			}
			S += f[k] * M;
		}
		fAnd[k] = S;
	}
	return fAnd;
}
double* ParallelOR( int procCount, double *f, int K )
{
	double *fOr = new double[K];
	double S, M, res;
	for( int k = 0; k < K; k++ )
	{
		S = 0;
		for( int n = 1; n <= procCount; n++ )
		{
			M = 1;
			for( int i = 1; i <= procCount; i++ )
			{
				if( n == i )
					res = 1;
				else
				{
					res = 0;
					for( int ki = 0; ki < K; ki++ )
						if( (i < n && ki > k) || (i >= n && ki >= k) )
							res += f[ki];
				}
				M *= res;
			}
			S += f[k] * M;
		}
		fOr[k] = S;
	}
	return fOr;
}
double* Consecutive( int procCount, double **f, int k )
{
	double S;
	double *f2 = new double[k * procCount];
	double *f1 = new double[k * procCount];
	for( int i = 0; i < k * procCount; i++ )
		f1[i] = f2[i] = 0;
	for( int i = 0; i < k; i++ )
		f1[i] = f[0][i];
	for( int j = 2; j <= procCount; j++ ){
		for( int l = 1;l <= k * j; l++ ){
			S =	0;
			for( int m = 1; m <= k * procCount; m++ )
				if( m < l && m <= k * (j - 1) && l - m <= k )
					S += f1[m - 1] * f[j - 1][l - m - 1];
			f2[l - 1] = S;
		}
		for( int i = 0; i < k * j; i++ )
			f1[i] = f2[i];
	}
	delete f2;
	return f1;
}
/// Other stuff
void InitMyMatrix( double** m, int dim )
{
	if( dim < 8 ) return;
	for( int i = 0; i < dim; i++ )
		for( int j = 0; j < dim; j++ )
			m[i][j] = 0;
	m[0][1] = m[2][3] = m[3][4] = m[4][5] = m[7][7] = 1;
	m[1][2] = m[1][3] = m[5][4] = m[5][6] = 0.5;
	m[6][4] = m[6][7] = 0.5;
  //0{0.0,   1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  //1{0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0},
  //2{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
  //3{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
  //4{0.0, 0.0, 0.0, 0.0, 0.0,   1, 0.0, 0.0},
  //5{0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0},
  //6{0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5},
  //7{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}
}
/// EOF