#include "curand_kernel.h"
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "Kalibrace.h"
#include <stdio.h>
#include <time.h>
#include <math.h>

__global__ void curand_setup(curandState * state, unsigned long seed)
{
	int id = threadIdx.x;
    curand_init ( seed, id, 0, &state[id] );
}

__global__ void cuda_fitness(double* data,Bod* Body,double* vysledek)
{
	int k = threadIdx.x;
	
}

__global__ void ElitaKernel(double** data, double* vysledek, int* int_parametry)
{
	int Maxi = 0;
	for (int i=1;i<int_parametry[1];i++)
		if (data[i][int_parametry[0]]<data[Maxi][int_parametry[0]])
			Maxi = i;
	for (int i=0;i<int_parametry[0]+1;i++)
		vysledek[i]=data[Maxi][i];
}
__global__ void StandartDeviationKernel(double** data, int* int_parametry, double* vysledek)
{
	double prum = 0;
	for (int i=0;i<int_parametry[1];i++)
		prum += data[i][int_parametry[0]];
	prum = prum/int_parametry[1];
	double suma = 0;
	for (int i=0;i<int_parametry[1];i++)
		suma += (data[i][int_parametry[0]]-prum)*(data[i][int_parametry[0]]-prum);
	*vysledek = sqrt(suma/(double)(int_parametry[1]-1));
}
__global__ void ReprodukceKernel(double** data, double** vysledek, double* horni_mez, double* spodni_mez, double* double_parametry, int* int_parametry, Bod* Body, curandState* state)
{
		int X[3];
		do
		for (int k=0;k<3;k++)
		{
			int q = abs((int)curand(&state[threadIdx.x])-1);
			X[k]=(int)((((double)int_parametry[1]-1)*(double)q-1)/2147483648.0);
			if ((X[k]>=threadIdx.x)||(X[k]<0))
				X[k]++;
		}
		while((X[0]==X[1])||(X[0]==X[2])||X[1]==X[2]);
		for (int k=0;k<int_parametry[0];k++)
		{
			int q = abs((int)curand(&state[threadIdx.x])-1);
			if (((double)q/2147483648.0)<double_parametry[1])
			{
				vysledek[threadIdx.x][k]=double_parametry[0]*(data[X[0]][k]-data[X[1]][k])+data[X[2]][k];
				if ((vysledek[threadIdx.x][k]<spodni_mez[k])||(vysledek[threadIdx.x][k]>horni_mez[k]))
					vysledek[threadIdx.x][k]=((horni_mez[k]-spodni_mez[k])*(double)abs((int)curand(&state[threadIdx.x]))/(2147483648)+spodni_mez[k]);
			}
			else
				vysledek[threadIdx.x][k]=data[threadIdx.x][k];
		}
		double odchylka = 0;
		double SinAlfa = sin(vysledek[threadIdx.x][2]);
		double CosAlfa = cos(vysledek[threadIdx.x][2]);
		double SinBeta = sin(vysledek[threadIdx.x][3]);
		double CosBeta = cos(vysledek[threadIdx.x][3]);
		double SinGama = sin(vysledek[threadIdx.x][4]);
		double CosGama = cos(vysledek[threadIdx.x][4]);
		for (int k=0;k<int_parametry[2];k++)
		{
			if ((Body[k].ZobrX>0)&&(Body[k].ZobrY>0))
			{
				double x = vysledek[threadIdx.x][0]*(CosGama*CosBeta*Body[k].x+CosGama*SinAlfa*SinBeta*Body[k].y+CosGama*CosAlfa*SinBeta*Body[k].z-SinGama*CosAlfa*Body[k].y+SinAlfa*SinGama*Body[k].z+vysledek[threadIdx.x][5])/(-SinBeta*Body[k].x+SinAlfa*CosBeta*Body[k].y+CosAlfa*CosBeta*Body[k].z+vysledek[threadIdx.x][7])+int_parametry[3];
				double y = vysledek[threadIdx.x][1]*(SinGama*CosBeta*Body[k].x+SinGama*SinAlfa*SinBeta*Body[k].y+SinGama*CosAlfa*SinBeta*Body[k].z+CosAlfa*CosGama*Body[k].y-SinAlfa*CosGama*Body[k].z+vysledek[threadIdx.x][6])/(-SinBeta*Body[k].x+SinAlfa*CosBeta*Body[k].y+CosAlfa*CosBeta*Body[k].z+vysledek[threadIdx.x][7])+int_parametry[4];
				odchylka += ((Body[k].ZobrX-x)*(Body[k].ZobrX-x)+(Body[k].ZobrY-y)*(Body[k].ZobrY-y));
			}
		}
		odchylka = sqrt(odchylka);
		if (odchylka<data[threadIdx.x][int_parametry[0]])
			vysledek[threadIdx.x][int_parametry[0]]=odchylka;
		else
			for (int k=0;k<int_parametry[0]+1;k++)
				vysledek[threadIdx.x][k]=data[threadIdx.x][k];
}


