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

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

この問題は,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つの象限についてそれぞれ計算すればよい.

本日数時間、線分・直線間の距離を求める問題にチャレンジしてみたけど、結果は散々だった。自分の考えを一度整理しないとね。

ACM問題はこれ、#10709 Intersection is Not that Easy。線分や直線間の距離を求める問題なので、中学校程度の数学を勉強したひとなら、誰でも簡単に計算できるはずだが、私にはまだ解けない。

混乱を避けるため、用語についてまずその定義をはっきりしておこう。

線分 (Line Segment) とは、XY平面上の両点を結む線のこと。両端点の間にしか存在しない。なお、両端点は異なるもの(同一点でないという意味)でないといけない。

直線 (Line) とは、XY平面上の両点を通る線のこと。直線上、端点となるものは存在せず、平面上無限に両側に伸びていく。なお、直線の通る両点は異なるものでないといけない。

線分間の距離、または、直線間の距離、または、直線・線分間の距離とは、 線分(または直線)上の任意の点から、もう片方の線分(または直線)上の任意の点までの直線距離のなかで、もっとも値の小さいもの。

The distance between 2 lines (or line-segments) is the minimal value of distance between any 2 points, p1 and p2, where p1 ranges over all points in the first line (or line-segment), and p2 ranges over all points in the second line (or line-segment).

ちなみに書くと、直線間の距離は2次元平面においてはあまり意味のないものかもしれない。なぜなら、2直線が平行していなければ、必ず交差するので、両平行直線のケース以外には、2直線間の距離は常にゼロだからだ。

しかし、線分のことが出てくると、平面においても、線分間の距離が現実的になってくる。というよりも、現実の世界ではここでいう直線というものは存在しないので、すべては線分と考えるほうがむしろ正しいだろう。

さて、自分もだいぶ混乱してたので、場合分けして、それぞれの距離を整理していく。

1.両直線間の距離

両直線が平行していなければ、必ず交差するので、距離はゼロになる。両直線の平行条件とは、つぎの通り。(点の座標は図の通り)。

(式1)  (y2-y1)*(x4-x3) = (y4-y3)*(x2-x1)   (両直線が平行である必要十分条件)

計算メモ⇒ 座標値が整数の場合には、計算誤差は一切発生しない。ただし、整数同士の乗算になるので、計算結果のオーバーフローには注意すべし。座標値が32ビット整数の場合では、64ビット整数変数(Gnu Cでは long long int)を使っても危険性あり。

さて、両直線が平行していれば、点(x1, y1) から相手直線に垂線を下ろすことができる。点(x1, y1) と垂線の足 (x, y) との距離が平行直線間の距離になるので、平行直線間の距離はつぎの計算式で求められる。(^2 は2乗の意。sqrt はルート(平方根)計算。)

(式2)  Us = (x1-x3)*(x4-x3)+(y1-y3)*(y4-y3)
(式3)  Ub = (x4-x3)^2+(y4-y3)^2

(式4)  x = x3 + (x4-x3) * Us/Ub
(式5)  y = y3 + (y4-y3) * Us/Ub

(式6)  距離 dist = sqrt((x1-x)^2 + (y1-y)^2)

計算メモ⇒ 式2 において、Us の値が0であれば、式4・式5では、x = x3, y = y3 になる。また、Us = Ubであれば、x = x4, y = y4 になる。
 計算誤差のことを検討すると、式4と式5では除算が入って来るので、左辺の x と y は、一般的に整数ではなくなり、式6は誤差が伴うものになり、注意すべし。なお、直線という定義に満たすものであれば、式3がゼロになることはありえない、安心して式4と式5が使えるだろう。

では、最後に擬似コードを使って、まとめておこう。

 if 式1が偽 then
   distance = 0    /* 両直線は平行せず */
 else
   distance = 式6  /* 両直線は平行 */
 end.

2.線分間の距離
 とてもややしくなり、数個の数式で終わる話しではない。

2.1 両線分が交差するケース

ここでいう交差とは、向きの異なる両線分が、両端点の間(両端点を含まない)のところで交差することをいう。線分は一部または全部が重なることや、端点がもう片方の線分上にあるようなケースは含まれない。

両線分が交差する場合は、両線分間の距離がゼロになる。両線分が交差する条件とは以下の通り。

