Trial Division

October 25, 2016

It was a gorgeous autumn weekend: brisk mornings, warm sunny afternoons, and absolutely no reason to sit at a computer screen indoors. And even though the programming on my big project at work is finished, there are still details to be looked after, so I had no time there, either. Thus, another recycled exercise today, another one that I’ve seen on the beginning-programmer message boards recently:

Write a program that determines the prime factors of a number n.

You should write at the level of a first-year programmer, nothing sophisticated. An appropriate algorithm is trial division.

Your task is to write a program to factor integers. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

Posted by programmingpraxis
Filed in Exercises
2 Comments »

2 Responses to “Trial Division”

  1. Alex B said
    October 25, 2016 at 11:18 AM

    Very basic Python 3 solution that divides out by 2 and then by each odd number until the reduced n becomes less than the square of the factor to be tested.
    I used Python’s builtin divmod instead of separate modulo and division.

    def factors_td(n):
    '''Return a list of the prime factors of n'''
    factors = []

    while n % 2 == 0:
        factors.append(2)
        n //= 2

    f = 3
    while n > f*f:
        q, r = divmod(n, f)
        if r == 0:
            factors.append(f)
            n = q
        else:
            f += 2
    if n != 1:
        factors.append(n)
    return factors
  2. matthew said
    October 25, 2016 at 8:33 PM

    Here’s another Python solution, to make things more interesting, it factorizes Gaussian integers (ie. integral complex numbers):

    import heapq

# Usual complex division but only for integers
def div(z,w):
    (a,b),(c,d) = z,w
    e = a*c + b*d
    f = b*c - a*d
    k = c*c + d*d
    if e%k == 0 and f%k == 0: return (e//k,f//k)

# norm(a*b) = norm(a)*norm(b)
def norm(z):
    (a,b) = z
    return a*a + b*b

def factorize(z):
    h, factors = [],[]
    # Heap of all pairs, ordered by norm
    def addentry(w): heapq.heappush(h,(norm(w),w))
    def extend(z):
        (a,b) = z
        addentry((a,b+1))
        if b == 0: addentry((a+1,b))
    extend((1,0))
    while True:
        (n,w) = heapq.heappop(h)
        if n*n > norm(z): break
        extend(w)
        while True:
            quot = div(z,w)
            if quot is None: break
            factors.append(w)
            z = quot
    if z != (1,0): factors.append(z)
    return factors

for a in range(1,5):
    for b in range(0,5):
        z = (a,b)
        print("%s: %s"%(z,factorize(z)))

    (1, 0): []
(1, 1): [(1, 1)]
(1, 2): [(1, 2)]
(1, 3): [(1, 1), (2, 1)]
(1, 4): [(1, 4)]
(2, 0): [(1, 1), (1, 1), (0, -1)]
(2, 1): [(2, 1)]
(2, 2): [(1, 1), (1, 1), (1, 1), (0, -1)]
(2, 3): [(2, 3)]
(2, 4): [(1, 1), (1, 1), (2, -1)]
(3, 0): [(3, 0)]
(3, 1): [(1, 1), (2, -1)]
(3, 2): [(3, 2)]
(3, 3): [(1, 1), (3, 0)]
(3, 4): [(2, 1), (2, 1)]
(4, 0): [(1, 1), (1, 1), (1, 1), (1, 1), (-1, 0)]
(4, 1): [(4, 1)]
(4, 2): [(1, 1), (1, 1), (1, -2)]
(4, 3): [(1, 2), (1, 2), (0, -1)]
(4, 4): [(1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (-1, 0)]

Leave a comment