#include "MUtils.h"


void Memms() { // ==> should be done after MMX and before FPU calculations <==
#ifdef ARCH_X86
  emms();
#endif
}


int MExists(const char* FileName) {
  FILE* File;
  if ((File=fopen(FileName,"r"))== NULL) return 0;
  fclose(File);
  return 1;
}
const char* MExistingFileName(const char* PrefixSuffix, int Number) {
  static char FileName[1024];
  if (PrefixSuffix== NULL || !strlen(PrefixSuffix))
    return NULL;
  int i, iMax= strlen(PrefixSuffix);
  for (i=0; i< iMax; i++)
    if (PrefixSuffix[i]=='*')
      break;
  if (i== iMax)
    return NULL;
  sprintf(FileName,"%s", PrefixSuffix);
  sprintf(FileName+i,"%d%s\0", int(Number), PrefixSuffix+i+1);
  if (MExists(FileName)) return FileName;
  sprintf(FileName+i,"0%d%s\0", int(Number), PrefixSuffix+i+1);
  if (MExists(FileName)) return FileName;
  sprintf(FileName+i,"00%d%s\0", int(Number), PrefixSuffix+i+1);
  if (MExists(FileName)) return FileName;
  sprintf(FileName+i,"000%d%s\0", int(Number), PrefixSuffix+i+1);
  if (MExists(FileName)) return FileName;
  sprintf(FileName+i,"0000%d%s\0", int(Number), PrefixSuffix+i+1);
  if (MExists(FileName)) return FileName;
  return NULL;
}


#include <sys/time.h>

double MTime() {
  struct timeval tv;
  gettimeofday(&tv, (struct timezone *)NULL);
  return tv.tv_sec + tv.tv_usec/1000000.0;
}

double MTime2()
{
	struct tms buffer;
	(void) times(&buffer);
	double NbClkTck=sysconf(_SC_CLK_TCK);
	return double(buffer.tms_utime/NbClkTck);
}





int MCpuid() { // This function is from the GNU mpeg2dec source.
  // Copyright (C) 1999-2001 Aaron Holtzman <aholtzma@ess.engr.uvic.ca>
#ifdef ARCH_X86
  uint32_t eax, ebx, ecx, edx;
  int AMD;
  uint32_t caps;

#define cpuid(op,eax,ebx,ecx,edx)       \
    asm ("cpuid"                        \
         : "=a" (eax),                  \
           "=b" (ebx),                  \
           "=c" (ecx),                  \
           "=d" (edx)                   \
         : "a" (op)                     \
         : "cc")

  asm ("pushfl\n\t"
       "popl %0\n\t"
       "movl %0,%1\n\t"
       "xorl $0x200000,%0\n\t"
       "pushl %0\n\t"
       "popfl\n\t"
       "pushfl\n\t"
       "popl %0"
       : "=a" (eax),
       "=b" (ebx)
       :
       : "cc");

  if (eax == ebx)             /* no cpuid */
    return 0;
  cpuid (0x00000000, eax, ebx, ecx, edx);
  if (!eax)                   /* vendor string only */
    return 0;
  AMD = (ebx == 0x68747541) && (ecx == 0x444d4163) && (edx == 0x69746e65);
  cpuid (0x00000001, eax, ebx, ecx, edx);
  if (! (edx & 0x00800000))   /* no MMX */
    return 0;
  caps = M_MM_ACCEL_X86_MMX;
  if (edx & 0x02000000)       /* SSE - identical to AMD MMX extensions */
    caps = M_MM_ACCEL_X86_MMX | M_MM_ACCEL_X86_MMXEXT;
  cpuid (0x80000000, eax, ebx, ecx, edx);
  if (eax < 0x80000001)       /* no extended capabilities */
    return caps;
  cpuid (0x80000001, eax, ebx, ecx, edx);
  if (edx & 0x80000000)
    caps |= M_MM_ACCEL_X86_3DNOW;
  if (AMD && (edx & 0x00400000))      /* AMD MMX extensions */
    caps |= M_MM_ACCEL_X86_MMXEXT;
  return caps;
#endif
  return 0;
};

