Add premutations module

ext/garside/init.hpp

@@ -28,6 +28,9 @@
 
 extern Gomu::Type* type_ArtinWordA;
 extern Gomu::Type* type_DualWordA;
+extern Gomu::Type* type_monoid_family;
+extern Gomu::Type* type_word;
+
 
 //****************
 //* MonoidFamily *

ext/permutations/Makefile

@@ -0,0 +1,15 @@
+CPP 	= g++ -g --std=c++11 -march=corei7 -Wno-return-local-addr -fPIC -rdynamic -fmax-errors=1 -I/usr/local/include
+LDFLAGS = #-L/usr/local/lib -lgmpxx -lgmp -lflint
+MOD 	= ../permutations.so
+
+all: $(MOD)
+
+
+permutations.o: ../permutations/permutations.cpp
+	$(CPP) -c $< -o $@
+
+$(MOD): init.cpp init.hpp permutations.o
+	$(CPP) -shared  $(LDFLAGS) $^ -o $@
+
+clean:
+	-$(RM) $(MOD) *.o

ext/permutations/coxeter.hpp

@@ -0,0 +1,26 @@
+/**
+ * This file is part of Gomu.
+ *
+ *  Copyright 2016 by Jean Fromentin <jean.fromentin@math.cnrs.fr>
+ *
+ * Gomu is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Gomu is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Gomu. If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+#ifndef COXETER_HPP
+#define COXETER_HPP
+
+enum CoxeterType{CoxeterA,CoxeterB,CoxeterD};
+
+
+#endif

ext/permutations/init.cpp

@@ -0,0 +1,57 @@
+#include "init.hpp"
+
+//******************
+//* Global objects *
+//******************
+
+Gomu::Type* type_PermutationA;
+Gomu::Type* type_PermutationB;
+Gomu::Type* type_PermutationD;
+Gomu::Type* type_PermutationEnumeratorA;
+Gomu::Type* type_PermutationEnumeratorB;
+Gomu::Type* type_PermutationEnumeratorD;
+
+//*************************
+//* Extension inilisation *
+//*************************
+
+extern "C"{
+  Gomu::Module::Type types[]={
+    {"PermutationA",Permutations::dispPerm<CoxeterA>,Permutations::delPerm<CoxeterA>,Permutations::copyPerm<CoxeterA>,Permutations::cmpPerm<CoxeterA>,&type_PermutationA},
+    {"PermutationB",Permutations::dispPerm<CoxeterB>,Permutations::delPerm<CoxeterB>,Permutations::copyPerm<CoxeterB>,Permutations::cmpPerm<CoxeterB>,&type_PermutationB},
+    {"PermutationD",Permutations::dispPerm<CoxeterD>,Permutations::delPerm<CoxeterD>,Permutations::copyPerm<CoxeterD>,Permutations::cmpPerm<CoxeterD>,&type_PermutationD},
+    {"PermutationEnumeratorA",Permutations::dispPermE<CoxeterA>,Permutations::delPermE<CoxeterA>,Permutations::copyPermE<CoxeterA>,Permutations::cmpPermE<CoxeterA>,&type_PermutationEnumeratorA},
+    {"PermutationEnumeratorB",Permutations::dispPermE<CoxeterB>,Permutations::delPermE<CoxeterB>,Permutations::copyPermE<CoxeterB>,Permutations::cmpPermE<CoxeterB>,&type_PermutationEnumeratorB},
+    {"PermutationEnumeratorD",Permutations::dispPermE<CoxeterD>,Permutations::delPermE<CoxeterD>,Permutations::copyPermE<CoxeterD>,Permutations::cmpPermE<CoxeterD>,&type_PermutationEnumeratorD},
+   TYPE_SENTINEL
+  };
+
+  //--- Functions ---//
+  Gomu::Module::Function functions[]={
+    {"PermutationA","permutationA",{"Integer"},(void*)Permutations::intToPerm<CoxeterA>},
+    {"PermutationB","permutationB",{"Integer"},(void*)Permutations::intToPerm<CoxeterB>},
+    {"PermutationD","permutationD",{"Integer"},(void*)Permutations::intToPerm<CoxeterD>},
+    {"PermutationEnumeratorA","permutationEnumeratorA",{"Integer"},(void*)Permutations::intToPermE<CoxeterA>},
+    {"PermutationEnumeratorB","permutationEnumeratorB",{"Integer"},(void*)Permutations::intToPermE<CoxeterB>},
+    {"PermutationEnumeratorD","permutationEnumeratorD",{"Integer"},(void*)Permutations::intToPermE<CoxeterD>},
+    FUNC_SENTINEL
+  };
+
+  //--- Member functions ---//
+  Gomu::Module::Function member_functions[]={
+    //--- PermutationEnumertor ---//
+    {"PermutationA","get",{"PermutationEnumeratorA"},(void*)Permutations::PE_get<CoxeterA>},
+    {"PermutationB","get",{"PermutationEnumeratorB"},(void*)Permutations::PE_get<CoxeterB>},
+    {"PermutationD","get",{"PermutationEnumeratorD"},(void*)Permutations::PE_get<CoxeterD>},
+    {"Boolean","next",{"PermutationEnumeratorA"},(void*)Permutations::PE_next<CoxeterA>},
+    {"Boolean","next",{"PermutationEnumeratorB"},(void*)Permutations::PE_next<CoxeterB>},
+    {"Boolean","next",{"PermutationEnumeratorD"},(void*)Permutations::PE_next<CoxeterD>},
+    {"Void","reset",{"PermutationEnumeratorA"},(void*)Permutations::PE_reset<CoxeterA>},
+    {"Void","reset",{"PermutationEnumeratorB"},(void*)Permutations::PE_reset<CoxeterB>},
+    {"Void","reset",{"PermutationEnumeratorD"},(void*)Permutations::PE_reset<CoxeterD>},
+    {"Integer","size",{"PermutationEnumeratorA"},(void*)Permutations::PE_size<CoxeterA>},
+    {"Integer","size",{"PermutationEnumeratorB"},(void*)Permutations::PE_size<CoxeterB>},
+    {"Integer","size",{"PermutationEnumeratorD"},(void*)Permutations::PE_size<CoxeterD>},
+    FUNC_SENTINEL
+  };
+}

ext/permutations/init.hpp

@@ -0,0 +1,99 @@
+/**
+ * This file is part of Gomu.
+ *
+ *  Copyright 2016 by Jean Fromentin <jean.fromentin@math.cnrs.fr>
+ *
+ * Gomu is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Gomu is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Gomu. If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+#include "../../module.hpp"
+#include "permutations.hpp"
+
+using namespace std;
+
+//*****************
+//* Global object *
+//*****************
+
+extern Gomu::Type* type_PermutationA;
+extern Gomu::Type* type_PermutationB;
+extern Gomu::Type* type_PermutationD;
+extern Gomu::Type* type_PermutationEnumeratorA;
+extern Gomu::Type* type_PermutationEnumeratorB;
+extern Gomu::Type* type_PermutationEnumeratorD;
+
+
+namespace Permutations{
+  template<CoxeterType T> string dispPerm(void*);
+  template<CoxeterType T> void delPerm(void*);
+  template<CoxeterType T> void* copyPerm(void*);
+  template<CoxeterType T> int cmpPerm(void*,void*);
+
+  template<CoxeterType T> string dispPermE(void*);
+  template<CoxeterType T> void delPermE(void*);
+  template<CoxeterType T> void* copyPermE(void*);
+  template<CoxeterType T> int cmpPermE(void*,void*);
+
+
+  //--- Constructors ---//
+  template<CoxeterType T> void* intToPerm(void*);
+  template<CoxeterType T> void* intToPermE(void*);
+
+  //--- Permutation functions ---//
+
+  //--- Permutation enumerator functions ---//
+  template<CoxeterType T> void* PE_get(void*);
+  template<CoxeterType T> void* PE_next(void*);
+  template<CoxeterType T> void* PE_reset(void*);
+  template<CoxeterType T> void* PE_size(void*);
+}
+
+
+namespace Permutations{
+  template<CoxeterType T> inline string dispPerm(void* P){return ((Permutation<T>*)P)->display();}
+  template<CoxeterType T> inline void delPerm(void* v){delete ((Permutation<T>*)v);}
+  template<CoxeterType T> inline void* copyPerm(void* v){return new Permutation<T>(*((Permutation<T>*)v));}
+  template<CoxeterType T> inline int cmpPerm(void* v1,void* v2){return ((Permutation<T>*)v1)->cmp(*(Permutation<T>*)v2);}
+
+  template<CoxeterType T> inline string dispPermE(void* P){return ((PermutationEnumerator<T>*)P)->display();}
+  template<CoxeterType T> inline void delPermE(void* v){delete ((PermutationEnumerator<T>*)v);}
+  template<CoxeterType T> inline void* copyPermE(void* v){return new PermutationEnumerator<T>(*((PermutationEnumerator<T>*)v));}
+  template<CoxeterType T> inline int cmpPermE(void* v1,void* v2){return ((PermutationEnumerator<T>*)v1)->cmp(*(PermutationEnumerator<T>*)v2);}
+
+  //--- Constructors ---//
+  template<CoxeterType T>
+  inline void*
+  intToPerm(void* z){
+    slong n=Gomu::get_slong(z);
+    if(n<0 or n>15) RuntimeError("A postive integer < 16 is needed");
+    return new Permutation<T>(n);
+  }
+
+  template<CoxeterType T>
+  inline void*
+  intToPermE(void* z){
+    slong n=Gomu::get_slong(z);
+    if(n<0 or n>15) RuntimeError("A postive integer < 16 is needed");
+    return new PermutationEnumerator<T>(n);
+  }
+
+  //***********************************
+  //* PermutationEnumerator functions *
+  //***********************************
+  template<CoxeterType T> inline void* PE_get(void* PE){return new Permutation<T>(((PermutationEnumerator<T>*)PE)->get());}
+  template<CoxeterType T> inline void* PE_next(void* PE){return Gomu::to_boolean(((PermutationEnumerator<T>*)PE)->next());}
+  template<CoxeterType T> inline void* PE_reset(void* PE){((PermutationEnumerator<T>*)PE)->reset();return nullptr;}
+  template<CoxeterType T> inline void* PE_size(void* PE){return Gomu::to_integer(((PermutationEnumerator<T>*)PE)->size());}
+  
+}

ext/permutations/permutations.cpp

@@ -0,0 +1,145 @@
+/**
+ * This file is part of Gomu.
+ *
+ *  Copyright 2016 by Jean Fromentin <jean.fromentin@math.cnrs.fr>
+ *
+ * Gomu is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Gomu is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Gomu. If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+
+#include "permutations.hpp"
+
+
+PermutationEnumerator<CoxeterA>::PermutationEnumerator(uint64 r):n(r),sigma(r){
+}
+
+void
+PermutationEnumerator<CoxeterA>::reset(){
+  sigma=Permutation<CoxeterA>(n);
+}
+
+const Permutation<CoxeterA>&
+PermutationEnumerator<CoxeterA>::get() const{
+  return sigma;
+}
+
+bool
+PermutationEnumerator<CoxeterA>::next(){
+  for(int first=n-1;first>0;--first){
+    if(sigma.win[first]<sigma.win[first+1]){
+      int j=n;
+      while(sigma.win[j]<sigma.win[first]) --j;
+      swap(sigma.win[j],sigma.win[first]);
+      int i=first+1;
+      j=n;
+      while(i<j) swap(sigma.win[i++],sigma.win[j--]);
+      return true;
+    }
+  }
+  return false;
+}
+
+PermutationEnumerator<CoxeterB>::PermutationEnumerator(uint64 r):PermutationEnumerator<CoxeterA>(r),sigmaB(r),sign(0){
+}
+
+void
+PermutationEnumerator<CoxeterB>::reset(){
+  PermutationEnumerator<CoxeterA>::reset();
+  sigmaB=Permutation<CoxeterB>(n);
+  sign=0;
+}
+
+const Permutation<CoxeterB>&
+PermutationEnumerator<CoxeterB>::get() const{
+  return sigmaB;
+}
+
+bool
+PermutationEnumerator<CoxeterB>::next(){
+  ++sign;
+  if(sign==(1<<n)){
+    sign=0;
+    if(not PermutationEnumerator<CoxeterA>::next()) return false;
+    sigmaB.m128=sigma.m128;
+  }
+  else{
+    uint64 mask=(1<<(n-1));
+    for(uint64 i=1;i<=n;++i){
+      sigmaB.win[i]=(mask&sign)?-sigma.win[i]:sigma.win[i];
+      mask>>=1;
+    }
+  }
+  return true;
+};
+
+PermutationEnumerator<CoxeterD>::PermutationEnumerator(uint64 r):PermutationEnumerator<CoxeterA>(r),sigmaD(r),sign(0){
+}
+
+void
+PermutationEnumerator<CoxeterD>::reset(){
+  PermutationEnumerator<CoxeterA>::reset();
+  sigmaD=Permutation<CoxeterD>(n);
+  sign=0;
+}
+
+const Permutation<CoxeterD>&
+PermutationEnumerator<CoxeterD>::get() const{
+  return sigmaD;
+}
+
+bool
+PermutationEnumerator<CoxeterD>::next(){
+  ++sign;
+  if(sign==(1<<n)){
+    sign=0;
+    if(not PermutationEnumerator<CoxeterA>::next()) return false;
+    sigmaD.m128=sigma.m128;
+  }
+  else{
+    if(__builtin_popcountll(sign)%2==0){
+      uint64 mask=1;
+      for(uint64 i=1;i<=n;++i){
+	sigmaD.win[i]=(mask&sign)?-sigma.win[i]:sigma.win[i];
+	mask<<=1;
+      }
+    }
+    else return next();
+  }
+  return true;
+};
+
+
+
+int
+PermutationEnumerator<CoxeterA>::cmp(const PermutationEnumerator<CoxeterA>& E) const{
+  if(n==E.n) return 0;
+  if(n<E.n) return 1;
+  return -1;
+}
+
+ int
+PermutationEnumerator<CoxeterB>::cmp(const PermutationEnumerator<CoxeterB>& E) const{
+  if(n==E.n) return 0;
+  if(n<E.n) return 1;
+  return -1;
+}
+
+int
+PermutationEnumerator<CoxeterD>::cmp(const PermutationEnumerator<CoxeterD>& E) const{
+  if(n==E.n) return 0;
+  if(n<E.n) return 1;
+  return -1;
+}
+
+

ext/permutations/permutations.hpp

@@ -0,0 +1,251 @@
+/**
+ * This file is part of Gomu.
+ *
+ *  Copyright 2016 by Jean Fromentin <jean.fromentin@math.cnrs.fr>
+ *
+ * Gomu is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Gomu is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Gomu. If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+
+#ifndef PERMUTATIONS_HPP
+#define PERMUTATIONS_HPP
+
+#include <cassert>
+#include <cstdint>
+#include <iostream>
+#include "smmintrin.h"
+#include "coxeter.hpp"
+
+using namespace std;
+
+typedef int8_t int8;
+typedef uint32_t uint32;
+typedef uint64_t uint64;
+typedef int64_t int64;
+
+static const __m128i id128=_mm_setr_epi8(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15);
+static const __m128i shmask=_mm_setr_epi8(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,15);
+
+void swap(int8& a,int8& b);
+uint64 fact(uint64 n);
+
+template<CoxeterType T> class Permutation{
+public:
+  uint64 n;
+  static Permutation Ptemp;
+  union{
+    __m128i m128;
+    int8 win[16];
+  };
+  Permutation(uint64=0);
+  void setInverse(Permutation&) const;
+  const Permutation&  getInverse() const;
+  uint64 descentInteger() const;
+  int cmp(const Permutation&) const;
+  string display() const;
+  template<CoxeterType U> friend ostream& operator<<(ostream&,const Permutation<U>&);
+};
+
+
+template<CoxeterType T> class PermutationEnumerator;
+
+  
+template<> class PermutationEnumerator<CoxeterA>{
+protected:
+  Permutation<CoxeterA> sigma;
+  uint64 n;
+public:
+  PermutationEnumerator(uint64);
+  void reset();
+  bool next();
+  size_t size() const;
+  const Permutation<CoxeterA>& get() const;
+  int cmp(const PermutationEnumerator<CoxeterA>& E) const;
+  string display() const;
+  friend ostream& operator<<(ostream&,const PermutationEnumerator<CoxeterA>&);
+};
+
+template<> class PermutationEnumerator<CoxeterB>:public PermutationEnumerator<CoxeterA>{
+protected:
+  Permutation<CoxeterB> sigmaB;
+  uint64 sign;
+public:
+  PermutationEnumerator(uint64);
+  void reset();
+  bool next();
+  size_t size() const;
+  const Permutation<CoxeterB>& get() const;
+  int cmp(const PermutationEnumerator<CoxeterB>& E) const;
+  string display() const;
+  friend ostream& operator<<(ostream&,const PermutationEnumerator<CoxeterB>&);
+};
+
+template<> class PermutationEnumerator<CoxeterD>:public PermutationEnumerator<CoxeterA>{
+protected:
+  Permutation<CoxeterD> sigmaD;
+  uint64 sign;
+public:
+  PermutationEnumerator(uint64);
+  void reset();
+  bool next();
+  size_t size() const;
+  const Permutation<CoxeterD>& get() const;
+  int cmp(const PermutationEnumerator<CoxeterD>& E) const;
+  string display() const;
+  friend ostream& operator<<(ostream&,const PermutationEnumerator<CoxeterD>&);
+};
+
+
+typedef Permutation<CoxeterA> PermutationA;
+typedef Permutation<CoxeterB> PermutationB;
+typedef Permutation<CoxeterD> PermutationD;
+
+inline void
+swap(int8& a,int8& b){
+  int8 t=a;
+  a=b;
+  b=t;
+}
+
+inline uint64
+fact(uint64 n){
+  if(n<2) return 1;
+  return n*fact(n-1);
+}
+
+//***************
+//* Permutation *
+//***************
+
+template<CoxeterType T> Permutation<T> Permutation<T>::Ptemp;
+
+template<CoxeterType T> inline
+Permutation<T>::Permutation(uint64 r):n(r),m128(id128){
+}
+
+template<CoxeterType T> void
+Permutation<T>::setInverse(Permutation<T>& inv) const{
+  inv.n=n;
+  for(uint64 i=1;i<=n;++i){
+    int8 v=win[i];
+    if(v>0) inv.win[v]=i;
+    else inv.win[-v]=-i;
+  }
+}
+
+template<CoxeterType T> const Permutation<T>&
+Permutation<T>::getInverse() const{
+  setInverse(Ptemp);
+  return Ptemp;
+}
+
+template<> inline uint64
+Permutation<CoxeterA>::descentInteger() const{
+  __m128i temp=_mm_shuffle_epi8(m128,shmask);
+  temp=_mm_cmpgt_epi8(m128,temp);
+  return _mm_movemask_epi8(temp);
+}
+
+template<> inline uint64
+Permutation<CoxeterB>::descentInteger() const{
+  __m128i temp=_mm_shuffle_epi8(m128,shmask);
+  temp=_mm_cmpgt_epi8(m128,temp);
+  return _mm_movemask_epi8(temp);
+}
+
+template<> inline uint64
+Permutation<CoxeterD>::descentInteger() const{
+  __m128i temp=_mm_shuffle_epi8(m128,shmask);
+  temp=_mm_cmpgt_epi8(m128,temp);
+  uint64 res=_mm_movemask_epi8(temp);
+  if(win[1]+win[2]<0) res|=1L;
+  else res&=(~1L);
+  return res;
+}
+
+template<CoxeterType T> int
+Permutation<T>::cmp(const Permutation& P) const{
+  if(n!=P.n){
+    if(n<P.n) return 1;
+    return -1;
+  }
+  for(uint64 i=1;i<=n;++i){
+    if(win[i]!=P.win[i]){
+      if(win[i]<P.win[i]) return 1;
+      return -1;
+    }
+  }
+  return 0;
+}
+
+
+template<CoxeterType T> string
+Permutation<T>::display() const{
+  string str;
+  if(n==0) return "()";
+  str="(\033[34m"+to_string((int64)win[1])+"\033[0m";
+  for(uint64 i=2;i<=n;++i){
+    str+=",\033[34m"+to_string((int64)win[i])+"\033[0m";
+  }
+  str+=')';
+  return str;
+}
+
+template<CoxeterType T> inline ostream&
+operator<<(ostream& os,const Permutation<T>& sigma){
+  return os<<sigma.display();
+}
+
+//
+// PermutationEnumerator
+//
+
+inline size_t
+PermutationEnumerator<CoxeterA>::size() const{
+  return fact(n);
+}
+
+inline size_t
+PermutationEnumerator<CoxeterB>::size() const{
+  return (1<<n)*fact(n);
+}
+
+inline size_t
+PermutationEnumerator<CoxeterD>::size() const{
+  if(n==0) return 1;
+  return (1<<(n-1))*fact(n);
+}
+
+inline string
+PermutationEnumerator<CoxeterA>::display() const{
+  return "Enumerator for permutations of type A and rank <= "+to_string(n);
+}
+
+inline string
+PermutationEnumerator<CoxeterB>::display() const{
+  return "Enumerator for permutations of type B and rank <= "+to_string(n);
+}
+
+inline string
+PermutationEnumerator<CoxeterD>::display() const{
+  return "Enumerator for permutations of type D and rank <= "+to_string(n);
+}
+
+template<CoxeterType T> inline ostream&
+operator<<(ostream& os,const PermutationEnumerator<T>& E){
+  return os<<E.display();
+}
+
+
+#endif

kernel.hpp

@@ -285,7 +285,6 @@ namespace Gomu{
   //------
   // Type
   //------
-
   
   inline string
   type_disp(void* T){