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.
Programming Praxis 
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 factorsHere’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)))