#include "boulanger.h"


//renvoie la matrice de la transformée de Fourier discrète
void fourierDiscrete (int dimension, 
						double complex **matrice) //[dimension][dimension]
{
	int i = 0, j = 0;
	for (i = 0 ; i < dimension ; i++)
	{
		for (j = 0 ; j < dimension ; j++)
		{
			matrice[i][j] = cexp(- 2 * I * M_PI * i * j / dimension) / sqrt(dimension);
		}
	}	
}	// fin void fourierDiscrete(int dimension, double complex matrice[dimension][dimension])


//renvoie la matrice de la transformée de Fourier inverse
void fourierInverse (int dimension, 
					double complex **matrice) //[dimension][dimension]
{
	int i = 0, j = 0;
	for (i = 0 ; i < dimension ; i++)
	{
		for (j = 0 ; j < dimension ; j++)
		{
			matrice[i][j] = cexp(2 * I * M_PI * i * j / dimension) / sqrt(dimension);
		}
	}	
}	// fin void fourierInverse(int dimension, double complex matrice[dimension][dimension])

double coeffBoulanger (double alpha, double x)
{
	if (floor (x) != x)
	{
		return sin (M_PI * alpha * x) / sin (M_PI * x);
	}
	else
	{
		return alpha ;
	}		
}	// fin double coeffBoulanger(double alpha, double x)


void boulangerFerme1 (int dimension, 
						int nombreRectangles, 
						int largeursBlocs[nombreRectangles], 
						double complex matrice[dimension][dimension])
{
	int i = 0 , j = 0 , k = 0;
	double *temp = {0}; 
	
	temp = malloc ((nombreRectangles + 1) * sizeof(double));
	temp[0] = 0;
	for (i = 1 ; i <= nombreRectangles ; i++)
	{
		temp[i] = sommeTableau(largeursBlocs, i);
	}
	for (i = 0 ; i < nombreRectangles ; i++)
	{
		for (j = temp[i] ; j < temp[i + 1] ; j++)
		{
			for (k = 0 ; k < dimension ; k++)
			{
				double x = (double) (k) / dimension - (double)((j - temp[i])) / largeursBlocs[i];
				matrice[j][k] = cexp(2 * I * M_PI * (double)(k) * temp[i] / (double)(dimension)) * cexp (I * M_PI * (largeursBlocs[i] - 1) * x) * coeffBoulanger(largeursBlocs[i], x) / sqrt(largeursBlocs[i] * dimension);
			}	
		}	
	}
	free(temp);	
}	// fin void boulangerFerme1


// Retourne la matrice du boulanger fermé en fonction de la dimension , du nombre de divisions et de la largeurs des rectangles (mis dans un tableau de taille nombreRectangles)

void boulangerFerme (int dimension, 
					int nombreRectangles, 
					int largeursBlocs[nombreRectangles], 
					double complex matrice[dimension][dimension])
{
	int i = 0 , j = 0 , k = 0 , l = 0;
	double complex** fourier = malloc (dimension * sizeof(double complex)); // matrice qui servira à écrire la Fourier inverse
			for (l = 0 ; l < dimension ; l++)
			{
				fourier[l] = malloc (dimension * sizeof(double complex));
			}
	double complex** blocs = malloc (dimension * sizeof(double complex)); // matrice qui servira à écrire la matrice diagonale par blocs
			for (l = 0 ; l < dimension ; l++)
			{
				blocs[l] = malloc (dimension * sizeof(double complex));
			}
	char TRANSA = 'N' , TRANSB = 'N'; // arguments de zgemm
	double complex alpha = 1, beta = 0; // arguments de zgemm
   	
	if (sommeTableau(largeursBlocs, nombreRectangles) == dimension)
	{
		for (i = 0 ; i < nombreRectangles ; i++) // cette boucle crée une matrice diagonale par blocs dont les blocs sont des matrices de Fourier
		{
			double complex** temp = malloc (largeursBlocs[i] * sizeof(double complex)); // bloc i de la matrice
			for (l = 0 ; l < largeursBlocs[i] ; l++)
			{
				temp[l] = malloc (largeursBlocs[i] * sizeof(double complex));
			} 
			for (j = sommeTableau(largeursBlocs, i) ; j < sommeTableau(largeursBlocs, i + 1) ; j++)
			{
				for (k = 0 ; k < dimension ; k++)
				{
					if (k >= sommeTableau(largeursBlocs, i) && k < sommeTableau(largeursBlocs, i + 1))
					{
						fourierDiscrete (largeursBlocs[i], temp);
						blocs[j][k] = temp[j - sommeTableau(largeursBlocs, i)][k - sommeTableau(largeursBlocs, i)];
					}
					else
					{
						blocs[j][k] = 0;
					}
					printf("%d, %d, %f, %f \n", j, k, creal(blocs[j][k]), cimag(blocs[j][k]));	
				}
			}
			for(l = 0 ; l < largeursBlocs[i] ; l++)
			{
    			free(temp[l]);
			}
			free(temp);	
		}
		
		fourierInverse (dimension, fourier);

		zgemm_ (&TRANSA, &TRANSB, &dimension, &dimension, &dimension, &alpha, &(fourier[0][0]), &dimension, &(blocs[0][0]), &dimension, &beta, &(matrice[0][0]), &dimension);
		
		for(l = 0 ; l < dimension ; l++)
			{
    			free(fourier[l]);
			}
			free(fourier);
		for(l = 0 ; l < dimension ; l++)
			{
    			free(blocs[l]);
			}
			free(blocs);
	}
	else
	{
		printf("Problème de dimensions \n");
	}
}	// fin void boulangerFerme

