自分自身を除く約数の和が自分自身に等しい自然数は完全数と呼ばれる(e.g. 6 = 1 + 2 + 3)が、そうでないものについて、不足数、過剰数という概念がある。
ある自然数を2つの過剰数の和で表すことを考える。例えば12は最小の過剰数(1 + 2 + 3 + 4 + 6 = 16)なので、2つの過剰数で表せる最小の自然数は24だ。
ところで、このような自然数は無限にある。それどころか、28123より大きな自然数は全て2つの過剰数の和として表すことができることが知られている。
では、過剰数の和として表すことができない自然数の総和を求めよ、というのがこの問題。
最初bruteforceでやって約80秒だった。forum覗いて少し良いやり方を見つけたものの、劇的な改善はせず38秒。そしてやり方はわかるもののなぜこれで正しい答えが出るのかも理解しきれていない。計算時間については、他の言語の様子を見ると普通にネストしたfor回して200msとかになってるし、まあ仕方ないかなという気持ちもある。
ところで、28123という上限は20161まで下げられる。Wikipediaにも書いてある(過剰数 - Wikipedia)。ループがネストしているので計算量は2/3より大きく減らせて、計算時間は20秒になった。
以下解答。
# limit = 20161 limit = 28123 abdnumlist <- (1:limit)[sapply(1:limit, function(n) sum((1:n)[n %% (1:n) == 0]) - n) > 1:limit] sumofadb <- function(n){ low = 1 high = length(abdnumlist) sum = abdnumlist[low] + abdnumlist[high] while(low < high && sum != n){ if(sum > n){ high <- high - 1 } else { low <- low + 1 } sum = abdnumlist[low] + abdnumlist[high] } return(sum == n) } sum((1:limit)[!sapply(1:limit, sumofadb)]) end = Sys.time()
追記
python3
python3でやったら28123までで34秒、20161までで19秒でRと大差なし(書き方悪いのかもしれないけど)。sympy使わないともっと遅い。
# coding: utf-8 import time import sympy as sym limit = 20161 #limit = 28123 #limit = 1000 start = time.time() def is_adbnum(n): if(sum(sym.divisors(n)) > 2*n): return True else: return False adbnumlist = [i for i in range(1, limit+1) if is_adbnum(i)] def is_sumofadb(n): low = 0 high = len(adbnumlist)-1 s = adbnumlist[low] + adbnumlist[high] while low < high and s != n: if s > n: high -= 1 else: low += 1 s = adbnumlist[low] + adbnumlist[high] return s == n print(sum([i for i in range(1, limit+1) if not is_sumofadb(i)])) print(time.time() - start)
julia
試しにjuliaで書いてみたけど28123までで33秒、20161までで17秒。書き方が悪いのかもと思って色々試してみたけど改善できず。定数にconst宣言つけるだけで28123までで5秒になった。処理を関数にすると早いという話も聞いたが、今回の例ではconst宣言が最も効果があった。
const limit = 28123 # limit = 20161 function divsum(n) return sum([i for i in 1:(n-1) if n % i == 0]) end function make_abdnumlist(lim) return [i for i in 1:lim if divsum(i) > i] end const abdnumlist = make_abdnumlist(limit) function is_not_sumobadb(n)::Bool low = 1 high = length(abdnumlist) s = abdnumlist[low] + abdnumlist[high] while low < high && s != n if s > n high -= 1 else low += 1 end s = abdnumlist[low] + abdnumlist[high] end return s != n end function solve_problem(lim) sum([i for i in 1:lim if is_not_sumobadb(i)]) end @show solve_problem(limit)
C++
C++で書いてみたら28123までで1.8秒、20161までで0.9秒だった…
#include <cstdio> using namespace std; const int LIMIT = 28123; int divsum(int n){ int result = 0; for (int i = 1; i < n; i++) { if (n % i == 0){ result += i; } } return result; } bool is_not_sumofadb(int n, int *adbnumlist, int index){ int low = 0; int high = index; int s = adbnumlist[low] + adbnumlist[high-1]; while(low < high & s != n){ if(s > n){ high--; } else { low++; } s = adbnumlist[low] + adbnumlist[high]; } return s != n; } int main(){ int n = 12; int ans = 0; int adbnumlist[LIMIT]; int index = 0; for(int i = 1; i <= LIMIT; i++){ if(divsum(i) > i){ adbnumlist[index] = i; index++; } } for(int i = 1; i <= LIMIT; i++){ if(is_not_sumofadb(i, adbnumlist, index)){ ans += i; } } printf("%d\n", ans); }
さらに追記
コメントで教えていただいた方法に少し手を加えたら0.5秒程度にまで短縮できた。
limit <- 28123 #temp[i]がiの約数の総和になるようにする temp <- rep(0, limit) for(i in 1:limit){ temp[(1:(limit %/% i)) * i] <- temp[(1:(limit %/% i)) * i] + i } abdnumlist <- (1:limit)[temp > 2 * (1:limit)] sieve = logical(limit*2) # ちょうど2つの過剰数で表せる自然数のindexをTRUEに for(i in seq_along(abdnumlist)){ sieve[abdnumlist[1:i] + abdnumlist[i]] <- TRUE } sieve <- sieve[1:limit] sum((1:limit)[!sieve])
