From: Jannis Harder <jannis@...>
Date: 2005-12-05T03:22:23+09:00
Subject: Re: [QUIZ] Weird Numbers (#57) Solution

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