Repository: sympy/sympy
Pull request: Return sympy integers instead of python integers
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234". See
https://github.com/blog/1506-closing-issues-via-pull-requests . Please also
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #18338


#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
*   functions
      *   Sympified output numbers for `factorial` and `binomial` used with `Mod`.
<!-- END RELEASE NOTES -->


Problem statement:
Clean up gotchas in _eval_Mod
Things like
```python
Mod(factorial(4, evaluate=False), 3)
Mod(factorial(5, evaluate=False), 6)
```
shouldn't return python integers, but sympy integers.

I see `_eval_Mod` in `factorial` and `binomial` have such pitfalls.
When you see `return 0`, it should be changed to `return S.Zero`.

Hints from issue discussion:
Shouldn't it be automatically sympified? 
I don’t think that it is.
I see a lot of _eval does not sympify, and I thought that this was left on purpose that post-sympify can slow down things slightly for stuff which can be resolved in a more careful review.
I'm working on this issue.

Commit message: Add more tests to improve coverage

Patch to review:
diff --git a/sympy/functions/combinatorial/factorials.py b/sympy/functions/combinatorial/factorials.py
index 69378c4a0c26..ed77c4343a13 100644
--- a/sympy/functions/combinatorial/factorials.py
+++ b/sympy/functions/combinatorial/factorials.py
@@ -209,7 +209,7 @@ def _eval_Mod(self, q):
             aq = abs(q)
             d = aq - n
             if d.is_nonpositive:
-                return 0
+                return S.Zero
             else:
                 isprime = aq.is_prime
                 if d == 1:
@@ -218,9 +218,9 @@ def _eval_Mod(self, q):
                     # its inverse (if n > 4 is a composite number, then
                     # (n-1)! = 0 mod n)
                     if isprime:
-                        return -1 % q
+                        return S(-1 % q)
                     elif isprime is False and (aq - 6).is_nonnegative:
-                        return 0
+                        return S.Zero
                 elif n.is_Integer and q.is_Integer:
                     n, d, aq = map(int, (n, d, aq))
                     if isprime and (d - 1 < n):
@@ -231,7 +231,7 @@ def _eval_Mod(self, q):
                     else:
                         fc = self._facmod(n, aq)
 
-                    return Integer(fc % q)
+                    return S(fc % q)
 
     def _eval_rewrite_as_gamma(self, n, **kwargs):
         from sympy import gamma
@@ -931,14 +931,14 @@ def _eval_Mod(self, q):
 
             # handle negative integers k or n
             if k < 0:
-                return 0
+                return S.Zero
             if n < 0:
                 n = -n + k - 1
                 res = -1 if k%2 else 1
 
             # non negative integers k and n
             if k > n:
-                return 0
+                return S.Zero
 
             isprime = aq.is_prime
             aq = int(aq)
@@ -997,7 +997,7 @@ def _eval_Mod(self, q):
                             res *= pow(prime, exp, aq)
                             res %= aq
 
-            return Integer(res % q)
+            return S(res % q)
 
     def _eval_expand_func(self, **hints):
         """
diff --git a/sympy/functions/combinatorial/tests/test_comb_factorials.py b/sympy/functions/combinatorial/tests/test_comb_factorials.py
index 9e45a69a7bef..d5473087e335 100644
--- a/sympy/functions/combinatorial/tests/test_comb_factorials.py
+++ b/sympy/functions/combinatorial/tests/test_comb_factorials.py
@@ -232,6 +232,8 @@ def test_factorial_Mod():
     assert Mod(factorial(q - 1800, evaluate=False), q) == 905504050
     assert Mod(factorial(153, evaluate=False), r) == Mod(factorial(153), r)
     assert Mod(factorial(255, evaluate=False), s) == Mod(factorial(255), s)
+    assert Mod(factorial(4, evaluate=False), 3) == S.Zero
+    assert Mod(factorial(5, evaluate=False), 6) == S.Zero
 
 
 def test_factorial_diff():
@@ -497,6 +499,9 @@ def test_binomial_Mod_slow():
     assert Mod(binomial(9734, 451, evaluate=False), q) == Mod(binomial(9734, 451), q)
     assert Mod(binomial(-10733, 4459, evaluate=False), q) == Mod(binomial(-10733, 4459), q)
     assert Mod(binomial(-15733, 4458, evaluate=False), q) == Mod(binomial(-15733, 4458), q)
+    assert Mod(binomial(23, 38), q) == S.Zero
+    assert Mod(binomial(4556, 38635874), q) == S.Zero
+    assert Mod(binomial(456, 789), q) == S.Zero
 
     # binomial factorize
     assert Mod(binomial(753, 119, evaluate=False), r) == Mod(binomial(753, 119), r)


Code review comments:
diff_hunk: @@ -497,6 +499,9 @@ def test_binomial_Mod_slow():
     assert Mod(binomial(9734, 451, evaluate=False), q) == Mod(binomial(9734, 451), q)
     assert Mod(binomial(-10733, 4459, evaluate=False), q) == Mod(binomial(-10733, 4459), q)
     assert Mod(binomial(-15733, 4458, evaluate=False), q) == Mod(binomial(-15733, 4458), q)
+    assert Mod(binomial(23, 38), q) == S.Zero
+    assert Mod(binomial(4556, 38635874), q) == S.Zero
+    assert Mod(binomial(456, 789), q) == S.Zero
path: sympy/functions/combinatorial/tests/test_comb_factorials.py
text: @user1:
These still doesn't increase the coverage for the diff for now.
What about changing these to
```python3
    assert Mod(binomial(23, -38, evaluate=False), q) is S.Zero
    assert Mod(binomial(23, 38, evaluate=False), q) is S.Zero
```