That’s currently true, but while discussing of the PEP draft — this constraints was removed. I think, this is a good entry point. I can add arguments for removing this to the “Rationale” section, does it make sense for you?
Edit:
Strictly speaking, we don’t have even accurate descriptions of function signatures in docs, see Update function signatures to use `*` and `/` as needed · Issue #131885 · python/cpython · GitHub.
“Suitable functions must be computed exactly” is sufficient for me. As the link provided by @skirpichev explains, it is possible to exactly compute non-integer results that would still qualify. The simplest example I can think of is a function to compute the list of prime factors of an integer.
Note that initially math.factorial() accepted floats. math.gcd() was the first function which only accepted integer-like arguments. This was one of reasons to add a separate module – so its factorial() would be more type strict.
Finally, accepting floats in math.factorial() was deprecacated and disabled. But it took some time. With imath.factorial() you would have a type-strict function much earlier.
It did but in a strange way:
>>> math.factorial(5)
120
>>> math.factorial(5.0)
120
>>> math.factorial(5.1)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ValueError: factorial() only accepts integral values
That kind of coercion of integer-valued floats to int is not really in keeping with Python conventions in stdlib or in third party packages.
I think that the math module should have a factorial function that explicitly coerces to float and computes factorial for real integer and noninteger values like the scipy one:
>>> math.factorial(5)
120.0
>>> math.factorial(5.0)
120.0
>>> math.factorial(5.1)
142.45194406567867
This matches the other functions in the math module and is more useful in many situations than the current math.factorial function.
@skirpichev, wouldn’t adding a historical overview help justify the PEP?
First there was only factorial() which accepted floats, like all other math functions. Then gcd() which accepted only integers was moved from fractions, because fractions was not expected place for it, and math at that time looked the best option. Then addition of imath was discussed, to fix factorial() and to add a proper place for future integer functions. Then isqrt(), comb() and perm() was added, and addition of imath was proposed second time, so all new functions would go directly in imath, without littering the math namespace. Then something else happened.
I agree that it should, but only if we have a sane home (like imath) for the exact (bigint) factorial on integer arguments too.
This one would have the issue of getting weird results if a user does from math import int. Is there any current example of a module with a name that could shadow a builtin?
There is idlelib.format, as well as _pyrepl.input. I agree that math.int would be a lot more likely to lead to confusion, though.
“Implicit” is no good here, and “explicit” is a charitable reading. As others have noted, things like “number of ways” do not imply int.It can be read that way, but it’s not what it says. comb(4, 2) -> 6.0 would also give the “number of ways”, as would comb(4, 2) returning 6+0j or Fraction(6) or …
When int is the truth, that’s what the docs should say. “Integer” is a mathematical concept; int is a concrete Python class.
Similarly, things like " the largest positive integer" are also sloppy. That doesn’t explicitly say int either. As the docs say:
Note that all floating-point numbers of sufficiently large magnitude are exact integers
thus making it explicitly ambiguous what the docs mean when they say “integer”. The only unambiguous way to communicate “the Python int class” is to explicitly name it.
Accepting floats and coercing them to int is not like any other math functions. There is a huge difference between an approximate floating point calculation with bounded cost vs an exact calculation with unbounded cost regardless of what types are accepted and coerced as inputs. The nature of the internal calculation and the cost vs output accuracy guarantees are the reason for choosing to use say scipy.factorial(n) or mpmath.factorial(n) over math.factorial(n):
In [5]: %time scipy.special.factorial(10_000_000)
CPU times: user 117 μs, sys: 11 μs, total: 128 μs
Wall time: 133 μs
Out[5]: np.float64(inf)
In [6]: %time mpmath.factorial(10_000_000)
CPU times: user 1.28 ms, sys: 0 ns, total: 1.28 ms
Wall time: 1.3 ms
Out[6]: mpf('1.2024234005159035e+65657059')
In [7]: %time out = math.factorial(10_000_000)
CPU times: user 8min 48s, sys: 1.16 s, total: 8min 49s
Wall time: 8min 47s
In [8]: out.bit_length() / (8*1024**2)
Out[8]: 26.000503301620483
The output from scipy here is guaranteed to be a 64 bit float so fixed size but that means we get inf in 0.1 ms. The output from mpmath has unbounded exponent but is computed to fixed precision mantissa and gives us an accurate but approximate value in 1.3 ms. The output from math.factorial is a 26 megabyte integer that took almost 10 minutes to compute on this computer. If I increase n then math.factorial will happily consume all available memory and crash the entire computer but first it will lock up the CPU for longer than I am prepared to wait for this demonstration.
Before these integer functions were added the math module only had functions like the scipy one that you can think of as being O(1) computational cost without getting into trouble. Coercing a float into an int and then performing unbounded computation and producing outputs of unbounded size is a completely different kind of operation regardless of whether floats are allowed as input. There are good reasons for using something like scipy.factorial by default over something like math.factorial in most cases.
While mpmath handles the case above efficiently it can still get you into trouble if you go one level up in size because it does not bound the exponent:
In [12]: %time out = mpmath.factorial(mpmath.factorial(10_000))
CPU times: user 71.3 ms, sys: 3.99 ms, total: 75.3 ms
Wall time: 73.6 ms
In [13]: %time out = mpmath.factorial(mpmath.factorial(100_000))
CPU times: user 1.13 s, sys: 0 ns, total: 1.13 s
Wall time: 1.13 s
In [14]: %time out = mpmath.factorial(mpmath.factorial(1_000_000))
^C^C---------------------------------------------------------------------------
KeyboardInterrupt
There is a huge difference between an operation that is bounded in time and space and an operation that is not and there are not good reasons to mix up the kind of exact calculation that math.factorial does with approximate fixed precision floating point calculations.
An imath module can come with a big banner clearly stating that these operations are unbounded. Ideally it would also provide contexts that can be used to bound them making them safe to use for potentially large inputs:
from imath import context
ctx = context(maxbits=1000)
ctx.factorial(1000) # OverflowError
Please see update pr, where some historical details were mentioned: PEP 791: address PEP review comments by skirpichev · Pull Request #4430 · python/peps · GitHub
Just to note that while SciPy defaults to the quicker floating-point implementation, it’s possible to get exact results (for integer inputs) by passing exact=True; this extends to array-shaped inputs too (with some logic to avoid redundant calculations between elements).
Should one think about good interoperability with gmpy2, which in turn provides interfaces with de facto standard fast arbitrary precision integers library gmp, here?
The mpz type provides __index__() dunder method — so, module functions will work with the mpz per default.
All current functions, except perm(), have gmpy2’s counterparts with compatible interfaces. Probably we should keep that in mind for future: if the gmpy2 already has function, that proposed for the stdlib — we should chose same interface for the stdlib. With this, people e.g. in SymPy or mpmath can easily replace stdlib’s functions by gmpy2’s and use a different integer type.
Should I add something about this in the PEP? I doubt.
It won’t hurt to engage gmpy2’s authors/maintainers here. I’ve opened PEP 791 (imath stdlib proposal) · Issue #561 · aleaxit/gmpy · GitHub to inform them about this PEP.
IMHO you should mention gmpy2 in
otherwise module-scope-and-possible-extensions section reads as if Python interface to GMP is some sort of unexplored terrain.
For what it’s worth; NumPy’s integer scalar types and 0-d arrays also __index__().
Also, the Array API.
There are many types that define __index__, it is not really productive to enumerate them all here.
Indeed. The __index__() method is documented as a standard way to say “it’s an integer type”: “Called to implement operator.index(), and whenever Python needs to losslessly convert the numeric object to an integer object (such as in slicing, or in the built-in bin(), hex() and oct() functions). Presence of this method indicates that the numeric object is an integer type. Must return an integer.”