Repository: sympy/sympy
Pull request: ntheory: divisor_sigma symbolic case recognition
Fixes: #16123

#### Brief description of what is fixed or changed
**Previously**, the only case which was recognized symbolically by divisor_sigma() function was: `divisor_sigma(a) -> a + 1.`

**Now**, it also recognizes:
`divisor_sigma(a**j*b**k)` where` a and b are prime` and` j and k are positive integers`
`divisor_sigma(a**j*b**k)-> (a**(j + 1) - 1)*(b**(k + 1) - 1)/((a - 1)*(b - 1))`


#### Other Comments
Regression Tests have been added


#### Release Notes


<!-- BEGIN RELEASE NOTES -->
* ntheory
  * Made divisor_sigma() symbolically recognize the case: `divisor_sigma(a**j*b**k)`
<!-- END RELEASE NOTES -->


Problem statement:
symbolic divisor_sigma
`divisor_sigma(a**k*b**j, 1)` where `a` and `b` are prime and `k` and `j` are positive integers is `(a**(k+1)-1)/(a-1)*(b**(j + 1) -1)/(b - 1)`. Right now, the only case that is recognized symbolically is `divisor_sigma(a) -> a + 1`.
```python
def eval(cls, n):
    ...
    if n.is_integer:
        if n.is_Mul or n.is_Pow:  # prime is already handled
            p = []
            for b, e in (_.as_base_exp() for _ in Mul.make_args(n)):
                if b.is_prime and e.is_positive:
                    p.append((b**(e + 1) - 1)/(b - 1))
                else:
                    return
            return Mul(*p)

>>> divisor_sigma(2**3*3**5)
5460
>>> a,b = var('a b',prime=True)
>>> j,k=var('j k', integer=True, positive=True)
>>> eq = a**j*b**k
>>> eq.subs({a:2,j:3,b:3,k:5})
5460
```

Hints from issue discussion:
I'm taking this.
@smichr can you please crosscheck the result of above substitution will be `1944` and for the current master I'm able to evaluate it.
sorry, I misunderstood the purpose of divisor_sigma(). I will do the required changes.

```
>>> divisors(1944)
[1, 2, 3, 4, 6, 8, 9, 12, 18, 24, 27, 36, 54, 72, 81, 108, 162, 216, 243, 324, 486, 648, 972, 1944]
>>> sum(divisors(1944))
5460
```

Commit message: Made the function recognise any value of k

Patch to review:
diff --git a/sympy/ntheory/factor_.py b/sympy/ntheory/factor_.py
index 504c1e1ae674..b1811013d6e1 100644
--- a/sympy/ntheory/factor_.py
+++ b/sympy/ntheory/factor_.py
@@ -2039,8 +2039,20 @@ class divisor_sigma(Function):
     def eval(cls, n, k=1):
         n = sympify(n)
         k = sympify(k)
+
         if n.is_prime:
             return 1 + n**k
+
+        if n.is_integer and not n.is_Integer:
+            p = []
+            for b, e in (_.as_base_exp() for _ in Mul.make_args(n)):
+                if b.is_prime and e.is_positive:
+                    p.append((b**(k*(e + 1)) - 1)/(b**k - 1) if k != 0
+                    else e + 1)
+                else:
+                    return
+            return Mul(*p)
+
         if n.is_Integer:
             if n <= 0:
                 raise ValueError("n must be a positive integer")
diff --git a/sympy/ntheory/tests/test_factor_.py b/sympy/ntheory/tests/test_factor_.py
index afea0d38eb04..ca9bd0e32544 100644
--- a/sympy/ntheory/tests/test_factor_.py
+++ b/sympy/ntheory/tests/test_factor_.py
@@ -402,6 +402,12 @@ def test_divisor_sigma():
     assert divisor_sigma(23450, 2) == 730747500
     assert divisor_sigma(23450, 3) == 14666785333344
 
+    a = Symbol("a", prime=True)
+    b = Symbol("b", prime=True)
+    j = Symbol("j", integer=True, positive=True)
+    k = Symbol("k", integer=True, positive=True)
+    assert divisor_sigma(a**j*b**k) == (a**(j + 1) - 1)*(b**(k + 1) - 1)/((a - 1)*(b - 1))
+
     m = Symbol("m", integer=True)
     k = Symbol("k", integer=True)
     assert divisor_sigma(m)


Code review comments:
diff_hunk: @@ -2039,8 +2039,20 @@ class divisor_sigma(Function):
     def eval(cls, n, k=1):
         n = sympify(n)
         k = sympify(k)
+
         if n.is_prime:
             return 1 + n**k
+
+        if n.is_integer and not n.is_Integer:
+            p = []
+            for b, e in (_.as_base_exp() for _ in Mul.make_args(n)):
+                if b.is_prime and e.is_positive:
path: sympy/ntheory/factor_.py
text: @user1:
If this fails then `None` will be returned? Is this expected? Any examples which can cause this to fail?

@author:
No that same expression will be returned and this case is possible:

```
In[1]: k, m, n = symbols('k m n', integer=True)
In[2]: divisor_sigma(k*m)
Out[2]: divisor_sigma(k*m)
```