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

#include <mpi.h>

#define DIM 8

/*****************************************************************************/
/*
   	Вычисление факториала.
    x - .
    функция возвращает факториал числа x.	(x!)
*/
int fact( int x )
{
	int i, res = 1;
	for( i = 2; i <= x; i++ )
		res *= i;
	return res;
}
/*
   	Вычисление числа сочетаний C sup n sub m.
    n - .
	m - .
    функция возвращает число сочетаний из n по m.
*/
int CNM( int n, int m )
{ return fact( n ) / (fact( m ) * fact( n - m )); }
/*****************************************************************************/
/*
   	Печать матрицы на экран.
    M - исходная матрица.
    функция не имеет возвращаемого значения.
*/
void MPrint( double M[DIM][DIM] )
{
	int i, j;
	printf( "  Matrix\n    " );
	for( i = 0; i < DIM; i++ )
	{
		printf( "    " );
		for( j = 0; j < DIM; j++ )
			printf( "%6.3f ", M[i][j] );
		printf( "\n" );
	}
}
/*
    Копирование матрицы Src в матрицу Dest.				(Dest = Src)
    функция не имеет возвращаемого значения. Результат возвращается
      косвенно через аргумент Dest.
*/
void MCopy( double Dest[DIM][DIM], double Src[DIM][DIM] )
{
	int i, j;
	for( i = 0; i < DIM; i++ )
		for( j = 0; j < DIM; j++ )
			Dest[i][j] = Src[i][j];
}
/*
    произведение квадратных матриц размерности DIM.		(DIM = 8)
    A, B - матрицы, произведение которых вычисляется.	(A * B)
    C - матрица, в которую записывается результат.		(C = A * B)
    функция не имеет возвращаемого значения. Результат возвращается
      косвенно через аргумент C.
*/
void MMul( double A[DIM][DIM], double B[DIM][DIM], double C[DIM][DIM] )
{
	double tmpMatr[DIM][DIM];
	double sum;
	int i, j, k;
	for( i = 0; i < DIM; i++ )
		for( j = 0; j < DIM; j++ )
		{
			sum = 0;
			for( k = 0; k < DIM; k++ )
				sum += A[i][k] * B[k][j];
			tmpMatr[i][j] = sum;
		}
	MCopy( C, tmpMatr );
}
/*
    матричное возведение в степень	(C = A ^ power)
    A - исходная матрица
    power -	степень	(должна быть больше 1 )
    C - результат
    функция не имеет возвращаемого значения. Результат возвращается
      косвенно через аргумент C.
*/
void MPow( double A[DIM][DIM], int power, double C[DIM][DIM] )
{
	int i, repeatTimes = 100;
	double tmpMatr[DIM][DIM];
	while( repeatTimes-- )
	{
		MCopy( tmpMatr, A );
		MCopy( C, A );
		for( i = 1; i < power; i++ )
			MMul( C, tmpMatr, C );
	}
}
/*****************************************************************************/
/*
    Поиск числа K.
    M - исходная матрица
    delta - погрешность вычислений
    Функция возвращает число K.
*/
int FindK( double M[DIM][DIM], double delta )
{
	double tmpMatr[DIM][DIM];
	int k = 1, numprocs, myid;
	MPI_Comm_size( MPI_COMM_WORLD, &numprocs );
	MPI_Comm_rank( MPI_COMM_WORLD, &myid );
	k += myid;
	do
	{
		k += numprocs;
		MPow( M, k, tmpMatr );
	} while( (1 - tmpMatr[0][DIM - 1]) >= delta );
	return k;
}
/*
   	Функция вычисления плотности распределния
    M - исходная матрица
    K - размерность плотности
    f - указатель на массив, содержащий плотность распределения
    функция не имеет возвращаемого значения. Результат возвращается
      косвенно через аргумент f.
*/
void FindF( double M[DIM][DIM], int K, double *f )
{
	double m1[DIM][DIM], m2[DIM][DIM];
	int i, numprocs, myid;
	MPI_Comm_size( MPI_COMM_WORLD, &numprocs );
	MPI_Comm_rank( MPI_COMM_WORLD, &myid );
	for( i = 1 + myid; i <= K; i += numprocs )
	{
		MPow( M, i, m1 );
		MPow( M, i - 1, m2 );
		f[i - 1] = m1[0][DIM - 1] - m2[0][DIM - 1];
	}
}
/*
   	функция распечатывает на экран значения плотности распределния.
    f - указатель на массив со значениями плотности распределния.
    K - длинна массива.
    функция не имеет возвращаемого значения.
*/
void FPrint( double *f, int K )
{
	int i;
	printf( "  f.Size = %i\n  {\n", K );
	for( i = 0; i < K; i++ )
		printf( "    f[%i] = \t%e\n", i + 1, f[i] );
	printf( "  }\n" );
}
/*
    Проверка суммы значений плотности распределения.
    f - указатель на массив значений плотности распределения.
    K - размерность массива.
    delta - допустимая погрешность.
   	Функция не имеет возвращаемого значения.
*/
void CheckSum( double *f, int K, double delta )
{
	double sum = 0;
	int i;
	for( i = 0; i < K; i++ )
		sum += f[i];
	printf( "sum f[i] = %f\n", sum );
}
/*
    Вычисления математического ожидания.
    f - указатель на массив значений плотности распределения.
    K - размерность массива.
    Функция возвращает значение математического ожидания.
*/
double FindM( double *f, int K )
{
	double res = 0;
	int i;
	for( i = 0; i < K; i++ )
		res += i * f[i];
	return res;
}
/*
    Вычисление дисперсии.
    f - указатель на массив значений плотности распределения.
    K - размерность массива.
   	Функция возвращает значение дисперсии.
*/
double FindD( double *f, int K )
{
	double res = 0, M = FindM( f, K );
	int i;
	for( i = 0; i < K; i++ )
		res += pow( i - M, 2 ) * f[i];
	return res;
}
/*****************************************************************************/
/*
    объединение группы процессов по AND.
    N - количество процессов для объединения.
    f - указатель на массив с
	    плотностями распределения исходных процессов.
	fAnd - указатель на массив с результирующей плотностью распределения.
	K - длинна плотностей распределения f и fAnd.
    функция не имеет возвращаемого значения. Результат возвращается
      косвенно через аргумент fAnd.
*/
void And( int N, double *f, double *fAnd, int K )
{
	double S, M, res;
	int k, n, i, ki;
	for( k = 0; k < K; k++ )
	{
		S = 0;
		for( n = 1; n <= N; n++ )
		{
			M = 1;
			for( i = 1; i <= N; i++ )
			{
				res = 0;
				if( n == i )
					res = 1;
				else
					for( 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;
	}
}
/*
    объединение группы процессов по OR.
    N - количество процессов для объединения.
    f - указатель на массив с
	    плотностями распределения исходных процессов.
	fOr - указатель на массив с результирующей плотностью распределения.
	K - длинна плотностей распределения f и fOr.
    функция не имеет возвращаемого значения. Результат возвращается
      косвенно через аргумент fMN.
*/
void Or( int N, double *f, double *fOr, int K )
{
	double S, M, res;
	int k, n, i, ki;
	for( k = 0; k < K; k++ )
	{
		S = 0;
		for( n = 1; n <= N; n++ )
		{
			M = 1;
			for( i = 1; i <= N; i++ )
			{
				res = 0;
				if( n == i )
					res = 1;
				else
					for( 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;
	}
}
/*
    объединение группы процессов по M из N.
    N - количество процессов для объединения.
    f - указатель на массив с
	    плотностями распределения исходных процессов.
	fMN - указатель на массив с результирующей плотностью распределения.
	K - длинна плотностей распределения f и fMN.
    функция не имеет возвращаемого значения. Результат возвращается
      косвенно через аргумент fMN.
*/
void MfromN( int N, int M, double *f, double *fMN, int K )
{
	int Cnm = CNM( N, M );
	double *fAnd = (double*)malloc( sizeof( double ) * K );
	And( M, f, fAnd, K );
	Or( Cnm, fAnd, fMN, K );
	free( fAnd );
}
/*
    последовательное объединение группы процессов.
    N - количество процессов для объединения.
    f - указатель на массивы с
	    плотностями распределения исходных процессов.
	fP - указатель на массив с результирующей плотностью распределения.
	K - длинна плотностей распределения f[i] и fP.
    функция не имеет возвращаемого значения. Результат возвращается
      косвенно через аргумент fP.
*/
void Posl( int N, double *f[], double *fP, int K )
{
	double S;
	double *f1 = (double*)malloc( sizeof( double ) * K * N );
	int i, j, l, m;
	for( i = 0; i < K * N; i++ )
		fP[i] = f1[i] = 0;
	for( i = 0; i < K; i++ )
		fP[i] = f[0][i];
	for( j = 2; j <= N; j++ ){
		for( l = 1;l <= K * j; l++ ){
			S =	0;
			for( m = 1; m <= K * N; m++ )
				if( m < l && m <= K * (j - 1) && l - m <= K )
					S += fP[m - 1] * f[j - 1][l - m - 1];
			f1[l - 1] = S;
		}
		for( i = 0; i < K * j; i++ )
			fP[i] = f1[i];
	}
	free( f1 );
}
/*****************************************************************************/
/* Entry point */
int main( int argc, char *argv, char *envp )
{
	double delta = 0.00001;
	double MyMatrix[DIM][DIM] = {
		{0,   1,   0,   0,   0,   0,   0,   0},
		{0,   0,   0.5, 0.5, 0,   0,   0,   0},
		{0,   0,   0.5, 0.5, 0,   0,   0,   0},
		{0.5, 0,   0,   0,   0.5, 0,   0,   0},
		{0,   0,   0,   0,   0,   1,   0,   0},
		{0,   0,   0,   0,   0.5, 0,   0.5, 0},
		{0,   0,   0,   0,   0.5, 0,   0,   0.5},
		{0,   0,   0,   0,   0,   0,   0,   1}
	};
	int K, K1, m = 3, i = 5, l = 7;
	double M, D;
	double *f, *fAll[3], *fRes;
	double start, end;
	int numprocs, myid;
	MPI_Init( &argc, &argv );
	MPI_Comm_size( MPI_COMM_WORLD, &numprocs );
	MPI_Comm_rank( MPI_COMM_WORLD, &myid );
	if( myid == 0 )
	{
		printf( "proc count = %i\n", numprocs );
		printf( "Used " );
		MPrint( MyMatrix );
		printf( "Matrix K = " );
		start = MPI_Wtime( );
	}
	K = FindK( MyMatrix, delta );
	MPI_Reduce( &K, &K1, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
	K = K1;
	MPI_Bcast( &K, 1, MPI_INT, 0, MPI_COMM_WORLD );
	/* Вычисление плотности распределения для процесса 0 */
	f = (double*)calloc( K, sizeof( double ) );
	fAll[0] = (double*)calloc( K, sizeof( double ) );
	FindF( MyMatrix, K, f );
	MPI_Reduce( f, fAll[0], K, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
	if( myid == 0 )
	{
		/* Вычисление плотности распределения для объединения
		   группы процессов с 1 по I по функции M из N */
		fAll[1] = (double*)malloc( sizeof( double ) * K );
		MfromN( i, m, fAll[0], fAll[1], K );
		/* Вычисление плотности распределения для объединения
		   группы процессов с I + 1 по I + L по функции AND */
		fAll[2] = (double*)malloc( sizeof( double ) * K );
		Or( l, fAll[0], fAll[2], K );
		/* Вычисление результирующей плотности распределения
		   последовательное объединение. */
		fRes = (double*)malloc( sizeof( double ) * 3 * K );
		Posl( 3, fAll, fRes, K );
		K *= 3;
		end = MPI_Wtime( );
		printf( "%i\n==============================\n", K1 );
		/* вывод плотности распределения */
		FPrint( fRes, K );
		/* проверка суммы */
		CheckSum( fRes, K, delta );
		/* вычисление математического ожидания */
		M = FindM( fRes, K );
		/* вычисление дисперсии */
		D = FindD( fRes, K );
		printf( "M = %f\nD = %f\n time = %f s\n", M, D, end - start );
		fflush( stdout );
		free( fAll[1] );
		free( fAll[2] );
		free( fRes );
	}
	free( f );
	free( fAll[0] );
	//_getch( );
	/* очищаем динамическую память */
	MPI_Finalize( );
	return 0;
}