divisors should return Integers

[asmeurer]

See https://www.johndcook.com/blog/2019/03/14/irreducible-polynomials/. He has to use sum(...)//n instead of sum(...)/n because the computation produces floats. But if divisors returned Integer types, this wouldn't happen.

In general, I think all number theory functions should return SymPy types. For those that return numbers, this is a consequence of #16244, but even divisors, which doesn't make sense to be symbolic, should as well in my opinion. That way the results of them combine the most naturally in larger expressions.

Apparently this sort of thing has been discussed before. See #6981. So we should maybe reopen the discussion here before changing things. factorint is another example.

[smichr]

There is a large performance hit when using Integer instead of int. There is also a hit for using a function, so when ntheory functions call one another there will be redundant calls to as_int since these functions are all user-facing.

One way around this would be to make all existing user-facing functions that call another ntheory function shells that basically just do type checking and then call a _foo form of the function (which calls a _foo2 function).

But is it unreasonable to expect that a user has to know what data type they were using? Before truedivision was the norm we had to worry about ints truncating to ints...now we have to worry above ints casting to floats.

It would be interesting to see how much things slow down if all output becomes Integer. That shouldn't be to hard to find all returns and make sure the the return value is wrapped in Integer.

[asmeurer]

This would already be an implicit consequence of #16244 for functions like prime that return integers. That issue is a good idea anyway because we can do some symbolic computations on the ntheory functions.

Splitting out _foo methods that do the calculation without type checks and return int seems fine to me. For the functions that return integers, which are turned into classes by #16244, we can make it a staticmethod on the class. But I'm also not so sure it would affect performance that much, unless there is somewhere where one of these functions is being called many times.

But is it unreasonable to expect that a user has to know what data type they were using? Before truedivision was the norm we had to worry about ints truncating to ints...now we have to worry above ints casting to floats.

Once you've used SymPy enough you start to get a feel for when you don't need to worry about / because one of the arguments comes from a SymPy function, and so it "absorbs" its sympyness. SymPy functions that don't return integers break this. Since int/int is the only place where things go wrong if you don't map to SymPy types first, people usually don't think to bother converting inputs to SymPy types first. If we can minimize the places where this happens by never returning int, that is even better.

[smichr]

So I just tried a decorator and applied it to all publicly facing funtions in factor_ and my fan is runing like crazy as I type....

def iout(func):
    from sympy.utilities.iterables import is_sequence
    from sympy.core.numbers import Integer

    def do(*args, **kwargs):
        out = func(*args, **kwargs)

        def ido(s):
            if is_sequence(s):
                return type(s)([ido(i) for i in s])
            elif type(s) is dict:
                return dict([(ido(k), ido(v)) for (k, v) in s.items()])
            elif type(s) is int:
                return Integer(s)
            else:
                return s  #  strings, Rational, Expr, etc...

        return ido(out)

    return do

The least destructive and time consuming method would be to make a function like that part of sympy.functions and if a user is concerned about the int issue they can just wrap the output in iout.

e.g.

>>> iout(factorint(2019))
{3: 1, 673: 1}
>>> _[3].is_Integer
True

[jksuom]

number theory functions should return SymPy types

I don't think that changes like return int(n) -> return Integer(n) would affect performance. I would be +1 for those. But then I am not so sure about what to do with functions returning lists or dicts.

[smichr]

would affect performance

decorating all factor_ functions with iout did not allow test_comb_factorial to finish in a reasonable amount of time.

diff --git a/sympy/ntheory/factor_.py b/sympy/ntheory/factor_.py
index 266fc30..7ea122a 100644
--- a/sympy/ntheory/factor_.py
+++ b/sympy/ntheory/factor_.py
@@ -20,6 +20,27 @@
 from .generate import sieve, primerange, nextprime
 
 
