Knuth's Power Tree Method

指数の計算

a^nを素朴に計算するとn-1回の乗算が必要になる.例えばa^5なら

  1. a^2 = a * a
  2. a^3 = a^2 * a
  3. a^4 = a^3 * a
  4. a^5 = a^4 * a

しかし,少し工夫することによって乗算の回数は節約できる.例えばa^15は5回の乗算で計算できる.

  1. a^2 = a * a
  2. a^4 = a^2 * a^2
  3. a^5 = a^4 * a
  4. a^10 = a^5 * a^5
  5. a^15 = a^10 * a^5

ただ,この乗算パターンは単純な方法によっては生成できそうにない.

Binary Method

乗算の回数を節約するための単純な方法というのは例えば指数を倍々に増やしていき,目的の値を超えそうになったとき調整するやり方だ.この方法はbinary methodと呼ばれる.

  1. a^2 = a * a
  2. a^4 = a^2 * a^2
  3. a^8 = a^4 * a^4
  4. a^12 = a^8 * a^4
  5. a^14 = a^12 * a^2
  6. a^15 = a^14 * a

この乗算パターンを生成するやり方は分かりやすいが,最良ではない.この例では本来5回で済むところ6回の乗算を必要としている.
また,さらに問題を複雑にするのは,時として商を求める操作が計算量の節約に役立つためである.例えばa^31をa^32をaで割る事により得る場合,5回の乗算と1回の割り算で値を得られる.a=0の例外に対する処理がコスト0で,商を得る操作にかかるコストが乗算のコストの2倍未満であれば,a^31をより高速に計算できることになる.しかし,商を得る操作にかかるコストは大抵のマシンで高く付き,さらにa=0の例外に対するコストがかかる.指数がそれほど大きくないのであれば乗算のみを考慮する方が有用である.
Binary methodによりa^nを計算する場合,floor(log(n, base=2))*1回までは指数を倍々に増やせる.例えばa^15ならばfloor(log(n, base=2))は3で,1 -> 2 -> 4 -> 8という具合に3回まで倍々に指数を増やせる.4回目はa^16になるので目的の値を超える.
"調整"が必要な回数はnを2進数表現したときに含まれる1の個数から1を引いたものに等しい.例えば15の場合2進数表現は1111である.Binary methodでは8=1000(2)まで増やしたあと,0100(2)=4,0010(2)=2,0001(2)=1の順に指数を増やすので2進数表現に含まれる1の個数-1回の調整が必要となる.
a^nを計算する際,floor(log(n, base=2))=lambda(n),nの二進数表現に含まれる1の個数をv(n)とすると,binary methodでa^nを求めるのに必要な乗算の回数はlambda(n) + v(n) - 1となる.
また,lambda(n),v(n)については次の規則によって生成できる.

  • lambda(1) = 0, lambda(2n) = lambda(2n + 1) = lambda(n) + 1
  • v(1) = 1, v(2n) = v(n), v(2n + 1) = v(n) + 1

先に述べたように,binary methodは最良の乗算パターンを生成するとは限らない.しかしbinary methodについて考察することで少なくとも次のことがいえる.
I(n)をa^nを計算するのに必要な最小の乗算回数とする.このとき,

  • I(n) ≦ lambda(n) + v(n) - 1
  • I(n) ≧ ceiling(n) ((ceiling(n)は天井関数(x以上の最小の整数を求める)))
  • I(2^A) = A
  • I(2^A + 2^B) = A+1 (A > B)

である.
I(2n)=I(n)+1は多くの場合成立するが,絶対ではない.例えばI(191)=I(382)=11であり,このようにI(n)=I(2n)となるようなnは他に701, 743, 1111などがある.

The Power Tree Method

