diff --git a/permutation.cc b/permutation.cc
index bf8889a..386ac32 100644
--- a/permutation.cc
+++ b/permutation.cc
@@ -161,6 +161,316 @@ for(int i=n-1; i>=1; --i)
 }
 
 
+template <typename T>
+bool NRPerm<T>::next(void)
+{
+this->copyonwrite();
+int n=this->size();
+// Find longest non-increasing suffix and the pivot
+int i = n;
+while (i > 1 && (*this)[i-1] >= (*this)[i]) --i;
+if (i<=1) return false;
+T piv=(*this)[i-1];
+
+// Find rightmost element greater than the pivot
+int j = n;
+while ((*this)[j] <= piv) --j;
+// Now the value array[j] will become the new pivot, Assertion: j >= i
+    
+// Swap the pivot with j
+(*this)[i - 1] = (*this)[j];
+(*this)[j] = piv;
+    
+// Reverse the suffix
+j = n;
+while (i < j) {
+        T temp = (*this)[i];
+        (*this)[i] = (*this)[j];
+        (*this)[j] = temp;
+        i++;
+        j--;
+    }
+
+return true;
+}
+
+//Algorithm L from Knuth's volume 4
+template <typename T>
+PERM_RANK_TYPE NRPerm<T>::generate_all(void (*callback)(const NRPerm<T>&), int select)
+{
+int n=this->size();
+NRVec_from1<T> c(0,n);
+NRVec_from1<T> d(1,n);
+int j,k,s;
+T q;
+T t;
+
+this->identity();
+PERM_RANK_TYPE sumperm=0;
+
+p2:
+        sumperm++;
+        if(!select || (select&1) == (sumperm&1)) (*callback)(*this);
+        j=n; s=0;
+p4:
+        q=c[j]+d[j];
+        if(q<0) goto p7;
+        if(q==j) goto p6;
+
+        t=(*this)[j-c[j]+s]; (*this)[j-c[j]+s]=(*this)[j-q+s]; (*this)[j-q+s]=t;
+        c[j]=q;
+        goto p2;
+p6:
+        if(j==1) goto end;
+        s++;
+p7:
+        d[j]= -d[j];
+        j--;
+        goto p4;
+
+end:
+return select? sumperm/2:sumperm;
+}
+
+
+template <typename T>
+static T _n;
+
+template <typename T>
+static void (*_callback)(const NRPerm<T>& ); 
+
+PERM_RANK_TYPE _sumperm;
+
+template <typename T>
+NRPerm<T> *_perm;
+
+template <typename T>
+void permg(T n)
+{
+if(n<= _n<T>)
+        {
+        for(T i=1;i<= _n<T>;i++)
+                {
+                if(!(*_perm<T>)[i]) //place not occupied
+			{
+                	(*_perm<T>)[i]=n;
+                	permg(n+1);
+                	(*_perm<T>)[i]=0;
+			}
+                }
+        }
+else
+        {
+        _sumperm++;
+	(*_callback<T>)(*_perm<T>);
+        }
+}
+
+
+template <typename T>
+PERM_RANK_TYPE NRPerm<T>::generate_all2(void (*callback)(const NRPerm<T>&))
+{
+this->copyonwrite();
+this->clear();
+_n<T> = this->size();
+_callback<T> =callback;
+_sumperm=0;
+_perm<T> = this;
+permg(1);
+return _sumperm;
+}
+
+
+template <typename T>
+PERM_RANK_TYPE NRPerm<T>::generate_all_lex(void (*callback)(const NRPerm<T>&))
+{
+PERM_RANK_TYPE np=0;
+this->identity();
+do{
+++np;
+(*callback)(*this);
+}while(this->next());
+return np;
+}
+
+template <typename T>
+PERM_RANK_TYPE NRPerm<T>::rank() const
+{
+int c,i,k;
+PERM_RANK_TYPE r;
+int n= this->size();
+
+r=0;
+for (k=1; k<=n; ++k)
+        {
+        T l=(*this)[k];
+        c=0;
+        for(i=k+1;i<=n;++i) if((*this)[i]<l) ++c;
+        r+= c;
+	i=n-k;
+        if(i) r*= i;
+        }
+return r;
+}
+
+
+template <typename T>
+NRVec_from1<T> NRPerm<T>::inversions(const int type, PERM_RANK_TYPE *prank) const
+{
+PERM_RANK_TYPE s=0;
+int n=this->size();
+int j,k;
+T l;
+
+NRVec_from1<T> i(n);
+i.clear();
+
+switch(type) {
+case 3:  /*number of elements right from p[j] smaller < p[j]*/
+for(j=1;j<n;++j) /*number of elements right from j smaller <j*/
+        {
+        l=(*this)[j];
+        for(k=n;k>j;--k)
+                {
+                if((*this)[k]<l) ++i[j];
+                }
+        }
+        break;
+case 2: /*number of elements left from p[j] bigger >p[j]*/
+for(j=2;j<=n;++j) /*number of elements left from j bigger >j*/
+        {
+        l=(*this)[j];
+        for(k=1;k<j;++k)
+                {
+                if((*this)[k]>l) ++i[j];
+                }
+        }
+        break;
+case 1:
+for(j=1;j<=n;++j) /*number of elements right from j smaller <j*/
+        {
+        for(k=n;k>=1;--k)
+                {
+                if((*this)[k]==j) break;
+                if((*this)[k]<j) ++i[j];
+                }
+        }
+        break;
+case 0:
+for(j=1;j<=n;++j) /*number of elements left from j bigger >j*/
+        {
+        for(k=1;k<=n;++k)
+                {
+                if((*this)[k]==j) break;
+                if((*this)[k]>j) ++i[j];
+                }
+        }
+        break;
+default: laerror("illegal type in inversions");
+
+if(prank)
+	{
+	if(type!=3) laerror("rank can be computed from inversions only for type 3");
+	l=1;
+	for(j=1;j<n;++j)
+                {
+                l*= j;
+                *prank += l*i[n-j];
+                }
+	}
+}
+
+return i;
+}
+
+template <typename T>
+NRPerm<T>::NRPerm(const int type, const NRVec_from1<T> &i)
+{
+int n=i.size();
+this->resize(n);
+
+int k,l;
+T j;
+
+switch(type) {
+case 2:
+case 1:
+        for(l=0,j=1; j<=n; ++j,++l)
+                {
+                /*shift left and place*/
+                for(k=n-l+1; k<=n-i[j]; ++k) (*this)[k-1]=(*this)[k];
+                (*this)[n-i[j]]=j;
+                }
+        break;
+case 3:
+case 0:
+        for(l=0,j=n; j>0; --j,++l)
+                {
+                /*shift right and place*/
+                for(k=l; k>=i[j]+1; --k) (*this)[k+1]=(*this)[k];
+                (*this)[i[j]+1]=j;
+                }
+        break;
+default: laerror("illegal type in nrperm from inversions");
+}
+
+if(type>=2) (*this) = this->inverse();
+}
+
+
+
+template <typename T>
+NRPerm<T>::NRPerm(const int n, const PERM_RANK_TYPE rank)
+{
+this->resize(n);
+NRVec_from1<T> inv(n) ;
+#ifdef DEBUG
+if(rank>=factorial(n)) laerror("illegal rank for this n");
+#endif
+inv[n]=0;
+int k;
+PERM_RANK_TYPE r=rank;
+for(k=n-1; k>=0; --k)
+        {
+        PERM_RANK_TYPE t;
+        t=factorial(k);
+        inv[n-k]=r/t;
+        r=r%t;
+        }
+
+*this =  NRPerm(3,inv);
+}
+
+#define MAXFACT 20
+PERM_RANK_TYPE factorial(const int n)
+{
+static int ntop=20;
+static  PERM_RANK_TYPE a[MAXFACT+1]={1,1,2,6,24,120,720,5040,
+40320,
+362880,
+3628800,
+39916800ULL,
+479001600ULL,
+6227020800ULL,
+87178291200ULL,
+1307674368000ULL,
+20922789888000ULL,
+355687428096000ULL,
+6402373705728000ULL,
+121645100408832000ULL,
+2432902008176640000ULL};
+
+int j;
+if (n < 0) laerror("negative argument of factorial");
+if (n > MAXFACT) laerror("overflow in factorial");
+        while (ntop<n) {
+                j=ntop++;
+                a[ntop]=a[j]*ntop;
+                }
+return a[n];
+}
+
+
 ////////////////////////////////////////////////////////
 
 template <typename T>
