From ce9749898c887f96980d1deb897c7e9110e76ad5 Mon Sep 17 00:00:00 2001
From: mhostetter <matthostetter@gmail.com>
Date: Sun, 18 Dec 2022 16:08:53 -0500
Subject: [PATCH] Add check in `perfect_power()` to see if n is a prime power
 for small primes

See #443
---
 src/galois/_prime.py | 25 +++++++++++++++++++------
 1 file changed, 19 insertions(+), 6 deletions(-)

diff --git a/src/galois/_prime.py b/src/galois/_prime.py
index 98583ee68..e3cd8107f 100644
--- a/src/galois/_prime.py
+++ b/src/galois/_prime.py
@@ -887,7 +887,7 @@ def factors(n: int) -> tuple[list[int], list[int]]:
     if not n > 1:
         raise ValueError(f"Argument 'n' must be greater than 1, not {n}.")
 
-    # Step 0: Check if n is in the prime factors database
+    # Step 0: Check if n is in the prime factors database.
     try:
         p, e, n = PrimeFactorsDatabase().fetch(n)
         if n == 1:
@@ -897,13 +897,13 @@ def factors(n: int) -> tuple[list[int], list[int]]:
     except LookupError:
         p, e = [], []
 
-    # Step 1
+    # Step 1: Test is n is prime.
     if is_prime(n):
         p.append(n)
         e.append(1)
         return p, e
 
-    # Step 2
+    # Step 2: Test if n is a perfect power. The base may be composite.
     base, exponent = perfect_power(n)
     if base != n:
         pp, ee = factors(base)
@@ -912,12 +912,12 @@ def factors(n: int) -> tuple[list[int], list[int]]:
         e.extend(ee)
         return p, e
 
-    # Step 3
+    # Step 3: Perform trial division up to medium-sized primes.
     pp, ee, n = trial_division(n, 10_000_000)
     p.extend(pp)
     e.extend(ee)
 
-    # Step 4
+    # Step 4: Use Pollard's rho algorithm to find a non-trivial factor.
     while n > 1 and not is_prime(n):
         c = 1
         while True:
@@ -1035,6 +1035,18 @@ def _perfect_power(n: int) -> tuple[int, int]:
 
     abs_n = abs(n)
 
+    # Test if n is a prime power of small primes. This saves a very costly calculation below for n like 2^571.
+    # See https://github.com/mhostetter/galois/issues/443. There may be a bug remaining in perfect_power() still
+    # left to resolve.
+    for prime in primes(1_000):
+        exponent = ilog(abs_n, prime)
+        if prime**exponent == abs_n:
+            if n < 0:
+                if exponent % 2 == 0:
+                    return n, 1
+                return -prime, exponent
+            return prime, exponent
+
     for k in primes(ilog(abs_n, 2)):
         x = iroot(abs_n, k)
         if x**k == abs_n:
@@ -1051,7 +1063,8 @@ def _perfect_power(n: int) -> tuple[int, int]:
                     else:
                         return -c, e
 
-                # Failed to convert the even exponent to and odd one, therefore there is no real factorization of this negative integer
+                # Failed to convert the even exponent to and odd one, therefore there is no real factorization of
+                # this negative integer
                 return n, 1
 
             return c, e