class CUDA_Kalibrace : public Kalibrace
{
private:
	double** devp_data;
	double** dev_data;
	double** devp_vysledek;
	double** dev_vysledek;
	double* dev_meze_h;
	double* dev_meze_d;
	double* dev_parametry_double; // [0]..F, [1]..CR
	int* dev_parametry_int; // [0]..D, [1]..NP, [2]..count, [3]..Cx, [4]..Cy
	Bod* dev_body; 
	curandState* dev_rand;

public:
	CUDA_Kalibrace() : Kalibrace()
	{	cudaError_t cudaStatus;
		
		cudaStatus = cudaSetDevice(0);
		

		dev_data = new double*[NP];
		dev_vysledek = new double*[NP];
		for (int i=0;i<NP;i++)
		{
			cudaStatus =cudaMalloc((void**)&dev_data[i],(D+1)*sizeof(double));
			cudaStatus =cudaMalloc((void**)&dev_vysledek[i],(D+1)*sizeof(double));
		}
		cudaStatus =cudaMalloc((void**)&devp_data,NP*sizeof(double*));
		cudaStatus =cudaMalloc((void**)&devp_vysledek,NP*sizeof(double*));
		cudaStatus =cudaMalloc((void**)&dev_meze_h,D*sizeof(double));
		cudaStatus =cudaMalloc((void**)&dev_meze_d,D*sizeof(double));
		cudaStatus =cudaMemcpy(devp_data,dev_data,NP*sizeof(double*),cudaMemcpyHostToDevice);
		cudaStatus =cudaMemcpy(devp_vysledek,dev_vysledek,NP*sizeof(double*),cudaMemcpyHostToDevice);
		cudaStatus =cudaMalloc((void**)&dev_body,Item.GetCount()*sizeof(Bod));
		cudaStatus =cudaMalloc((void**)&dev_parametry_double,2*sizeof(double));
		cudaStatus =cudaMalloc((void**)&dev_parametry_int,5*sizeof(int));
		cudaStatus =cudaMalloc((void**)&dev_rand, NP*sizeof( curandState ) );
		curand_setup <<< 1, NP >>> ( dev_rand, time(NULL) );

	}

	~CUDA_Kalibrace()
	{
		for (int i=0;i<NP;i++)
		{
			cudaFree(dev_data[i]);
			cudaFree(dev_vysledek[i]);
		}
		cudaFree(devp_data);
		cudaFree(devp_vysledek);
		cudaFree(dev_meze_h);
		cudaFree(dev_meze_d);
		cudaFree(dev_parametry_double);
		cudaFree(dev_parametry_int);
		cudaFree(dev_body);
		cudaFree(dev_rand);
		delete(dev_data);
		delete(dev_vysledek);
	}

