Python port of the factor command from GNU Coreutils. Supports both Python 2 and 3.

factor.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Teal Dulcet

# Run: python3 -OO factor.py <numbers(s)>...
# Run: python2 -OO factor.py <numbers(s)>...
# Run: pypy3 -OO factor.py <numbers(s)>...
# Run: pypy -OO factor.py <numbers(s)>...

# python3 -X dev factor.py {2..100}
# diff <(factor {2..100}) <(python3 -X dev factor.py {2..100})

from __future__ import division, print_function, unicode_literals

import optparse
import sys
from collections import Counter

try:
    from math import gcd
except ImportError:
    from fractions import gcd

try:
    # Python 3.8+
    from math import isqrt
except ImportError:
    from math import sqrt

    def isqrt(x):
        return int(sqrt(x))


exponents = ("⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹", "⁻")  # 0-9, -

# debugging for developers.  Enables devmsg().
dev_debug = False

# Prove primality or run probabilistic tests.
flag_prove_primality = True

# Number of Miller-Rabin tests to run when not proving primality.
MR_REPS = 25


def primes(limit):
    if not limit & 1:
        limit -= 1
    size = (limit - 1) // 2
    sieve = bytearray((1,)) * size
    for i in range(isqrt(size) + 1):
        if sieve[i]:
            p = 3 + 2 * i
            j = (p * p - 3) // 2
            # sieve[j : size : p] = bytes(len(range(j, size, p)))
            sieve[j: size: p] = bytearray(len(range(j, size, p)))

    # nprimes = sieve.count(bytearray((1,)))
    nprimes = sum(sieve)
    primes = bytearray(nprimes)
    # magic = array('Q', [0] * nprimes)

    p = 2
    j = 0
    for i in range(size):
        if sieve[i]:
            ap = 3 + 2 * i
            primes[j] = ap - p
            # magic[j] = sys.maxsize // ap + 1
            j += 1
            p = ap

    p = limit
    is_prime = False
    while not is_prime:
        p += 2
        is_prime = True
        i = 0
        ap = 2
        while is_prime and ap * ap <= p:
            # if (magic[i] * p) % sys.maxsize < magic[i]:
            if not p % ap:
                is_prime = False
                break
            ap += primes[i]
            i += 1

    return primes, p


# primes_diff, FIRST_OMITTED_PRIME = primes(5000)
# primes_diff, FIRST_OMITTED_PRIME = primes(50000)
primes_diff, FIRST_OMITTED_PRIME = primes(1 << 16)  # 2^16 = 65536
# primes_diff, FIRST_OMITTED_PRIME = primes(500000)

# print(repr(tuple(primes_diff)), FIRST_OMITTED_PRIME)


PRIMES_PTAB_ENTRIES = len(primes_diff)
# assert PRIMES_PTAB_ENTRIES == 669 - 1
# assert FIRST_OMITTED_PRIME == 5003 # 4999
# assert PRIMES_PTAB_ENTRIES == 5133 - 1
# assert FIRST_OMITTED_PRIME == 50021  # 49999
assert PRIMES_PTAB_ENTRIES == 6542 - 1
assert FIRST_OMITTED_PRIME == 65537  # 65521
# assert PRIMES_PTAB_ENTRIES == 41538 - 1
# assert FIRST_OMITTED_PRIME == 500009 # 499979

# Output number as exponent


def outputexponent(number):
    base = 10
    anumber = abs(number)
    astr = ""

    while True:
        anumber, index = divmod(anumber, base)
        astr = exponents[index] + astr

        if anumber <= 0:
            break

    if number < 0:
        astr = exponents[-1] + astr

    return astr


def factor_using_division(t, factors):
    if dev_debug:
        print("[trial division] ", end="", file=sys.stderr)

    p = 0
    while not t & 1:
        t >>= 1
        p += 1
    if p:
        factors.update({2: p})

    p = 3
    i = 1
    while i <= PRIMES_PTAB_ENTRIES:
        if t % p:
            if i < PRIMES_PTAB_ENTRIES:
                p += primes_diff[i]
            i += 1
            if t < p * p:
                break
        else:
            t //= p
            factors.update([p])

    return t


def millerrabin(n, nm1, x, y, q, k):
    y = pow(x, q, n)

    if y in {1, nm1}:
        return True

    for _ in range(1, k):
        y = pow(y, 2, n)

        if y == nm1:
            return True
        if y == 1:
            return False
    return False


