Repository: sympy/sympy

Issue:
Incorrect symbolic residue summation
This summation should clearly be a positive number:
```python
In [2]: s = Sum(1/(n**2+1), (n, 1, oo))

In [3]: s
Out[3]: 
  ∞         
 ____       
 ╲          
  ╲     1   
   ╲  ──────
   ╱   2    
  ╱   n  + 1
 ╱          
 ‾‾‾‾       
n = 1       

In [4]: s.doit() # <--- incorrect
Out[4]: -1/2
```
The correct result seems to be:
```python
In [5]: (pi*coth(pi)-1)/2
Out[5]: 
  1   π⋅coth(π)
- ─ + ─────────
  2       2    

In [6]: N(_)
Out[6]: 1.07667404746858

In [7]: N(s)
Out[7]: 1.07667404746858
```
The bug is in `eval_sum_residue` somewhere.

Discussion:

Patch:
diff --git a/sympy/concrete/summations.py b/sympy/concrete/summations.py
index 3c175f18b8..830785990c 100644
--- a/sympy/concrete/summations.py
+++ b/sympy/concrete/summations.py
@@ -1475,11 +1475,15 @@ def get_shift(denom):
         shift = - b / a / n
         return shift
 
+    #Need a dummy symbol with no assumptions set for get_residue_factor
+    z = Dummy('z')
+
     def get_residue_factor(numer, denom, alternating):
+        residue_factor = (numer.as_expr() / denom.as_expr()).subs(i, z)
         if not alternating:
-            residue_factor = (numer.as_expr() / denom.as_expr()) * cot(S.Pi * i)
+            residue_factor *= cot(S.Pi * z)
         else:
-            residue_factor = (numer.as_expr() / denom.as_expr()) * csc(S.Pi * i)
+            residue_factor *= csc(S.Pi * z)
         return residue_factor
 
     # We don't know how to deal with symbolic constants in summand
@@ -1520,7 +1524,7 @@ def get_residue_factor(numer, denom, alternating):
             return None
 
         residue_factor = get_residue_factor(numer, denom, alternating)
-        residues = [residue(residue_factor, i, root) for root in nonint_roots]
+        residues = [residue(residue_factor, z, root) for root in nonint_roots]
         return -S.Pi * sum(residues)
 
     if not (a.is_finite and b is S.Infinity):
@@ -1568,7 +1572,7 @@ def get_residue_factor(numer, denom, alternating):
             return None
 
     residue_factor = get_residue_factor(numer, denom, alternating)
-    residues = [residue(residue_factor, i, root) for root in int_roots + nonint_roots]
+    residues = [residue(residue_factor, z, root) for root in int_roots + nonint_roots]
     full_sum = -S.Pi * sum(residues)
 
     if not int_roots:


Test patch:
diff --git a/sympy/concrete/tests/test_sums_products.py b/sympy/concrete/tests/test_sums_products.py
index db7f4120e9..c0c6869957 100644
--- a/sympy/concrete/tests/test_sums_products.py
+++ b/sympy/concrete/tests/test_sums_products.py
@@ -319,6 +319,8 @@ def test_geometric_sums():
             (1, Eq(exp(-2*I*pi*(k - q)/n), 1)), (0, True)
     )
 
+    #Issue 23491
+    assert Sum(1/(n**2 + 1), (n, 1, oo)).doit() == S(-1)/2 + pi/(2*tanh(pi))
 
 def test_harmonic_sums():
     assert summation(1/k, (k, 0, n)) == Sum(1/k, (k, 0, n))