// trou est un tableau contenant des 0 et des 1 : 0 si le rectangle i est un trou, 1 sinon. Il faudra par la suite mettre dedans le nombre par lequel on divise la largeur du cutoff.
void boulangerOuvert1 (int dimension, 
						int nombreRectangles, 
						int largeursBlocs[nombreRectangles], 
						int trou[nombreRectangles], 
						double complex matrice[dimension][dimension])
{
	int i = 0 , j = 0 , k = 0 ;
	double *temp = {0} ; 
	
	temp = malloc ((nombreRectangles + 1) * sizeof(double)) ;
	temp[0] = 0 ;
	for (i = 1 ; i <= nombreRectangles ; i++)
	{
		temp[i] = sommeTableau(largeursBlocs, i);
	}
	
	for (i = 0 ; i < nombreRectangles ; i++)
	{
		if (trou[i] != 0)
		{
			for (j = temp[i] ; j < temp[i + 1] ; j++ )
			{
				for (k = 0 ; k < dimension ; k++ )
				{
					double x = (double) (k) / dimension - (double)((j - temp[i])) / largeursBlocs[i] ;
					matrice[j][k] = cexp(2 * I * M_PI * j * temp[i] / dimension) * cexp (I * M_PI * (largeursBlocs[i] - 1) * x) * coeffBoulanger(largeursBlocs[i], x) / sqrt(largeursBlocs[i] * dimension);
				}	
			}	
		}
	}
	free(temp) ;	
}	// fin void boulangerOuvert1



double cutoff (double delta, 
				double lambda, 
				double x) //delta : largeur à gauche, lambda : largeur à droite, x : variable
{
	if (x <= 0 || x >=1)
	{
		return 0;	
	}		
	else if (0 < x && x < delta)
	{
		return exp(1 / pow(delta, 2) - 1 / pow(x, 2)) * (1 - exp(-1 / pow(x - delta, 2)));	
	}
	else if (delta < x && x < lambda)
	{
		return 1;	
	}
	else if (lambda < x && x < 1)
	{
		return exp(-1 / pow(x - 1, 2) + 1 / pow(lambda - 1, 2)) * (1 - exp(-1 / pow(x - lambda, 2)));	
	}
	else
	{
		return 0;
	}	// fin double cutoff
}	

int correlation (double x, 
					int tailleTableau, 
					double tableau[])
{
	int i = 0;
	double nombre = 0;
	for (i = 0 ; i < tailleTableau ; i++)
	{
		if (tableau[i] <= x)
		{
			nombre = nombre + 1;
		}
		else
		{
			nombre = nombre;
		}		
	}
	return nombre;		
}	// fin int correlation

int rectMin (int nombreRectangles, double trou[nombreRectangles])
{
	int resultat = 0, i = 0;
	while (trou[i] == 0)
	{
		resultat++ ;
		i++ ;
	}	
	return resultat;	
}	// fin int rectMin

int rectMax (int nombreRectangles, double trou[])
{
	int resultat = nombreRectangles - 1, i = nombreRectangles - 1;
	while (trou[i] == 0)
	{
		resultat-- ;
		i--;		
	}		
	return resultat;	
}	// fin int rectMax	

