diff --git a/permutation.cc b/permutation.cc
index 03600cb..4978202 100644
--- a/permutation.cc
+++ b/permutation.cc
@@ -21,6 +21,7 @@
 #include <stdio.h>
 #include <string.h>
 #include <list>
+#include "qsort.h"
 
 
 namespace LA {
@@ -246,6 +247,39 @@ end:
 return select? sumperm/2:sumperm;
 }
 
+//Algorithm L2 from Knuth's volume 4 for multiset permutations
+//input must be initialized with (repeated) numbers in any order
+template <typename T>
+PERM_RANK_TYPE NRPerm<T>::generate_all_multi(void (*callback)(const NRPerm<T>&))
+{
+PERM_RANK_TYPE sumperm=0;
+this->copyonwrite();
+myheapsort(this->size(),&(*this)[1],1);
+while(1)
+        {
+        int j,k,l;
+        T t;
+
+        sumperm++;
+        (*callback)(*this);
+        j=this->size()-1;
+        while(j>0 && (*this)[j]>=(*this)[j+1]) j--;
+        if(j==0) break; /*finished*/
+        l=this->size();
+        while((*this)[j]>=(*this)[l]) l--;
+        t=(*this)[j];(*this)[j]=(*this)[l];(*this)[l]=t;
+        k=j+1;
+        l=this->size();
+        while(k<l)
+                {
+                t=(*this)[k];(*this)[k]=(*this)[l];(*this)[l]=t;
+                k++;l--;
+                }
+        }
+return sumperm;
+}
+
+
 
 template <typename T>
 static T _n;
@@ -253,13 +287,13 @@ static T _n;
 template <typename T>
 static void (*_callback)(const NRPerm<T>& ); 
 
-PERM_RANK_TYPE _sumperm;
+static PERM_RANK_TYPE _sumperm;
 
 template <typename T>
-NRPerm<T> *_perm;
+static NRPerm<T> *_perm;
 
 template <typename T>
