How to wrap C code in Python (with a few fixed bugs)

_chi2.c

#include <Python.h>
#include <numpy/arrayobject.h>
#include "chi2.h"

/* Docstrings */
static char module_docstring[] =
    "This module provides an interface for calculating chi-squared using C.";
static char chi2_docstring[] =
    "Calculate the chi-squared of some data given a model.";

/* Available functions */
static PyObject *chi2_chi2(PyObject *self, PyObject *args);

/* Module specification */
static PyMethodDef module_methods[] = {
    {"chi2", chi2_chi2, METH_VARARGS, chi2_docstring},
    {NULL, NULL, 0, NULL}
};

/* Initialize the module */
PyMODINIT_FUNC init_chi2(void)
{
    PyObject *m = Py_InitModule3("_chi2", module_methods, module_docstring);
    if (m == NULL)
        return;

    /* Load `numpy` functionality. */
    import_array();
}

static PyObject *chi2_chi2(PyObject *self, PyObject *args)
{
    double m, b;
    PyObject *x_obj, *y_obj, *yerr_obj;
    PyObject *x_array, *y_array, *yerr_array;
    int N;
    double *x, *y, *yerr;
    double value;
    PyObject *ret;

    /* Parse the input tuple */
    if (!PyArg_ParseTuple(args, "ddOOO", &m, &b, &x_obj, &y_obj,
                                         &yerr_obj))
        return NULL;

    /* Interpret the input objects as numpy arrays. */
    x_array = PyArray_FROM_OTF(x_obj, NPY_DOUBLE, NPY_IN_ARRAY);
    y_array = PyArray_FROM_OTF(y_obj, NPY_DOUBLE, NPY_IN_ARRAY);
    yerr_array = PyArray_FROM_OTF(yerr_obj, NPY_DOUBLE,
                                            NPY_IN_ARRAY);

    /* If that didn't work, throw an exception. */
    if (x_array == NULL || y_array == NULL || yerr_array == NULL) {
        Py_XDECREF(x_array);
        Py_XDECREF(y_array);
        Py_XDECREF(yerr_array);
        return NULL;
    }

    /* How many data points are there? */
    N = (int)PyArray_DIM(x_array, 0);

    /* Get pointers to the data as C-types. */
    x    = (double*)PyArray_DATA(x_array);
    y    = (double*)PyArray_DATA(y_array);
    yerr = (double*)PyArray_DATA(yerr_array);

    /* Call the external C function to compute the chi-squared. */
    value = chi2(m, b, x, y, yerr, N);

    /* Clean up. */
    Py_DECREF(x_array);
    Py_DECREF(y_array);
    Py_DECREF(yerr_array);

    if (value < 0.0) {
        PyErr_SetString(PyExc_RuntimeError,
                    "Chi-squared returned an impossible value.");
        return NULL;
    }

    /* Build the output tuple */
    ret = Py_BuildValue("d", value);
    return ret;
}

chi2.c

#include "chi2.h"

double chi2(double m, double b, double *x, double *y, double *yerr, int N) {
    int n;
    double result = 0.0, diff;

    for (n = 0; n < N; n++) {
        diff = (y[n] - (m * x[n] + b)) / yerr[n];
        result += diff * diff;
    }

    return result;
}

chi2.h

double chi2(double m, double b, double *x, double *y, double *yerr, int N);

setup.py

try:
    from setuptools import setup
    from setuptools import Extension
except ImportError as E:
    print("{}. Using distutils instead.".format(E))
    from distutils.core import setup
    from distutils.extension import Extension

#from distutils.core import setup, Extension
import numpy.distutils.misc_util

c_ext = Extension("_chi2", ["_chi2.c", "chi2.c"])

setup(
    ext_modules=[c_ext],
    include_dirs=numpy.distutils.misc_util.get_numpy_include_dirs(),
)

test.py

import _chi2
result = _chi2.chi2(
    2.0, 1.0,
    [-1.0, 4.2, 30.6],
    [-1.5, 8.0, 63.0],
    [1.0, 1.5, 0.6]
)
print("Result: {} (expected: 2.89888)".format(result))