格子点とは座標値が整数の点のことだ。原点から見える格子点の数はどれぐらいだろうか。

上の図では,黒丸(●)は見える格子点,白丸(○)は見えない格子点を表す。無限までに広がる平面の場合,その比率は案外簡単に求められるが,ある限られた領域ではなかなか高速に計算できないかもしれない.

この問題は,UVa問題集 10214 Trees in a Wood. に集められている.有限領域は |x| <= a, |y| <= b,ただし 1 <= a <= 2000,1 <= b <= 2000000 という. さて,簡単に証明できるので省略するが,x軸とy軸上以外の格子点 P(x, y) が原点から見えるための必要十分条件とは,x と y は互いに素(つまり,x と y の最大公約数は 1) である.

例えば,上図において,点(2, 2) では座標値 2と2 は互いに素でないので,原点からは見えない(事実,原点との間に点(1,1) が存在している).

こういう理屈で,数学的には,有限領域内のすべての格子点について,その座標値 (x, y) に対し,x と y が互いに素であるかどうかを判定しカウントしておけば,本問題の答えが出てくるので,小学生レベルでも解ける問題といえよう.

しかし,プログラム的に(格好よく言えば,アルゴリズム的に),10秒以内に答えを出すのはそう簡単ではない.ふつうのひとよりも頭が悪いとは思っていないが,私も解けるまで数ヶ月かかってしまった.

<ヒント>
 本問題の説明文(英語)にヒントとして書かれているのが,つぎのオイラー関数の存在.

オイラー関数 φ(n) とは 1 <= a < n である整数 a のうち,n と互いに素 となる a の個数のことである.

例として,1 <= a < 10 のうち,10 と互いに素である整数は 1, 3, 7, 9 の4つだけなので,φ(10)=4 となる.

このとき、このオイラー関数を次の式で表せるという次の定理がある。

nを素因数分解したものを n=p1a1p2a2…pak とすると、

  φ(n) = n(1-1/p1)(1-1/p2)…(1-1/p)

例えば、n = 10 なら 10 = 2*5より、φ(10)=10*(1-1/2)(1-1/5) = 4

さて,大事なことだが,オイラー関数をここでつぎのように拡張する.

拡張オイラー関数 φ(m, n) とは 1 <= a < m である整数 a のうち,n と互いに素 となる a の個数のことである.なお,拡張オイラー関数とは自分が勝手につけた名前.

例.φ(100, 3) = 集合 (1, 2, 4, 5, 7, 8, 10, 11, 13,14, 16, 17, …, 95, 97, 98, 100) の要素数 = 100-33 = 67

また,拡張オイラー関数の計算には似たような実用的で高速な方法がある.

nを素因数分解したものを n=p1a1p2a2p3a3p4a4…pak とすると、

  φ(m, n) = m – [m/p1] – [m/p2] -[m/p3] – [m/p4] – … – [m/p]
      +[m/(p1 p2)] + [m/(p1 p3)] + [m/(p1 p4)] + … + [m/(pk-1 pk)]
      -[m/(p1 p2 p3)] + [m/(p1 p2 p4)] + …+[m/(pk-2 pk-1 pk)]
      +[m/(p1 p2 p3 p4)] …

という形になる.ただし,[a/b] は a/bの商の部分を表す.つまり,C言語でいう整数同士の割り算のことだ,

上の展開式はわかりにくいものだが,本質的には

  φ(m, n) = m(1-1/p1)(1-1/p2)…(1-1/p)

の展開式と考えてよい.商の部分を正確に表現しようとすると訳のわからない式になってしまう.
 
例.φ(19, 10) について計算してみる.10 = 2*5 なので,
  φ(19, 10) = 19 – [19/2] – [19/5] + [19/(2*5)] = 19 – 9 – 3 + 1 = 8

これは,1~19のうち,10と互いに素である整数が 1, 3, 7, 9, 11, 13, 17, 19 の8個であることと一致する.

<実際の計算アルゴリズム>
 つぎの手順でプログラムをつくる.

1.整数 2 ~ 2000 について,それぞれを素因数分解し,テーブルにしまう.最小素数4つの積は  2*3*5*7 = 210 だが,そのつぎの5番目の素数 11 を入れると 2000 を越えるので,本問題では素因数の個数は最大 4 まで.
2.有限領域の大きさを表す入力値 a と b を得る.
3.見える格子点の累計変数を Kとし, K = b とする.(← x = 1 上の見える格子点の数)
4.整数 n = 2 ~ a のそれぞれについて,拡張オイラー関数 φ(b, n) を素因数テーブルを利用して計算し,累計変数に足していく.つまり K = K + φ(b, n)
5.ここまでの計算により第1象限内の見える格子点の結果が得られたので,最後に,領域全体の見える格子点の数に直す.第2項の +4 は x軸上の2点とy軸上の2点を意味する.
  K = 4*K + 4
6.Kを出力して上記の手順2に戻り,つぎの計算を行う.

参考のため,プログラムを下に置いておく.
C Source Solution for UVa P10214 Online Judge

<感想と反省>
 オイラー関数φ(n) の計算定理が洗練すぎて,定理に昇華させるまでのプロセスがなかなか見えず,結果的に数ヶ月も悩んでいた訳だ.しかし,

 x=1という直線上の見える格子点の数は?
 x=2という直線上の見える格子点の数は?
 x=3という直線上の見える格子点の数は?
 x=6という直線上の見える格子点の数は?

というふうに小学生レベルに戻して考えれば,拡張オイラー関数に早くたどり着いていたのかもしれない.

<残る課題>
 x軸の大きさが2000までなので,本問題は結果的にそれなりに高速に計算できた訳だ.事実,上記の解答プログラムをUVaサーバに判定してもらったら,計算時間は 0.006秒 と 報告 されており,現時点ではもっとも高速なプログラムとなった.しかし,x軸も 2000000 までに広げるならば,本方法でも遅すぎて話にならない.しかし,これ以上の速いアルゴリズムが存在するかどうかは不明.

なお,原点ではなく,ある特定の位置から見える周りの格子点の数も同様な方法で計算できる.なぜなら,その特定の位置を原点にすればいいからだ.つまり平行移動して,4つの象限についてそれぞれ計算すればよい.

Comments are closed.

Post Navigation