diff --git a/sumpy/expansion/level_to_order.py b/sumpy/expansion/level_to_order.py
index cc8440fd6989edd33e14169e4e8c066ed794e1ed..5d00c5b01642a5c821ce0d0237a9a7d1313227aa 100644
--- a/sumpy/expansion/level_to_order.py
+++ b/sumpy/expansion/level_to_order.py
@@ -25,6 +25,8 @@ __doc__ = """
 .. autoclass:: SimpleExpansionOrderFinder
 """
 
+import math
+
 import numpy as np
 
 
@@ -143,7 +145,6 @@ class SimpleExpansionOrderFinder:
                     * helmholtz_k
                     / (2*float(np.pi)))
 
-            from math import factorial
             helm_order = 1
             helm_rec_error = self.err_const_helmholtz * factor
             while True:
@@ -152,7 +153,7 @@ class SimpleExpansionOrderFinder:
                 if helm_order < 4:
                     # this may overflow for large orders
                     helm_error_direct = (
-                            1/factorial(helm_order+1)
+                            1/math.factorial(helm_order+1)
                             * self.err_const_helmholtz
                             * factor**(helm_order+1))
                     assert (abs(helm_rec_error - helm_error_direct)
diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py
index f760998dbd98284fd082453ade3e36d17ae9bfb1..72d91c01de07422de8c9cbc01c3f0c5a19acef04 100644
--- a/sumpy/expansion/local.py
+++ b/sumpy/expansion/local.py
@@ -20,17 +20,21 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-import sumpy.symbolic as sym
-from sumpy.tools import add_to_sac, fft
-
-from sumpy.expansion import (
-    ExpansionBase, VolumeTaylorExpansion, LinearPDEConformingVolumeTaylorExpansion)
-
-from sumpy.tools import mi_increment_axis, matvec_toeplitz_upper_triangular
-from pytools import single_valued
+import math
 from typing import Tuple, Any
+
 import pymbolic
 import loopy as lp
+from pytools import single_valued
+
+import sumpy.symbolic as sym
+from sumpy.expansion import (
+        ExpansionBase,
+        VolumeTaylorExpansion,
+        LinearPDEConformingVolumeTaylorExpansion)
+from sumpy.tools import (
+        add_to_sac, fft,
+        mi_increment_axis, matvec_toeplitz_upper_triangular)
 
 import logging
 logger = logging.getLogger(__name__)
@@ -236,9 +240,8 @@ class LineTaylorLocalExpansion(LocalExpansionBase):
 
     def evaluate(self, tgt_kernel, coeffs, bvec, rscale, sac=None):
         # no point in heeding rscale here--just ignore it
-        from pytools import factorial
         return sym.Add(*(
-                coeffs[self.get_storage_index(i)] / factorial(i)
+                coeffs[self.get_storage_index(i)] / math.factorial(i)
                 for i in self.get_coefficient_identifiers()))
 
 # }}}
@@ -561,7 +564,6 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
         # not coming from a Taylor multipole: expand via derivatives
         rscale_ratio = add_to_sac(sac, tgt_rscale/src_rscale)
 
-        from math import factorial
         src_wrangler = src_expansion.expansion_terms_wrangler
         src_coeffs = (
             src_wrangler.get_full_kernel_derivatives_from_stored(
@@ -649,7 +651,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
                         if src_mi in src_mi_to_index:
                             cur_dim_output_coeffs[out_i] += (dvec[d]/src_rscale)**q \
                                 * cur_dim_input_coeffs[src_mi_to_index[src_mi]] \
-                                / factorial(q)
+                                / math.factorial(q)
                 # Y at the end of the iteration becomes the source coefficients
                 # for the next iteration
                 cur_dim_input_coeffs = cur_dim_output_coeffs
diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py
index 4af5cec332e925cdd05a8d8bfc2331fe35a0b7db..fd5a3b1d1ac3ca8fb1f05bb48b9d845a0bf4ff97 100644
--- a/sumpy/expansion/multipole.py
+++ b/sumpy/expansion/multipole.py
@@ -20,12 +20,12 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-import sumpy.symbolic as sym  # noqa
+import math
 
+import sumpy.symbolic as sym
 from sumpy.expansion import (
     ExpansionBase, VolumeTaylorExpansion, LinearPDEConformingVolumeTaylorExpansion)
-from pytools import factorial
-from sumpy.tools import mi_set_axis, add_to_sac
+from sumpy.tools import mi_set_axis, add_to_sac, mi_power, mi_factorial
 
 import logging
 logger = logging.getLogger(__name__)
@@ -58,7 +58,6 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
         coefficients, compressing and then summing.
         """
         from sumpy.kernel import KernelWrapper
-        from sumpy.tools import mi_power, mi_factorial
 
         if not self.use_rscale:
             rscale = 1
@@ -124,8 +123,6 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
                     type(self).__name__,
                     self.order)
 
-        from sumpy.tools import mi_factorial
-
         src_mi_to_index = {mi: i for i, mi in enumerate(
             src_expansion.get_coefficient_identifiers())}
 
@@ -284,7 +281,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
                         contrib = cur_dim_input_coeffs[tgt_mi_to_index[input_mi]]
                         for n, k, dist in zip(tgt_mi, input_mi, dvec):
                             assert n >= k
-                            contrib /= factorial(n-k)
+                            contrib /= math.factorial(n-k)
                             contrib *= \
                                 sym.UnevaluatedExpr(dist/tgt_rscale)**(n-k)
 
diff --git a/sumpy/tools.py b/sumpy/tools.py
index b5943051bc5f85943fba22a5a33d3fc2ea76a1a4..14152d706b4f9f49ed61fadb1e5a42463f902350 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -59,10 +59,10 @@ def add_mi(mi1, mi2):
 
 
 def mi_factorial(mi):
-    from pytools import factorial
+    import math
     result = 1
     for mi_i in mi:
-        result *= factorial(mi_i)
+        result *= math.factorial(mi_i)
     return result