両線分交差の必要十分条件 =式11
(式7)   t1 = (y1-y3)*(x3-x4)-(x1-x3)*(y3-y4)
(式8)   t2 = (y2-y3)*(x3-x4)-(x2-x3)*(y3-y4)
(式9)   t3 = (y3-y1)*(x1-x2)-(x3-x1)*(y1-y2)
(式10)   t4 = (y4-y1)*(x1-x2)-(x4-x1)*(y1-y2)
(式11)   t1*t2 < 0 かつ t3*t4 < 0

計算メモ⇒ 座標値が整数であれば、計算誤差は一切発生しない。しかし、計算結果のオーバーフローに注意すべし。

2.2 両線分が互いに重なったり、タッチするケース

基本的には、式7~式9の計算結果のいずれかがゼロになることが必要だが、昔に書いた記事(点が線分上にあるかどうかの判定法)で正確に判定してくれる。

つまり、次の流れで判定を行う。
  1.点(x1, y1) が線分(x3, y3)-(x4, y4) 上にあるかどうかを判定する。
    線分上にあると判定すれば、6へ。
  2.点(x2, y2) が線分(x3, y3)-(x4, y4) 上にあるかどうかを判定する。
    線分上にあると判定すれば、6へ。
  3.点(x3, y3) が線分(x1, y1)-(x2, y2) 上にあるかどうかを判定する。
    線分上にあると判定すれば、6へ。
  4.点(x4, y4) が線分(x1, y1)-(x2, y2) 上にあるかどうかを判定する。
    線分上にあると判定すれば、6へ。
  5.本ケース(両線分が互いに重なったり、タッチするケース)に相当しないと判定し、7へ。
  6.本ケースに相当すると判定。7へ。
  7.判定終了。

両線分が互いに重なったり、タッチする本ケースでは、両線分間の距離はゼロである。

2.3 両線分が離れているケース
 さらに場合わけが必要かな。苦しい。。。

2.3.1 垂線の足が線分上に存在するケース
 線分の両端点のどれかが、もう片方の線分に垂線を下ろしたとき、その交点(垂線の足)が相手線分上にあるケース。

図にある青い線は垂線のこと。点線で描かれた垂線はその足が相手線分の外側になってしまうので、最短距離の計算には使えない。

垂線の足4本すべてが相手線分の外側にあるケースはつぎの 2.3.2 で考える。

左側に図示したケースのように、2端点ペア各々の距離よりも、点(x1, y1)と垂線の足(x, y) との距離のほうが短いので、垂線の長さが線分間の距離になる。しかし、右側ケースのように、必ずしも垂線の長さが最短ではない!

ということで、本ケースでは、線分間の距離とは以下の値のなかの最小値である。

 垂線の長さ。垂線の足が相手線分上にある場合はすべて算出しないといけない。最大4本。
 端点間、つまり、(x1,y1)と(x3,y3)ペア、(x2,y2)と(x3,y3)ペア、(x1,y1)と(x4,y4) ペア、(x2, y2)と(x4,y4)ペアの4ペア、の距離。

垂線の足が相手線分上にあるかどうかの判定は、上記の式2と式3を利用すればできる。
つまり、
       0≦Us≦Ub      ⇔  垂線の足が線分上にある
   Us < 0 または Us > Ub ⇔  垂線の足が線分の外側にある

また、点(x1,y2)から垂線の足 (x, y) までのの距離は上記の式6で計算できる。

2.3.2 垂線の足が線分の外側に存在するケース
 線分の両端点のどれもが、もう片方の線分に垂線を下ろしたとき、その交点(垂線の足)が相手線分の外側にあるケース。

本ケースでは、線分間の距離とは端点間(4ペア、以下)の距離のなかからの最小値である。

 端点間の距離、つまり、(x1,y1)と(x3,y3)ペア、(x2,y2)と(x3,y3)ペア、(x1,y1)と(x4, y4)ペア、(x2,y2)と(x4,y4)ペアの4ペア間の距離。

3.直線と線分間の距離
 上の2を理解しておけば、直線と線分との間の距離は問題なく計算できるはず。つまり、線分と直線とは交差するかしないか、線分の両端点からの垂線の長さはどうかで距離が判る。詳しい説明は省略していいだろう。

最後に、C言語ソースプログラムを載せておく。ACM問題を解答するものだが、重要な部分はそのまま利用できるだろう。

p10709.c

○ ○ ○

上の文章を書いている内に、自分がどこで過ちを犯したのか、気づいたので、プログラムのほうはすんなり書けた。文章に纏めることが、私にとって最善の整理法なのかもしれないね。

点 P(x,y) が線分 (x1, y1)-(x2, y2) 上にあるかどうかを判定する関数。

