diff --git a/qsort.h b/qsort.h
index 6f7f078..4f25b75 100644
--- a/qsort.h
+++ b/qsort.h
@@ -1,7 +1,6 @@
 #ifndef _QSORT_H
 #define _QSORT_H
-//general quicksort suitable if we do not want to create extra class for the member type but prefer to encapsulate it into compare and swap soutines
-//returns parity of the permutation
+//quicksort, returns parity of the permutation
 
 template<typename INDEX, typename COMPAR>
 int genqsort(INDEX l, INDEX r,COMPAR (*cmp)(const INDEX, const INDEX), void (*swap)(const INDEX,const INDEX)) 
@@ -12,7 +11,7 @@ int parity=0;
 if(r<=l) return parity; //1 element
 if(cmp(r,l)<0) {parity^=1; swap(l,r);}
 if(r-l==1) return parity; //2 elements and preparation for median
-piv= (l+r)/2; //pivoting by median of 3 - safer 
+piv= l+(r-l)/2; //pivoting by median of 3 - safer 
 if(cmp(piv,l)<0) {parity^=1; swap(l,piv);} //and change the pivot element implicitly
 if(cmp(r,piv)<0) {parity^=1; swap(r,piv);} //and change the pivot element implicitly
 if(r-l==2) return parity; //in the case of 3 elements we are finished too
@@ -43,19 +42,21 @@ return parity;
 }
 
 
+
 //for SORTABLE classes which provide LA_sort_traits<SORTABLE,INDEX,type>::compare and swap member functions
