feat: add dqagi

NumbaQuadpack/driver.py

@@ -9,11 +9,12 @@
 
 rootdir = os.path.dirname(os.path.realpath(__file__))+'/'
 libquadpack = ct.CDLL(rootdir+'libcquadpack.so')
+
+# # # # dqags
 dqags_ = libquadpack.dqags
 dqags_.argtypes = [ct.c_void_p, ct.c_double, ct.c_double, ct.c_double, ct.c_double, \
                  ct.c_void_p, ct.c_void_p, ct.c_void_p, ct.c_void_p,]
 dqags_.restype = ct.c_double
-
 @njit
 def dqags(funcptr, a, b, data = np.array([0.0], np.float64), epsabs = 1.49e-08, epsrel = 1.49e-08):
     abserr = np.array(0.0,np.float64)
@@ -29,3 +30,34 @@ def dqags(funcptr, a, b, data = np.array([0.0], np.float64), epsabs = 1.49e-08,
         success = False
 
     return sol, abserr.item(), success
+
+# # # # dqagi
+dqagi_ = libquadpack.dqagi
+dqagi_.argtypes = [ct.c_void_p, ct.c_double, ct.c_void_p, ct.c_double, ct.c_double, \
+                 ct.c_void_p, ct.c_void_p, ct.c_void_p, ct.c_void_p,]
+dqagi_.restype = ct.c_double
+@njit
+def dqagi(funcptr, a, b, data = np.array([0.0], np.float64), epsabs = 1.49e-08, epsrel = 1.49e-08):
+    '''
+    a - optional finite bound on integral.
+
+    b - specifies range of integration as follows:
+        b = -1 -- range is from -infinity to bound,
+        b =  1 -- range is from bound to +infinity,
+        b =  2 -- range is from -infinity to +infinity,
+                (bound is immaterial in this case).
+    '''
+    abserr = np.array(0.0,np.float64)
+    neval = np.array(0,np.int32)
+    ier = np.array(0,np.int32)
+
+    sol = dqagi_(funcptr, a, b, epsabs, epsrel, \
+                 abserr.ctypes.data, neval.ctypes.data, \
+                 ier.ctypes.data, data.ctypes.data)
+
+    success = True
+    if ier != 0:
+        success = False
+
+    return sol, abserr.item(), success
+

README.md

@@ -59,3 +59,26 @@ def test():
     return sol
 sol = test() # this works!
 ```
+
+Integration over (semi-) infinite intervals.
+```python
+from NumbaQuadpack import dqagi
+
+@nb.njit(nogil=True)
+def exp_test(x):
+    return np.exp(-x**2)
+
+@nb.cfunc(quadpack_sig)
+def f(x, data_):
+    c, d, e = nb.carray(data_, (3,))
+    return exp_test(x) + c + d * e
+funcptr = f.address
+
+@nb.njit()
+def test():
+    data = np.array([1.0, 1.0, 2.0],np.float64)
+    sol, abserr, success = dqagi(funcptr, 0, 2, data = data)
+    return sol, abserr, success
+sol, abserr, success = test()
+sol, abserr, success
+```