diff --git a/c_ext/Makefile b/c_ext/Makefile
new file mode 100644
index 0000000..71c3eb6
--- /dev/null
+++ b/c_ext/Makefile
@@ -0,0 +1,5 @@
+# You must specify your env variable LSAPE_DIR
+#LSAPE_DIR=/home/bgauzere/Téléchargements/lsape/include/
+
+liblsap.so:lsap.cpp
+	g++ -fPIC -I/home/bgauzere/Téléchargements/lsape/include/ -shared lsap.cpp -o liblsap.so -O3 -I$(LSAPE_DIR)
diff --git a/c_ext/lsap.cpp b/c_ext/lsap.cpp
new file mode 100644
index 0000000..712a092
--- /dev/null
+++ b/c_ext/lsap.cpp
@@ -0,0 +1,43 @@
+/*
+Python wrapper
+*/
+
+#include "hungarian-lsape.hh"
+#include "hungarian-lsap.hh"
+
+#include <cstdio>
+
+extern "C" int lsap(double * C, const int nm, long * rho, long * varrho){
+  double * u = new double[nm];
+  double * v = new double[nm];
+
+  int * rho_int = new int[nm];
+  int * varrho_int = new int[nm];
+
+  hungarianLSAP(C,nm,nm,rho_int,u,v,varrho_int);
+  //Find a better way to do
+  for (int i =0;i<nm;i++){
+    rho[i] = (long)(rho_int[i]);
+    varrho[i] = (long)(varrho_int[i]);
+  }  
+  return 0;
+}
+
+
+
+extern "C" int * lsape(double * C, const int n, const int m, long * rho, long * varrho){
+  double * u = new double[n];
+  double * v = new double[m];
+
+  int * rho_int = new int[n];
+  int * varrho_int = new int[m];
+
+  hungarianLSAPE(C,n,m,rho_int,varrho_int,u,v);
+  for (int i =0;i<n;i++)
+    rho[i] = (long)(rho_int[i]);
+
+  for (int i =0;i<m;i++)
+    varrho[i] = (long)(varrho_int[i]);
+  
+  return 0;
+}
diff --git a/c_ext/lsape_binders.py b/c_ext/lsape_binders.py
new file mode 100644
index 0000000..5efd0ff
--- /dev/null
+++ b/c_ext/lsape_binders.py
@@ -0,0 +1,23 @@
+import numpy as np
+import ctypes as c
+from ctypes import cdll
+import os.path
+
+def lsap_solverHG(C):
+    ''' Binding for lsape hungarian solver '''
+
+    nm = C.shape[0]
+    dll_name = 'liblsap.so'
+    lib = cdll.LoadLibrary(os.path.abspath(
+        os.path.join(os.path.dirname(__file__), dll_name)))
+    lib.lsap.restype = c.c_int
+    rho = np.zeros((nm, 1), int)
+    varrho = np.zeros((nm, 1), int)
+    C[C == np.inf] = 10000
+
+    lib.lsap(c.c_void_p(C.transpose().ctypes.data),
+             c.c_int(nm),
+             c.c_void_p(rho.ctypes.data),
+             c.c_void_p(varrho.ctypes.data))
+
+    return np.array(range(0, nm)), np.array([c.c_int(i).value for i in varrho])