Cython quadrature problem : big array element access

Question

My problem is that I want to to a 3d quadrature on some function. The objective is to compute electronic-integrals. The electron-electron potential integrals being the problem now. I am using 20 Gauss-Chebyshev integration points in each direction of the second kind. Meaning 20³=8000 points. But when I use this function in my test, my result is nan :

# My test file that is ran by pytest
from qmodelling.basis.primitive_gaussian import PrimitiveGaussian, gaussian_integral, integrate_3d

# No need to precise the Test Class

def test_integrate_3d_unit_cube_volume(self):
        I = integrate_3d(lambda x : 1.0)
        assert isclose(I, 8) # OK since the function returns 1 anyway

    def test_integrate_3d_test_1(self):
        I = integrate_3d(lambda x : x[0]) # NAN
        assert isclose(I, 0.0)

I feel like the way I store my quadrature points in array views (double[:] or double[:]) is not optimal. Or maybe the problem is memory leaks. If I can have any piece of advise from any Cython user for integral quadratures I would appreciate it.

The code is below :

In the file primitive_gaussian.pyx

from qmodelling.integral.quadrature.primitive_gaussian_quadrature cimport get_Chebyshev_01, get_Chebyshev_3d
data_3d = get_Chebyshev_3d()
cdef Chebyshev_weights_3d = data_3d[0] # double[:]
cdef Chebyshev_abscissa_3d = data_3d[1]# double[:, :]

#######################
# quadrature test functions 
#######################
cpdef double integrate_3d(fun):

    I = 0.0
    for k in range(ngauss_3d):
        rk = Chebyshev_abscissa_3d[k,:]
        wk = Chebyshev_weights_3d[k]

        # apply arctanh to each elt of double[:]
        Rk = array_atanh(rk)

        # prod is like np.prod but for double[:] and power double apply
        # exponents element-wise
        # not really the interesting part since it a variable substitution to
        #integrate of R3 with points on [-1,1]³
        # the line of code below is bug-free
        I += wk * fun(Rk) / (-prod(add_double(power_double(rk,2.0), - 1.0)))

    return I

In file primitive_gaussian_quadrature.pyx

import numpy as np
from cython cimport tuple
from qmodelling.constants cimport dtype
from qmodelling.integral.quadrature.get_quadrature import get_quadrature_points

ngauss = 20
ngauss_3d = ngauss**3

cdef cnp.ndarray[dtype, ndim=2] get_Chebyshev_01():
    # pts on [-1,1]
    Chebyshev_quadrature_points = get_quadrature_points(
        n=ngauss, quadrature_type="gauss_chebyshev_2"
    )
    
    # pts on [0,1]
    Chebyshev_01 = Chebyshev_quadrature_points
    Chebyshev_01[:, 1] = (
        Chebyshev_01[:, 1] + 1.0              
    ) / 2.0
    # sum weights = 1
    Chebyshev_01[:, 0] /= 2.0 

    return Chebyshev_01                                          

cdef tuple get_Chebyshev_3d():
    # pts on [-1,1] 
    Chebyshev_quadrature_points = get_quadrature_points(
        n=ngauss, quadrature_type="gauss_chebyshev_2" 
    )

    # Generating points for R³ integration with variable substitution on [-1,1]³
    weights_1d = Chebyshev_quadrature_points[:, 0]                        
    weights_3d = np.meshgrid(weights_1d, weights_1d, weights_1d, indexing="ij")
    cdef double[:] Chebyshev_weights_3d = (weights_3d[0] * weights_3d[1] * weights_3d[2]).flatten("A")

    abscissa_1d = Chebyshev_quadrature_points[:, 1]
    abscissa_3d = np.meshgrid(abscissa_1d, abscissa_1d, abscissa_1d, indexing="ij")     
    cdef double[:,:] Chebyshev_abscissa_3d = np.stack([p.flatten("A") for p in abscissa_3d], axis=1).reshape(
       -1, 3
    )
        
    return Chebyshev_weights_3d, Chebyshev_abscissa_3d