NumPy Factorial: Why np.factorial Fails and What to Use

There is no NumPy factorial function. np.factorial() doesn’t exist and never has, which is why you’re here after an AttributeError.

import math
math.factorial(5)                       # 120, one number, exact

from scipy.special import factorial
factorial(np.array([3, 4, 5]))          # array([ 6., 24., 120.])

Use math.factorial for a single integer and scipy.special.factorial for a whole array. The rest of this page is the detail that decides between them.

Versions used below: NumPy 2.5.3, SciPy 1.18.1, Python 3.12.5.

Why does np.factorial raise an AttributeError?

NumPy has never shipped a factorial. The function lives in the standard library and in SciPy, not in the array package:

import numpy as np

print("numpy version:", np.__version__)
print("does np.factorial exist?", hasattr(np, "factorial"))

np.factorial(5)

What Python says:

numpy version: 2.5.3
does np.factorial exist? False
Traceback (most recent call last):
  File "C:\pyguides\numpy_factorial_error.py", line 6, in <module>
    np.factorial(5)
    ^^^^^^^^^^^^
  File "C:\pyguides\venv\Lib\site-packages\numpy\__init__.py", line 782, in __getattr__
    raise AttributeError(f"module {__name__!r} has no attribute {attr!r}")
AttributeError: module 'numpy' has no attribute 'factorial'

The second half of the confusion is np.math.factorial, which used to work by accident. np.math was a passthrough to the standard library module and it was removed in NumPy 2.0:

import numpy as np

# np.math was a passthrough to the standard library and was removed in NumPy 2.0
print("does np.math exist?", hasattr(np, "math"))

np.math.factorial(5)

On NumPy 2.x:

does np.math exist? False
Traceback (most recent call last):
  File "C:\pyguides\runs\factorial\ex_npmath.py", line 6, in <module>
    np.math.factorial(5)
    ^^^^^^^
  File "C:\pyguides\venv\Lib\site-packages\numpy\__init__.py", line 782, in __getattr__
    raise AttributeError(f"module {__name__!r} has no attribute {attr!r}")
AttributeError: module 'numpy' has no attribute 'math'. Did you mean: 'emath'?

If an old tutorial or a Stack Overflow answer tells you to use np.math.factorial, it was written before that removal. Import math yourself instead.

math.factorial for a single number

For one integer this is the right answer, and it stays exact no matter how large the result gets:

import math

print("math.factorial(5) :", math.factorial(5))
print("math.factorial(20):", math.factorial(20))

# exact, however large it gets: Python integers do not overflow
big = math.factorial(50)
print("math.factorial(50):", big)
print("digits            :", len(str(big)))

Output:

math.factorial(5) : 120
math.factorial(20): 2432902008176640000
math.factorial(50): 30414093201713378043612608166064768844377641568960512000000000000
digits            : 65

Python integers have no fixed width, so they don’t overflow. That’s the property you give up the moment you move to floats.

scipy.special.factorial for a NumPy array

When you have an array of values, SciPy applies the factorial element by element:

import numpy as np
from scipy.special import factorial

values = np.array([0, 1, 2, 3, 4, 5, 6])

print("scipy.special.factorial:", factorial(values))
print("dtype                  :", factorial(values).dtype)

# exact integers instead, if the numbers stay small
print("exact=True             :", factorial(values, exact=True))
print("dtype                  :", factorial(values, exact=True).dtype)

Output:

scipy.special.factorial: [  1.   1.   2.   6.  24. 120. 720.]
dtype                  : float64
exact=True             : [  1   1   2   6  24 120 720]
dtype                  : int32

Note the dtype. The default returns float64, which is fine for plotting or probability work and wrong the moment you need exact integers.

exact=True switches to Python integers, so the dtype becomes object once the values grow and nothing overflows. You pay for it in speed on large arrays.

Where scipy.special.factorial stops being exact

This is the part that catches people, and it has a definite boundary:

import math
import numpy as np
from scipy.special import factorial

for n in (21, 22, 23, 30):
    exact = math.factorial(n)
    approx = factorial(n)              # returns a float64
    print(f"{n}!  exact {exact}")
    print(f"     float {approx:.0f}   match: {int(approx) == exact}")

Output:

21!  exact 51090942171709440000
     float 51090942171709440000   match: True
22!  exact 1124000727777607680000
     float 1124000727777607680000   match: True
23!  exact 25852016738884976640000
     float 25852016738884978212864   match: False
30!  exact 265252859812191058636308480000000
     float 265252859812191104246398737973248   match: False

A float64 keeps 53 bits of mantissa. Factorials pick up a lot of factors of two along the way, so they stay exact further than digit-counting suggests — all the way to 22!.

23! is the first that needs more precision than a double has. It still prints in full, which is what makes this easy to miss.

If the result feeds into a count, an index or a comparison, use exact=True or math.factorial. The same reasoning applies to any NumPy dtype choice.

Factorials with np.prod and np.cumprod

A factorial is just a product, so NumPy can do it without any special function:

import numpy as np

n = 10

# a factorial is just the product of 1..n
print("np.prod   :", np.prod(np.arange(1, n + 1)))

# cumprod gives every factorial up to n in one pass
print("np.cumprod:", np.cumprod(np.arange(1, n + 1)))

# watch the dtype: the default integer overflows silently past 20!
print("21! via np.prod:", np.prod(np.arange(1, 22, dtype=np.int64)))
print("21! exact      :", 51090942171709440000)

Output:

np.prod   : 3628800
np.cumprod: [      1       2       6      24     120     720    5040   40320  362880
 3628800]
21! via np.prod: -4249290049419214848
21! exact      : 51090942171709440000

np.cumprod is genuinely useful when you want 1! through n! in one pass, for a series expansion for instance.

The overflow is the catch. A NumPy integer array has a fixed width, so 21! silently wraps rather than raising, which is exactly the failure Python’s own integers avoid.

math.perm and math.comb instead of dividing factorials

Plenty of factorial searches are really permutation or combination problems. Python 3.8 added both, and they avoid building the huge intermediate numbers:

import math

# what people often actually want
print("math.perm(10, 3):", math.perm(10, 3), "  ordered arrangements")
print("math.comb(10, 3):", math.comb(10, 3), "   unordered combinations")

# the long way round, which overflows and is slower
print("by hand         :", math.factorial(10) // math.factorial(7))
print("same answer     :", math.perm(10, 3) == math.factorial(10) // math.factorial(7))

Output:

math.perm(10, 3): 720   ordered arrangements
math.comb(10, 3): 120    unordered combinations
by hand         : 720
same answer     : True

math.perm(n, k) counts ordered arrangements and math.comb(n, k) counts unordered ones. Both are faster and safer than dividing one factorial by another.

More NumPy and Python maths guides:

• NumPy data types
• np.add.at for repeated indices
• Create a 2D array with NumPy
• Round numbers to 2 decimal places
• SciPy root finding

Frequently asked questions