BHH定理を応用した四分円の中の格子点の個数の概算

0.参考文献

[1]山本芳嗣・久保幹雄著(1997年2月20日初版第1刷) 巡回セールスマン問題への招待 株式会社 朝倉書店, 26

[2]https://manabitimes.jp/math/661 ピタゴラス数の求め方とその証明

1.BHH定理とは

BHH定理とは、Beardwood, Halton, Hammersleyの3名の研究者の頭文字をとって名付けられた定理です。定理の内容は以下の通りです。

「面積がSの長方形の中に、ランダムにn個の地点を用意する。全ての地点を経由する最短の巡回路の長さlは、β√(nS)に収束する。そのβは約0.72、あるいは約0.765だろう。」[1]

2.四分円の中の格子点の個数

円の方程式は、x^2+y^2=r^2で表せる。(x,y)はほとんど有理数(稠密性より)、あるいは実数であり、反対に格子点は有限個しかなく、rによって概算できそうです。

ここで重要なのは、円の方程式がピタゴラスの定理を応用して作られることと、そのピタゴラスの定理を成り立たせるピタゴラス数であることです。

3.原子ピタゴラス数の生成法

a^2+b^2=c^2を満たす原子ピタゴラス数の組み(a,b,c)を考えるとき、以下の性質が成り立ちます。

自然数m,n(m>n)が存在して、

a=m^2-n^2, b=2mn, c=m^2+n^2

これは、実際に代入すれば恒等式になることからわかります。

4. 原子ピタゴラス数の、四分円の中の格子点の個数問題への応用

まず、r=c, c-b=t^2を考えます。

この時、mとnの関係は、

m=n+t

となります。

この仮定の下、(a,b,c)を求めます。

a=t(2n+t), b=2n(n+t), c=2n^2+2nt+t^2

ここで対称性から、a<=bの範囲まで考えたら良いので、

t(2n+t)<=2n(n+t)

→t<=n・sqrt(2)

すなわち、nを決定し、tを1からn・sqrt(2)まで変化させる中での(a,b,c)の組みがわかれば、四分円の中の格子点の個数が数えやすくなります。

例えば、r=2c=10の場合、

半径が決まる→端点が決まる→2個確定します。

c=5→(n,t)=(1,1)

2a=6, 2b=8

つまり、(1+6)*(1+8)+(1+6)が確定しました。

2b=7,9の場合をそれぞれ別に数えると、

2b=7→8つ

2b=9→5つ

以上から、

2+(1+6)*(1+8)+(1+6)+8+5+5=90

となります。

実は、kc=r, k=const.の場合、cはn,tの2次のオーダーで大きくなるためこの方法でかなり多くの場合について適応することができますが、それでも場合分けをする必要があります。

5.BHH定理の、四分円の中の格子点の個数問題への応用(その1 応用の仕方)

そこで、個数Nと半径rを用いた公式N=N(r)なるものを見つけたいと思います。そこで、BHH定理を応用したいと思います。

まず、BHH定理で必要なのは、面積S、格子点の個数N、半径r、最短経路の長さlです。幸い、S、lは次のように表せます。

S=πr^2/4

l=a+b・sqrt(2), 0<=b<=2, a+b+1=N(※1)

(※1) 1つを除く格子点全てに、片手があると考える。原点から最も遠い(rだけ離れている)2点からそれぞれ出発して、格子点に触れていくことを考える。この時、ほとんどは距離1だけ手を伸ばせば届く。

ところが、偶奇の関係から、2つのルートがうまく結びつくためには斜めに手を伸ばす必要が出てくる。

この時、sqrt(2)だけ手を伸ばす必要がある。さらに、偶奇の関係から時々、最遠の2点からではなく、一方をその隣の点から出発させた方が最短経路になる場合がある、この際、使用しなかった1点を結ぶために、さらに斜めに手を伸ばす必要が出てくるのです。

以上から、BHH定理のβを未知数として、

2(N-1)=β・sqrt(N・πr^2), b=0

2(N-3+sqrt(2))=β・sqrt(N・πr^2), b=2

と表せます。このβがrについての関数であることを、数値実験を通して近似曲線を用いて示せたら、N=N(r)の公式がわかるようになります。

ところで一度、βについて解いてみましょう。

b=0→β=2(N-1)/sqrt(N・πr^2)

b=2→β=2(N-3+sqrt(2))/sqrt(N・πr^2)

両者の差は、

(2-sqrt(2))/sqrt(N・πr^2)

となり、rが大きくなるほど無視できるようになります。

近似曲線用に、両者の平均のβについて、rの関数として表してみます。

6.BHH定理の、四分円の中の格子点の個数問題への応用(その2 βの近似曲線)

実際にr=12までの格子点の個数と、βの計算値を表に示してみました。1に近いことから、今回の方法でそこそこの精度で個数が求まるだろうことが期待できますね。

実際に、横軸をr、縦軸をβとしてグラフ化してみました。

これを見ると、ar^(-b)+cの形状をしていることが推測できます。matlabで以下のコマンドをしました。

f=fit(r,β,'power2')

これにより、0.524r^(-1.068)+1.036に近似していることがわかりました。

近似曲線との比較をしました。

そこそこの近似ができていますが、所々でずれがあるようですね。そこで、β/(近似曲線の値)のグラフを作りました。

どうやら、周期関数のように見えます。広く一般に、周期関数はフーリエ級数で表せることが知られています。(より一般には、全ての関数がフーリエ級数で表せるとされています。)そこで、フーリエ級数による近似を考えます。(このyの値をg_1とします。)

とりあえず、a_0+b_0cos(t_0・r)+c_0sin(t_0・r)として考えます。matlabで以下のコマンドをしました。

f2=fit(r,g_1,'fourier1')

これにより、

1.1013+0.0389cos(0.8194・r)-0.04892sin(0.8194・r)

に近似していることがわかりました。

g_1とf2の比較をしてみました。

また、f*f2の値と比較してみました。

まずまずの近似ですね。同様にして、β/(近似曲線2の値)のグラフを作ってみました。

これもまた、周期関数です。同様にフーリエ級数による近似を考えます。以下、同様に繰り返すことで、より精度の良い近似ができると考えられます。しかし今回は、次で最後にします。(このyの値をg_2とします。)

a_1+b_1cos(t_1・r)+c_1sin(t_1・r)として以下のコマンドをしました。

f3=fit(circle.VarName9,circle3.VarName20,'fourier1');

結果、

1.004-0.0336cos(1.76・r)-0.0336sin(1.76・r)

を得られました。

g_2とf3の比較をしてみました。

また、f*f2*f3の値と比較してみました。

かなり良い近似ができているようです。

ちなみにβ/(近似曲線3の値)のグラフを作ってみました。

(1±0.04)倍以内には収まっているようなので、良いかと思われます。

以上から、βは以下のようにrの関数で表せます。

7. 個数nを半径rで表す

matlabで、以下のコマンドを打ちました。

S=solve(β==2*(n-2+sqrt(2)/2)/(r*sqrt(n*pi)),n,"ReturnConditions",true)

これを解いて、以下の値が返ってきました。ただし、n>0の仮定をしています。

val =

((pi*β^2*r^2 - 8*2^(1/2) + 32)^(1/2)/4 + (β*r*pi^(1/2))/4)^2

βに代入します！

はい。がんばりました。