Like I've said in a previous post, I have a thing for prime numbers. Factoring falls along that same line. After I while I thought it would be interesting to write some code to return the prime factors of an input number, mostly because I wanted to see how fast my old machine would able to do it.
I decided on
Pollard's rho algorithm because it is easy to code and a whole lot faster than trial division. This would also need a primality (or compositeness) test, so since I had already coded one for my prime number generator, I stuck with
Rabin-Miller. There are some trial divisions in there, but it's just to catch the small factors.
The hardest part of this, for me, was coming up with a way for Python to handle the list of intermediate factors. The way it's coded works, but I am sure it could be done better. That actually applies for the rest of it as well. An earlier version of this tore though the
8th Fermat number in a reasonable amount of time, so this works for me.
Maybe I'll update this with a different factoring algorithm and significantly cleaner code as a follow-up project.
Here's the code:
import math
import random
import fractions
################################################
def removerepeats(alist):
alist.sort()
for element in alist:
while alist.count(element)>1:
alist.remove(element)
return(alist)
################################################
def milrabprimality(totest):
if totest==1:
return(-1)
smallprimes=[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211]
if totest in smallprimes:
return(1)
for smtest in smallprimes:
if totest%smtest==0 and totest!=smtest:
return(-1)
s=totest-1
twos=0
while s%2==0: #Figures out the seed and powers of two
s=s//2
twos+=1
eventest=totest-1 #Make even -> eventest=(2^r)*s
trials=10 #Change if you want to
agen=0
a=[]
while agen<trials:
newa=random.randint(5,25)
if newa not in a:
a.append(newa)
agen+=1
trial1count=0
totprime=0
while trial1count<trials: #Miller-Rabin test happens here
primecount=0
part1=0
#Part 1
res1=pow(a[trial1count],s,totest)
if res1==1:
primecount=1
part1=1
#Part2
twopow=twos-1
itercount=0
if part1==0: #Only do this loop if part 1 was unsuccessful
while twopow>=(0):
res2=pow(a[trial1count],s*(2**twopow),totest)
twopow-=1
itercount+=1
if res2==(totest-1):
primecount+=1
break
if primecount!=0:
totprime+=1
else:
return(-1) #This ends it early since it has to pass for all a
trial1count+=1 #increments the random base
if totprime==trial1count:
return(1)
################################################
def pollardwrapper(number):
factors=[]
nonprimefactors=[]
factors.append(number)
if milrabprimality(number)==1: #Just return it if prime
return(factors)
for factor in factors:
if milrabprimality(factor)!=1:
result1=pollardsrho(factor)
result2=factor//result1
factors.append(result1)
factors.append(result2)
nonprimefactors.append(factor)
for element in nonprimefactors:
factors.remove(element)
factors.sort()
factors=removerepeats(factors)
return(factors)
################################################
def pollardsrho(number):
if number==1:
return(1)
if number==2:
return(2)
if number==3:
return(3)
if number==4:
return(2)
prime=milrabprimality(number)
if prime==1:
return(number)
#Algorithm from the Crandall and Pomerance book
a=random.randint(1,number-3)
s=random.randint(0,number-1)
u=s
v=s
g=1
while g==1:
u=((u**2)+a)%number
v=((v**2)+a)%number
v=((v**2)+a)%number
g=fractions.gcd((u-v),number) #If g!=1 then g is a factor (unless it is equal to number).
if g==number: #The bad seed case
a=random.randint(1,number-3)
s=random.randint(0,number-1)
u=s
v=s
g=1
if g!=1:
return(g) #End of algorithm
################################################
tofactor=int(input('\nWhat number do you want to factor? '))
pollardfactors=pollardwrapper(tofactor)
print('The prime factors of ',tofactor,'are',pollardfactors)
.jpg)
.jpg)
.jpg)