diff --git a/permutation.h b/permutation.h
index db3258e..d75a84e 100644
--- a/permutation.h
+++ b/permutation.h
@@ -23,6 +23,8 @@
 #include "la_traits.h"
 #include "vec.h"
 
+typedef   unsigned long long PERM_RANK_TYPE;
+
 //permutations are always numbered from 1; offset is employed when applied to vectors and matrices
 
 namespace LA {
@@ -39,7 +41,6 @@ public:
 	NRPerm(): NRVec_from1<T>() {};
 	NRPerm(const int n) : NRVec_from1<T>(n) {};
 	NRPerm(const NRVec_from1<T> &rhs): NRVec_from1<T>(rhs) {};
-	NRPerm(const T &a, const int n): NRVec_from1<T>(a, n) {};
         NRPerm(const T *a, const int n): NRVec_from1<T>(a, n) {};
 	explicit NRPerm(const CyclePerm<T> &rhs, const int n=0); 
 
@@ -52,15 +53,17 @@ public:
 	NRPerm conjugate_by(const NRPerm q) const; //q^-1 p q
 	int parity() const;
 	void randomize(void); //uniformly random by Fisher-Yates shuffle
-
-
-	//TODO:
-	//@@@permgener
-	//@@@next permutation
-	//@@@lex rank
-	//@@@inversion tables
+	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_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 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
+	explicit NRPerm(const int n, const PERM_RANK_TYPE rank); //compute permutation from its rank
 };
 