	void Start_Evoluce()
	{
		cudaError_t cudaStatus;
		Diferencialni::Start_Evoluce();
		
		for(int i=0;i<NP;i++)
			cudaStatus =cudaMemcpy(dev_data[i],data[i],(D+1)*sizeof(double),cudaMemcpyHostToDevice);
		double* Help = new double[D];
		for (int i=0;i<D;i++)
			Help[i]=Meze[i][0];
		cudaStatus =cudaMemcpy(dev_meze_d,Help,D*sizeof(double),cudaMemcpyHostToDevice);
		for (int i=0;i<D;i++)
			Help[i]=Meze[i][1];
		cudaStatus =cudaMemcpy(dev_meze_h,Help,D*sizeof(double),cudaMemcpyHostToDevice);
		delete(Help);
		Bod* Help2 = new Bod[Item.GetCount()];
		for (int i=0;i<Item.GetCount();i++)
			Help2[i]=Item.GetBod(i);
		cudaStatus =cudaMemcpy(dev_body,Help2,(Item.GetCount())*sizeof(Bod),cudaMemcpyHostToDevice);
		delete(Help2);
		cudaStatus =cudaMemcpy(dev_parametry_double, &F,sizeof(double),cudaMemcpyHostToDevice);
		cudaStatus =cudaMemcpy(&dev_parametry_double[1], &CR,sizeof(double),cudaMemcpyHostToDevice);
		cudaStatus =cudaMemcpy(dev_parametry_int, &D,sizeof(int),cudaMemcpyHostToDevice);
		cudaStatus =cudaMemcpy(&dev_parametry_int[1], &NP, sizeof(int),cudaMemcpyHostToDevice);
		int l = Item.GetCount();
		cudaStatus =cudaMemcpy(&dev_parametry_int[2], &l, sizeof(int),cudaMemcpyHostToDevice);
		l = Q.GetCx();
		cudaStatus =cudaMemcpy(&dev_parametry_int[3], &l, sizeof(int),cudaMemcpyHostToDevice);
		l = Q.GetCy();
		cudaStatus =cudaMemcpy(&dev_parametry_int[4], &l, sizeof(int),cudaMemcpyHostToDevice);
	}

	virtual void Reprodukce()
	{
		ReprodukceKernel<<<1,NP>>>(devp_data,devp_vysledek,dev_meze_h,dev_meze_d,dev_parametry_double,dev_parametry_int,dev_body,dev_rand);
		double** help = devp_data;
		devp_data = devp_vysledek;
		devp_vysledek = help;
	}

	virtual void Elita(double* hodnota)
	{
		cudaError_t cudaStatus;
		double* dev_elita;
		cudaStatus =cudaMalloc((void**)&dev_elita,(D+1)*sizeof(double));
		ElitaKernel<<<1,1>>>(devp_data,dev_elita,dev_parametry_int);
		cudaStatus =cudaMemcpy(hodnota,dev_elita,(D+1)*sizeof(double),cudaMemcpyDeviceToHost);
		cudaFree(dev_elita);
		
	}

	virtual double HodnotaElita()
	{
		double* dev_elita;
		double hodnota;
		cudaMalloc((void**)&dev_elita,(D+1)*sizeof(double));
		ElitaKernel<<<1,1>>>(devp_data,dev_elita,dev_parametry_int);
		cudaMemcpy(&hodnota,&dev_elita[D],sizeof(double),cudaMemcpyDeviceToHost);
		cudaFree(dev_elita);
		return hodnota;
	}

	void SetNP(int A)
	{
		Kalibrace::SetNP(A);
		cudaFree(dev_rand);
		cudaMalloc((void**)&dev_rand, A*sizeof( curandState ) );
		curand_setup <<< 1, A >>> ( dev_rand, time(NULL) );

	}

	void AddObjectPoint(double X,double Y, double Z)
	{
		Kalibrace::AddObjectPoint(X,Y,Z);
		cudaFree(dev_body);
		cudaMalloc((void**)&dev_body,Item.GetCount()*sizeof(Bod));
	}

	void BlockAsObject(double a,double b, double c)
	{
		Kalibrace::BlockAsObject(a,b,c);
		cudaFree(dev_body);
		cudaMalloc((void**)&dev_body,Item.GetCount()*sizeof(Bod));
	}

	virtual double StandartDeviation()
	{
		double hodnota;
		double* dev_deviation;
		cudaMalloc((void**)&dev_deviation,sizeof(double));
		StandartDeviationKernel<<<1,1>>>(devp_data,dev_parametry_int,dev_deviation);
		cudaMemcpy(&hodnota,dev_deviation,sizeof(double),cudaMemcpyDeviceToHost);
		cudaFree(dev_deviation);
		return hodnota;
	}

};


