MAINT: Lattice sum eta estimation

src/treams/lattice/_esum.pyx

@@ -51,10 +51,8 @@ cdef number_t _check_eta(number_t eta, number_t k, double *a, long ds, long dl)
     else:
         return NAN
     cdef double absres = abs(res)
-    if absres < 0.1:
-        res *= 0.1 / absres
-    elif absres > 10:
-        res *= 10 / absres
+    if absres < 0.125:
+        res *= 0.125 / absres
     return res * sqrtd(2 * pi)
 
 
@@ -83,11 +81,11 @@ cdef double complex _redincgamma(double n, double complex z) nogil:
     cdef long twicen = <long>(2 * n)
     if twicen != 2 * n:
         raise ValueError("l is not integer of half-integer")
-    cdef double singularity = 0.5e-3  # Value of the singularity: smaller is stronger
+    cdef double singularity = 1e-7  # Value of the singularity: smaller is stronger
     cdef double complex res
     if creal(z) * creal(z) + cimag(z) * cimag(z) < 1e-12:
         if twicen == 0:
-            return -<double>EULER - 2 * log(singularity) + 0.5j * pi
+            return -<double>EULER - log(singularity) + 0.5j * pi
         if twicen > 0 and (twicen + 2) % 4 != 0:
             res = tgamma(n) * (1 + exp(-1j * pi * n)) / (2 * singularity)
         else:
