add cquad_wrapper and adapt makefile

Makefile

-CFLAGS = -lm -g
-
-SRC = $(wildcard *.c)
-OBJ = $(SRC:.c=.o)
+CFLAGS = -g -O2 -I/usr/include/python2.7
 
 .PHONY: run_test
 run_test: test
 	./test
 
-test: $(OBJ)
-	gcc $(CFLAGS) $(OBJ) -o test
+test: test.o tests.c cquad.o cquad_const.o
+	gcc $(CFLAGS) $^ -lm -o test
+
+cquad_wrapper.so: cquad_wrapper.o cquad.o cquad_const.o
+	gcc -shared -pthread -fPIC -fwrapv -fno-strict-aliasing $^ -lm -lstdc++ -o cquad_wrapper.so
+
+%.o: %.c
+	gcc $(CFLAGS) -fPIC -c -Wall -Wextra -o $@ $<
+
+%.o: %.cpp
+	g++ $(CFLAGS) -fPIC -c -Wall -Wextra -o $@ $<
 
-%.o:    %.c
-	gcc $(CFLAGS) -c -Wall -Wextra -o $@ $<
+%.cpp: %.pyx
+	cython --cplus $<
 
 .PHONY: clean
 clean:
-	rm -f $(OBJ)
+	rm -f *.o *.so *.cpp
 	rm -f test
 
 .DELETE_ON_ERROR:

cquad_wrapper.pyx

+import kwant
+import numpy as np
+
+cimport cython
+cimport cython.operator as co
+cimport libcpp.vector
+cimport libcpp.map
+
+
+cdef extern from "gsl_integration.h":
+    ctypedef struct gsl_function:
+        double (* function) (double x, void* params)
+        void* params
+
+cdef extern from "gsl_integration.h":
+
+    # Data of a single interval
+    ctypedef struct gsl_integration_cquad_ival:
+        double a, b
+        double c[64]
+        double fx[33]
+        double igral, err
+        int depth, rdepth, ndiv
+
+    # The workspace is just a collection of intervals 
+    ctypedef struct gsl_integration_cquad_workspace:
+        size_t size
+        gsl_integration_cquad_ival *ivals
+        size_t *heap
+        int nivals
+        
+    gsl_integration_cquad_workspace* gsl_integration_cquad_workspace_alloc(size_t n)
+
+    void  gsl_integration_cquad_workspace_free(gsl_integration_cquad_workspace* w)
+
+    int  gsl_integration_cquad(gsl_function* f, double a, double b, double epsabs, double epsrel, gsl_integration_cquad_workspace* workspace, double* result, double* abserr, size_t* nevals)
+
+
+cdef gsl_integration_cquad_workspace* workspace
+
+# wrap function and cache
+cdef struct Entry:
+    double complex evaluated
+
+cdef libcpp.map.map[double, Entry] cache
+
+cdef double real_f(double e, void* ignore):
+    return func(e).real
+
+cdef double imag_f(double e, void* ignore):
+    return func(e).imag
+
+cdef double complex func(double x):
+    cdef Entry* cached
+    cdef libcpp.map.map[double, Entry].iterator cached_iter
+   
+    cached_iter = cache.find(x)
+        
+    if cached_iter == cache.end():
+
+        # Create a default-constructed entry at energy E.
+        cached = &cache[x]
+        cached[0].evaluated = py_func(x)
+
+    else:
+        cached = &co.dereference(cached_iter).second
+
+    return cached[0].evaluated
+
+cdef double complex complex_quad(double a, double b,
+                                 double epsabs, double epsrel, int limit):
+    cdef double real_res, imag_res, error
+
+    cdef gsl_function func
+    func.params = NULL
+
+    func.function = real_f
+    if gsl_integration_cquad(&func, a, b, epsabs, epsrel, workspace,
+                             &real_res, &error, NULL) != 0:
+        raise RuntimeError()
+
+    func.function = imag_f
+    if gsl_integration_cquad(&func, a, b, epsabs, epsrel, workspace,
+                             &imag_res, &error, NULL) != 0:
+        raise RuntimeError()
+
+    return complex(real_res, imag_res)
+
+
+def cquad(funcin, a, b, args=[], epsabs=1e-6, epsrel=1e-6, limit=1000, full_output=0):
+    global cache, workspace, py_func
+    py_func = lambda x: funcin(x, *args)
+    cache.clear()
+    
+    try:
+        workspace = gsl_integration_cquad_workspace_alloc(limit)
+        if workspace == NULL:
+            raise RuntimeError()
+
+        G = complex_quad(a, b, epsabs, epsrel, limit)
+        size = co.dereference(workspace).size
+        dicodico = [co.dereference(workspace).ivals[i] for i in xrange(size)]
+        heap = [co.dereference(workspace).ivals[i] for i in xrange(size)]
+        nivals = co.dereference(workspace).nivals
+        
+    finally:
+        gsl_integration_cquad_workspace_free(workspace)
+
+    if full_output:
+        return G, dicodico, heap, nivals, cache
+    else:
+        return G

gsl_integration.h

@@ -21,6 +21,11 @@
 #define __GSL_INTEGRATION_H__
 #include <stdlib.h>
 
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
 enum {
     GSL_SUCCESS = 0,
     GSL_EDOM = 1,               /* input domain error, e.g sqrt(-1) */
@@ -70,4 +75,8 @@ int gsl_integration_cquad(const gsl_function * f, double a, double b,
                           gsl_integration_cquad_workspace * ws,
                           double *result, double *abserr, size_t * nevals);
 
+#ifdef __cplusplus
+}
+#endif
+
 #endif /* __GSL_INTEGRATION_H__ */