Cython¶

Christian Elsasser (che@physik.uzh.ch)

The content of this lecture might be distributed under CC by-sa.

An introductionnary example - Fibonacci numbers¶

Pure Python¶

In [1]:
def fib0(n):
    a,b=1,1
    for i in range(n):
        a,b=a+b,a
    return a
In [2]:
fib0(5)
Out[2]:
13
In [3]:
fib0(6)
Out[3]:
21
In [4]:
%timeit fib0(10000)

1.86 ms ± 46.9 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Cython: Compile unmodified¶

In [5]:
import cython
In [6]:
%load_ext Cython
In [7]:
%%cython
def helloWorld():
    print("Hello, World!")
In [8]:
helloWorld
Out[8]:
<function _cython_magic_118d8ca079ebf334845d8489efb17859.helloWorld>
In [9]:
helloWorld()

Hello, World!
In [10]:
%%cython
def fib1(n):
    a,b=1,1
    for i in range(n):
        a,b=a+b,a
    return a
In [11]:
fib1
Out[11]:
<function _cython_magic_21f003b0a40db19747c9753af9912dc5.fib1>
In [12]:
fib0
Out[12]:
<function __main__.fib0(n)>
In [13]:
fib1(5)
Out[13]:
13
In [14]:
%timeit fib1(10000)

1.7 ms ± 17.4 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Cython: Declaring Types¶

In [15]:
%%cython
def fib2(int n):
    cdef int a,b,i
    a,b=1,1
    for i in range(n):
        a,b=a+b,a
    return a
In [16]:
fib2
Out[16]:
<function _cython_magic_ca13a30b54ad3711d7532bab9c6ca7ab.fib2>
In [17]:
fib2(5)
Out[17]:
13
In [18]:
%timeit fib2(10000)

2.38 µs ± 64.2 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Impressive, but¶

In [19]:
print("python result: ", fib0(500))
print("cython result: ", fib2(500))

python result:  365014740723634211012237077906479355996081581501455497852747829366800199361550174096573645929019489792751
cython result:  -1583205649
In [20]:
print("python result: ", fib0(40))
print("cython result: ", fib2(40))

python result:  267914296
cython result:  267914296
In [21]:
%timeit fib0(40)
%timeit fib2(40)

1.67 µs ± 29.5 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
47.7 ns ± 0.538 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)

Use better data types¶

In [22]:
%%cython
def fib3(int n):
    cdef long a,b
    cdef int i
    a,b=1,1
    for i in range(n):
        a,b=a+b,a
    return a
In [23]:
print("python result: ", fib0(50))
print("cython result: ", fib3(50))

python result:  32951280099
cython result:  32951280099
In [24]:
%timeit fib3(40)

46.5 ns ± 1.76 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)

Conclusion: The speedup is real, but we pay with new potential errors.

Second Example - Numerical integration¶

start with declaring types¶

In [25]:
%%cython -a

from math import sin,exp

def f0(double x):
    return sin(x)*exp(-x)

def integrate0(double a,double b,int N):
    cdef double dx,x,s
    cdef int _
    dx = (b-a)/N
    x = a+dx/2
    s = 0.0
    for _ in range(N):
        s+=f0(x)
        x += dx
    return s*dx
Out[25]:
Cython: _cython_magic_2585458f5223e3779c320efdcd5fa44b.pyx

Generated by Cython 0.29.30

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it. 

 01: 
