2021-11-20

CYTHON-GSL

Python extension modules can be written in Cython a compiled language rather similar to Python itself, and capable to be significantly faster than standard Python code. Moreover, existing C and C++ libraries can be easily interfaced within a Cython module as shown in this note.

CythonGSL, for instance, provides a Cython interface for the GNU Scientific Library (GSL). The GSL can be installed on Windows from the MINGW64 package manager as follows:

   $ pacman -S mingw-w64-x86_64-gsl 

On Windows, I prefer to use a Python instance that supports the MSVC compiler, as described in one of my previous notes, where I have installed CythonGSL with PIP.

   $ pip install cythongsl

Find the roots of a second order polynomial

A root (solution) of a quadratic equation is the value of x for which the second order univariate polynomial a x²+b x + c is equal to zero, a,b,c are termed the coefficients.

A graphical representation will make it easier to understand this simple problem. Let's take for example the following quadratic equation 2x² + 4 x - 4 = 0, which has two non-trivial real solutions x≅2.7 and x≅0.8, the two points where our parabola intersects the x-axis (see the graph below).

Cython module

The following Cython code allows calling the external GSL function gsl_poly_complex_solve_quadratic from a Python module.

    # filename: quadratic.pyx
    from cython_gsl cimport *
    cdef double a
    cdef double b
    cdef double c

    def c_quadratic(a,b,c):
        cdef gsl_complex z0
        cdef gsl_complex z1
        cdef int n
        n = gsl_poly_complex_solve_quadratic(a, b, c, &z0, &z1)
        print '\n roots %d' % n
        print "\nz0 = %f + %f*i" % (GSL_REAL(z0),GSL_IMAG(z0))
        print "\nz1 = %f + %f*i" % (GSL_REAL(z1),GSL_IMAG(z1))
    

The C-compiler and the GSL libraries belong to a different set of instances, in order to allow the C-compiler linking the required static libraries I have moved MINGW64 gsl.a and gslcblas.a files into a new directory and renamed them as gsl.lib and gslcblas.lib.

In the setup.py, I specified the paths for the MINGW64's header files and the previous two .lib files.

    # filename:setup.py
    from distutils.core import setup
    from Cython.Distutils import Extension
    from Cython.Distutils import build_ext
    import cython_gsl

    setup(
        include_dirs = [
            cython_gsl.get_include(),
            "C:\\...\\mingw64\\include"
            ],
         cmdclass = {'build_ext': build_ext},
        ext_modules = [Extension("quadratic",
                                 ["quadratic.pyx"],
                                 libraries=cython_gsl.get_libraries(),
                                 library_dirs=[cython_gsl.get_library_dir(),
                                     "C:\\...\\cython_gsl\\lib"],
                                 include_dirs=[cython_gsl.get_cython_include_dir()])]
        )
    

Now, I can compile this Cython module with:

   $ python setup.py build

To import the new module in a Python REPL, I just need the following script:

    # filename: libpath.py
    import sys
    sys.path.append('./build/lib.win-amd64-3.9')
    import quadratic