permutation: in-place inversion

permutation.cc

@@ -72,6 +72,48 @@ for(T i=1; i<=n; ++i) if(used[i]!=1) return 0;
 return 1;
 }
 
+
+
+
+//in-place o(n^2) cf. arxiv 1901.01926v2
+//better alg. exists but more complicated
+template <typename T>
+void NRPerm<T>::reverse_cycle(int start)
+{
+int cur=(*this)[start];
+int prev=start;
+while(cur!=start)
+	{
+	int next=(*this)[cur];
+	(*this)[cur]=prev;
+	prev=cur;
+	cur=next;
+	}
+(*this)[start]=prev;
+}
+
+
+template <typename T>
+int NRPerm<T>::cycle_leader(int start) const
+{
+int cur=(*this)[start];
+int smallest=start;
+while(cur!=start)
+	{
+	if(cur<smallest) smallest=cur;
+	cur=(*this)[cur];
+	}
+return smallest;
+}
+
+
+template <typename T>
+void NRPerm<T>::inverseme()
+{
+for(int i=1; i<=size(); ++i) if(cycle_leader(i)==i) reverse_cycle(i);
+}
+
+
 template <typename T>
 NRPerm<T> NRPerm<T>::inverse() const
 {

permutation.h

@@ -64,6 +64,9 @@ public:
 	bool is_identity() const;
 	CompressedPartition<T> cycles() const {return CyclePerm<T>(*this).cycles(size());};
 	NRPerm inverse() const;
+	void reverse_cycle(int i); //aux for inverseme
+	int cycle_leader(int i) const; //aux for inverseme
+	void inverseme(); //in-place
 	NRPerm reverse() const; //backward order
 	NRPerm operator&(const NRPerm &rhs) const; //concatenate the permutations this,rhs, renumbering rhs (not commutative)
 	NRPerm operator|(const NRPerm &rhs) const; //concatenate the permutations rhs,this, renumbering rhs (not commutative)

t.cc

@@ -4788,7 +4788,7 @@ cout<<tt;
 cout <<"Error = "<<(t-tt).norm()<<endl;
 }
 
-if(1)
+if(0)
 {
 Tensor<double> tt;
 tt.numpy_read(argv[1]);
@@ -4798,4 +4798,15 @@ cout<<"part\n"<<tt.subtensor1(1);
 cout<<"part\n"<<tt.subtensor1(2);
 }
 
+if(1)
+{
+int n;
+cin >>n;
+NRPerm<int> p(n); p.randomize();
+NRPerm<int> q=p.inverse();
+p.inverseme();
+if(p!=q) laerror("inverseme failed");
+else cout <<p<<"OK\n";
+}
+
 }//main