Repository: more-itertools/more-itertools

Issue:
Consider making factor() more powerful
The current `factor()` from the itertools recipes uses trial division. That is reasonably effective up to about `10 ** 15` though the underlying `sieve()` eats memory like crazy.  

There are [other algorithms](https://stackoverflow.com/a/2274520/424499) that handle much larger inputs and use very little memory.  Here is a first draft using [Pollard's Rho Algorithm](https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm). It doesn't begin to bog down until `10 ** 30`.  

Further improvements are possible (listed in referenced article). I just wanted to get a sketch up and running to see if you would be interested in something like this.

**Main code**:

```
from more_itertools import factor as factor_small, is_prime
from math import gcd

def pollard(n):
    'Return a non-trivial factor of a composite n.'
    for b in range(1, n - 2):
        x = y = 2
        d = 1
        while d == 1:
            x = (x * x + b) % n
            y = (y * y + b) % n
            y = (y * y + b) % n
            d = gcd(x - y, n)
        if d != n:
            return d
    raise ValueError(f'prime or under 5')

def factor(n):
    todo = [n] if n > 1 else []
    while todo:
        n = todo.pop()
        if is_prime(n):
            yield n
        elif n < 1000:
            yield from factor_small(n)
        else:
            fact = pollard(n)
            todo += (n // fact, fact)
```

**Demonstration and tests**:

```
from collections import Counter
from random import seed, randrange
from math import prod

for n in range(100_000):
    assert Counter(factor(n)) == Counter(factor_small(n))

seed(8675309)
for i in range(20):
    n = randrange(10 ** 25)
    factors = list(factor(n))
    assert prod(factors) == n and all(map(is_prime, factors))
    print(f'{n:,}: {factors}')
```

**Output**:

```
9,869,602,325,083,211,630,749,747: [4136608051, 2385916722928897]
4,573,804,804,328,949,139,154,601: [3, 3, 7027, 36997, 1954783451491231]
3,985,828,998,088,783,078,306,188: [3, 3, 2, 2, 4999, 2411, 9186198265711447]
1,213,036,987,311,438,672,884,787: [3, 3, 11, 37, 13, 61, 15227, 27425186502559]
320,097,820,909,839,163,017,286: [2, 2063, 80143757, 968018810273]
2,056,565,487,725,687,980,272,202: [2, 2131, 266089, 1813435764650639]
6,491,522,496,841,641,382,228,938: [3, 7, 3, 2, 71329, 722287145707965947]
2,603,976,484,230,244,518,758,553: [3, 107, 1756663807, 4617887749799]
6,500,694,336,153,090,155,592,949: [6500694336153090155592949]
3,903,368,629,694,986,387,124,384: [2, 2, 2, 2, 2, 13, 19553, 479880205980464633]
7,630,023,212,177,997,234,899,609: [814519661, 9367512630463069]
4,284,851,646,165,198,181,579,572: [3, 2, 2, 11, 3049, 27308417, 389859308837]
2,076,514,068,333,688,913,860,647: [3, 29, 10137187, 2354497139782363]
6,125,479,605,228,782,433,229,343: [7, 109473107, 7993456466277907]
7,147,764,838,882,517,029,314,012: [3, 3, 3, 2, 2, 43, 103, 30013, 20873, 23853193909]
3,336,096,809,997,208,438,930,133: [7, 11, 3467, 88499, 424913, 332319401]
5,499,158,913,066,851,896,481,729: [83, 7537, 8790623147599316299]
5,988,544,417,700,083,510,486,399: [5988544417700083510486399]
1,647,403,876,764,101,672,307,088: [2, 2, 2, 2, 19, 23, 109471, 13571009, 158594251]
1,755,395,861,949,707,706,066,732: [3, 2, 2, 26463804541, 5527662822221]
```

**Open question**: 

Would it be worth giving up laziness in order to sort the output?

Discussion:
@bbayles: Hmm, I'm conflicted on this. It's better to have a more efficient algorithm, but this changes the order of output as well as the laziness property. I could add it as an additional function, but then why have the original aside from the fact that it was there first?

Perhaps I should explain that the existing recipe uses trial division and link to a well-regarded [alternative implementation](https://docs.sympy.org/latest/modules/ntheory.html#sympy.ntheory.factor_.pollard_rho)?@rhettinger: > but this changes the order of output as well as the laziness property. 

Small clarification:  It either loses output ordering or laziness _but not both_.  The above code is lazy but doesn't sort the output. 

> I could add it as an additional function

That would solve the 100% compatibility problem. However, given a choice between `factor()` and `much_better_factor()`, I don't think I would ever use the first one. Unlike most itertools, the current `factor()` code is unexpectedly memory intensive. Eliminating the memory requirement is a huge win.  Also the speed-up can be enormous — I get a 35x improvement for `1637927 * 4514317`, 95x for `54493 * 51563 * 63337`, and 249x for `4789 * 2999 * 8243 * 6217`.  And when you consider the vastly larger input range, the second function would always win.

> Hmm, I'm conflicted on this.

After more thought, I prefer keeping the ordered output. We really want deterministic output that doesn't change as the algorithms get better. Here are some options:

* Forgo laziness and return an iterator over the sorted factors. The rationale being that most use cases need all the factors.
* Let the user choose by extending API:  `factor(n, *, ordered=True)` or `factor(n, *, lazy=False)`
* Two separate functions. I'm -1 on this for the reason noted above.

I'm glad you're in the driver seat on this one. You have a great track record for making the right decisions when trade-offs need to be made. As always, I 100% support whatever you decide. You've been an excellent steward for this project.



    @rhettinger: Here's what the first option looks like:

```
def factor(n):
    todo = [n] if n > 1 else []
    factors = []
    while todo:
        n = todo.pop()
        if is_prime(n):
            factors += [n]
        elif n < 1000:
            factors += factor_small(n)
        else:
            fact = pollard(n)
            todo += (n // fact, fact)
    return iter(sorted(factors))
```

Solution patch:
diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml
index 01935a6..f2d6956 100644
--- a/.github/workflows/python-app.yml
+++ b/.github/workflows/python-app.yml
@@ -8,7 +8,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12", "3.13", "pypy-3.10"]
+        python-version: ["3.9", "3.10", "3.11", "3.12", "3.13", "pypy-3.10"]
 
     steps:
     - uses: actions/checkout@v4
@@ -22,16 +22,16 @@ jobs:
     - name: Run tests
       run: make coverage
     - name: Static checks
-      if: "matrix.python-version == '3.8'"
+      if: "matrix.python-version == '3.9'"
       run: make check
     - name: Build docs with sphinx
-      if: "matrix.python-version == '3.8'"
+      if: "matrix.python-version == '3.9'"
       run: make docs
     - name: Build packages
-      if: "matrix.python-version == '3.8'"
+      if: "matrix.python-version == '3.9'"
       run: make package
     - name: Upload packages
-      if: "matrix.python-version == '3.8'"
+      if: "matrix.python-version == '3.9'"
       uses: actions/upload-artifact@v4
       with:
         name: more-itertools-packages
diff --git a/docs/conf.py b/docs/conf.py
index 1e28765..11d08fc 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -10,9 +10,9 @@
 # All configuration values have a default; values that are commented out
 # serve to show the default.
 
-import sys, os
+import sys
+import os
 
-import sphinx_rtd_theme
 
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
diff --git a/more_itertools/more.py b/more_itertools/more.py
index 0e05986..9cd0140 100755
--- a/more_itertools/more.py
+++ b/more_itertools/more.py
@@ -218,8 +218,8 @@ def first(iterable, default=_marker):
         return item
     if default is _marker:
         raise ValueError(
-            'first() was called on an empty iterable, and no '
-            'default value was provided.'
+            'first() was called on an empty iterable, '
+            'and no default value was provided.'
         )
     return default
 
@@ -240,15 +240,14 @@ def last(iterable, default=_marker):
         if isinstance(iterable, Sequence):
             return iterable[-1]
         # Work around https://bugs.python.org/issue38525
-        elif hasattr(iterable, '__reversed__') and (hexversion != 0x030800F0):
+        if hasattr(iterable, '__reversed__'):
             return next(reversed(iterable))
-        else:
-            return deque(iterable, maxlen=1)[-1]
+        return deque(iterable, maxlen=1)[-1]
     except (IndexError, TypeError, StopIteration):
         if default is _marker:
             raise ValueError(
-                'last() was called on an empty iterable, and no default was '
-                'provided.'
+                'last() was called on an empty iterable, '
+                'and no default value was provided.'
             )
         return default
 
diff --git a/more_itertools/more.pyi b/more_itertools/more.pyi
index f8543b4..b33dcff 100644
--- a/more_itertools/more.pyi
+++ b/more_itertools/more.pyi
@@ -12,7 +12,6 @@ from typing import (
     ContextManager,
     Generic,
     Hashable,
-    Mapping,
     Iterable,
     Iterator,
     Mapping,
diff --git a/more_itertools/recipes.py b/more_itertools/recipes.py
index cf61fc6..ce8c2cf 100644
--- a/more_itertools/recipes.py
+++ b/more_itertools/recipes.py
@@ -1009,12 +1009,8 @@ def matmul(m1, m2):
     return batched(starmap(_sumprod, product(m1, transpose(m2))), n)
 
 
-def factor(n):
-    """Yield the prime factors of n.
-
-    >>> list(factor(360))
-    [2, 2, 2, 3, 3, 5]
-    """
+def _factor_trial(n):
+    # Factor using trial division, used for small inputs
     for prime in sieve(math.isqrt(n) + 1):
         while not n % prime:
             yield prime
@@ -1025,6 +1021,44 @@ def factor(n):
         yield n
 
 
+def _factor_pollard(n):
+    # Return a factor of n using Pollard's rho algorithm
+    for b in range(1, n - 2):
+        x = y = 2
+        d = 1
+        while d == 1:
+            x = (x * x + b) % n
+            y = (y * y + b) % n
+            y = (y * y + b) % n
+            d = math.gcd(x - y, n)
+        if d != n:
+            return d
+    raise ValueError('prime or under 5')
+
+
+def factor(n):
+    """Yield the prime factors of n.
+
+    >>> list(factor(360))
+    [2, 2, 2, 3, 3, 5]
+
+    This function uses trial division for when *n* is less than 1000 and
+    Pollard's rho algorithm for larger inputs.
+    """
+    todo = [n] if n > 1 else []
+    factors = []
+    while todo:
+        n = todo.pop()
+        if is_prime(n):
+            factors += [n]
+        elif n < 1000:
+            factors += _factor_trial(n)
+        else:
+            fact = _factor_pollard(n)
+            todo += (n // fact, fact)
+    return iter(sorted(factors))
+
+
 def polynomial_eval(coefficients, x):
     """Evaluate a polynomial at a specific value.
 
diff --git a/more_itertools/recipes.pyi b/more_itertools/recipes.pyi
index 1f73af9..63f641a 100644
--- a/more_itertools/recipes.pyi
+++ b/more_itertools/recipes.pyi
@@ -129,6 +129,8 @@ def reshape(
     matrix: Iterable[Iterable[_T]], cols: int
 ) -> Iterator[tuple[_T, ...]]: ...
 def matmul(m1: Sequence[_T], m2: Sequence[_T]) -> Iterator[tuple[_T]]: ...
+def _factor_trial(n: int) -> Iterator[int]: ...
+def _factor_pollard(n: int) -> int: ...
 def factor(n: int) -> Iterator[int]: ...
 def polynomial_eval(coefficients: Sequence[_T], x: _U) -> _U: ...
 def sum_of_squares(it: Iterable[_T]) -> _T: ...
diff --git a/pyproject.toml b/pyproject.toml
index 682f20e..b6a9be4 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -6,7 +6,7 @@ build-backend = "flit_core.buildapi"
 name = "more-itertools"
 authors = [{name = "Erik Rose", email = "erikrose@grinchcentral.com"}]
 readme = "README.rst"
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 license = {file = "LICENSE"}
 keywords = [
     "itertools",
@@ -24,7 +24,6 @@ classifiers = [
     "Natural Language :: English",
     "License :: OSI Approved :: MIT License",
     "Programming Language :: Python :: 3",
-    "Programming Language :: Python :: 3.8",
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
@@ -39,11 +38,12 @@ dynamic = ["version", "description"]
 
 [project.urls]
 Homepage = "https://github.com/more-itertools/more-itertools"
+Documentation = "https://more-itertools.readthedocs.io/en/stable/"
 
 [tool.flit.module]
 name = "more_itertools"
 
 [tool.black]
 line-length = 79
-target-version = ['py38']
+target-version = ["py39"]
 skip-string-normalization = true
diff --git a/requirements/development.txt b/requirements/development.txt
index a879a9e..03c2f8f 100644
--- a/requirements/development.txt
+++ b/requirements/development.txt
@@ -2,7 +2,7 @@ black
 coverage
 flake8
 flit
-mypy
+mypy==1.13
 setuptools
 sphinx
 sphinx_rtd_theme


Test patch:
diff --git a/tests/test_recipes.py b/tests/test_recipes.py
index 6f8bfaf..c31d2b1 100644
--- a/tests/test_recipes.py
+++ b/tests/test_recipes.py
@@ -1153,8 +1153,12 @@ class FactorTests(TestCase):
             (6, [2, 3]),
             (360, [2, 2, 2, 3, 3, 5]),
             (128_884_753_939, [128_884_753_939]),
-            (999953 * 999983, [999953, 999983]),
-            (909_909_090_909, [3, 3, 7, 13, 13, 751, 113797]),
+            (999_953 * 999_983, [999_953, 999_983]),
+            (909_909_090_909, [3, 3, 7, 13, 13, 751, 1_137_97]),
+            (
+                1_647_403_876_764_101_672_307_088,
+                [2, 2, 2, 2, 19, 23, 109471, 13571009, 158594251],
+            ),
         ):
             with self.subTest(n=n):
                 actual = list(mi.factor(n))