φの和
こちらの記事に（https://codeiq.jp/magazine/2014/08/13926/）
φの和を高速に求める方法が書いてあるのだが、
大事な部分だけを取り出してみた。
以下において、
n = Σφ(d) (d は n の約数) …①
という性質を用いる。
f(n) = Σφ(k) (k は 1〜n)、
g(n) = n * (n + 1) / 2
とおく。
例えば、
Σf([6 / k]) (k は 1〜6)
= f([6 / 1]) + f([6 / 2]) + f([6 / 3]) + f([6 / 4]) + f([6 / 5]) + f([6 / 6])
= φ(1) + φ(2) + φ(3) + φ(4) + φ(5) + φ(6)
+ φ(1) + φ(2) + φ(3)
+ φ(1) + φ(2)
+ φ(1)
+ φ(1)
+ φ(1)
= φ(1)
+ φ(1) + φ(2)
+ φ(1) + φ(3)
+ φ(1) + φ(2) + φ(4)
+ φ(1) + φ(5)
+ φ(1) + φ(2) + φ(3) + φ(6)
= 1 + 2 + 3 + 4 + 5 + 6 (∵①)
= g(6)
一般に、
g(n) = Σf([n / k]) (k は 1〜n) …②
が言える。
②より、
f（n) = g(n) - Σf([n / k]) (k は 2〜n)
また、f（1) = g(1) より、
f は g で表される。
これを実装すると以下のとおりになる。
@memo = Hash.new
def f(n)
return 1 if n == 1
s = n * (n + 1) / 2
(2..n).each{|k|
@memo[n / k] = f(n / k) if @memo[n / k] == nil
s -= @memo[n / k]
}
s
end
p f(7777777)
k が大きくなると、
f([n / k])が同じ値になることが多いので、
次のようにまとめると速くなる。
@memo = Hash.new
def f(n)
return 1 if n == 1
s = n * (n + 1) / 2
sq = Math.sqrt(n).to_i
# 前半部分は普通に計算
(2..n / sq).each{|k|
@memo[n / k] = f(n / k) if @memo[n / k] == nil
s -= @memo[n / k]
}
# 後半部分は同じ値になる部分をまとめて計算
(1..sq - 1).each{|k|
@memo[k] = f(k) if @memo[k] == nil
s -= @memo[k] * (n / k - n / (k + 1))
}
s
end
p f(7777777)
高速化できたが念のため、
オンライン整数列大辞典の
A002088（http://oeis.org/A002088/list）
と比較し、答え合わせしてみる。
高速化できたが念のため、
オンライン整数列大辞典の
A002088（http://oeis.org/A002088/list）
と比較し、答え合わせしてみる。
@memo = Hash.new
def f(n)
return 1 if n == 1
s = n * (n + 1) / 2
sq = Math.sqrt(n).to_i
# 前半部分は普通に計算
(2..n / sq).each{|k|
@memo[n / k] = f(n / k) if @memo[n / k] == nil
s -= @memo[n / k]
}
# 後半部分は同じ値になる部分をまとめて計算
(1..sq - 1).each{|k|
@memo[k] = f(k) if @memo[k] == nil
s -= @memo[k] * (n / k - n / (k + 1))
}
s
end
ary = (1..56).map{|i| f(i)}.unshift(0)
# OEIS A002088のデータ
ary0 =
[0,1,2,4,6,10,12,18,22,28,32,42,46,58,64,72,80,96,
102,120,128,140,150,172,180,200,212,230,242,270,
278,308,324,344,360,384,396,432,450,474,490,530,
542,584,604,628,650,696,712,754,774,806,830,882,
900,940,964]
# 一致の確認
p ary == ary0
0 件のコメント:
コメントを投稿
注: コメントを投稿できるのは、このブログのメンバーだけです。