-void permg(T n)
+static void permg(T n)
 {
 if(n<= _n<T>)
         {
@@ -307,6 +341,81 @@ do{
 return np;
 }
 
+template <typename T>
+static T _n2;
+
+template <typename T>
+static void (*_callback2)(const NRPerm<T>& );
+
+static PERM_RANK_TYPE _sumperm2;
+
+template <typename T>
+static NRPerm<T> *_perm2;
+
+template <typename T>
+static const T *pclasses2;
+
+template <typename T>
+static T *_pperm;
+
+static int sameclass;
+
+template <typename T>
+static void permg2(T n) //place number n in a free box in all possible ways and according to restrictions
+{
+T i;
+if(n<=_n2<T>)
+        {
+        T c=pclasses2<T>[n];
+        for(i=1;i<=_n2<T>;i++)
+                {
+                if(_pperm<T>[i]>=1) goto skipme; /*place already occupied*/
+                if (sameclass)
+                        {
+                        if(sameclass==1 && c!=pclasses2<T>[i]) goto skipme;
+                        if(sameclass== -1 && c==pclasses2<T>[i] ) goto skipme;
+                        if(sameclass== -2) /*allow only permutations which leave elements of each class sorted*/
+                                {
+                                T j,k;
+                                for(j=i+1; j<=_n2<T>; j++)
+                                        if((k=_pperm<T>[j])>=1)
+                                                 if(/* automatically fulfilled k<n &&*/ c==pclasses2<T>[k]) goto skipme;
+                                }
+                        }
+                /*put it there and next number*/
+                _pperm<T>[i]=n;
+                permg2(n+1);
+                _pperm<T>[i]=0;
+                skipme:;
+                }
+        }
+else
+        {
+        _sumperm2++;
+	(*_callback2<T>)(*_perm2<T>);
+        }
+}
+
+
+
+template <typename T>
+PERM_RANK_TYPE NRPerm<T>::generate_restricted(void (*callback)(const NRPerm<T>&), const NRVec_from1<T> &classes, int restriction_type)
+{
+this->copyonwrite();
+this->clear();
+_n2<T> = this->size();
+_callback2<T> =callback;
+_sumperm2=0;
+_perm2<T> = this;
+_pperm<T> = &(*this)[1]-1;
+pclasses2<T> = &classes[1]-1;
+sameclass=restriction_type;
+permg2<T>(1);
+return _sumperm2;
+}
+
+
+
 template <typename T>
 PERM_RANK_TYPE NRPerm<T>::rank() const
 {
diff --git a/permutation.h b/permutation.h
index 5b6583b..ab76be9 100644
--- a/permutation.h
+++ b/permutation.h
@@ -58,12 +58,14 @@ public:
 	NRPerm operator*(const NRPerm q) const; //q is rhs and applied first, this applied second
 	NRPerm operator*(const CyclePerm<T> &r) const;
 	NRPerm conjugate_by(const NRPerm q) const; //q^-1 p q
-	int parity() const;
+	int parity() const; //returns +/- 1
 	void randomize(void); //uniformly random by Fisher-Yates shuffle
 	bool next(); //generate next permutation in lex order
-	PERM_RANK_TYPE generate_all(void (*callback)(const NRPerm<T>&), int parity_select=0); //Algorithm from Knuth's vol.4, efficient but not in lex order!
+	PERM_RANK_TYPE generate_all(void (*callback)(const NRPerm<T>&), int parity_select=0); //Algorithm L from Knuth's vol.4, efficient but not in lex order!
+	PERM_RANK_TYPE generate_all_multi(void (*callback)(const NRPerm<T>&));  //Algorithm L2 from Knuth's vol.4, for multiset (repeated numbers, not really permutations)
 	PERM_RANK_TYPE generate_all2(void (*callback)(const NRPerm<T>&)); //recursive method, also not lexicographic
 	PERM_RANK_TYPE generate_all_lex(void (*callback)(const NRPerm<T>&)); //generate in lex order using next()
+	PERM_RANK_TYPE generate_restricted(void (*callback)(const NRPerm<T>&), const NRVec_from1<T> &classes, int restriction_type=0);
 	PERM_RANK_TYPE rank() const; //counted from 0 to n!-1
 	NRVec_from1<T> inversions(const int type, PERM_RANK_TYPE *prank=NULL) const; //inversion tables
 	explicit NRPerm(const int type, const NRVec_from1<T> &inversions); //compute permutation from inversions
diff --git a/qsort.h b/qsort.h
index b42f914..da91a4f 100644
--- a/qsort.h
+++ b/qsort.h
@@ -243,5 +243,43 @@ else
 return parity;
 }
 
+//heapsort
+//DATA must have assignmend and comparison operators implemented
+template<typename INDEX, typename DATA>
+void myheapsort(INDEX n, DATA *ra,int typ=1)
+{
+        INDEX l,j,ir,i;
+        DATA rra;
+	ra--; //below indexed from 1
+
+        l=(n >> 1)+1;
+        ir=n;
+        for (;;) {
+                if (l > 1) {
+                        rra=ra[--l];
+                } else {
+                        rra=ra[ir];
+                        ra[ir]=ra[1];
+                        if (--ir == 1) {
+                                ra[1]=rra;
+                                return;
+                        }
+                }
+                i=l;
+                j=l << 1;
+                while (j <= ir) {
+                        if (j < ir && (typ>0?ra[j] < ra[j+1]:ra[j]>ra[j+1])) ++j;
+                        if (typ>0?rra < ra[j]:rra>ra[j]) {
+                                ra[i]=ra[j];
+                                j += (i=j);
+                        }
+                        else j=ir+1;
+                }
+                ra[i]=rra;
+        }
+}
+
+
+
 }//namespace
 #endif
diff --git a/t.cc b/t.cc
index a93a2a5..cd8ac55 100644
--- a/t.cc
+++ b/t.cc
@@ -75,6 +75,19 @@ for(int i=0; i<4; ++i)
 	}
 }
 
+void printme0(const NRPerm<int> &p)
+{
+cout<<p;
+}
+
+void printme1(const NRPerm<int> &p)
+{
+cout <<p.parity()<<' ';
+cout<<p;
+//cout<<p.inverse();
+}
+
+
 static int unitary_n;
 static PERM_RANK_TYPE space_dim;
 
@@ -2187,6 +2200,33 @@ int tot=p.generate_all_lex(printme);
 cout <<"generated "<<tot<<endl;
 }
 
+if(0)
+{
+int n;
+cin >>n;
+NRPerm<int> p(n);
+int seed;
+int f=open("/dev/random",O_RDONLY);
+if(sizeof(int)!=read(f,&seed,sizeof(int))) laerror("cannot read /dev/random");
+close(f);
+srand(seed);
+for(int i=1; i<=n; ++i) p[i]=1+RANDINT32()%(n-1);
+cout <<"Initial = "<<p<<endl;
+int tot=p.generate_all_multi(printme0);
+cout <<"generated "<<tot<<endl;
+}
+
+if(1)
+{
+int n; cin >>n;
+NRPerm<int> p(n); 
+NRVec_from1<int> classes; cin>>classes;
+if(classes.size()!=p.size()) laerror("sizes do not match");
+int tot=p.generate_restricted(printme1,classes,-2);
+cout <<"generated "<<tot<<endl;
+}
+
+
 if(0)
 {
 int n,unitary_n;;
@@ -2969,7 +3009,7 @@ cout <<c%ir<<endl;
 cout <<cc<<endl;
 }
 
-if(1)
+if(0)
 {
 int seed;
 int f=open("/dev/random",O_RDONLY);