例えばa^3はa^1 * a^2により計算できる.a^3とそれまでに計算されている値との1回の乗算を加えて計算できる値は,a^4, a^5, a^6の3種類である.このように,a^nとa^nを計算するのに使用された"材料"のリストから,a^nにプラス1回の乗算で計算できる値のリストを求められる.この関係は木構造を作る.
a^1からはa^2のみが生成される.a^2からはa^3とa^4が生成される.a^3からa^4, a^5, a^6が,a^4からa^5, a^6, a^8が生成される.このように,階層が一回増えると,前の階層の葉の数×(前の階層の葉の数+1)だけ葉が増える.この方法では階層が少し深くなるだけで葉の数が膨大になり,大量のメモリを消費することになる.
そこで,nをメモした木に既に含まれる値が葉の候補として生成された場合,その値は葉に加えない事にしよう.このようにして生成された木をpower treeと呼ぶ.power treeを生成するアルゴリズムは次のようになる.

  1. 現在最も深い階層の葉を調査リストに加える
  2. そのリストの葉に対し小さいものから順に
    1. 根からその葉に対するパスに含まれる全てのノードの値を求める(根と注目している葉も含める)
    2. 注目している葉の値 + 今求めたノードの値を加える葉の候補とする
    3. 全ての候補に対し
      1. もし候補の値が木の中に含まれていないのであれば葉として加え,そうでなければ加えない

詳しくはThe Art of Computer Programming (2) 日本語版 Seminumerical algorithms Ascii Addison Wesley programming series*2
こうして作られた木の深さが葉の値を生成するのに必要な乗算の回数を表す.
この方法はヒューリスティックな方法である.多くの場合最良の乗算回数とパターンを求められるが,例えばn=77, 154, 233などが例外である*3
Rによる実装を次に示す.

library(igraph)
## Knuth's power tree method
## ノードeに次の葉を追加
nextLeaves <- function(g, e){
  nextl <- V(g)[unlist(get.shortest.paths(g,
                                          from=V(g)[V(g)$name==1],
                                          to=V(g)[V(g)$name==e]))]$name + e
  pass <- logical(length(nextl))
  for(i in nextl){
    if(is.element(i, V(g)$name)) next
    pass[nextl==i] <- TRUE
    g <- add.vertices(g, 1)
    V(g)[is.na(V(g)$name)]$name <- i
    g <- add.edges(g, c(V(g)[V(g)$name==e], V(g)[V(g)$name==i]))
  }
  return(list(g, nextl[pass]))
}

## get power tree
knuth.power.tree <- function(limit){
  g <- graph(c(0, 1))
  V(g)$name <- 1:2
  nextlist <- 2
  cost <- numeric(limit)
  i <- 1
  while(!is.element(limit, V(g))){
    cost[nextlist] <- i
    nextlist.new <- numeric(0)
    for(n in nextlist){
      tmp <- nextLeaves(g, n)
      g <- tmp[[1]]
      nextlist.new <- c(nextlist.new, tmp[[2]])
    }
    nextlist <- nextlist.new
    i <- i + 1
  }
  return(g)
}

たとえばn^15までのパスを含むpower treeを求めたい場合次のように入力する.

g <- knuth.power.tree(15)
plot(g, vertex.label=V(g)$name, layout=layout.reingold.tilford)

これにより得られるpower treeは次のようなものになる.

根から目的の値を持つ葉までのエッジの個数が必要な乗算の回数となる.n^15の場合パスは1→2→3→5→10→15であり,5回の乗算が必要であることが分かる.

*1:log(n, base=2)は2を底とするnの対数を表し,floor(x)は床関数(x以下の最大の整数を求める)を表す.ついでにRのコードにもなっている.

*2:ぼくもってない…今度買おう。

*3:n=5000程度までのpower treeを使用した場合に計算される乗算回数と最良の乗算回数のリストがhttp://www.research.att.com/~njas/sequences/a003313.txtにある.

Problem 122

さっぱり解ける気がしないので色々調べている過程でKnuthのpower tree methodというものの存在を知った.それでどうにか理解したものの残念なことにこの問題を単純なpower treeで解くことはできない.
結局Evaluation of Powers | COME ON CODE ONにあるアルゴリズムを使ったのだが,δ_nを計算するときに使っているリストがどのようにして求められているのか,なぜこのリストが必要なのかを理解できていない.とりあえず30msで答えは出たけど今ひとつという感じは否めない.

続きを読む

Problem 124

nの素因数の積をrad(n)とする.例えばrad(504) = 2*3*7 = 42である.
ある範囲について求めたrad(n)を,rad(n)を第1キー,nを第2キーとしてソートする.このときソート結果のk番目にあるnの値をE(k)とする.例えば1≦n≦10についてrad(n)を求めこの条件でソートした場合,E(4)=8,E(6)=9となる.
1≦n≦100000についてソートしたときのE(10000)を求めよ.
何のひねりもないただのbrute forceだけど解けてしまった.16.5秒かかったけど.これは単にgmpのfactorizeが優秀なだけではないか.

続きを読む