#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>

#include "nbtrixx.h"

namespace SOPHYA {

// attention, Numerical recipes prend des tableaux de 1 a n on remet
// de 0 a n-1 en decramentant le pointeur du tableau d'entree
 
//-------------------------------------------------------------
#define SWAP_NBTRI(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50
template <class T>
void TabSort(uint_4 n, T* arr_c)
{
	uint_4 i,ir=n,j,k,l=1,*istack;
	int jstack=0;
	T a,temp;
        T *arr = arr_c - 1; // MODIF

	//istack=lvector(1,NSTACK);
        istack = new uint_4[NSTACK+1]; // MODIF
	for (;;) {
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				a=arr[j];
				for (i=j-1;i>=l;i--) {
					if (arr[i] <= a) break;
					arr[i+1]=arr[i];
				}
				arr[i+1]=a;
			}
			if (jstack == 0) break;
			ir=istack[jstack--];
			l=istack[jstack--];
		} else {
			k=(l+ir) >> 1;
			SWAP_NBTRI(arr[k],arr[l+1])
			if (arr[l] > arr[ir]) {
				SWAP_NBTRI(arr[l],arr[ir])
			}
			if (arr[l+1] > arr[ir]) {
				SWAP_NBTRI(arr[l+1],arr[ir])
			}
			if (arr[l] > arr[l+1]) {
				SWAP_NBTRI(arr[l],arr[l+1])
			}
			i=l+1;
			j=ir;
			a=arr[l+1];
			for (;;) {
				do i++; while (arr[i] < a);
				do j--; while (arr[j] > a);
				if (j < i) break;
				SWAP_NBTRI(arr[i],arr[j]);
			}
			arr[l+1]=arr[j];
			arr[j]=a;
			jstack += 2;
			if (jstack > NSTACK) throw("NSTACK too small in sort.");
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
	//free_lvector(istack,1,NSTACK);
        delete [] istack;  // MODIF
}
#undef M
#undef NSTACK
#undef SWAP_NBTRI

//-------------------------------------------------------------
#define SWAP_NBTRI(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50
template <class T>
void TabSort2(uint_4 n, T *arr_c, T *brr_c)
{
	uint_4 i,ir=n,j,k,l=1,*istack;
	int jstack=0;
	T a,b,temp;
        T *arr = arr_c - 1; // MODIF
        T *brr = brr_c - 1; // MODIF

	//istack=lvector(1,NSTACK);
        istack = new uint_4[NSTACK+1]; // MODIF
	for (;;) {
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				a=arr[j];
				b=brr[j];
				for (i=j-1;i>=l;i--) {
					if (arr[i] <= a) break;
					arr[i+1]=arr[i];
					brr[i+1]=brr[i];
				}
				arr[i+1]=a;
				brr[i+1]=b;
			}
			if (!jstack) {
			        //free_lvector(istack,1,NSTACK);
			        delete [] istack;  // MODIF
				return;
			}
			ir=istack[jstack];
			l=istack[jstack-1];
			jstack -= 2;
		} else {
			k=(l+ir) >> 1;
			SWAP_NBTRI(arr[k],arr[l+1])
			SWAP_NBTRI(brr[k],brr[l+1])
			if (arr[l] > arr[ir]) {
				SWAP_NBTRI(arr[l],arr[ir])
				SWAP_NBTRI(brr[l],brr[ir])
			}
			if (arr[l+1] > arr[ir]) {
				SWAP_NBTRI(arr[l+1],arr[ir])
				SWAP_NBTRI(brr[l+1],brr[ir])
			}
			if (arr[l] > arr[l+1]) {
				SWAP_NBTRI(arr[l],arr[l+1])
				SWAP_NBTRI(brr[l],brr[l+1])
			}
			i=l+1;
			j=ir;
			a=arr[l+1];
			b=brr[l+1];
			for (;;) {
				do i++; while (arr[i] < a);
				do j--; while (arr[j] > a);
				if (j < i) break;
				SWAP_NBTRI(arr[i],arr[j])
				SWAP_NBTRI(brr[i],brr[j])
			}
			arr[l+1]=arr[j];
			arr[j]=a;
			brr[l+1]=brr[j];
			brr[j]=b;
			jstack += 2;
			if (jstack > NSTACK) throw("NSTACK too small in sort2.");
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
}
#undef M
#undef NSTACK
#undef SWAP_NBTRI

//-------------------------------------------------------------
#define SWAP_NBTRI(a,b) itemp=(a);(a)=(b);(b)=itemp;
#define M 7
#define NSTACK 50
template <class T>
void TabSortInd(uint_4 n, T *arr_c, uint_4 *indx_c)
{
	uint_4 i,indxt,ir=n,itemp,j,k,l=1;
	int jstack=0,*istack;
	T a;
        T *arr = arr_c - 1; // MODIF
        uint_4 *indx = indx_c - 1; // MODIF

	//istack=ivector(1,NSTACK);
        istack = new int[NSTACK+1]; // MODIF
	for (j=1;j<=n;j++) indx[j]=j;
	for (;;) {
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				indxt=indx[j];
				a=arr[indxt];
				for (i=j-1;i>=l;i--) {
					if (arr[indx[i]] <= a) break;
					indx[i+1]=indx[i];
				}
				indx[i+1]=indxt;
			}
			if (jstack == 0) break;
			ir=istack[jstack--];
			l=istack[jstack--];
		} else {
			k=(l+ir) >> 1;
			SWAP_NBTRI(indx[k],indx[l+1]);
			if (arr[indx[l]] > arr[indx[ir]]) {
				SWAP_NBTRI(indx[l],indx[ir])
			}
			if (arr[indx[l+1]] > arr[indx[ir]]) {
				SWAP_NBTRI(indx[l+1],indx[ir])
			}
			if (arr[indx[l]] > arr[indx[l+1]]) {
				SWAP_NBTRI(indx[l],indx[l+1])
			}
			i=l+1;
			j=ir;
			indxt=indx[l+1];
			a=arr[indxt];
			for (;;) {
				do i++; while (arr[indx[i]] < a);
				do j--; while (arr[indx[j]] > a);
				if (j < i) break;
				SWAP_NBTRI(indx[i],indx[j])
			}
			indx[l+1]=indx[j];
			indx[j]=indxt;
			jstack += 2;
			if (jstack > NSTACK) throw("NSTACK too small in indexx.");
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
	//free_ivector(istack,1,NSTACK);
        delete [] istack;  // MODIF
        /* pour avoir un retour d'indice entre 0 et n-1 */
	for(j=0;j<n;j++) indx_c[j]--;  // MODIF
}
#undef M
#undef NSTACK
#undef SWAP_NBTRI

//-------------------------------------------------------------
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template TabSort<int_2>
#pragma define_template TabSort<uint_2>
#pragma define_template TabSort<int_4>
#pragma define_template TabSort<uint_4>
#pragma define_template TabSort<int_8>
#pragma define_template TabSort<uint_8>
#pragma define_template TabSort<r_4>
#pragma define_template TabSort<r_8>

#pragma define_template TabSort2<int_2>
#pragma define_template TabSort2<uint_2>
#pragma define_template TabSort2<int_4>
#pragma define_template TabSort2<uint_4>
#pragma define_template TabSort2<int_8>
#pragma define_template TabSort2<uint_8>
#pragma define_template TabSort2<r_4>
#pragma define_template TabSort2<r_8>

#pragma define_template TabSortInd<int_2>
#pragma define_template TabSortInd<uint_2>
#pragma define_template TabSortInd<int_4>
#pragma define_template TabSortInd<uint_4>
#pragma define_template TabSortInd<int_8>
#pragma define_template TabSortInd<uint_8>
#pragma define_template TabSortInd<r_4>
#pragma define_template TabSortInd<r_8>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template void TabSort(uint_4, int_2*);
template void TabSort(uint_4, uint_2*);
template void TabSort(uint_4, int_4*);
template void TabSort(uint_4, uint_4*);
template void TabSort(uint_4, int_8*);
template void TabSort(uint_4, uint_8*);
template void TabSort(uint_4, r_4*);
template void TabSort(uint_4, r_8*);

template void TabSort2(uint_4, int_2*, int_2*);
template void TabSort2(uint_4, uint_2*, uint_2*);
template void TabSort2(uint_4, int_4*, int_4*);
template void TabSort2(uint_4, uint_4*, uint_4*);
template void TabSort2(uint_4, int_8*, int_8*);
template void TabSort2(uint_4, uint_8*, uint_8*);
template void TabSort2(uint_4, r_4*, r_4*);
template void TabSort2(uint_4, r_8*, r_8*);

template void TabSortInd(uint_4, int_2*, uint_4*);
template void TabSortInd(uint_4, uint_2*, uint_4*);
template void TabSortInd(uint_4, int_4*, uint_4*);
template void TabSortInd(uint_4, uint_4*, uint_4*);
template void TabSortInd(uint_4, int_8*, uint_4*);
template void TabSortInd(uint_4, uint_8*, uint_4*);
template void TabSortInd(uint_4, r_4*, uint_4*);
template void TabSortInd(uint_4, r_8*, uint_4*);
#endif

} // FIN namespace SOPHYA