void boulangerOuvertLisse (int dimension, 
							int nombreRectangles, 
							int largeursBlocs[nombreRectangles], 
							double trou[nombreRectangles], 
							double complex **matrice)
{
	int i = 0, j = 0, k = 0;
	int m = rectMin (nombreRectangles, trou), M = rectMax (nombreRectangles, trou); // premier et dermier rectangles qui ne sont pas dans le trou
	int *abscisse = malloc ((nombreRectangles + 1) * sizeof(int));
	abscisse[0] = 0;
	for (i = 1 ; i <= nombreRectangles ; i++)
	{
		abscisse[i] = (int)(sommeTableau(largeursBlocs, i));
	}
		// largeurs de montée des cutoffs
	double deltaq = 0, lambdaq = 0; 
	double deltap = floor((double)(abscisse[m]) / dimension + ((double)(largeursBlocs[m] * abscisse[m]) / pow(dimension, 2)) / (1 - (double)(largeursBlocs[m]) / dimension));
	double lambdap = ceil((double)(abscisse[M]) / dimension + ((double)(largeursBlocs[M] * abscisse[M]) / pow(dimension, 2)) / (1 - (double)(largeursBlocs[M]) / dimension));
		
	double complex** blocs = malloc (dimension * sizeof(double complex)); // matrice qui servira à écrire la matrice diagonale par blocs
	for (i = 0 ; i < dimension ; i++)
	{
		blocs[i] = malloc (dimension * sizeof(double complex));
	}
	
	char TRANSA = 'N', TRANSB = 'N'; // arguments de zgemm
	double complex alpha = 1, beta = 0; // arguments de zgemm

	if (abscisse[nombreRectangles] == dimension)
	{
		for (i = 0 ; i < nombreRectangles ; i++) // cette boucle crée une matrice diagonale par blocs dont les blocs sont des matrices de Fourier
		{
			if (trou[i] != 0)
			{
				deltaq = floor(((double)(largeursBlocs[i] * abscisse[m]) / pow(dimension, 2)) / (1 - (double)(largeursBlocs[m]) / dimension));
				lambdaq = ceil(((double)(largeursBlocs[i] * abscisse[M]) / pow(dimension, 2)) / (1 - (double)(largeursBlocs[M]) / dimension));
							
				double complex** temp = malloc (largeursBlocs[i] * sizeof(double complex)); // bloc i de la matrice
				for (j = 0 ; j < largeursBlocs[i] ; j++)
				{
					temp[j] = malloc (largeursBlocs[i] * sizeof(double complex));
				} 
				fourierDiscrete (largeursBlocs[i], temp);
				for (j = abscisse[i] ; j < abscisse[i + 1] ; j++)
				{
					for (k = 0 ; k < dimension ; k++)
					{
						if (abscisse[i] <= k && k < abscisse[i + 1])
						{
							blocs[j][k] = temp[j - abscisse[i]][k - abscisse[i]] * cutoff(deltaq, lambdaq, (double)(j - abscisse[i]) / dimension) * cutoff(deltap, lambdap, (double)(k) / dimension);
						}
						else
						{
							blocs[j][k] = 0;
						}		
					}
				//free(temp[j - abscisse[i]]);
				}
				for(j = 0 ; j < largeursBlocs[i] ; j++)
				{
    				free(temp[j]);
				}
				free(temp);
			}
		}
		free(abscisse);
		double complex** fourier = malloc (dimension * sizeof(double complex)); // matrice qui servira à écrire la Fourier inverse
		for (i = 0 ; i < dimension ; i++)
		{
			fourier[i] = malloc (dimension * sizeof(double complex));
		}
		fourierInverse (dimension, fourier);
		zgemm_ (&TRANSA, &TRANSB, &dimension, &dimension, &dimension, &alpha, &(fourier[0][0]), &dimension, &(blocs[0][0]), &dimension, &beta, &(matrice[0][0]), &dimension);
		
		for(i = 0 ; i < dimension ; i++)
			{
    			free(fourier[i]);
			}
			free(fourier);
		for(i = 0 ; i < dimension ; i++)
			{
    			free(blocs[i]);
			}
			free(blocs);
	}
	else
	{
		printf("Problème de dimensions \n");
		exit(1);
	}
}	// fin void boulangerOuvertLisse