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

long double **multMatrix(long double **matrixA, long double **matrixB, int dim)
{
	int i, j, elements;

	long double **returnValue;

	elements = 0;

	returnValue = (long double**)malloc(sizeof (long double *) * dim);
	for (i=0; i < dim; i++)
	{
		returnValue[i] = (long double *)malloc(sizeof(long double)* dim);
		for (j=0; j < dim; j++)
		{
			returnValue[i][j] = matrixA[i][j] * matrixB[i][j];
			elements++;
		}
	}
	return returnValue;
}

void matrixOperations(int dim, int dummy)
{
	long double **matA;
	long double **matB;
	long double **matC;
	long double average=0.0;


	int i = 0;
	int j = 0;
	int weight = 0;

	matA = (long double**)malloc(sizeof (long double *) * dim);
	matB = (long double**)malloc(sizeof (long double *) * dim);

	for (i = 0; i < dim; i++)
	{
		matA[i] = (long double *)malloc(sizeof(long double) * dim);
		matB[i] = (long double *)malloc(sizeof(long double) * dim);
		for (j = 0; j < dim; j++)
		{
			matA[i][j] = (i+j)/(dummy * M_PI);
			matB[i][j] = (i+j)/(M_PI);
		}
	}

	matC = multMatrix(matA, matB, dim);
	for (i = 0; i < dim; i++)
	{
		for (j = 0; j < dim; j++)
		{
			average = (average*weight + matC[i][j])/(weight + 1);
		}
		free(matA[i]);
		free(matB[i]);
		free(matC[i]);
	}
	printf("Average matrix value: %Lf\n", average);
	free(matA);
	free(matB);
	free(matC);
}

void normalFpOperations()
{
	long double value, power, squareRoot, sine;
	long double aValue = 0.0;
	long double aPower = 0.0;
	long double aSquareRoot = 0.0;
	long double aSine = 0.0;

	int i,j, weight;

	weight = 0;

	for (i=0; i<500; i++)
	{
		for (j=0; j<500; j++)
		{
			value = (i+j)/M_PI;
			squareRoot = sqrt(value);
			power = pow(value, M_PI);
			sine = sin(value);
			aValue = (aValue*weight + value)/(weight + 1);
			aPower = (aPower*weight + power)/(weight + 1);
			aSquareRoot = (aSquareRoot*weight + squareRoot)/(weight + 1);
			aSine = (aSine*weight + sine)/(weight + 1);
			weight++;
		}
	}

	printf("Number of calculations of each type: %d\n", weight);
	printf("Average value: %Lf\n", aValue);
	printf("Average power: %Lf\n", aPower);
	printf("Average square root: %Lf\n", aSquareRoot);
	printf("Average sine: %Lf\n", aSine);
}

int main(int argc, char *argv[]) {
	clock_t startTime, endTime;


	startTime = clock();
	printf("FP operations...\n");
	normalFpOperations();
	endTime = clock();
	printf("FP operations done in %f seconds\n",(double)(endTime - startTime)/CLOCKS_PER_SEC);
	startTime = clock();
	printf("Matrix operations, 1000x1000 elements\n");
	matrixOperations(1000, 1);
	matrixOperations(1000, 2);
	matrixOperations(1000, 3);
	matrixOperations(1000, 4);
	matrixOperations(1000, 5);
	endTime = clock();
	printf("Matrix operations done in %f seconds\n",(double)(endTime - startTime)/CLOCKS_PER_SEC);
	printf("Done.\n");
	return EXIT_SUCCESS;
}