semi working code
import random
import fractions
def primality(number, iterations): #converted from http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Algorithm_and_running_time
n = number #number to test for primality
k = iterations #corresponds to line 'k, a parameter that determines the accuracy of the test
n_=n-1 #n-1 is substituted by n_, for later use when i need to write n - 1 as 2**s * d,
s=0
d=0
if n == 3 or n == 2:
return True
if n % 2 == 0:
return False
while n_ % 2 == 0: #writes n - 1 aka n_ as 2**s * d by dividing by two and checking if n_ is still even
s = s + 1
n_= n_/ 2
d = n_
while k != 0: #corresponds to line in pseudocode | Loop: repeat k times:
k = k - 1
a = random.randint(2, n-2) #picks random integer in range 2, n-2
x = pow(a,d,n) # sets x to a^d mod n
if x == 1 or x == n - 1: # corresponds to line in pseudocode| if x = 1 or x = n - 1 then do next Loop
continue
for r in xrange (1,s): #corresponds to line in pseudocode | for r = 1...s-1
x = pow(x, 2, n) #sets x to x^2 mod n
if x == 1: #corresponds to line | if x = 1 then return composite
return False
elif x == n - 1: #corresponds to line | if x = n - 1 then do next Loop
break
else: #corrresponds to line | return composite
return False
return True #corresponds to line | return probably prime
def pollard_brent(number): #returns a factor of an integer: code from http://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/#viewSource
if number % 2 == 0:
return 2
y = random.randint(1,number-1)
c = random.randint(1,number-1)
m = random.randint(1,number-1)
g = 1
r = 1
q = 1
while g == 1:
x = y
for i in range(r):
y = (pow(y,2,number)+c) % number
k = 0
while k < r and g == 1:
ys = y
for i in range(min(m, r-k)):
y = (pow(y,2,number)+c) % number
q = q*(abs(x-y)) % number
g = fractions.gcd(q, number)
k = k + m
r = r * 2
if g == number:
while True:
ys = (pow(ys, 2, number)+c) % number
g = fractions.gcd(abs(x-ys), number)
if g > 1:
break
return g
def prime_factorization(number, iterations): #returns list of prime factors of a number
n = number
primefactors = []
while primality(number, iterations) == False: # if number is not a prime number
if number == 1: # and if number is not 1
break
factor = pollard_brent(number) #factor the number
if factor != number and primality(factor, 10) == True: #if the factor is composite or equal to number then try again.
primefactors.append(int(factor)) # if factor is composite and not equal to number then add it to the list of prime factors
number = number / factor #set number to number/factor
if number != 1 and n != number: #if number is not 1 or n, then
primefactors.append(int(number)) #append it
return primefactors #return array with all the prime factors
def divisorlist(num): #returns list of proper divisors when given prime factors
#all the divisors of a number can be generated by multiplying the prime factors of the number in all possible ways
divisors=[]
f = lambda l: reduce(lambda z, x: z + [y + [x] for y in z], l, [[]]) #returns the set of all subsets of its argument http://wiki.python.org/moin/Powerful%20Python%20One-Liners
for list in f(prime_factorization(num,10)): #for all lists in the set of subsets that can be created from the list of the prime factors of num
n=1 #set n to one
for int in list: #then multiply all the numbers in each list
n = n * int
divisors.append(n) #and that is one possible divisor of num
if num in divisors: #remove num if it exists because only proper divisors needed
divisors.remove(num)
divisors_=[]
for int in divisors: #remove duplicates
if int not in divisors_:
divisors_.append(int)
divisors = divisors_
return divisors #return all the proper divisors
def subset(num): #figure out if a number is semiperfect or not
semiperfect = False
f = lambda l: reduce(lambda z, x: z + [y + [x] for y in z], l, [[]]) #creates all possible subsets,
for list in f(divisorlist(num)): #of the list of proper divisors
if sum(list) == num: #if the sum any list is equal to the number
semiperfect = True #then the number is semiperfect
return semiperfect
#weird number is both NOT semiperfect and abundant
for num in xrange(2,1000000):
if sum(divisorlist(num)) > num: #if the sum of proper divisors of num is larger than num then the number is abundant
if subset(num) is False: # if the num is NOT semiperfect subset(num) should be False
print num #if both of the above are true then the number is weird and print it.