diff --git a/setup.py b/setup.py
index 9ad58995a366e3d9f8a37ec702c7bb27d4ea9513..33d9b0169bf0f5642a5e8ecaa1867209b421e519 100644
--- a/setup.py
+++ b/setup.py
@@ -45,6 +45,7 @@ setup(
           "xbob/math/histogram.cpp",
           "xbob/math/linsolve.cpp",
           "xbob/math/pavx.cpp",
+          "xbob/math/norminv.cpp",
           "xbob/math/main.cpp",
           ],
         packages = packages,
diff --git a/xbob/math/main.cpp b/xbob/math/main.cpp
index a4bb4bbaacc4c97d596764d27a0c99a1e40bde32..04d3d00e5bdf79ee5e6dcd979132ca946669b616 100644
--- a/xbob/math/main.cpp
+++ b/xbob/math/main.cpp
@@ -16,6 +16,7 @@
 #include "histogram.h"
 #include "linsolve.h"
 #include "pavx.h"
+#include "norminv.h"
 
 PyDoc_STRVAR(s_histogram_intersection_str, "histogram_intersection");
 PyDoc_STRVAR(s_histogram_intersection_doc,
@@ -239,6 +240,29 @@ The width and height arrays are returned. The width array is a 64-bit\n\
 of the returned tuple) is a 64-bit **float** 1D array of the same size.\n\
 ");
 
+PyDoc_STRVAR(s_norminv_str, "norminv");
+PyDoc_STRVAR(s_norminv_doc,
+"norminv(p, mu, sigma) -> scalar\n\
+\n\
+Computes the inverse normal cumulative distribution for a probability\n\
+``p``, given a distribution with mean ``mu`` and standard deviation\n\
+``sigma``. The value ``p`` must lie in the range [0,1].\n\
+\n\
+Reference: `<http://home.online.no/~pjacklam/notes/invnorm/>`_\n\
+");
+
+PyDoc_STRVAR(s_normsinv_str, "normsinv");
+PyDoc_STRVAR(s_normsinv_doc,
+"normsinv(p) -> scalar\n\
+\n\
+Computes the inverse normal cumulative distribution for a probability\n\
+``p``, given a distribution with mean 0.0 and standard deviation 1.0.\n\
+It is equivalent as calling :py:func:`norminv(p, 0, 1)`. The value\n\
+``p`` must lie in the range [0,1].\n\
+\n\
+Reference: `<http://home.online.no/~pjacklam/notes/invnorm/>`_\n\
+");
+
 static PyMethodDef module_methods[] = {
     {
       s_histogram_intersection_str,
@@ -318,6 +342,18 @@ static PyMethodDef module_methods[] = {
       METH_VARARGS|METH_KEYWORDS,
       s_pavx_width_height_doc
     },
+    {
+      s_norminv_str,
+      (PyCFunction)py_norminv,
+      METH_VARARGS|METH_KEYWORDS,
+      s_norminv_doc
+    },
+    {
+      s_normsinv_str,
+      (PyCFunction)py_normsinv,
+      METH_VARARGS|METH_KEYWORDS,
+      s_normsinv_doc
+    },
     {0}  /* Sentinel */
 };
 
diff --git a/xbob/math/norminv.cpp b/xbob/math/norminv.cpp
index e97c328bcfa65ff109b395f60e86321875c63cc5..2d0cc0a331cb5bc6f39c6ba4ab22345f66faaa60 100644
--- a/xbob/math/norminv.cpp
+++ b/xbob/math/norminv.cpp
@@ -4,34 +4,59 @@
  * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
  *
  * @brief Binds the inverse normal cumulative distribution into python
- *
- * Copyright (C) 2011-2013 Idiap Research Institute, Martigny, Switzerland
- * 
- * This program 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, version 3 of the License.
- * 
- * This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <boost/python.hpp>
-
+#include "pavx.h"
+#include <xbob.blitz/cppapi.h>
 #include <bob/math/norminv.h>
 
-using namespace boost::python;
+PyObject* py_norminv (PyObject*, PyObject* args, PyObject* kwds) {
+
+  /* Parses input arguments in a single shot */
+  static const char* const_kwlist[] = { "p", "mu", "sigma", 0 /* Sentinel */ };
+  static char** kwlist = const_cast<char**>(const_kwlist);
+
+  double p = 0.;
+  double mu = 0.;
+  double sigma = 0.;
 
-static const char* NORMSINV_DOC = "Compute the inverse normal cumulative distribution for a probability p, given a distribution with zero mean and and unit variance.\nReference: http://home.online.no/~pjacklam/notes/invnorm/";
-static const char* NORMINV_DOC = "Compute the inverse normal cumulative distribution for a probability p, given a distribution with mean mu and standard deviation sigma.\nReference: http://home.online.no/~pjacklam/notes/invnorm/";
+  if (!PyArg_ParseTupleAndKeywords(args, kwds, "ddd", kwlist, &p, &mu, &sigma)) 
+    return 0;
+
+  try {
+    return PyBlitzArrayCxx_FromCScalar(bob::math::norminv(p, mu, sigma));
+  }
+  catch (std::exception& e) {
+    PyErr_SetString(PyExc_RuntimeError, e.what());
+  }
+  catch (...) {
+    PyErr_SetString(PyExc_RuntimeError, "norminv failed: unknown exception caught");
+  }
+
+  return 0;
 
-void bind_math_norminv()
-{
-  def("normsinv", &bob::math::normsinv, (arg("p")), NORMSINV_DOC);
-  def("norminv", &bob::math::norminv, (arg("p"), arg("mu"), arg("sigma")), NORMINV_DOC);
 }
 
+PyObject* py_normsinv (PyObject*, PyObject* args, PyObject* kwds) {
+
+  /* Parses input arguments in a single shot */
+  static const char* const_kwlist[] = { "p", 0 /* Sentinel */ };
+  static char** kwlist = const_cast<char**>(const_kwlist);
+
+  double p = 0.;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwds, "d", kwlist, &p)) return 0;
+
+  try {
+    return PyBlitzArrayCxx_FromCScalar(bob::math::normsinv(p));
+  }
+  catch (std::exception& e) {
+    PyErr_SetString(PyExc_RuntimeError, e.what());
+  }
+  catch (...) {
+    PyErr_SetString(PyExc_RuntimeError, "normsinv failed: unknown exception caught");
+  }
+
+  return 0;
+
+}
diff --git a/xbob/math/norminv.h b/xbob/math/norminv.h
new file mode 100644
index 0000000000000000000000000000000000000000..e740827585cb67416c091b286639d6d619154cc4
--- /dev/null
+++ b/xbob/math/norminv.h
@@ -0,0 +1,11 @@
+/**
+ * @author Andre Anjos <andre.anjos@idiap.ch>
+ * @date Thu  5 Dec 12:01:57 2013 
+ *
+ * @brief Declaration of components relevant for main.cpp
+ */
+
+#include <Python.h>
+
+PyObject* py_norminv(PyObject*, PyObject* args, PyObject* kwds);
+PyObject* py_normsinv(PyObject*, PyObject* args, PyObject* kwds);