+02: from math import sin,exp
  __pyx_t_1 = PyList_New(2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_INCREF(__pyx_n_s_sin);
  __Pyx_GIVEREF(__pyx_n_s_sin);
  PyList_SET_ITEM(__pyx_t_1, 0, __pyx_n_s_sin);
  __Pyx_INCREF(__pyx_n_s_exp);
  __Pyx_GIVEREF(__pyx_n_s_exp);
  PyList_SET_ITEM(__pyx_t_1, 1, __pyx_n_s_exp);
  __pyx_t_2 = __Pyx_Import(__pyx_n_s_math, __pyx_t_1, 0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_ImportFrom(__pyx_t_2, __pyx_n_s_sin); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_sin, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_ImportFrom(__pyx_t_2, __pyx_n_s_exp); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_exp, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;

 03: 
+04: def f0(double x):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_1f0(PyObject *__pyx_self, PyObject *__pyx_arg_x); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_1f0 = {"f0", (PyCFunction)__pyx_pw_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_1f0, METH_O, 0};
static PyObject *__pyx_pw_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_1f0(PyObject *__pyx_self, PyObject *__pyx_arg_x) {
  double __pyx_v_x;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("f0 (wrapper)", 0);
  assert(__pyx_arg_x); {
    __pyx_v_x = __pyx_PyFloat_AsDouble(__pyx_arg_x); if (unlikely((__pyx_v_x == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 4, __pyx_L3_error)
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_2585458f5223e3779c320efdcd5fa44b.f0", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_f0(__pyx_self, ((double)__pyx_v_x));
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_f0(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_x) {
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("f0", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_AddTraceback("_cython_magic_2585458f5223e3779c320efdcd5fa44b.f0", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple_ = PyTuple_Pack(2, __pyx_n_s_x, __pyx_n_s_x); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 4, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple_);
  __Pyx_GIVEREF(__pyx_tuple_);
/* … */
  __pyx_t_2 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_1f0, NULL, __pyx_n_s_cython_magic_2585458f5223e3779c); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 4, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_f0, __pyx_t_2) < 0) __PYX_ERR(0, 4, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_codeobj__2 = (PyObject*)__Pyx_PyCode_New(1, 0, 2, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple_, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_che_cache_ipython_cython, __pyx_n_s_f0, 4, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__2)) __PYX_ERR(0, 4, __pyx_L1_error)

+05:     return sin(x)*exp(-x)
  __Pyx_XDECREF(__pyx_r);
  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_sin); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = PyFloat_FromDouble(__pyx_v_x); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_2))) {
    __pyx_t_4 = PyMethod_GET_SELF(__pyx_t_2);
    if (likely(__pyx_t_4)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2);
      __Pyx_INCREF(__pyx_t_4);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_2, function);
    }
  }
  __pyx_t_1 = (__pyx_t_4) ? __Pyx_PyObject_Call2Args(__pyx_t_2, __pyx_t_4, __pyx_t_3) : __Pyx_PyObject_CallOneArg(__pyx_t_2, __pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_exp); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = PyFloat_FromDouble((-__pyx_v_x)); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_3))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_3);
    if (likely(__pyx_t_5)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
      __Pyx_INCREF(__pyx_t_5);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_3, function);
    }
  }
  __pyx_t_2 = (__pyx_t_5) ? __Pyx_PyObject_Call2Args(__pyx_t_3, __pyx_t_5, __pyx_t_4) : __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = PyNumber_Multiply(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_r = __pyx_t_3;
  __pyx_t_3 = 0;
  goto __pyx_L0;

 06: 
+07: def integrate0(double a,double b,int N):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_3integrate0(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_3integrate0 = {"integrate0", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_3integrate0, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_3integrate0(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
  double __pyx_v_a;
  double __pyx_v_b;
  int __pyx_v_N;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("integrate0 (wrapper)", 0);
  {
    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_a,&__pyx_n_s_b,&__pyx_n_s_N,0};
    PyObject* values[3] = {0,0,0};
    if (unlikely(__pyx_kwds)) {
      Py_ssize_t kw_args;
      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
      switch (pos_args) {
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        CYTHON_FALLTHROUGH;
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        CYTHON_FALLTHROUGH;
        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = PyDict_Size(__pyx_kwds);
      switch (pos_args) {
        case  0:
        if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_a)) != 0)) kw_args--;
        else goto __pyx_L5_argtuple_error;
        CYTHON_FALLTHROUGH;
        case  1:
        if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_b)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("integrate0", 1, 3, 3, 1); __PYX_ERR(0, 7, __pyx_L3_error)
        }
        CYTHON_FALLTHROUGH;
        case  2:
        if (likely((values[2] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_N)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("integrate0", 1, 3, 3, 2); __PYX_ERR(0, 7, __pyx_L3_error)
        }
      }
      if (unlikely(kw_args > 0)) {
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "integrate0") < 0)) __PYX_ERR(0, 7, __pyx_L3_error)
      }
    } else if (PyTuple_GET_SIZE(__pyx_args) != 3) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
      values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
      values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
    }
    __pyx_v_a = __pyx_PyFloat_AsDouble(values[0]); if (unlikely((__pyx_v_a == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
    __pyx_v_b = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_b == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
    __pyx_v_N = __Pyx_PyInt_As_int(values[2]); if (unlikely((__pyx_v_N == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("integrate0", 1, 3, 3, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 7, __pyx_L3_error)
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_2585458f5223e3779c320efdcd5fa44b.integrate0", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_2integrate0(__pyx_self, __pyx_v_a, __pyx_v_b, __pyx_v_N);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_2integrate0(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_a, double __pyx_v_b, int __pyx_v_N) {
  double __pyx_v_dx;
  double __pyx_v_x;
  double __pyx_v_s;
  CYTHON_UNUSED int __pyx_v__;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("integrate0", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_XDECREF(__pyx_t_6);
  __Pyx_XDECREF(__pyx_t_7);
  __Pyx_XDECREF(__pyx_t_8);
  __Pyx_XDECREF(__pyx_t_9);
  __Pyx_AddTraceback("_cython_magic_2585458f5223e3779c320efdcd5fa44b.integrate0", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple__4 = PyTuple_Pack(7, __pyx_n_s_a, __pyx_n_s_b, __pyx_n_s_N, __pyx_n_s_dx, __pyx_n_s_x, __pyx_n_s_s, __pyx_n_s__3); if (unlikely(!__pyx_tuple__4)) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__4);
  __Pyx_GIVEREF(__pyx_tuple__4);
/* … */
  __pyx_t_2 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_2585458f5223e3779c320efdcd5fa44b_3integrate0, NULL, __pyx_n_s_cython_magic_2585458f5223e3779c); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_integrate0, __pyx_t_2) < 0) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;

 08:     cdef double dx,x,s
 09:     cdef int _
+10:     dx = (b-a)/N
  __pyx_t_1 = (__pyx_v_b - __pyx_v_a);
  if (unlikely(__pyx_v_N == 0)) {
    PyErr_SetString(PyExc_ZeroDivisionError, "float division");
    __PYX_ERR(0, 10, __pyx_L1_error)
  }
  __pyx_v_dx = (__pyx_t_1 / ((double)__pyx_v_N));

+11:     x = a+dx/2
  __pyx_v_x = (__pyx_v_a + (__pyx_v_dx / 2.0));

+12:     s = 0.0
  __pyx_v_s = 0.0;

+13:     for _ in range(N):
  __pyx_t_2 = __pyx_v_N;
  __pyx_t_3 = __pyx_t_2;
  for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
    __pyx_v__ = __pyx_t_4;

+14:         s+=f0(x)
    __pyx_t_5 = PyFloat_FromDouble(__pyx_v_s); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 14, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_GetModuleGlobalName(__pyx_t_7, __pyx_n_s_f0); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 14, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_7);
    __pyx_t_8 = PyFloat_FromDouble(__pyx_v_x); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 14, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_8);
    __pyx_t_9 = NULL;
    if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_7))) {
      __pyx_t_9 = PyMethod_GET_SELF(__pyx_t_7);
      if (likely(__pyx_t_9)) {
        PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_7);
        __Pyx_INCREF(__pyx_t_9);
        __Pyx_INCREF(function);
        __Pyx_DECREF_SET(__pyx_t_7, function);
      }
    }
    __pyx_t_6 = (__pyx_t_9) ? __Pyx_PyObject_Call2Args(__pyx_t_7, __pyx_t_9, __pyx_t_8) : __Pyx_PyObject_CallOneArg(__pyx_t_7, __pyx_t_8);
    __Pyx_XDECREF(__pyx_t_9); __pyx_t_9 = 0;
    __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
    if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 14, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_6);
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
    __pyx_t_7 = PyNumber_InPlaceAdd(__pyx_t_5, __pyx_t_6); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 14, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_7);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
    __pyx_t_1 = __pyx_PyFloat_AsDouble(__pyx_t_7); if (unlikely((__pyx_t_1 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 14, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
    __pyx_v_s = __pyx_t_1;

+15:         x += dx
    __pyx_v_x = (__pyx_v_x + __pyx_v_dx);
  }

+16:     return s*dx
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_7 = PyFloat_FromDouble((__pyx_v_s * __pyx_v_dx)); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 16, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_7);
  __pyx_r = __pyx_t_7;
  __pyx_t_7 = 0;
  goto __pyx_L0;
In [26]:
from math import pi
In [27]:
integrate0(0.0,pi,1000000)
Out[27]:
0.5216069591343814
In [28]:
%timeit integrate0(0.0,pi,1000000)

94.1 ms ± 3.14 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

turn f into a c-function¶

In [40]:
%%cython -a

from math import sin,exp

cdef double f1(double x):
    return sin(x)*exp(-x)

def integrate1(double a,double b,int N):
    cdef double dx,s
    cdef int i
    dx = (b-a)/N
    s = 0.0
    for i in range(N):
        s+=f1(a+(i+0.5)*dx)
    return s*dx
Out[40]:
Cython: _cython_magic_fdd88313f639ed94b7566857525b5d73.pyx

Generated by Cython 0.29.30

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it. 

 01: 
+02: from math import sin,exp
  __pyx_t_1 = PyList_New(2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_INCREF(__pyx_n_s_sin);
  __Pyx_GIVEREF(__pyx_n_s_sin);
  PyList_SET_ITEM(__pyx_t_1, 0, __pyx_n_s_sin);
  __Pyx_INCREF(__pyx_n_s_exp);
  __Pyx_GIVEREF(__pyx_n_s_exp);
  PyList_SET_ITEM(__pyx_t_1, 1, __pyx_n_s_exp);
  __pyx_t_2 = __Pyx_Import(__pyx_n_s_math, __pyx_t_1, 0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_ImportFrom(__pyx_t_2, __pyx_n_s_sin); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_sin, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_ImportFrom(__pyx_t_2, __pyx_n_s_exp); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_exp, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;

 03: 
+04: cdef double f1(double x):
static double __pyx_f_46_cython_magic_fdd88313f639ed94b7566857525b5d73_f1(double __pyx_v_x) {
  double __pyx_r;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("f1", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_WriteUnraisable("_cython_magic_fdd88313f639ed94b7566857525b5d73.f1", __pyx_clineno, __pyx_lineno, __pyx_filename, 1, 0);
  __pyx_r = 0;
  __pyx_L0:;
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

+05:     return sin(x)*exp(-x)
  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_sin); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = PyFloat_FromDouble(__pyx_v_x); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_2))) {
    __pyx_t_4 = PyMethod_GET_SELF(__pyx_t_2);
    if (likely(__pyx_t_4)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2);
      __Pyx_INCREF(__pyx_t_4);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_2, function);
    }
  }
  __pyx_t_1 = (__pyx_t_4) ? __Pyx_PyObject_Call2Args(__pyx_t_2, __pyx_t_4, __pyx_t_3) : __Pyx_PyObject_CallOneArg(__pyx_t_2, __pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_exp); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = PyFloat_FromDouble((-__pyx_v_x)); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_3))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_3);
    if (likely(__pyx_t_5)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
      __Pyx_INCREF(__pyx_t_5);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_3, function);
    }
  }
  __pyx_t_2 = (__pyx_t_5) ? __Pyx_PyObject_Call2Args(__pyx_t_3, __pyx_t_5, __pyx_t_4) : __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = PyNumber_Multiply(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_t_3); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_r = __pyx_t_6;
  goto __pyx_L0;

 06: 