ここでいう線分とは、両端点をもつ有限長さの直線のこと。線分が実質点であっても、垂直でも水平でも、すべて正しく判定される。また、整数座標ならば、計算誤差は出ない。

/*
  引数: 点の座標 x, y
      線分の座標 (x1,y1) - (x2, y2)
  関数値: 点が線分上にある場合は1、線分上にない場合は 0
*/
int isPointOnSegment(int x, int y, int x1, int y1, int x2, int y2)
{
    int d;
    if (x1 > x2) {
        d = x1, x1 = x2, x2 = d;
        d = y1, y1 = y2, y2 = d;
    }
    return x1 <= x && x <= x2 &&
        ((y1 <= y2 && y1 <= y && y <= y2) ||
        (y1 > y2 && y2 <= y && y <= y1))
        && (y-y1)*(x2-x1) == (y2-y1)*(x-x1);
}

二点P1、P2 を通る 直線L があり、また他に二点 A、B があるとする。

A と B が直線の同一側、あるいは、反対側にあるかどうかの判定方法。

各点の座標: P(x1, y1)、 P2(x2, y2)、 A(xA, yA)、 B(xB, yB) とする。

P1、P2を通る直線Lの方程式は

   (y-y1) / (x-x1) = (y2-y1) / (x2-x1)

であり、分母を払うと、

   (x2-x1) * (y-y1) + (y2-y1) * (x1-x) = 0

となる。x と y に、A点とB点をそれぞれ代入し、左辺の値だけを計算する。

   SA = (x2-x1) * (yA-y1) + (y2-y1) * (x1xA)
   SB = (x2-x1) * (yB-y1) + (y2-y1) * (x1xB)

得られた SAとSBが違う符号であれば、直線の反対側にあり、同符号であれば、同一側にある。

         SA * SB &lt 0 ⇔ 反対側
         SA * SB &gt 0 ⇔ 同一側

直線が水平であろうと、垂直でろうと、気にする必要は全くない。なお、SA または SB は値がゼロの時には、その点が直線上にあることになる。

本判定法は整数同士の減算と乗算だけでできるので、誤差はなく、高速である。

全ての頂点に対する最短経路問題(All-Pairs Shortest Paths Problem、APSP問題)。つまり、いくつかの市があり、市と市との間の直接距離が分かったとき、ある市から他の市までの最短距離を求める問題。

市の代わりに駅を考え、距離の代わりに電車代を考えると、この問題は駅から駅までの最低乗り換え電車運賃を求める問題となる。したがって、応用範囲は実に広い。最近のカーナビもその一例だろう。

<関数名>
  apsp — 全ての頂点に対する最短経路問題

<形式>
  void apsp(int *shortest, int *distance, int numOfCity);

<引数>
  shortest  (出力)最短距離が入っている配列
  distance  (入力)市と市との直接距離
  numOfCity (入力)市の数

<関数値>
  なし

<注意事項>
  引数 shortest、distance は numOfCity X numOfCity の2次元配列のつもりで利用すること。市に番号 1 ~ n を割りふったとき、
   distance[0]         市1から市1まで(自分自身へ)の距離
   distance[1]         市1から市2までの距離
   distance[numOfCity-1] 市1から市nまでの距離
   distance[numOfCity]   市2から市1までの距離
   distance[numOfCity+1] 市2から市2(自分自身へ)の距離

となる。つまり、numOfCity分の配列要素ごとに、ある市から他のすべての市への距離がdistance配列に入っている。
 なお、自分自身への距離は常に0であり、市間の直接経路がないときは、無限大のつもりで大きい値を使うとよい。