@@ -533,7 +531,7 @@ cdef double complex realsumcw2d(long l, number_t k, double *kpar, double *a, dou
     cdef double mone = -1, zero = 0
     cdef int dim = 2, inc = 1
     cdef char c = b'N'
-    for i in range(start, 50):
+    for i in range(start, 20):
         pprev = prev
         prev = realsum
         point[0] = -i
@@ -571,7 +569,7 @@ cdef double complex recsumcw2d(long l, number_t k, double *kpar, double *a, doub
     cdef char c = b'N'
     cdef double b[4]
     _misc.recvec2(a, a + 2, b , b + 2)
-    for i in range(50):
+    for i in range(20):
         pprev = prev
         prev = recsum
         point[0] = -i
@@ -619,7 +617,7 @@ cdef double complex realsumsw2d(long l, long m, number_t k, double *kpar, double
     cdef double mone = -1, zero = 0
     cdef int dim = 2, inc = 1
     cdef char c = b'N'
-    for i in range(start, 50):
+    for i in range(start, 20):
         pprev = prev
         prev = realsum
         point[0] = -i
@@ -672,7 +670,7 @@ cdef double complex recsumsw2d(long l, long m, number_t k, double *kpar, double
     cdef char c = b'N'
     cdef double b[4]
     _misc.recvec2(a, a + 2, b , b + 2)
-    for i in range(50):
+    for i in range(20):
         pprev = prev
         prev = recsum
         point[0] = -i
@@ -841,7 +839,7 @@ cdef double complex realsumsw3d(long l, long m, number_t k, double *kpar, double
     cdef double mone = -1, zero = 0
     cdef int dim = 3, inc = 1
     cdef char c = b'N'
-    for i in range(start, 50):
+    for i in range(start, 10):
         pprev = prev
         prev = realsum
         point[0] = -i
@@ -884,7 +882,7 @@ cdef double complex recsumsw3d(long l, long m, number_t k, double *kpar, double
     cdef double one = 1
     cdef int dim = 3, inc = 1
     cdef char c = b'N'
-    for i in range(50):
+    for i in range(10):
         pprev = prev
         prev = recsum
         point[0] = -i
@@ -1051,7 +1049,7 @@ cdef double complex realsumsw2d_shift(long l, long m, number_t k, double *kpar,
     cdef long point[3]
     cdef int dim = 2, inc = 1
     cdef char c = b'N'
-    for i in range(50):
+    for i in range(20):
         pprev = prev
         prev = realsum
         point[0] = -i
@@ -1116,7 +1114,7 @@ cdef double complex recsumsw2d_shift(long l, long m, number_t k, double *kpar, d
     cdef char c = b'N'
     cdef double b[4]
     _misc.recvec2(a, a + 2, b, b + 2)
-    for i in range(50):
+    for i in range(20):
         pprev = prev
         prev = recsum
         point[0] = -i
@@ -1291,7 +1289,7 @@ cdef double complex recsumsw1d_shift(long l, long m, number_t k, double kpar, do
     """
     if l < 0:
         return NAN
-    eta = _check_eta(eta, k, &a, 3, 2)
+    eta = _check_eta(eta, k, &a, 3, 1)
     if fabs(r[0]) < 1e-100 and fabs(r[1]) < 1e-100:
         if m == 0:
             return recsumsw1d(l, k, kpar, a, r[2], eta)

tests/integration/test_compare_direct_sum.py

@@ -40,7 +40,7 @@ def test_singular(self, l, r):  # noqa: E741
         kpar = 0
         assert isclose(
             la.lsumsw1d(l, k, kpar, a, r, 0),
-            np.sum(la.dsumsw1d(l, k, kpar, a, r, np.arange(1_100_000))),
+            np.sum(la.dsumsw1d(l, k, kpar, a, r, np.arange(1_300_000))),
             rel_tol=0.05,
             abs_tol=EPSSQ,
         )
@@ -68,7 +68,7 @@ def test_singular(self, l, r):  # noqa: E741
         kpar = 0
         assert isclose(
             la.lsumcw1d(l, k, kpar, a, r, 0),
-            np.sum(la.dsumcw1d(l, k, kpar, a, r, np.arange(1_600_000))),
+            np.sum(la.dsumcw1d(l, k, kpar, a, r, np.arange(2_000_000))),
             rel_tol=0.05,
             abs_tol=EPSSQ,
         )
@@ -453,14 +453,14 @@ def test_singular(self, lm, r):
         if lm == (2, 0) and np.array_equal(r, [0.2, 0.2, 1.1]):
             assert isclose(
                 la.lsumsw1d_shift(*lm, k, kpar, a, r, 0),
-                la.dsumsw1d_shift(*lm, k, kpar, a, r, np.arange(2_200_000)).sum(),
+                la.dsumsw1d_shift(*lm, k, kpar, a, r, np.arange(700_000)).sum(),
                 rel_tol=1e-1,
                 abs_tol=EPS,
             )
         else:
             assert isclose(
                 la.lsumsw1d_shift(*lm, k, kpar, a, r, 0),
-                la.dsumsw1d_shift(*lm, k, kpar, a, r, np.arange(800_000)).sum(),
+                la.dsumsw1d_shift(*lm, k, kpar, a, r, np.arange(1_100_000)).sum(),
                 rel_tol=1e-1,
                 abs_tol=EPS,
             )

tests/unit/test_lattice_module.py

@@ -260,7 +260,7 @@ def test_singular(self):
         a = np.eye(2)
         assert isclose(
             la.lsumsw2d_shift(6, 0, 2 * np.pi, [0, 0], a, [0.2, 0.2, 0.1], 0),
-            110.69823590056764 - 138.74329901305353j,
+            554145.857367349-554173.9024304616j,
         )
 
     def test_z_zero(self):
@@ -307,7 +307,7 @@ def test2(self):
     def test_singular(self):
         assert isclose(
             la.lsumsw1d_shift(6, 0, 2 * np.pi, 0, 1, [0.2, 0.2, 0.1], 0),
-            -0.3552853450806117 - 25.183882160999745j,
+            -0.3552862673579646-24.943884092634878j,
         )
 
     def test_xy_zero(self):
@@ -364,7 +364,7 @@ def test2(self):
     def test_singular(self):
         assert isclose(
             la.lsumcw1d_shift(6, 2 * np.pi, 0, 1, [0.2, 0.1], 0),
-            -224.72932005925895 + 675.8601828886733j,
+            -1743317.7346371387+1743768.8654999665j,
             1e-5,
         )