def prime_p(n):
    is_prime = False
    tmp = 0
    factors = Counter()

    if n <= 1:
        return False

    # We have already casted out small primes.
    if n < FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME:
        return True

    # Precomputation for Miller-Rabin.
    nm1 = n - 1

    # Find q and k, where q is odd and n = 1 + 2**k * q.
    k = 0
    q = nm1
    while not q & 1:
        q >>= 1
        k += 1

    a = 2

    # Perform a Miller-Rabin test, finds most composites quickly.
    if not millerrabin(n, nm1, a, tmp, q, k):
        return False

    if flag_prove_primality:
        # Factor n-1 for Lucas.
        tmp = nm1
        factor(tmp, factors)

    # Loop until Lucas proves our number prime, or Miller-Rabin proves our
    # number composite.
    for r in range(PRIMES_PTAB_ENTRIES):
        if flag_prove_primality:
            is_prime = True
            # for p, e in factors.items():
            for p in factors:
                tmp = pow(a, nm1 // p, n)
                is_prime = tmp != 1

                if not is_prime:
                    break
        else:
            # After enough Miller-Rabin runs, be content.
            is_prime = r == MR_REPS - 1

        if is_prime:
            return is_prime

        a += primes_diff[r]  # Establish new base.

        if not millerrabin(n, nm1, a, tmp, q, k):
            return False

    print("Lucas prime test failure.  This should not happen", file=sys.stderr)
    sys.exit(1)


def factor_using_pollard_rho(n, a, factors):
    x = 2
    z = 2
    y = 2
    P = 1
    t = 0

    if dev_debug:
        print("[pollard-rho ({0})] ".format(a), end="", file=sys.stderr)

    k = 1
    l = 1

    while n != 1:
        assert a < n
        while True:
            factor_found = False
            while True:
                x = ((x * x) % n) + a

                P = (P * (z - x)) % n

                if k % 32 == 1:
                    if gcd(P, n) != 1:
                        factor_found = True
                        break
                    y = x

                k -= 1
                if k == 0:
                    break

            if factor_found:
                break

            z = x
            k = l
            l *= 2
            for _ in range(k):
                x = ((x * x) % n) + a

            y = x

        while True:
            y = ((y * y) % n) + a
            t = gcd(z - y, n)

            if t != 1:
                break

        n //= t  # divide by t, before t is overwritten

        if not prime_p(t):
            if dev_debug:
                print("[composite factor--restarting pollard-rho] ",
                      end="", file=sys.stderr)
            t = factor_using_pollard_rho(t, a + 1, factors)
        else:
            factors.update([t])

        if prime_p(n):
            factors.update([n])
            break

        x %= n
        z %= n
        y %= n

    return n

# Use Pollard-rho to compute the prime factors of
# arbitrary-precision T, and put the results in FACTORS.


def factor(t, factors):
    if t:
        assert t >= 2
        t = factor_using_division(t, factors)

        if t != 1:
            assert t >= 2
            if dev_debug:
                print("[is number prime?] ", end="", file=sys.stderr)
            if prime_p(t):
                factors.update([t])
            else:
                t = factor_using_pollard_rho(t, 1, factors)


def outputfactors(number):
    if number < 1:
        print("Error: Number must be > 0", file=sys.stderr)
        return ""

    # strm = ""
    counts = Counter()
    factor(number, counts)

    if options.print_exponents:
        strm = ("{0}^{1}".format(prime, exponent) if exponent >
                1 else prime for prime, exponent in sorted(counts.items()))
    else:
        strm = sorted(counts.elements())

    return " ".join(map(str, strm))

    # for prime, exponent in sorted(counts.items()):
        # for j in range(exponent):
            # if strm:
                # strm += " × " if options.unicode else " * "
            # strm += str(prime)
            # if options.print_exponents and exponent > 1:
                # if options.unicode:
                    # strm += outputexponent(exponent)
                # else:
                    # strm += "^" + str(exponent)
                # break

    # return strm


def integers(token):
    try:
        n = int(token, options.frombase)
    except ValueError:
        print("Error: Invalid integer number: {0!r}".format(
            token), file=sys.stderr)
        return 1

    print("{0}: ".format(n), end="")

    if dev_debug:
        print("[using single-precision arithmetic] ", end="", file=sys.stderr)

    print(outputfactors(n))

    return 0


parser = optparse.OptionParser(
    version="%prog 1.0", description="Print the prime factors of each specified integer NUMBER. If none are specified on the command line, read them from standard input.")
parser.add_option("--from-base", dest="frombase", type="int",
                  default=0, help="Input in bases 2 - 36")
# parser.add_option("-u", "--unicode", action="store_true", dest="unicode", help="Unicode")
parser.add_option("--exponents", action="store_true", dest="print_exponents",
                  help="Print repeated factors in form p^e unless e is 1")
options, args = parser.parse_args()

if options.frombase and not 2 <= options.frombase <= 36:
    print("Error: Base must be 2 - 36.", file=sys.stderr)
    sys.exit(1)

if args:
    for arg in args:
        integers(arg)
else:
    for line in sys.stdin:
        for token in line.split():
            integers(token)