+def iout(func):
+    from sympy.utilities.iterables import is_sequence
+    from sympy.core.numbers import Integer
+
+    def do(*args, **kwargs):
+        out = func(*args, **kwargs)
+
+        def ido(s):
+            if is_sequence(s):
+                return type(s)([ido(i) for i in s])
+            elif type(s) is dict:
+                return dict([(ido(k), ido(v)) for (k, v) in s.items()])
+            elif type(s) is int:
+                return Integer(s)
+            else:
+                return s  #  strings, Rational, Expr, etc...
+
+        return ido(out)
+
+    return do
+
 # Note: This list should be updated whenever new Mersenne primes are found.
 # Refer: https://www.mersenne.org/
 MERSENNE_PRIME_EXPONENTS = (2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203,
@@ -32,6 +53,7 @@
     small_trailing[1<<j::1<<(j+1)] = [j] * (1<<(7-j))
 
 
+@iout
 def smoothness(n):
     """
     Return the B-smooth and B-power smooth values of n.
@@ -62,6 +84,7 @@ def smoothness(n):
     return max(facs), max(m**facs[m] for m in facs)
 
 
+@iout
 def smoothness_p(n, m=-1, power=0, visual=None):
     """
     Return a list of [m, (p, (M, sm(p + m), psm(p + m)))...]
@@ -165,6 +188,7 @@ def smoothness_p(n, m=-1, power=0, visual=None):
     return '\n'.join(lines)
 
 
+@iout
 def trailing(n):
     """Count the number of trailing zero digits in the binary
     representation of n, i.e. determine the largest power of 2
@@ -214,6 +238,7 @@ def trailing(n):
     return t
 
 
+@iout
 def multiplicity(p, n):
     """
     Find the greatest integer m such that p**m divides n.
@@ -282,6 +307,7 @@ def multiplicity(p, n):
     return m
 
 
+@iout
 def perfect_power(n, candidates=None, big=True, factor=True):
     """
     Return ``(b, e)`` such that ``n`` == ``b**e`` if ``n`` is a
@@ -386,6 +412,7 @@ def perfect_power(n, candidates=None, big=True, factor=True):
         return False
 
 
+@iout
 def pollard_rho(n, s=2, a=1, retries=5, seed=1234, max_steps=None, F=None):
     r"""
     Use Pollard's rho method to try to extract a nontrivial factor
@@ -501,6 +528,7 @@ def pollard_rho(n, s=2, a=1, retries=5, seed=1234, max_steps=None, F=None):
     return None
 
 
+@iout
 def pollard_pm1(n, B=10, a=2, retries=0, seed=1234):
     """
     Use Pollard's p-1 method to try to extract a nontrivial factor
@@ -823,6 +851,7 @@ def done(n, d):
     return done(n, d)
 
 
+@iout
 def factorint(n, limit=None, use_trial=True, use_rho=True, use_pm1=True,
               verbose=False, visual=None, multiple=False):
     r"""
@@ -1234,6 +1263,7 @@ def factorint(n, limit=None, use_trial=True, use_rho=True, use_pm1=True,
         low, high = high, high*2
 
 
+@iout
 def factorrat(rat, limit=None, use_trial=True, use_rho=True, use_pm1=True,
               verbose=False, visual=None, multiple=False):
     r"""
@@ -1298,6 +1328,7 @@ def factorrat(rat, limit=None, use_trial=True, use_rho=True, use_pm1=True,
 
 
 
+@iout
 def primefactors(n, limit=None, verbose=False):
     """Return a sorted list of n's prime factors, ignoring multiplicity
     and any composite factor that remains if the limit was set too low
@@ -1359,6 +1390,7 @@ def rec_gen(n=0):
         yield p
 
 
+@iout
 def divisors(n, generator=False):
     r"""
     Return all divisors of n sorted from 1..n by default.
@@ -1405,6 +1437,7 @@ def divisors(n, generator=False):
     return rv
 
 
+@iout
 def divisor_count(n, modulus=1):
     """
     Return the number of divisors of ``n``. If ``modulus`` is not 1 then only
@@ -1449,6 +1482,7 @@ def _udivisors(n):
         yield d
 
 
+@iout
 def udivisors(n, generator=False):
     r"""
     Return all unitary divisors of n sorted from 1..n by default.
@@ -1496,6 +1530,7 @@ def udivisors(n, generator=False):
     return rv
 
 
+@iout
 def udivisor_count(n):
     """
     Return the number of unitary divisors of ``n``.
@@ -1544,6 +1579,7 @@ def _antidivisors(n):
             yield d
 
 
+@iout
 def antidivisors(n, generator=False):
     r"""
     Return all antidivisors of n sorted from 1..n by default.
@@ -1582,6 +1618,7 @@ def antidivisors(n, generator=False):
     return rv
 
 
+@iout
 def antidivisor_count(n):
     """
     Return the number of antidivisors [1]_ of ``n``.
@@ -1670,6 +1707,7 @@ def _eval_is_integer(self):
         return fuzzy_and([self.args[0].is_integer, self.args[0].is_positive])
 
 
+@iout
 class reduced_totient(Function):
     r"""
     Calculate the Carmichael reduced totient function lambda(n)
@@ -1719,6 +1757,7 @@ def _eval_is_integer(self):
         return fuzzy_and([self.args[0].is_integer, self.args[0].is_positive])
 
 
+@iout
 class divisor_sigma(Function):
     r"""
     Calculate the divisor function `\sigma_k(n)` for positive integer n
@@ -1789,6 +1828,7 @@ def eval(cls, n, k=1):
                            else e + 1 for p, e in factorint(n).items()])
 
 
+@iout
 def core(n, t=2):
     r"""
     Calculate core(n, t) = `core_t(n)` of a positive integer n
@@ -1856,6 +1896,7 @@ def core(n, t=2):
         return y
 
 
+@iout
 def digits(n, b=10):
     """
     Return a list of the digits of n in base b. The first element in the list
@@ -2047,6 +2088,7 @@ def eval(cls, n):
                 return sum(factorint(n).values())
 
 
+@iout
 def mersenne_prime_exponent(nth):
     """Returns the exponent ``i`` for the nth Mersenne prime (which
     has the form `2^i - 1`).
@@ -2128,6 +2170,7 @@ def is_mersenne_prime(n):
     return b and r in MERSENNE_PRIME_EXPONENTS
 
 
+@iout
 def abundance(n):
     """Returns the difference between the sum of the positive
     proper divisors of a number and the number.