-template<int type, typename SORTABLE, typename INDEX>
-int memqsort(SORTABLE &object, INDEX l, INDEX r)
+//this allows to use it in general templates also for complex elements, for which comparison falls back to error
+template<int type, typename SORTABLE, typename INDEX, typename PERMINDEX>
+int memqsort(SORTABLE &object, PERMINDEX *perm, INDEX l, INDEX r)
 {
 INDEX i,j,piv;
 int parity=0;
 
 if(r<=l) return parity; //1 element
-if(LA_sort_traits<SORTABLE,INDEX,type>::compare(object,l,r)) {parity^=1; object.swap(l,r);}
+if(LA_sort_traits<SORTABLE,INDEX,type>::compare(object,l,r)) {parity^=1; object.swap(l,r); 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)/2; //pivoting by median of 3 - safer 
-if(LA_sort_traits<SORTABLE,INDEX,type>::compare(object,l,piv)) {parity^=1; object.swap(l,piv);} //and change the pivot element implicitly
-if(LA_sort_traits<SORTABLE,INDEX,type>::compare(object,piv,r)) {parity^=1; object.swap(r,piv);} //and change the pivot element implicitly
+piv= l+(r-l)/2; //pivoting by median of 3 - safer 
+if(LA_sort_traits<SORTABLE,INDEX,type>::compare(object,l,piv)) {parity^=1; object.swap(l,piv); 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>::compare(object,piv,r)) {parity^=1; object.swap(r,piv); 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
@@ -70,16 +71,96 @@ do{
   if(i<j)
 	{
 	// swap and keep track of position of pivoting element 
-	parity^=1; object.swap(i,j); 
+	parity^=1; object.swap(i,j);  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 ^=memqsort<type,SORTABLE,INDEX>(object,l,j); if(i<r) parity ^=memqsort<type,SORTABLE,INDEX>(object,i,r);}
+	{if(l<j) parity ^=memqsort<type,SORTABLE,INDEX,PERMINDEX>(object,perm,l,j); if(i<r) parity ^=memqsort<type,SORTABLE,INDEX,PERMINDEX>(object,perm,i,r);}
 else
-	{if(i<r) parity ^=memqsort<type,SORTABLE,INDEX>(object,i,r); if(l<j) parity ^=memqsort<type,SORTABLE,INDEX>(object,l,j);}
+	{if(i<r) parity ^=memqsort<type,SORTABLE,INDEX,PERMINDEX>(object,perm,i,r); if(l<j) parity ^=memqsort<type,SORTABLE,INDEX,PERMINDEX>(object,perm,l,j);}
+return parity;
+}
+
+
+template<typename S, typename PERMINDEX>
+int ptrqsortup(S *l, S *r, PERMINDEX *perm=NULL)
+{
+S *i,*j,*piv;
+int parity=0;
+
+if(r-l<=0) return parity; //1 element
+if(*l > *r) {parity^=1; {S tmp; tmp=*l; *l= *r; *r=tmp;} if(perm) {PERMINDEX tmp=*perm; *perm=perm[r-l]; perm[r-l]=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(*l>*piv) {parity^=1; {S tmp; tmp=*l; *l=*piv; *piv=tmp;} if(perm) {PERMINDEX tmp= *perm; *perm=perm[piv-l]; perm[piv-l]=tmp;}} //and change the pivot element implicitly
+if(*piv>*r) {parity^=1; {S tmp; tmp=*r; *r=*piv; *piv=tmp;} if(perm) {PERMINDEX tmp=perm[r-l]; perm[r-l]=perm[piv-l]; perm[piv-l]=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(*piv > *i++);
+  i--;
+  while(*j-- > *piv);
+  j++;
+  if(i<j)
+	{
+	// swap and keep track of position of pivoting element 
+	parity^=1; {S tmp; tmp=*i; *i=*j; *j=tmp;} if(perm) {PERMINDEX tmp=perm[i-l]; perm[i-l]=perm[j-l]; perm[j-l]=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 ^=ptrqsortup(l,j,perm); if(i<r) parity ^=ptrqsortup(i,r,perm+(i-l));}
+else
+	{if(i<r) parity ^=ptrqsortup(i,r,perm+(i-l)); if(l<j) parity ^=ptrqsortup(l,j,perm);}
+return parity;
+}
+
+
+template<typename S, typename PERMINDEX>
+int ptrqsortdown(S *l, S *r, PERMINDEX *perm=NULL)
+{
+S *i,*j,*piv;
+int parity=0;
+
+if(r-l<=0) return parity; //1 element
+if(*l < *r) {parity^=1; {S tmp; tmp=*l; *l= *r; *r=tmp;} if(perm) {PERMINDEX tmp=*perm; *perm=perm[r-l]; perm[r-l]=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(*l<*piv) {parity^=1; {S tmp; tmp=*l; *l=*piv; *piv=tmp;} if(perm) {PERMINDEX tmp= *perm; *perm=perm[piv-l]; perm[piv-l]=tmp;}} //and change the pivot element implicitly
+if(*piv<*r) {parity^=1; {S tmp; tmp=*r; *r=*piv; *piv=tmp;} if(perm) {PERMINDEX tmp=perm[r-l]; perm[r-l]=perm[piv-l]; perm[piv-l]=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(*piv < *i++);
+  i--;
+  while(*j-- < *piv);
+  j++;
+  if(i<j)
+	{
+	// swap and keep track of position of pivoting element 
+	parity^=1; {S tmp; tmp=*i; *i=*j; *j=tmp;} if(perm) {PERMINDEX tmp=perm[i-l]; perm[i-l]=perm[j-l]; perm[j-l]=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 ^=ptrqsortdown(l,j,perm); if(i<r) parity ^=ptrqsortdown(i,r,perm+(i-l));}
+else
+	{if(i<r) parity ^=ptrqsortdown(i,r,perm+(i-l)); if(l<j) parity ^=ptrqsortdown(l,j,perm);}
 return parity;
 }
 
diff --git a/test.cc b/test.cc
index ed9e14d..c2a58d5 100644
--- a/test.cc
+++ b/test.cc
@@ -1,4 +1,6 @@
 #include "bitvector.h"
+#include "qsort.h"
+
 
 int main(void)
 {
@@ -23,6 +25,10 @@ for(int i=0; i<10; ++i) t[i]=i;
 cout <<t;
 t.sort(1);
 cout <<t;
-t.sort(0);
+NRVec<int> u(10);
+for(int i=0; i<10;++i) u[i]=i;
+ptrqsortup(&t[0],&t[9],&u[0]);
+cout<<t <<"U= "<<u;
+ptrqsortdown<int,int>(&t[0],&t[9]);
 cout<<t;
 }
diff --git a/vec.cc b/vec.cc
index df827e0..24a14a1 100644
--- a/vec.cc
+++ b/vec.cc
@@ -494,10 +494,10 @@ NRVec< complex<double> >::operator|(const NRVec< complex<double> > &b) const
 
 
 template<typename T>
-int NRVec<T>::sort(int direction, int from, int to)
+int NRVec<T>::sort(int direction, int from, int to, int *perm)
 {
 copyonwrite();
 if(to == -1) to=nn-1;
-if(direction) return memqsort<1,NRVec<T>,int>(*this,from,to);
-else return memqsort<0,NRVec<T>,int>(*this,from,to);
+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);
 }
diff --git a/vec.h b/vec.h
index e4d38b8..ae72268 100644
--- a/vec.h
+++ b/vec.h
@@ -106,7 +106,7 @@ public:
 	bool bigger(int i, int j) const {return LA_traits<T>::bigger(v[i],v[j]);};
 	bool smaller(int i, int j) const {return LA_traits<T>::smaller(v[i],v[j]);};
 	void swap(int i, int j) {T tmp; tmp=v[i]; v[i]=v[j]; v[j]=tmp;};
-	int sort(int direction=0, int from=0, int to= -1); //sort, ascending by default, returns parity of permutation
+	int sort(int direction=0, int from=0, int to= -1, int *perm=NULL); //sort, ascending by default, returns parity of permutation
 };