diff --git a/miscprob.cc b/miscprob.cc
index 89f7623bb00db7e90f685607fcb253fee599551d..00b08a31b1d0357ba42e7e633bb1648136de4d1e 100644
--- a/miscprob.cc
+++ b/miscprob.cc
@@ -145,4 +145,33 @@ ReturnMatrix mvnrnd(const RowVector& mu, const SymmetricMatrix& covar, int nsamp
   return mvn.next(nsamp);
 }
 
+ReturnMatrix perms(const int n){
+  if(n<=1){
+    Matrix P(1,1);
+    P << n;
+    P.Release();
+    return P;
+  }
+  Matrix Q = perms(n-1);  // recursive calls
+  int m = Q.Nrows();
+  Matrix P(n*m,n);
+  for(int i=1;i<=m;i++){
+    P(i,1)=n;
+    for(int j=1;j<=Q.Ncols();j++)
+      P(i,j+1)=Q(i,j);
+  }
+  for(int i=n-1;i>=1;i--){
+    int jj=1;
+    for(int j=(n-i)*m+1;j<=(n-i+1)*m;j++){
+      P(j,1)=i;
+      for(int k=1;k<=n-1;k++){
+	P(j,k+1)= (Q(jj,k)==i) ? n : Q(jj,k);
+      }
+      jj++;
+    } 
+  }
+  P.Release();
+  return P;
+}
+
 }
diff --git a/miscprob.h b/miscprob.h
index 26c2e0494c018d02cd321829c1ed32f43ac5db4c..f0fb80bb7991a75170e1f490298f63ac6acc2124 100644
--- a/miscprob.h
+++ b/miscprob.h
@@ -41,6 +41,9 @@ namespace MISCMATHS {
 
   ReturnMatrix gammacdf(const RowVector& vals, const float mu = 0, const float var = 1);
 
+  // returns n! * n matrix of all possible permutations
+  ReturnMatrix perms(const int n);
+
   
   class Mvnormrandm
     {