+07: def integrate1(double a,double b,int N):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_fdd88313f639ed94b7566857525b5d73_1integrate1(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_fdd88313f639ed94b7566857525b5d73_1integrate1 = {"integrate1", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_46_cython_magic_fdd88313f639ed94b7566857525b5d73_1integrate1, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_fdd88313f639ed94b7566857525b5d73_1integrate1(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
  double __pyx_v_a;
  double __pyx_v_b;
  int __pyx_v_N;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("integrate1 (wrapper)", 0);
  {
    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_a,&__pyx_n_s_b,&__pyx_n_s_N,0};
    PyObject* values[3] = {0,0,0};
    if (unlikely(__pyx_kwds)) {
      Py_ssize_t kw_args;
      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
      switch (pos_args) {
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        CYTHON_FALLTHROUGH;
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        CYTHON_FALLTHROUGH;
        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = PyDict_Size(__pyx_kwds);
      switch (pos_args) {
        case  0:
        if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_a)) != 0)) kw_args--;
        else goto __pyx_L5_argtuple_error;
        CYTHON_FALLTHROUGH;
        case  1:
        if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_b)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("integrate1", 1, 3, 3, 1); __PYX_ERR(0, 7, __pyx_L3_error)
        }
        CYTHON_FALLTHROUGH;
        case  2:
        if (likely((values[2] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_N)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("integrate1", 1, 3, 3, 2); __PYX_ERR(0, 7, __pyx_L3_error)
        }
      }
      if (unlikely(kw_args > 0)) {
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "integrate1") < 0)) __PYX_ERR(0, 7, __pyx_L3_error)
      }
    } else if (PyTuple_GET_SIZE(__pyx_args) != 3) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
      values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
      values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
    }
    __pyx_v_a = __pyx_PyFloat_AsDouble(values[0]); if (unlikely((__pyx_v_a == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
    __pyx_v_b = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_b == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
    __pyx_v_N = __Pyx_PyInt_As_int(values[2]); if (unlikely((__pyx_v_N == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("integrate1", 1, 3, 3, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 7, __pyx_L3_error)
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_fdd88313f639ed94b7566857525b5d73.integrate1", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_46_cython_magic_fdd88313f639ed94b7566857525b5d73_integrate1(__pyx_self, __pyx_v_a, __pyx_v_b, __pyx_v_N);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_fdd88313f639ed94b7566857525b5d73_integrate1(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_a, double __pyx_v_b, int __pyx_v_N) {
  double __pyx_v_dx;
  double __pyx_v_s;
  int __pyx_v_i;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("integrate1", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_AddTraceback("_cython_magic_fdd88313f639ed94b7566857525b5d73.integrate1", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple_ = PyTuple_Pack(6, __pyx_n_s_a, __pyx_n_s_b, __pyx_n_s_N, __pyx_n_s_dx, __pyx_n_s_s, __pyx_n_s_i); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple_);
  __Pyx_GIVEREF(__pyx_tuple_);
/* … */
  __pyx_t_2 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_fdd88313f639ed94b7566857525b5d73_1integrate1, NULL, __pyx_n_s_cython_magic_fdd88313f639ed94b7); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_integrate1, __pyx_t_2) < 0) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;

 08:     cdef double dx,s
 09:     cdef int i
+10:     dx = (b-a)/N
  __pyx_t_1 = (__pyx_v_b - __pyx_v_a);
  if (unlikely(__pyx_v_N == 0)) {
    PyErr_SetString(PyExc_ZeroDivisionError, "float division");
    __PYX_ERR(0, 10, __pyx_L1_error)
  }
  __pyx_v_dx = (__pyx_t_1 / ((double)__pyx_v_N));

+11:     s = 0.0
  __pyx_v_s = 0.0;

+12:     for i in range(N):
  __pyx_t_2 = __pyx_v_N;
  __pyx_t_3 = __pyx_t_2;
  for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
    __pyx_v_i = __pyx_t_4;

+13:         s+=f1(a+(i+0.5)*dx)
    __pyx_v_s = (__pyx_v_s + __pyx_f_46_cython_magic_fdd88313f639ed94b7566857525b5d73_f1((__pyx_v_a + ((__pyx_v_i + 0.5) * __pyx_v_dx))));
  }

+14:     return s*dx
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_5 = PyFloat_FromDouble((__pyx_v_s * __pyx_v_dx)); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 14, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_r = __pyx_t_5;
  __pyx_t_5 = 0;
  goto __pyx_L0;
In [30]:
%timeit integrate1(0.0,pi,1000000)

68.2 ms ± 797 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

f is now not visible in Python anymore

In [31]:
f1(5)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [31], in <cell line: 1>()
----> 1 f1(5)

NameError: name 'f1' is not defined

use cpdef to keep f in Python

use math functions from c¶

In [32]:
%%cython -a

from libc.math cimport sin,exp

cpdef double f2(double x):
    return sin(x)*exp(-x)

def integrate2(double a,double b,int N):
    cdef double dx,s
    cdef int i
    dx = (b-a)/N
    s = 0.0
    for i in range(N):
        s+=f2(a+(i+0.5)*dx)
    return s*dx
Out[32]:
Cython: _cython_magic_f3a5d4af9f8f3f5c037181c3247ca121.pyx

Generated by Cython 0.29.30

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it. 

 01: 
 02: from libc.math cimport sin,exp
 03: 
+04: cpdef double f2(double x):
static PyObject *__pyx_pw_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_1f2(PyObject *__pyx_self, PyObject *__pyx_arg_x); /*proto*/
static double __pyx_f_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_f2(double __pyx_v_x, CYTHON_UNUSED int __pyx_skip_dispatch) {
  double __pyx_r;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("f2", 0);
/* … */
  /* function exit code */
  __pyx_L0:;
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_1f2(PyObject *__pyx_self, PyObject *__pyx_arg_x); /*proto*/
static PyObject *__pyx_pw_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_1f2(PyObject *__pyx_self, PyObject *__pyx_arg_x) {
  double __pyx_v_x;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("f2 (wrapper)", 0);
  assert(__pyx_arg_x); {
    __pyx_v_x = __pyx_PyFloat_AsDouble(__pyx_arg_x); if (unlikely((__pyx_v_x == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 4, __pyx_L3_error)
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121.f2", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_f2(__pyx_self, ((double)__pyx_v_x));
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_f2(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_x) {
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("f2", 0);
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_1 = PyFloat_FromDouble(__pyx_f_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_f2(__pyx_v_x, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 4, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_r = __pyx_t_1;
  __pyx_t_1 = 0;
  goto __pyx_L0;

  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_AddTraceback("_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121.f2", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

+05:     return sin(x)*exp(-x)
  __pyx_r = (sin(__pyx_v_x) * exp((-__pyx_v_x)));
  goto __pyx_L0;

 06: 
+07: def integrate2(double a,double b,int N):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_3integrate2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_3integrate2 = {"integrate2", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_3integrate2, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_3integrate2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
  double __pyx_v_a;
  double __pyx_v_b;
  int __pyx_v_N;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("integrate2 (wrapper)", 0);
  {
    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_a,&__pyx_n_s_b,&__pyx_n_s_N,0};
    PyObject* values[3] = {0,0,0};
    if (unlikely(__pyx_kwds)) {
      Py_ssize_t kw_args;
      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
      switch (pos_args) {
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        CYTHON_FALLTHROUGH;
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        CYTHON_FALLTHROUGH;
        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = PyDict_Size(__pyx_kwds);
      switch (pos_args) {
        case  0:
        if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_a)) != 0)) kw_args--;
        else goto __pyx_L5_argtuple_error;
        CYTHON_FALLTHROUGH;
        case  1:
        if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_b)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("integrate2", 1, 3, 3, 1); __PYX_ERR(0, 7, __pyx_L3_error)
        }
        CYTHON_FALLTHROUGH;
        case  2:
        if (likely((values[2] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_N)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("integrate2", 1, 3, 3, 2); __PYX_ERR(0, 7, __pyx_L3_error)
        }
      }
      if (unlikely(kw_args > 0)) {
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "integrate2") < 0)) __PYX_ERR(0, 7, __pyx_L3_error)
      }
    } else if (PyTuple_GET_SIZE(__pyx_args) != 3) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
      values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
      values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
    }
    __pyx_v_a = __pyx_PyFloat_AsDouble(values[0]); if (unlikely((__pyx_v_a == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
    __pyx_v_b = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_b == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
    __pyx_v_N = __Pyx_PyInt_As_int(values[2]); if (unlikely((__pyx_v_N == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("integrate2", 1, 3, 3, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 7, __pyx_L3_error)
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121.integrate2", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_2integrate2(__pyx_self, __pyx_v_a, __pyx_v_b, __pyx_v_N);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_2integrate2(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_a, double __pyx_v_b, int __pyx_v_N) {
  double __pyx_v_dx;
  double __pyx_v_s;
  int __pyx_v_i;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("integrate2", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_AddTraceback("_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121.integrate2", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple_ = PyTuple_Pack(6, __pyx_n_s_a, __pyx_n_s_b, __pyx_n_s_N, __pyx_n_s_dx, __pyx_n_s_s, __pyx_n_s_i); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple_);
  __Pyx_GIVEREF(__pyx_tuple_);
/* … */
  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_3integrate2, NULL, __pyx_n_s_cython_magic_f3a5d4af9f8f3f5c03); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_integrate2, __pyx_t_1) < 0) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 08:     cdef double dx,s
 09:     cdef int i
+10:     dx = (b-a)/N
  __pyx_t_1 = (__pyx_v_b - __pyx_v_a);
  if (unlikely(__pyx_v_N == 0)) {
    PyErr_SetString(PyExc_ZeroDivisionError, "float division");
    __PYX_ERR(0, 10, __pyx_L1_error)
  }
  __pyx_v_dx = (__pyx_t_1 / ((double)__pyx_v_N));

+11:     s = 0.0
  __pyx_v_s = 0.0;

+12:     for i in range(N):
  __pyx_t_2 = __pyx_v_N;
  __pyx_t_3 = __pyx_t_2;
  for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
    __pyx_v_i = __pyx_t_4;

+13:         s+=f2(a+(i+0.5)*dx)
    __pyx_v_s = (__pyx_v_s + __pyx_f_46_cython_magic_f3a5d4af9f8f3f5c037181c3247ca121_f2((__pyx_v_a + ((__pyx_v_i + 0.5) * __pyx_v_dx)), 0));
  }

+14:     return s*dx
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_5 = PyFloat_FromDouble((__pyx_v_s * __pyx_v_dx)); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 14, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_r = __pyx_t_5;
  __pyx_t_5 = 0;
  goto __pyx_L0;
In [33]:
%timeit integrate2(0.0,pi,1000000)

10.5 ms ± 232 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Line 10 is yellow, due to checking for ZeroDivisionError. We can deactivate this with

import cython

@cython.cdivision(True)
def integrate(...

But:

  • we would again open another door for errors
  • this line is only executed once

Was it worth it?¶

In [34]:
import numpy as np

def f3(x):
    return np.sin(x)*np.exp(-x)

def integrate3(a,b,N):
    dx = (b-a)/N
    x = (np.arange(N)+0.5)*dx+a
    return np.sum(f3(x))*dx
In [35]:
%timeit integrate3(0.0,pi,1000000)

20 ms ± 422 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Cython and STL containers¶

In [42]:
%%cython -+
#distutils: language = c++
from libcpp.vector cimport vector
def manipulateVector(vector[double] v):
    v.push_back(42.0) # Equivalent to append for list in python
    return v
In [43]:
manipulateVector([1,2,3,4])
Out[43]:
[1.0, 2.0, 3.0, 4.0, 42.0]

Cython and exceptions¶

In [44]:
%%cython
cpdef int raiseError():
    raise RuntimeError("A problem!")
    return 1
In [45]:
raiseError()
print("Still going strong!")

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
RuntimeError: A problem!
Exception ignored in: '_cython_magic_edc180d4d46526990c4f883521624a73.raiseError'
Traceback (most recent call last):
  File "/var/folders/z1/pxv05sfd3l59nrf3m04qfz9c0000gn/T/ipykernel_5673/3788160389.py", line 1, in <cell line: 1>
RuntimeError: A problem!

Still going strong!
In [46]:
%%cython
cpdef int raiseErrorBetter() except *:
    raise RuntimeError("A problem!")
    return 1
In [47]:
raiseErrorBetter()
print("Still going strong!")

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Input In [47], in <cell line: 1>()
----> 1 raiseErrorBetter()
      2 print("Still going strong!")

File _cython_magic_ff2685e2d9bb4fb65879ee8d14905368.pyx:1, in _cython_magic_ff2685e2d9bb4fb65879ee8d14905368.raiseErrorBetter()

File _cython_magic_ff2685e2d9bb4fb65879ee8d14905368.pyx:2, in _cython_magic_ff2685e2d9bb4fb65879ee8d14905368.raiseErrorBetter()

RuntimeError: A problem!
In [48]:
try:
    raiseErrorBetter()
except RuntimeError:
    print("Error caught")

Error caught
In [49]:
%%cython
cpdef int raiseException(int n) except -1:
    # Assuming C function returning an integer; if -1 --> Error appeared
    return n
In [50]:
raiseException(1)
Out[50]:
1
In [51]:
raiseException(-1)

---------------------------------------------------------------------------
SystemError                               Traceback (most recent call last)
Input In [51], in <cell line: 1>()
----> 1 raiseException(-1)

SystemError: <built-in function raiseException> returned NULL without setting an exception

Back to integration - Cython and Classes¶

In [52]:
%%cython

from libc.math cimport sin,exp

cdef class Integrand(object):
    cpdef double evaluate(self,double x):
        raise NotImplementedError()
        
cdef class SinExpFunction(Integrand):
    cpdef double evaluate(self,double x):
        return sin(x)*exp(-x)
    
def integrate(Integrand f, double a, double b, int N):
    cdef int i
    cdef double dx,s
    dx = (b-a)/N
    s = 0.0
    for i in range(N):
        s += f.evaluate(a+(i+0.5)*dx)
    return s*dx
In [53]:
f = SinExpFunction()
In [54]:
import math
In [55]:
%timeit integrate(f,0,math.pi,1000000)

13.8 ms ± 295 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [56]:
class Poly(Integrand):
    def evaluate(self,x):
        return 3*x-x**3
In [57]:
g = Poly()
integrate(g,0,math.pi,10000)
Out[57]:
-9.5478660351052
In [58]:
%timeit integrate(g,0,math.pi,1000000)

191 ms ± 8.65 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

C function wrapping¶

In [ ]:
# %load fastinvsqrt/fis.c
#include <stdio.h>

double fast_inv_sqrt( double number )
  {
      double y = number;
      double x2 = y * 0.5;
      long long i = *(long long*) &y;
      // The magic number is for doubles is from https://cs.uwaterloo.ca/~m32rober/rsqrt.pdf
      i = 0x5fe6eb50c7b537a9 - (i >> 1);
      y = *(double *) &i;
      y = y * (1.5 - (x2 * y * y));   // 1st iteration
      //      y  = y * ( 1.5 - ( x2 * y * y ) );   // 2nd iteration, this can be removed
      return y;
  }
In [60]:
%%cython -I.
# -I followed by include directory i.e. root directory where the files specified 
cdef extern from "fastinvsqrt/fis.c":
    double fast_inv_sqrt(double number)

def py_fis(number:double) -> double:
    return fast_inv_sqrt(number)

def norm_vector(values:list) -> list:
    length_squared = sum([x**2 for x in values])
    return [x*fast_inv_sqrt(length_squared) for x in values]
In [61]:
norm_vector([1,2,3])
Out[61]:
[0.26721425303629537, 0.5344285060725907, 0.8016427591088862]
In [62]:
py_fis(102)
Out[62]:
0.09885709754491949
In [ ]:
%%cython?