title: Cython boundscheck & nonecheck
body: <p>I struggle for several day with a script giving me unexpected results. 
  Today, I just realize that if I used a cython function with or without <code>boundscheck</code> and <code>nonecheck</code> decorators, I do not obtain the same results !</p>
  
  <p>Here is an example : </p>
  
  <pre><code>import numpy as np
  cimport numpy as np
  cimport cython
  
  cdef double[4] c
  c[0] = 0.1
  c[1] = 0.2
  c[2] = 0.3
  c[3] = 0.4
  
  
  def cp1(double[:,::1] u, double[:,::1] K, int ixmin, int ixmax, int izmin, int izmax):
      cpc1(u, K, ixmin, ixmax, izmin, izmax)
  
  
  def cp2(double[:,::1] u, double[:,::1] K, int ixmin, int ixmax, int izmin, int izmax):
      cpc2(u, K, ixmin, ixmax, izmin, izmax)
  
  
  @cython.boundscheck(False)
  @cython.nonecheck(False)
  cdef void cpc1(double[:,::1] u, double[:,::1] K, int ixmin, int ixmax, int izmin, int izmax) nogil:
  
      cdef Py_ssize_t ix, iz
      cdef double dpu, dmu
  
      for ix in range(ixmin+2, ixmax-1):
          for iz in range(izmin, izmax):
              dpu = c[0]*u[ix-1, iz] + c[1]*u[ix, iz] + c[2]*u[ix+1, iz]
              dmu = c[1]*u[ix-1, iz] + c[2]*u[ix, iz] + c[3]*u[ix+1, iz]
              K[ix, iz] = 0.5*dpu - 0.5*dmu
  
  @cython.boundscheck(True)
  @cython.nonecheck(True)
  cdef void cpc2(double[:,::1] u, double[:,::1] K, int ixmin, int ixmax, int izmin, int izmax) nogil:
  
      cdef Py_ssize_t ix, iz
      cdef double dpu, dmu
  
      for ix in range(ixmin+2, ixmax-1):
          for iz in range(izmin, izmax):
              dpu = c[0]*u[ix-1, iz] + c[1]*u[ix, iz] + c[2]*u[ix+1, iz]
              dmu = c[1]*u[ix-1, iz] + c[2]*u[ix, iz] + c[3]*u[ix+1, iz]
              K[ix, iz] = 0.5*dpu - 0.5*dmu
  </code></pre>
  
  <p>If I run these lines : </p>
  
  <pre><code>u = np.random.rand(256, 256)
  K1 = np.zeros_like(u)
  K2 = np.zeros_like(u)
  
  cp1(u, K1, 100, 150, 100, 150)
  cp2(u, K2, 100, 150, 100, 150)
  </code></pre>
  
  <p>the instruction <code>np.all(K1 == K2)</code> returns <code>False</code>. The difference between the two arrays is close to the machine precision (about 5e-17) but using this function thousand of times is enough to give me large differences on the final results.</p>
  
  <p>Now, if remove the <code>nogil</code> instructions in <code>cpc1</code> and <code>cpc2</code> and I replace <code>cdef double[4] c</code> by <code>c = np.zeros(4)</code>, both <code>cp1</code> and <code>cp2</code> functions return the same results. The problem is that I lose about 50% performance using <code>ndarray</code> instead of <code>c array</code>.</p>
  
  <p>In think the problem comes from the c array precision, but why the value of <code>boundscheck</code> and <code>nonecheck</code> have an impact on the results in this case (no access out of bounds, no none, ...)</p>
  
  <p>Is there a way to solve this ?</p>
  
  <p><strong>EDIT</strong></p>
  
  <p>As highlighted by ead, if I compile the code without <code>-03 --ffast-math -march=native</code>, both <code>cp1</code> and <code>cp2</code> return the same results ! But at the cost of doubling the execution time ! I can more or less understand why <code>O3</code> and <code>ffast-math</code> lead to unexpected results doing aggressive optimizations, but I don't undertand why <code>march=native</code> also <em>breaks</em> the code.</p>
  
  <p>Is there a way to preserve both performance and precision ?</p>