Here is my solution (it's optimized for statistics not for speed).
I don't use things like someone tried all odd numbers up to 10**18
because doesn't change the algorithm.. I'll use that in my speed
optimized version..
#############
=begin
References:
http://en.wikipedia.org/wiki/Divisor_function
http://en.wikipedia.org/wiki/Weird_number
Eric W. Weisstein. "Weird Number." From MathWorld--A Wolfram Web
Resource. http://mathworld.wolfram.com/WeirdNumber.html
Eric W. Weisstein. "Semiperfect Number." From MathWorld--A Wolfram
Web Resource. http://mathworld.wolfram.com/SemiperfectNumber.html
Eric W. Weisstein. "Perfect Number." From MathWorld--A Wolfram Web
Resource. http://mathworld.wolfram.com/PerfectNumber.html
Eric W. Weisstein. "Divisor Function." From MathWorld--A Wolfram Web
Resource. http://mathworld.wolfram.com/DivisorFunction.html
Eric W. Weisstein. "Mersenne Prime." From MathWorld--A Wolfram Web
Resource. http://mathworld.wolfram.com/MersennePrime.html
Eric W. Weisstein. "Abundance." From MathWorld--A Wolfram Web
Resource. http://mathworld.wolfram.com/Abundance.html
=end
class Integer
$prime_factors = Hash.new # we cache prime factors...
def prime_factors position = -1
if cached = $prime_factors[self] # cached?
return cached # yes
end
if self == 1 # we have 1 we are done
return $prime_factors[self]=[] # return no factors
elsif position<0 # we havn't reached factor 5 yet
if self&1 == 0 # test factor 2
return $prime_factors[self]=[2,*(self>>1).prime_factors]
elsif self%3 == 0 # and factor 3
return $prime_factors[self]=[3,*(self/3).prime_factors]
end
end
loop do
position+=6 # increment position by 6
if position*position > self # we have a prime number return it
return $prime_factors[self]=[self]
elsif (quo,rem = divmod(position)) and rem.zero? # test 6n-1
return $prime_factors[self]=[position,*quo.prime_factors
(position-6)]
elsif (quo,rem = divmod(position+2)) and rem.zero? # and 6n+1
return $prime_factors[self]=[position+2,*quo.prime_factors
(position-6)]
end
end
end
def distinct_prime_factors # rle encode the prime factors ;)
distinct_prime_fac = Hash.new{0} # setup the hash
prime_factors.each do |prime_factor| # fill it
distinct_prime_fac[prime_factor]+=1
end
distinct_prime_fac.to_a.sort_by{|(fac,count)|fac} # and return
it as sorted array
end
def divisors # get the divisors (not needed for divisor sum)
divs = [] # store divisors here
n = 1 # start with 1
loop do
break if n*n > self # done
if (qua,rem = divmod(n)) and rem.zero? # test for division
divs << qua # add divisors
divs << n
end
n+=1
end
divs.uniq.sort[0..-2] # we don't want self
end
def semi_perfect? deficient=false # final test
cached_abundance = abundance
return deficient if cached_abundance < 0 # deficient return the
argument
return true if cached_abundance == 0 # perfect => semi perfect too
possible_values = {0=>true} # store all possible values in a hash
divs = self.divisors # get the divisors
div_sum_left = divs.inject(0){|a,b|a+b} # get the divisor sum
pos_2 = div_sum_left - self # this is a possibility too
divs.reverse.each do |div| # for each divisor
possible_values.keys.each do |value| # and each possible value
if value+div_sum_left < self # check wether it can reach the
number with the divisors left
possible_values.delete(value) # if not delete the number
(huge speedup)
end
new_value = value+div # we create a new possible value
including the divisor
if new_value == self or new_value == pos_2 # if it is the
number it's semi perfect
return true
elsif new_value < self # if it's less than the number it
could be semi perfect
possible_values[new_value]=true # add it to the possiblities
end # if it's more than the value we can ignore it
end
div_sum_left-=div # the current divisor isn't left anymore
end
false # we found no way to compose the number using the divisors
end
def restricted_divisor_sum # uses the formular from wikipedia
distinct_prime_factors.map do |(fac,count)|
comm_sum = 1
comm_mul = 1
count.times do
comm_sum += (comm_mul*=fac)
end
comm_sum
end.inject(1){|a,b|a*b}-self
end
def perfect? # perfect numbers have the form 2**(n-1)*(2**n-1)
where n is prime (and small ;) )
return false if self&1 == 1 # it isn't known weather there are
odd perfect numbers.. but it's irrelevant for my algorithm
doubled = self * 2 # the perfect number is a triangular number
of the form (n*(n+1))/2
doubled_root = Math.sqrt(doubled).floor # if we double it and
take the floored square root we get n
return false unless doubled == doubled_root*(doubled_root+1) #
if we don't get n it isn't perfect
doubled_root_string = (doubled_root+1).to_s(2) # from the first
line we know n+1 has to be of the form 2**n
return false unless doubled_root_string.count('1')==1 # we check
that here
return false unless
(doubled_root_string.size-1).prime_factors.size == 1 # and n ha to be
a prime
return false unless self.abundance == 0 # if it passes all the
earlier test we check it using the abundance
true # and if it passes it's perfect
end
def abundance
self.restricted_divisor_sum-self
end
end
require 'benchmark'
max_num = Integer(ARGV.shift||'1000') rescue 1000 # read a number
from the command line
new_semi_perfects = [] # store semi perfect numbers that can't be
constructed using other semi perfect numbers
STDOUT.sync = true
puts "Weird Numbers Up To #{max_num}:"
#begin_stat
perfect_count = 0
composed_semi_perfect_count = 0
new_semi_perfect_count = 0
weird_count = 0
deficient_count = 0
record_nums = (1..80).map{|n|max_num*n/80}
record_nums_left = record_nums.dup
next_record = record_nums_left.shift
recorded_times = []
init_time = [Benchmark.times,Time.now]
min_odd = 10**17
#end_stat
(1..max_num).each do |test_num| # counting loop
if test_num == next_record #stat
recorded_times << [Benchmark.times,Time.now] #stat
next_record = record_nums_left.shift #stat
end #stat
if test_num.perfect? # it's perfect
new_semi_perfects << test_num
perfect_count += 1 #stat
next
end
do_next = false
new_semi_perfects.each do |semi_per|
if test_num % semi_per == 0 # is it possible to compose the
current number using a semi-perfect number?
do_next = true # yes
composed_semi_perfect_count += 1 #stat
break
end
end
next if do_next
# no
case test_num.semi_perfect? nil # we don't care about deficient
numbers
when true # but we care about semi perfects
new_semi_perfects << test_num
new_semi_perfect_count += 1 #stat
when false # and even more about abundand non semi perfects
puts test_num
weird_count += 1 #stat
else #stat
deficient_count += 1 #stat
end
end
#end
#begin_stat
final_time = [Benchmark.times,Time.now]
digit_length = max_num.to_s.size
def form_float(num)
"%12s" % ("%im%#{7}ss" % [(num/60).floor,("%.4f" % (num %60))]).tr
(" ","0")
end
def rel(x,y)
"%#{y.to_s.size}i (%6s%%)" % [x,"%3.2f" % (x/y.to_f*100)]
end
puts "Stats"
puts "Time:"
puts "- User: "+form_float(ut = final_time.first.utime -
init_time.first.utime)
puts "- System: "+form_float(st = final_time.first.stime -
init_time.first.stime)
puts "- U+S: "+form_float(ut+st)
puts "- Real: "+form_float(final_time.last - init_time.last)
puts
puts "Numbers:"
puts "- Total "+rel(max_num,max_num)
puts "- Weird "+rel(weird_count,max_num)
puts "- Perfect "+rel(perfect_count,max_num)
puts "- Deficient "+rel(deficient_count,max_num)
puts "- Abundand "+rel(abundand_count = max_num-perfect_count-
deficient_count,max_num)
puts "- Semi-Perfect "+rel(abundand_count-weird_count,max_num)
puts " (w/o perfects)"
puts ""
puts "- Passed 1st "+rel(max_num-perfect_count,max_num)
puts " (perfect test)"
puts "- Passed 2nd "+rel(new_semi_perfect_count+weird_count
+deficient_count,max_num)
puts " (composed test)"
puts "- Passed 3rd "+rel(new_semi_perfect_count+weird_count,max_num)
puts " (deficient test)"
puts "- Uncomposed "+rel(new_semi_perfects.size,max_num)
puts " (semi-perfects that arn't a multiply of another semi-perfect)"
puts
if recorded_times.size >= 80
puts "Graphs:"
puts
process_plot = []
real_plot = []
first_ustime = init_time.first.utime+init_time.first.stime
first_realtime = init_time.last
recorded_times.each do |(process_time,realtime)|
process_plot << process_time.utime+process_time.stime -
first_ustime
real_plot << realtime - first_realtime
end
max_process = process_plot.last
step_process = max_process / 22
max_real = real_plot.last
step_real = max_real / 22
def plot(plot_data,max,step)
22.downto(0) do |k|
res = ""
res = form_float(max) if k == 22
while res.size != plot_data.size
val = plot_data[res.size]
lower_range = k*step
res << ( val >= lower_range ? (val < lower_range+step ?
'#' : ':') : ' ' )
end
puts res
end
end
puts "Y: Realtime X: Number"
plot(real_plot,max_real,step_real)
puts "%-40s%40s"% [1,max_num]
puts "Y: System+Usertime X: Number"
plot(process_plot,max_process,step_process)
puts "%-40s%40s"% [1,max_num]
puts
end
#end_stat
__END__
--
Jannis Harder