int main()
{
	CUDA_Kalibrace X;
	Kalibrace K;
	X.SetCenterOfImage(400,300);
	K.SetCenterOfImage(400,300);

	K.BlockAsObject(6,3,8);
	K.SetCenterOfImage(200.5,150.5);
	K.SetDisplayedPoint(38,210,0);
	K.SetDisplayedPoint(237,273,1);
	K.SetDisplayedPoint(65,109,2);
	K.SetDisplayedPoint(267,171,3);
	K.SetDisplayedPoint(113,194,4);
	K.SetDisplayedPoint(287,247,5);
	K.SetDisplayedPoint(137,104,6);
	K.SetDisplayedPoint(313,159,7);

	X.BlockAsObject(6,3,8);
	
	X.SetCenterOfImage(200.5,150.5);
	X.SetDisplayedPoint(38,210,0);
	X.SetDisplayedPoint(237,273,1);
	X.SetDisplayedPoint(65,109,2);
	X.SetDisplayedPoint(267,171,3);
	X.SetDisplayedPoint(113,194,4);
	X.SetDisplayedPoint(287,247,5);
	X.SetDisplayedPoint(137,104,6);
	X.SetDisplayedPoint(313,159,7);

	/*X.AddObjectPoint(-7,1.5,-4);
	X.AddObjectPoint(-1,1.5,-4);
	X.AddObjectPoint(-7,-1.5,-4);
	X.AddObjectPoint(-1,-1.5,-4);
	X.AddObjectPoint(-7,1.5,4);
	X.AddObjectPoint(-1,1.5,4);
	X.AddObjectPoint(-7,-1.5,4);
	X.AddObjectPoint(-1,-1.5,4);
	X.AddObjectPoint(1,1.5,-4);
	X.AddObjectPoint(7,1.5,-4);
	X.AddObjectPoint(1,-1.5,-4);
	X.AddObjectPoint(7,-1.5,-4);
	X.AddObjectPoint(1,1.5,4);
	X.AddObjectPoint(7,1.5,4);
	X.AddObjectPoint(1,-1.5,4);
	X.AddObjectPoint(7,-1.5,4);
	X.SetDisplayedPoint(126,300,0);
	X.SetDisplayedPoint(318,361,1);
	X.SetDisplayedPoint(151,200,2);
	X.SetDisplayedPoint(349,261,3);
	X.SetDisplayedPoint(212,290,4);
	X.SetDisplayedPoint(381,343,5);
	X.SetDisplayedPoint(238,202,6);
	X.SetDisplayedPoint(409,254,7);
	X.SetDisplayedPoint(387,382,8);
	X.SetDisplayedPoint(591,447,9);
	X.SetDisplayedPoint(415,282,10);
	X.SetDisplayedPoint(620,345,11);
	X.SetDisplayedPoint(441,361,12);
	X.SetDisplayedPoint(620,419,13);
	X.SetDisplayedPoint(468,273,14);
	X.SetDisplayedPoint(647,329,15);
	K.AddObjectPoint(-7,1.5,-4);
	K.AddObjectPoint(-1,1.5,-4);
	K.AddObjectPoint(-7,-1.5,-4);
	K.AddObjectPoint(-1,-1.5,-4);
	K.AddObjectPoint(-7,1.5,4);
	K.AddObjectPoint(-1,1.5,4);
	K.AddObjectPoint(-7,-1.5,4);
	K.AddObjectPoint(-1,-1.5,4);
	K.AddObjectPoint(1,1.5,-4);
	K.AddObjectPoint(7,1.5,-4);
	K.AddObjectPoint(1,-1.5,-4);
	K.AddObjectPoint(7,-1.5,-4);
	K.AddObjectPoint(1,1.5,4);
	K.AddObjectPoint(7,1.5,4);
	K.AddObjectPoint(1,-1.5,4);
	K.AddObjectPoint(7,-1.5,4);
	K.SetDisplayedPoint(126,300,0);
	K.SetDisplayedPoint(318,361,1);
	K.SetDisplayedPoint(151,200,2);
	K.SetDisplayedPoint(349,261,3);
	K.SetDisplayedPoint(212,290,4);
	K.SetDisplayedPoint(381,343,5);
	K.SetDisplayedPoint(238,202,6);
	K.SetDisplayedPoint(409,254,7);
	K.SetDisplayedPoint(387,382,8);
	K.SetDisplayedPoint(591,447,9);
	K.SetDisplayedPoint(415,282,10);
	K.SetDisplayedPoint(620,345,11);
	K.SetDisplayedPoint(441,361,12);
	K.SetDisplayedPoint(620,419,13);
	K.SetDisplayedPoint(468,273,14);
	K.SetDisplayedPoint(647,329,15);*/
	
	/*X.AddObjectPoint(-7,3.5,-5);
	X.AddObjectPoint(-1,3.5,-5);
	X.AddObjectPoint(-7,0.5,-5);
	X.AddObjectPoint(-1,0.5,-5);
	X.AddObjectPoint(-7,3.5,3);
	X.AddObjectPoint(-1,3.5,3);
	X.AddObjectPoint(-7,0.5,3);
	X.AddObjectPoint(-1,0.5,3);
	X.AddObjectPoint(1,3.5,-3);
	X.AddObjectPoint(7,3.5,-3);
	X.AddObjectPoint(1,0.5,-3);
	X.AddObjectPoint(7,0.5,-3);
	X.AddObjectPoint(1,3.5,5);
	X.AddObjectPoint(7,3.5,5);
	X.AddObjectPoint(1,0.5,5);
	X.AddObjectPoint(7,0.5,5);
	X.AddObjectPoint(-7,-0.5,-5);
	X.AddObjectPoint(-1,-0.5,-5);
	X.AddObjectPoint(-7,-3.5,-5);
	X.AddObjectPoint(-1,-3.5,-5);
	X.AddObjectPoint(-7,-0.5,3);
	X.AddObjectPoint(-1,-0.5,3);
	X.AddObjectPoint(-7,-3.5,3);
	X.AddObjectPoint(-1,-3.5,3);
	X.AddObjectPoint(1,-0.5,-3);
	X.AddObjectPoint(7,-0.5,-3);
	X.AddObjectPoint(1,-3.5,-3);
	X.AddObjectPoint(7,-3.5,-3);
	X.AddObjectPoint(1,-0.5,5);
	X.AddObjectPoint(7,-0.5,5);
	X.AddObjectPoint(1,-3.5,5);
	X.AddObjectPoint(7,-3.5,5);
	X.SetDisplayedPoint(94,367,0);
	X.SetDisplayedPoint(289,431,1);
	X.SetDisplayedPoint(121,268,2);
	X.SetDisplayedPoint(318,329,3);
	X.SetDisplayedPoint(184,350,4);
	X.SetDisplayedPoint(356,406,5);
	X.SetDisplayedPoint(209,263,6);
	X.SetDisplayedPoint(382,316,7);
	X.SetDisplayedPoint(374,446,8);
	X.SetDisplayedPoint(573,509,9);
	X.SetDisplayedPoint(401,345,10);
	X.SetDisplayedPoint(604,410,11);
	X.SetDisplayedPoint(431,419,12);
	X.SetDisplayedPoint(607,475,13);
	X.SetDisplayedPoint(454,332,14);
	X.SetDisplayedPoint(632,387,15);
	X.SetDisplayedPoint(130,234,16);
	X.SetDisplayedPoint(327,295,17);
	X.SetDisplayedPoint(159,131,18);
	X.SetDisplayedPoint(357,192,19);
	X.SetDisplayedPoint(218,231,20);
	X.SetDisplayedPoint(392,285,21);
	X.SetDisplayedPoint(242,142,22);
	X.SetDisplayedPoint(418,195,23);
	X.SetDisplayedPoint(441,312,24);
	X.SetDisplayedPoint(613,376,25);
	X.SetDisplayedPoint(440,212,26);
	X.SetDisplayedPoint(644,274,27);
	X.SetDisplayedPoint(463,302,28);
	X.SetDisplayedPoint(642,357,29);
	X.SetDisplayedPoint(489,214,30);
	X.SetDisplayedPoint(669,267,31);

	K.AddObjectPoint(-7,3.5,-5);
	K.AddObjectPoint(-1,3.5,-5);
	K.AddObjectPoint(-7,0.5,-5);
	K.AddObjectPoint(-1,0.5,-5);
	K.AddObjectPoint(-7,3.5,3);
	K.AddObjectPoint(-1,3.5,3);
	K.AddObjectPoint(-7,0.5,3);
	K.AddObjectPoint(-1,0.5,3);
	K.AddObjectPoint(1,3.5,-3);
	K.AddObjectPoint(7,3.5,-3);
	K.AddObjectPoint(1,0.5,-3);
	K.AddObjectPoint(7,0.5,-3);
	K.AddObjectPoint(1,3.5,5);
	K.AddObjectPoint(7,3.5,5);
	K.AddObjectPoint(1,0.5,5);
	K.AddObjectPoint(7,0.5,5);
	K.AddObjectPoint(-7,-0.5,-5);
	K.AddObjectPoint(-1,-0.5,-5);
	K.AddObjectPoint(-7,-3.5,-5);
	K.AddObjectPoint(-1,-3.5,-5);
	K.AddObjectPoint(-7,-0.5,3);
	K.AddObjectPoint(-1,-0.5,3);
	K.AddObjectPoint(-7,-3.5,3);
	K.AddObjectPoint(-1,-3.5,3);
	K.AddObjectPoint(1,-0.5,-3);
	K.AddObjectPoint(7,-0.5,-3);
	K.AddObjectPoint(1,-3.5,-3);
	K.AddObjectPoint(7,-3.5,-3);
	K.AddObjectPoint(1,-0.5,5);
	K.AddObjectPoint(7,-0.5,5);
	K.AddObjectPoint(1,-3.5,5);
	K.AddObjectPoint(7,-3.5,5);
	K.SetDisplayedPoint(94,367,0);
	K.SetDisplayedPoint(289,431,1);
	K.SetDisplayedPoint(121,268,2);
	K.SetDisplayedPoint(318,329,3);
	K.SetDisplayedPoint(184,350,4);
	K.SetDisplayedPoint(356,406,5);
	K.SetDisplayedPoint(209,263,6);
	K.SetDisplayedPoint(382,316,7);
	K.SetDisplayedPoint(374,446,8);
	K.SetDisplayedPoint(573,509,9);
	K.SetDisplayedPoint(401,345,10);
	K.SetDisplayedPoint(604,410,11);
	K.SetDisplayedPoint(431,419,12);
	K.SetDisplayedPoint(607,475,13);
	K.SetDisplayedPoint(454,332,14);
	K.SetDisplayedPoint(632,387,15);
	K.SetDisplayedPoint(130,234,16);
	K.SetDisplayedPoint(327,295,17);
	K.SetDisplayedPoint(159,131,18);
	K.SetDisplayedPoint(357,192,19);
	K.SetDisplayedPoint(218,231,20);
	K.SetDisplayedPoint(392,285,21);
	K.SetDisplayedPoint(242,142,22);
	K.SetDisplayedPoint(418,195,23);
	K.SetDisplayedPoint(441,312,24);
	K.SetDisplayedPoint(613,376,25);
	K.SetDisplayedPoint(440,212,26);
	K.SetDisplayedPoint(644,274,27);
	K.SetDisplayedPoint(463,302,28);
	K.SetDisplayedPoint(642,357,29);
	K.SetDisplayedPoint(489,214,30);
	K.SetDisplayedPoint(669,267,31);
*/


	for (int k=0;k<10;k++)
	{
	
	
	X.SetCR(0.99);
	X.SetF(0.5);
	double Eli[9];
	clock_t t1,t2,t3,t4;
	X.Start_Evoluce();
	int i =0;
	double q = 500;
	t1 = clock();
	while((i++)<4000&&q>0.01)
	{
		X.Reprodukce();
		if (i%10==0)
		{
			q = X.StandartDeviation();
			//printf("\%f\t\%f\n",X.HodnotaElita(),q);
		}
	}
	t2 = clock();
	X.Elita(Eli);
	//printf("f/px= \%f\nf/py= \%f\nalfa= \%f\nbeta= \%f\ngama= \%f\nX= \%f\nY= \%f\nZ= \%f\nOdchylka= \%f\n",Eli[0],Eli[1],Eli[2],Eli[3],Eli[4],Eli[5],Eli[6],Eli[7],Eli[8]);
	printf("%f ms\t %d\t",1000*(double)(t2-t1)/CLOCKS_PER_SEC,i);
	
	
	K.SetCR(0.99);
	K.SetF(0.5);
	K.Start_Evoluce();
	i =0;
	q = 500;
	t3 = clock();
	while((i++)<4000&&q>0.01)
	{
		K.Reprodukce();
		if (i%10==0)
		{
			q = K.StandartDeviation();
			//printf("\%f\t\%f\n",X.HodnotaElita(),q);
		}
	}
	t4 = clock();
	printf("%f ms %d\n",1000*(double)(t4-t3)/CLOCKS_PER_SEC,i);
	}
	

	
	

    return 0;
}