+extern PERM_RANK_TYPE factorial(const int n);
 
 //permutations represented in the cycle format
 template <typename T>
@@ -71,6 +74,7 @@ public:
 
 	bool is_valid() const; //is it really a permutation
 	bool is_identity() const; //no cycles of length > 1
+	void identity() {this->resize(0);};
 	CyclePerm inverse() const; //reverse all cycles
 	int parity() const; //negative if having odd number of even-length cycles
 	T max() const {T m=0; for(int i=1; i<=this->size(); ++i) {T mm= (*this)[i].max(); if(mm>m) m=mm;} return m;}
diff --git a/t.cc b/t.cc
index 4ba4c42..2d5aa5a 100644
--- a/t.cc
+++ b/t.cc
@@ -53,6 +53,21 @@ inline int randind(const int n)
 
 complex<double> mycident (const complex<double>&x) {return x;}
 
+void printme(const NRPerm<int> &p)
+{
+PERM_RANK_TYPE rank=p.rank();
+int n=p.size();
+cout<<p.rank()<<" "<<p;
+NRPerm<int> qq(n,rank);
+if(qq!=p)  laerror("error in rank algorithm");
+for(int i=0; i<4; ++i)
+	{
+	NRVec_from1<int>  inv=p.inversions(i);
+	NRPerm<int> q(i,inv);
+	if(q!=p) laerror("error in inversions algorithm");
+	}
+}
+
 
 int main()
 {
@@ -2038,7 +2053,7 @@ cout<<vvv;
 cout<<"error "<<(v-vvv).norm()<<endl;
 }
 
-if(1)
+if(0)
 {
 int seed;
 int f=open("/dev/random",O_RDONLY);
@@ -2057,7 +2072,26 @@ NRVec<double> v4=vv.permuted(p,false);
 cout<<v<<vv;
 cout<<vvv<<v4<<p;
 cout <<"error "<<(vv-vvv).norm() <<" "<<(v-v4).norm()<<endl;
+}
 
+if(0)
+{
+int n;
+cin >>n;
+NRPerm<int> p(n);
+p.identity();
+do{
+cout <<p;
+}while(p.next());
+}
+
+if(1)
+{
+int n;
+cin >>n;
+NRPerm<int> p(n);
+int tot=p.generate_all_lex(printme);
+cout <<"generated "<<tot<<endl;
 }
 
 }