void MCheckCpu() {
  assert(sizeof(char)== 1);
  assert(sizeof(short)== 2);
  assert(sizeof(int)== 4);
  assert(sizeof(float)== 4);
  assert(sizeof(double)== 8);
  assert(sizeof(long long)== 8);
  assert((MCpuid()&(M_MM_ACCEL_X86_MMX|M_MM_ACCEL_X86_MMXEXT))==
	 (M_MM_ACCEL_X86_MMX|M_MM_ACCEL_X86_MMXEXT));
  WithMmxSseUse= 1;
}
int MWithMmxSseUse() { return WithMmxSseUse; };
void MSetWithMmxSseUse(int w) {WithMmxSseUse=w;}

//---------------------------------------------------------------------------------------------
// Fonctions stats de base
//---------------------------------------------------------------------------------------------
float CalcMoyenne(float *Tab,int n)
{
	float moy=0;
	for (int i=0;i<n;i++) moy+=Tab[i];
	return (moy/n);
}
float CalcEcartType(float *Tab,int n)
{
	float moy=CalcMoyenne(Tab,n);
	return CalcEcartType(Tab,moy,n);
}
float CalcEcartType(float *Tab,float Moy,int n)
{
	float ec=0;
	for (int i=0;i<n;i++) ec+=(Tab[i]-Moy)*(Tab[i]-Moy);
	ec=ec/n;
	return sqrt(ec);	
}
float CalcVariance(float *Tab,int n)
{
	float moy=CalcMoyenne(Tab,n);
	return CalcVariance(Tab,moy,n);
}
float CalcVariance(float *Tab,float Moy,int n)
{
	float ec=0;
	for (int i=0;i<n;i++) ec+=(Tab[i]-Moy)*(Tab[i]-Moy);
	ec=ec/n;
	return ec;
}

//Returns the kth smallest value in the array arr[1..n].
//The input array will be rearranged to have this value in location arr[k]
//with all smaller elements moved to arr[1..k-1] (in arbitrary order) and all larger elements in arr[k+1..n]
// (also in arbitrary order).
// (numerical recipies)
static float MedianSelect(unsigned long k, unsigned long n, float arr[])
{
	unsigned long i,ir,j,l,mid;
	float a,temp;
	l=1;
	ir=n;
	for (;;)
	{
		if (ir <= l+1)
		{
			//Active partition contains 1 or 2e lements.
			if (ir == l+1 && arr[ir] < arr[l]);
			{
				//Case of 2el ements.
				Swap(arr[l],arr[ir]);
			}
			return arr[k];
		}
		else
		{
			mid=(l+ir) >> 1;
			//Choose median of left, center, and right elements as partitioning element a. Also rearrange so that arr[l] d arr[l+1], arr[ir] e arr[l+1].
			Swap(arr[mid],arr[l+1]);
			if (arr[l] > arr[ir])
			{
				Swap(arr[l],arr[ir]);
			}
			if (arr[l+1] > arr[ir])
			{
				Swap(arr[l+1],arr[ir]);
			}
			if (arr[l] > arr[l+1])
			{
				Swap(arr[l],arr[l+1]);
			}
			i=l+1;
			// Initialize pointers for partitioning.
			j=ir;
			a=arr[l+1];
			//Partitioning element.
			for (;;)
			{
				//Beginning of innermost loop.
				do i++;while (arr[i] < a);
				//Scan up to  nd element > a.
				do j--;while (arr[j] > a);
				//Scan down to  nd element < a.
				if (j < i) break;
				//Pointers crossed. Partitioning complete.
				Swap(arr[i],arr[j]);
			} //End of innermost loop.
			arr[l+1]=arr[j]; //Insert partitioning element.
			arr[j]=a;
			if (j >= k) ir=j-1; //Keep active the partition that contains the kth element.
			if (j <= k) l=i;
		}
	}
}

float CalcMediane(float *Tab,int n)
{
  if (n%2!=0) n=n-1;
	float *Tab2=new float[n];
  for (int i=0;i<n;i++) Tab2[i]=Tab[i];
  int demi=n/2;
  float m=MedianSelect(demi,n,Tab2);  
	delete [] Tab2;
	return m;
 /* int i;
	list<float> liste;
	for (i=0;i<n;i++) liste.push_back(Tab[i]);
	liste.sort();
  i=0;
  list<float>::iterator it=liste.begin();
	while (i<n/2)
	{
		i++;
		++it;
	}
	return (*it); */
}