jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

implementation of stabilized quicksort

Browse Source
This commit is contained in:
Jiri Pittner 2023-07-30 11:00:41 +02:00
parent ec42999812
commit a439e0be94
4 changed files with 92 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
la_traits.h
View File

@ -145,14 +145,28 @@ template<typename T, typename I, int type> struct LA_sort_traits;
template<typename T, typename I>
struct LA_sort_traits<T,I,0>
{
static inline bool compare(T object, I i, I j) {return object.bigger(i,j);};
static inline bool compare(const T &object, I i, I j) {return object.bigger(i,j);};
static inline bool comparestable(const T &object, I i, I j, I *auxkey)
{
bool r= object.bigger(i,j);
if(!r && object[i]==object[j]) r= auxkey[i]>auxkey[j];
return r;
};
};
template<typename T, typename I>
struct LA_sort_traits<T,I,1>
{
static inline bool compare(T object, I i, I j) {return object.smaller(i,j);};
static inline bool compare(const T &object, I i, I j) {return object.smaller(i,j);};
static inline bool comparestable(const T &object, I i, I j, I *auxkey)
{
bool r= object.smaller(i,j);
if(!r && object[i]==object[j]) r= auxkey[i]<auxkey[j];
return r;
};
};
//we will need to treat char and unsigned char as numbers in << and >> I/O operators

60
qsort.h
View File

@ -104,6 +104,66 @@ return parity;
}
//stabilized version of quicksort employing auxiliary key of the original position
template<int type, typename SORTABLE, typename INDEX, typename PERMINDEX>
int memqsortaux(SORTABLE &object, PERMINDEX *perm, INDEX l, INDEX r, INDEX *aux)
{
INDEX i,j,piv;
int parity=0;
if(r<=l) return parity; //1 element
if(LA_sort_traits<SORTABLE,INDEX,type>::comparestable(object,l,r,aux)) {parity^=1; object.swap(l,r); {INDEX tmp=aux[l]; aux[l]=aux[r]; aux[r]=tmp;} if(perm) {PERMINDEX tmp=perm[l]; perm[l]=perm[r]; perm[r]=tmp;}}
if(r-l==1) return parity; //2 elements and preparation for median
piv= l+(r-l)/2; //pivoting by median of 3 - safer
if(LA_sort_traits<SORTABLE,INDEX,type>::comparestable(object,l,piv,aux)) {parity^=1; object.swap(l,piv); {INDEX tmp=aux[l]; aux[l]=aux[piv]; aux[piv]=tmp;} if(perm) {PERMINDEX tmp=perm[l]; perm[l]=perm[piv]; perm[piv]=tmp;}} //and change the pivot element implicitly
if(LA_sort_traits<SORTABLE,INDEX,type>::comparestable(object,piv,r,aux)) {parity^=1; object.swap(r,piv); {INDEX tmp=aux[r]; aux[r]=aux[piv]; aux[piv]=tmp;} if(perm) {PERMINDEX tmp=perm[r]; perm[r]=perm[piv]; perm[piv]=tmp;}} //and change the pivot element implicitly
if(r-l==2) return parity; //in the case of 3 elements we are finished too
//general case , l-th r-th already processed
i=l+1; j=r-1;
do{
//important sharp inequality - stops at sentinel element for efficiency
// this is inefficient if all keys are equal - unnecessary n log n swaps are done, but we assume that it is atypical input
while(LA_sort_traits<SORTABLE,INDEX,type>::comparestable(object,piv,i++,aux));
i--;
while(LA_sort_traits<SORTABLE,INDEX,type>::comparestable(object,j--,piv,aux));
j++;
if(i<j)
{
// swap and keep track of position of pivoting element
parity^=1;
object.swap(i,j);
{INDEX tmp =aux[i]; aux[i]=aux[j]; aux[j]=tmp;}
if(perm) {PERMINDEX tmp=perm[i]; perm[i]=perm[j]; perm[j]=tmp;}
if(i==piv) piv=j; else if(j==piv) piv=i;
}
if(i<=j) {i++; j--;}
}while(i<=j);
if(j-l < r-i) //because of the stack in bad case process first the shorter subarray
{if(l<j) parity ^=memqsortaux<type,SORTABLE,INDEX,PERMINDEX>(object,perm,l,j,aux); if(i<r) parity ^=memqsortaux<type,SORTABLE,INDEX,PERMINDEX>(object,perm,i,r,aux);}
else
{if(i<r) parity ^=memqsortaux<type,SORTABLE,INDEX,PERMINDEX>(object,perm,i,r,aux); if(l<j) parity ^=memqsortaux<type,SORTABLE,INDEX,PERMINDEX>(object,perm,l,j,aux);}
return parity;
}
template<int type, typename SORTABLE, typename INDEX, typename PERMINDEX>
int memqsortstable(SORTABLE &object, PERMINDEX *perm, INDEX l, INDEX r)
{
if(r<=l) return 0; //1 element
INDEX *aux= new INDEX[r-l+1];
if(aux==NULL) laerror("allocation failed in memqsortstable");
INDEX *auxshifted = aux-l;
for(INDEX i=l; i<=r; ++i) auxshifted[i]=i;
int parity=memqsortaux<type,SORTABLE,INDEX,PERMINDEX>(object,perm,l,r,auxshifted);
delete[] aux;
return parity;
}
template<typename S, typename PERMINDEX>
int ptrqsortup(S *l, S *r, PERMINDEX *perm=NULL)
{

2
t.cc
View File

@ -2758,7 +2758,7 @@ cout <<"Check adjacency of the clique: "<<adj.submatrix(clique)<<endl;
NRVec<int> cover = cliquecover(adj);
cout <<"Clique cover is "<<cover<<endl;
NRPerm<int> p(cover.size());
cover.sort(0,p);
cover.sort(0,p,true);
cout<<"permutation to disentabgle the cliques = "<<p<<endl;
NRSMat<char> adjperm = adj.permuted(p);
cout <<"resorted graph = "<<adjperm<<endl;

18
vec.h
View File

@ -451,8 +451,8 @@ public:
};
//! sort by default in ascending order and return the parity of corresponding permutation resulting to this order
int sort(int direction = 0, int from = 0, int to = -1, int *perm = NULL);
int sort(int direction, NRPerm<int> &perm);
int sort(int direction = 0, int from = 0, int to = -1, int *perm = NULL, bool stable=false);
int sort(int direction, NRPerm<int> &perm, bool stable=false);
//! apply given function to each element
NRVec& call_on_me(T (*_F)(const T &) ){
@ -518,21 +518,29 @@ namespace LA {
template<typename T>
int NRVec<T>::sort(int direction, int from, int to, int *perm) {
int NRVec<T>::sort(int direction, int from, int to, int *perm, bool stable) {
NOT_GPU(*this);
copyonwrite();
if(to == -1) to = nn - 1;
if(stable)
{
if(direction) return memqsortstable<1, NRVec<T>, int, int>(*this, perm, from, to);
else return memqsortstable<0, NRVec<T>, int, int>(*this, perm, from, to);
}
else
{
if(direction) return memqsort<1, NRVec<T>, int, int>(*this, perm, from, to);
else return memqsort<0, NRVec<T>, int, int>(*this, perm, from, to);
}
}
template<typename T>
int NRVec<T>::sort(int direction, NRPerm<int> &perm)
int NRVec<T>::sort(int direction, NRPerm<int> &perm, bool stable)
{
if(nn!=perm.size()) laerror("incompatible vector and permutation");
perm.identity();
int r=sort(direction,0,nn-1,&perm[1]);
int r=sort(direction,0,nn-1,&perm[1],stable);
return r;
}