Reply to thread

Message: [QUOTE="Jannis Harder, post: 4495105"] 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: [URL]http://en.wikipedia.org/wiki/Divisor_function[/URL] [URL]http://en.wikipedia.org/wiki/Weird_number[/URL] Eric W. Weisstein. "Weird Number." From MathWorld--A Wolfram Web Resource. [URL]http://mathworld.wolfram.com/WeirdNumber.html[/URL] Eric W. Weisstein. "Semiperfect Number." From MathWorld--A Wolfram Web Resource. [URL]http://mathworld.wolfram.com/SemiperfectNumber.html[/URL] Eric W. Weisstein. "Perfect Number." From MathWorld--A Wolfram Web Resource. [URL]http://mathworld.wolfram.com/PerfectNumber.html[/URL] Eric W. Weisstein. "Divisor Function." From MathWorld--A Wolfram Web Resource. [URL]http://mathworld.wolfram.com/DivisorFunction.html[/URL] Eric W. Weisstein. "Mersenne Prime." From MathWorld--A Wolfram Web Resource. [URL]http://mathworld.wolfram.com/MersennePrime.html[/URL] Eric W. Weisstein. "Abundance." From MathWorld--A Wolfram Web Resource. [URL]http://mathworld.wolfram.com/Abundance.html[/URL] =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__ [/QUOTE]

Verification