<用例>
  C言語プログラム(apsp-test.c
  使い方ですが、まずデータファイルを用意しておく。市の数をデータファイルの1行目に書く。2行目からは1行ごとに各市間の距離を書く。2市間の直接経路がないときは、無限大のつもりで大きな数字(たとえば 20000)を使う。

   3    <-- 市の数
   0    <-- 市1から市1(自分自身)への距離
   20   <-- 市1から市2への距離
   20000  <-- 市1から市3への距離(直接経路がないので、無限大)
   20   <-- 市2から市1への距離
   0    <-- 市2から市2(自分自身)への距離
   10   <-- 市2から市3への距離
   25   <-- 市3から市1への距離(片通行の道はある)
   10   <-- 市3から市2への距離
   0    <-- 市3から市3(自分自身)への距離

  でデータファイルを書く。つぎに、用例プログラムは
    apsp-test データファイル
  で実行すれば、各市間の最短距離が表示される。

<C言語関数>
  (apsp.c

<説明>
  本プログラムには Floyd のアルゴリズムを使っている。実行時間は市の数の3乗オーダ。

平面上、整数座標で与えられた線分L と凸多角形P がある。

線分とは、有限長さの直線のこと。

凸多角形とは,多角形内部の任意の2点を結ぶ線分がつねにその多角形に含まれることをいい、凸でない多角形は凹多角形という。なお、ここでいう内部とは、多角形の頂点自体、あるいはその隣り合う頂点を結ぶ線分を含むものとする。つまり、多角形の境界線も内部と見なす。

凸多角形の判定はこのようにできる。つまり、頂点と頂点を結ぶ任意の対角線は,その多角形の内部に含まれる場合、多角形が凸である。

本関数では、与えられた線分Lが、完全に凸多角形Pの内部にある場合は 1と返事、線分Lの一部でも凸多角形Pの外部にある場合は0と返事する。

線分Lは両端が整数座標で与えられるものとする。

多角形Pは頂点数Nが3以上であり、頂点の座標を表す整数の配列 Xk, Yk (ただし、 0<=k<=N-1) が与えられるものとする。

k番目の頂点 (xk, yk) は k+1番目の頂点 (xk+1, yk+1) と腺分でつなぐ。最後のN番目の頂点は最初の1番目の頂点と線分で結ぶ。多角形の線分同士は頂点でしかぶつからない。1つの頂点から線分3本以上出ることはないものとする。

<関数名>
  segmentInsidePolygon —- 線分Lが凸多角形Pの内部にあるかどうかを判定する

<形式>
  int segmentInsidePolygon(int x1, int y1, int x2, y2, int pn, int *px, int *py);

<引数>
  x1, y1, x2, y2  入力 線分Lの両端の座標 (x1,y1), (x2,y2)
  pn   入力 凸多角形Pの頂点数
  px, py 入力 凸多角形Pの頂点の座標配列

<関数値>
  1  線分Lが完全に、凸多角形Pの内部、あるいは、ちょうど凸多角形Pの線分上にある場合。
  0  線分Lが一部でも凸多角形の外部にある場合。

<注意事項>
  凸多角形になっているかどうか、本関数内ではチェックしない。
  
用例

<C言語関数本体>
  segmentInsidePolygon.c

<説明>
  関数の内部では、線分Lの両端が凸多角形Pにあるかどうかの判定を利用している。

平面上、整数座標で与えられた点と多角形がある。

その点が、多角形の内部、あるいは、線分(エッジ)にある場合は 1と返事、多角形の外部にある場合は0と返事する。

点は整数座標で与えられるものとする。

多角形は頂点数Nが3以上であり(当たり前)、頂点の座標を表す整数の配列 Xk, Yk (ただし、 0<=k<=N-1) が与えられるものとする。

k番目の頂点 (xk, yk) は k+1番目の頂点 (xk+1, yk+1) と腺分でつなぐ。最後のN番目の頂点は最初の1番目の頂点と線分で結ぶ。線分同士は頂点でしかぶつからない。1つの頂点から線分3本以上出ることはないものとする。

<関数名>
  insidePolygon —- 点が多角形の内部にあるかどうかを判定する

<形式>
  int insidePolygon(int x, int y, int pn, int *px, int *py);

<引数>
  x, y  入力 点の座標
  pn   入力 多角形の頂点数
  px, py 入力 多角形の頂点の座標配列

<関数値>
  1  点が多角形の内部、あるいは、ちょうど多角形の線分上にある。
  0  点が多角形の外部にある。

<注意事項>
  多角形の頂点数N は3以上でないといけないが、誤って入力された場合では以下のように処理する。
  N < 1  無条件に 関数値=0 とする。
  N = 1  点の座標と比較し、一致すれば 関数値=1、一致しなければ 関数値=0 とする。
  N = 2  点の座標と線分と比較し、線分上にあれば 関数値=1、線分上になければ 関数値=0 とする。

  関数内の定数 JUST_ONの値をほかのもの(例えば2)にすれば、多角形の線分や頂点に点が乗っかっている場合には、関数値が2として正しく検出される。

  本関数の計算量はO(N)。

用例

<C言語関数本体>
  insidePolygon.c

<説明>
  判定の方法ですが、与えられた点 P から右の水平方向に直線 L を引き、多角形との交点の数を数える。交点の数が偶数(下図1のような、交点のないケースを 0、つまり、偶数とする)なら外部、奇数なら内部と判定する。下図の1,2,4は交点が偶数個、3は奇数個の例。