平方根を使わずに距離を求める の履歴(No.10)
更新平面上の2点間の距離を平方根を使わずに計算する†
こちらの記事で、
http://d.hatena.ne.jp/nowokay/20120604
2点間の距離を平方根を使わず近似的に求める手法を知り、その精度について詳しく調べてみたのがこのページです。
最終的に、このようにすることで 約 1% の誤差で距離を近似できることが分かりました。
リンク先の内容を Mathematica で確かめる†
正確な距離関数はこう:
例えば点 と点 との間の距離を として求められる。
この関数をグラフにすると、
LANG:Mathematica distance[x_, y_] := Sqrt[x^2 + y^2] Plot3D[distance[x, y], {x, -1, 1}, {y, -1, 1}, ImageSize -> Large]
最大値、最小値を利用した近似関数のグラフはこう。
LANG:Mathematica distance2[x_, y_, a_, b_] := a Max[Abs[x], Abs[y]] + b Min[Abs[x], Abs[y]] Plot3D[distance2[x, y, 1007/1024, 441/1024], {x, -1, 1}, {y, -1, 1}, ImageSize -> Large]
プロットには、上記リンク先で紹介されていた係数 を使った。
近似値から真値を引いた誤差を表示すると:
LANG:Mathematica Plot3D[distance2[x, y, 1007/1024, 441/1024] - distance[x, y], {x, -1, 1}, {y, -1, 1}, ImageSize -> Large, PlotPoints -> 80]
誤差は -2% から +8% くらいまでの間に分布する。
どうやらこのパラメータは近似値が真の距離を大きく下回らないように選ばれているようです?
とはいえ、 近辺では少々下回っているという謎仕様。。。
リンク先のコメント欄で紹介されている max + (3/8) min という非常に単純化された処理でも 誤差は -2.5% から +7% の間に収まるので、ちょっとこれはいただけない。
パラメータの最適化を考える†
エラーの相対誤差の二乗を最小化してみる。
LANG:Mathematica error[a_, b_] := NIntegrate[ (distance2[x, y, a, b] - distance[x, y])^2/ distance[x, y]^2, {x, 10^(-5), 1}, {y, 10^(-5), 1}] ContourPlot[error[a/1024, b/1024], {a, 500, 1500}, {b, 0, 1000}, ImageSize -> Large]
横軸が 1024 a、縦軸が 1024 b で、青い部分が誤差が小さいところ。
最小のエラーを与えるパラメータを求める。
LANG:Mathematica FindMinimum[error[a/1024, b/1024], {a, 1000}, {b, 400}] Out[]= {0.000552491, {a -> 963.831, b -> 419.924}}
プロットしてみる。
LANG:Mathematica Plot3D[distance2[x, y, 964/1024, 420/1024] - distance[x, y], {x, -1, 1}, {y, -1, 1}, ImageSize -> Large, PlotPoints -> 80]
誤差は -6% から +2.5% の範囲に分布する。
誤差がゼロの上下に分散したため、上記のパラメータを使うより絶対値が小さくなっている。
ただし、「誤差の二乗の平均値」を最小にしたため、 の周辺の狭い範囲では大きなエラーが出るままだ。
最適値の周りでの誤差の変化を見てみる。
LANG:Mathematica ContourPlot[error[a/1024, b/1024], {a, 964 - 10, 964 + 10}, {b, 420 - 10, 420 + 10}, ImageSize -> Large]
左下から右上への成分と、左上から右下への成分とに分けられそう。
まず最小値付近を2次関数でフィッティングして、
LANG:Mathematica Fit[Flatten[Table[{a, b, error[(963.8306122312068` + a)/1024, (419.92395268555106` + b)/ 1024]}, {a, -1, 1}, {b, -1, 1}], 1], {a^2, a b, b^2, 1}, {a, b}] Out[]= 0.000552491 + 7.48994*10^-7 a^2 + 6.61035*10^-7 a b + 2.04662*10^-7 b^2
係数を成分に分解します。
LANG:Mathematica Eigensystem[{{7.48994, 6.61035/2}, {6.61035/2, 2.04662}}] // N Out[]= {{9.04982, 0.486736}, {{-0.904343, -0.426806}, {0.426806, -0.904343}}}
それぞれの成分に相当する関数をプロットすると、
固有値 9.0 の主成分は、
固有値 0.5 の副次成分は、
のような形をしている。
主成分のみでほぼ完全に距離関数を復元しているのに対して、
副次成分は
と
とで異なる符号を持ち、
一方を上げると一方が下がるために平均的な誤差への影響が小さいことが見て取れる。
誤差の最大値を最小化する†
誤差の評価を、
のようにして、誤差の最大値を最小化することを考える。
パラメータ a, b に対して誤差をプロットすると、
LANG:Mathematica error2[a_, b_] := Table[ With[{x = Cos[t], y = Sin[t]}, Abs[distance2[x, y, a, b] - distance[x, y]] ], {t, 0, Pi/4, Pi/4/1024} ] // Max ContourPlot[error2[a/1024, b/1024], {a, 500, 1500}, {b, 0, 1000}, ImageSize -> Large]
拡大して、
LANG:Mathematica ContourPlot[error2[a/1024, b/1024], {a, 984 - 10, 984 + 10}, {b, 407 - 10, 407 + 10}, ImageSize -> Large, PlotPoints -> 100] Plot3d[error2[a/1024, b/1024], {a, 984 - 10, 984 + 10}, {b, 407 - 10, 407 + 10}, ImageSize -> Large, PlotPoints -> 100]
最小値を求める(始め から探したところうまく見つからなかった)。
LANG:Mathematica FindMinimum[error2[a/1024, b/1024], {a, 984}, {b, 407}] Out[]= {0.0395663, {a -> 983.485, b -> 407.372}} Plot3D[distance2[x, y, 983/1024, 407/1024] - distance[x, y], {x, -1, 1}, {y, -1, 1}, ImageSize -> Large, PlotPoints -> 80]
という近似で誤差は % の範囲に収まった。
※ 下で見るように、 でもほぼ 4% に収まる
もうちょっとちゃんと考えてみる†
もともと、どうして min と max で距離が求まるか、というと、、
&math( d(x,y) &=\sqrt{x^2+y^2}\\ &=\sqrt{|x|^2+|y|^2}\\ &=\sqrt{|y|^2+|x|^2}\\ );
という関係から、
&math( \begin{cases} x'=\mathrm{max}(|x|,|y|)\\ y'=\mathrm{min}(|x|,|y|)\\ \end{cases} );
と置けば、
が成り立つところがキモなのだった。
この は、 かつ、 より、 第一象限の 以下の部分に存在する点となる。
LANG:mathematica Show[{ Plot[x, {x, 0, 10}, Filling -> Axis, AspectRatio -> 1, PlotRange -> {{-2, 2}, {-2, 2}}], ListPlot[{ {3, 1}, {3, -1}, {-3, 1}, {-3, -1}, {1, 3}, {-1, 3}, {1, -3}, {-1, -3} }/2] }]
図中の8個の点はすべて、網掛け領域内の (3/2,1/2) と同じだけ原点と離れている、というのがこの意味である。
で、最終的に max と min との一次関数で表そうというのは、 この領域の距離関数を平面で近似しようということになる。
LANG:mathematica Show[{ Plot3D[{Sqrt[x^2 + y^2]}, {x, -1, 1}, {y, -1, 1}, RegionFunction -> (#1 > 0 && #2 > 0 && #1 > #2 && #3 < 1 &), Mesh -> None], Plot3D[{Sqrt[x^2 + y^2]}, {x, -1, 1}, {y, -1, 1}, RegionFunction -> ( (#1 < 0 || #2 < 0 || #1 < #2) && #3 < 1 &), PlotStyle -> Directive[Opacity[0.2], Orange, Specularity[White, 50]], Mesh -> None] }] Plot3D[{Sqrt[x^2 + y^2]}, {x, 0, 1}, {y, 0, 0.71}, RegionFunction -> (#1 > 0 && #2 > 0 && #1 > #2 && #3 < 1 &), Mesh -> None]
正確な距離関数は本来円錐状なのだけれど、 1/8 だけ取り出してみるとたしかにかなり平面っぽい。
これを本当の平面で近似した雰囲気はこんな感じ。
LANG:mathematica Plot3D[ {Sqrt[x^2 + y^2], 123/128 x + 51/128 y}, {x, 0, 1}, {y, 0, 0.74}, RegionFunction -> (#1 > 0 && #2 > 0 && #1 > #2 && #3 < 1 &), Mesh -> None ]
本当はたわんでいる黄色のグラフをビシッと伸ばしてしまい青の三角形で近似する。
との交線をプロットすると、
LANG:mathematica Plot[{(1 - 123/128 x)/(51/128), Sqrt[1 - x^2], x}, {x, 0, 1.2}, RegionFunction -> (#1 >= #2 &), PlotRange -> {{0, 1.1}, {0, 0.75}}, AspectRatio -> Automatic, PlotStyle -> {Automatic, Automatic, Dotted} ]
正しいグラフは の 1/8 円弧になるけれど、 近似グラフは直線になる。
昔から粗い近似として使われていたとリンク先のコメントで書かれていた であればこうなる。
LANG:mathematica Plot[{(1 - x)/(3/8), Sqrt[1 - x^2], x}, {x, 0, 1.2}, RegionFunction -> (#1 >= #2 &), PlotRange -> {{0, 1.1}, {0, 0.75}}, AspectRatio -> Automatic, PlotStyle -> {Automatic, Automatic, Dotted} ]
ずれてると言えばずれてるけど、まあそこそあってる。
こんな近似式をどこから導いたんだろう、と疑問だったのだけれど、 こうしてちゃんと考えればまずまず納得できる。(でもやっぱり賢いなあ)
近似精度の評価†
近似の精度を評価することを考えると、
&math( \begin{cases} x=r\cos\theta\\ y=r\sin\theta \end{cases} );
と置けば、
&math( \tilde d(x,y;a,b) &=a\,\mathrm{max}(|x|,|y|)+b\,\mathrm{min}(|x|,|y|)\\ &=|r|\,\Big[a\,\mathrm{max}(|\cos\theta|,|\sin\theta|)+b\,\mathrm{min}(|\cos\theta|,|\sin\theta|)\Big]\\ );
となることから、相対的な近似精度を見るには に固定して の関数として扱えば良い。
そこで、
を、上記で求めた「誤差の最大値を最小化するパラメータ」に対してプロットすると、
LANG:mathematica d2[t_, a_, b_] := a Max[Abs[Cos[t Pi]], Abs[Sin[t Pi]]] + b Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]] Plot[{1, d2[t, 983/1024, 407/1024]}, {t, 0, 2}, ImageSize -> Large]
となる。(横軸は を単位とした方位角)
今度は正しいグラフが直線に、近似グラフが曲線になっている。 (さっきまでとグラフの色が変わっていてわかりにくいので気をつけて)
どうしてこんな関数になるかを理解するため、 と とを別々にプロットしてみると、
LANG:mathematica Plot[{Max[Abs[Cos[t Pi]], Abs[Sin[t Pi]]], Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]]}, {t, 0, 2}]
青が max、黄色が min。sin と cos の大きい方、小さい方がそれぞれ max, min になる。
これらの線形結合で最も 1 に近い関数を作ったのが上記の結果になる。
max も min も の範囲が反転しながら繰り返し現れるから、 誤差を最小化するには の範囲だけ考えれば良い。 これは、上で 1/8 の領域だけ考えればいいと書いたのと対応する。。
LANG:mathematica Plot[{Max[Abs[Cos[t Pi]], Abs[Sin[t Pi]]], Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]]}, {t, 0, 1/4}]
このプロットから、このパラメータに対する誤差の最大値は約 4% であることが分かる。
さらに誤差を小さくするには、例えば min の2倍の周期の三角波を
のようにして作って、
LANG:mathematica Plot[{ Max[Abs[Cos[t Pi]], Abs[Sin[t Pi]]], Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]], Abs[Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]]-49/128 d2[t,a/1024,b/1024]] }, {t, 0, 2}]
LANG:mathematica With[{a = 938, b = 384}, Plot[{1, d2[t, a/1024, b/1024] + 55/256 Abs[ Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]] - 49/128 d2[t, a/1024, b/1024]]) }, {t, 0, 1/4}, PlotRange -> {Full, {0.98, 1.02}}, GridLines -> Automatic] ]
とすると、最大誤差を 1% 程度に抑えられる。
ただ、計算量が増えてしまうのでもう少し何とかならないか?
計算量を減らす†
計算量を増加させずに精度を上げるには、 リンク先と同様に min と max の大小関係で分岐するのが良さそう。
LANG:mathematica Plot3D[Sqrt[x^2 + y^2], {x, -1, 1}, {y, -1, 1}, ImageSize -> Large, PlotPoints -> 100] Plot3D[Module[{max, min}, max = Max[Abs[x], Abs[y]]; min = Min[Abs[x], Abs[y]]; If[max < (5/2) min, 864 max + 569 min, 1016 max + 190 min]/1024 ], {x, -1, 1}, {y, -1, 1}, ImageSize -> Large, PlotPoints -> 100 ]
左が正確な距離関数、右が近似関数。
相対誤差を評価すると、
LANG:mathematica Plot[{1, Module[{x, y, max, min}, x = Cos[t Pi]; y = Sin[t Pi]; max = Max[x, y]; min = Min[x, y]; If[max < (5/2) min, 864 max + 569 min, 1016 max + 190 min]/1024 ]}, {t, 0, 1/4}, GridLines -> Automatic, PlotRange -> {0.98, 1.02} ]
横軸は を単位とする方位角、縦軸は相対誤差です。
max と min とを混ぜ合わせる簡単な処理で 1% 程度の精度で距離関数を再現できることが分かります。
コードにするとこんな感じ。
LANG:cpp int distance2(int x, int y) { int absx = abs(x); int absy = abs(y); int max = absx > absy ? absx : absy; int min = absx > absy ? absy : absx; int result; if( 2 * max < 5 * min ) { result = 864 * max + 569 * min; } else { result = 1016 * max + 190 * min; } return (result + 512) / 1024; // 四捨五入 }
掛け算を使わないで書くとこんな感じ。
LANG:cpp int distance3(int x, int y) { int absx = abs(x); int absy = abs(y); int max = absx > absy ? absx : absy; int min = absx > absy ? absy : absx; int result; if( (max << 1) < ((min << 2) + min) ) { // 2 * max < 5 * min // 864 * max + 569 * min // 864 = 0x360 = 0b001101100000 // 569 = 0x239 = 0b001000111001 result = ((max << 10) - (max << 7) - (max << 5)) + ((min << 9) + (min << 6) - (min << 3) + min); } else { // 1016 * max + 190 * min // 1016 = 0x3f8 = 0b001111111000 // 190 = 0xbe = 0b000010111110 result = ((max << 10) - (max << 3)) + ((min << 8) - (min << 6) - (min << 1)); } return (result + 512) >> 10; // 四捨五入 }
こちらで
以下のコードを試したところ、 の範囲で 相対誤差の最大値は 1.07576% でした。
LANG:cpp #include <math.h> #include <stdio.h> double distance(int x, int y) { return sqrt((double)x*x+(double)y*y); } int distance2(int x, int y) { int absx = abs(x); int absy = abs(y); int max = absx > absy ? absx : absy; int min = absx > absy ? absy : absx; int result; if( 2 * max < 5 * min ) { result = 864 * max + 569 * min; } else { result = 1016 * max + 190 * min; } return (result + 512) / 1024; } int main() { int x, y; double max = 0; for(x=1000; x<10000; x+=10) { for(y=1000; y<10000; y+=10) { double d = distance(x, y); double d2 = distance2(x, y); if (max < abs(d-d2) / d) max = abs(d-d2) / d; } } printf("%lg", max); return 0; }
ビット幅などいろいろ工夫の余地は残る†
リンク先のコメント欄にある、
LANG:c int distance(int x, int y) { x = abs(x); y = abs(y); int max = x>y ? x:y; int min = x>y ? y:x; return max + ((min + min + min) >> 3); // max + (3/8)*min }
これだけでも、誤差は -2.5% から 7% 程度に収まる。
分母を 32 にすれば -5% ~ 6% の範囲にできる。
max にも係数を付ければ分母が 128 の段階で誤差をほぼ 4% 以内に収められる。
区間を分けない限り、これがほぼ最良なのは、この周辺での誤差分布をプロットすれば分かる。
LANG:mathematica ContourPlot[ Max[Table[ Module[{x, y, max, min}, x = Cos[t Pi]; y = Sin[t Pi]; max = Max[x, y]; min = Min[x, y]; Abs[(a max + b min)/128 - 1] ], {t, 0, 1/4, 1/400} ]], {a, 123 - 2, 123 + 2}, {b, 51 - 2, 51 + 2}, PlotPoints -> 400, PlotLegends -> Automatic, GridLines -> Automatic, Method -> {"GridLinesInFront" -> True} ] FindMinimum[ Max[Table[ Module[{x, y, max, min}, x = Cos[t Pi]; y = Sin[t Pi]; max = Max[x, y]; min = Min[x, y]; Abs[(a max + b min)/128 - 1] ], {t, 0, 1/4, 1/400} ]], {{a, 123}, {b, 51}} ] Out[]= {0.0396404, {a -> 122.926, b -> 50.9694}}
最小値は (a,b)=(122.926, 50.9694) の時の 0.03964
どうして三角形状の谷になっているかというと・・・
Table[ Module[{a, b}, a = 122.92602775001741` + 3 Cos[tt]; b = 50.96935967667149` + 3 Sin[tt]; GraphicsRow[{ Show[%110, ListPlot[{{a, b}}, PlotStyle -> White]], Plot[{1, Module[{x, y, max, min}, x = Cos[t Pi]; y = Sin[t Pi]; max = Max[x, y]; min = Min[x, y]; (a max + b min)/128] }, {t, 0, 1/4}, PlotRange -> {0.9, 1.1} ]}, ImageSize -> Large ] ], {tt, 0, 2 Pi, Pi/16} ]; Export["animation.gif", %%, "DisplayDurations" -> 0.2]
パラメータを右上に動かすと中央部分、左や左下へ動かすと両端部分のエラーが増加する。
区間を2つに分ける場合、次の式なら分母を 512 にできて、精度は 1024 を分母とした上記のものより むしろほんの少しだけ良くなる。ただし、加減算の回数が少し増える。
LANG:cpp int d = 5 * max < 12 * min ? 430 * max + 287 * min : 507 * max + 101 * min; return (d+ 256) >> 9; // 5 = 0b0101 加減算 2 回 // 12 = 0b1100 加減算 2 回 // 430 = 0x1ae = 0b000110101110 加減算 3 回 // 287 = 0x11f = 0b000100011111 加減算 2 回 // 507 = 0x1f8 = 0b000111111000 加減算 1 回 // 101 = 0x065 = 0b000001100101 加減算 3 回
区間を3分割すれば 0.5% の精度が出る。
LANG:cpp int d = 4 * max > 15 * min ? 510 * max + 67 * min : 4 * max > 7 * min ? 476 * max + 195 * min : 409 * max + 312 * min ; return (d+ 256) >> 9;
区分法による最適化の方針†
区間を分けた場合、ベストな結果を与えるのは区間の両端の値が等しい時なので、 ここから a と b との関係が決まる。
つまり、区間 での最良値を与える a, b の間には、
が成り立つので、ここから のような線形関係を導ける。
例えば に対してなら、
LANG:mathematica Solve[d2[0, a, b] == d2[0.1, a, b], b][[1]] // N Out[]= {b -> 0.158384 a}
この線形関係が表す直線は上の三角形状の谷グラフの最小値から左下へ伸びる谷筋に相当する。
この直線上で a を動かした際の誤差は、両端での誤差と中央での誤差の大きい方で決まって、
LANG:mathematica maxerrorinsegment[t1_, t2_, a_] := Solve[d2[t1, a, b] == d2[t2, a, b], b][[1]] // N // Max[ Abs[d2[ t1 , a, b /. #1] - 1], Abs[d2[(t1 + t2)/2, a, b /. #1] - 1] ] & Plot[maxerrorinsegment[0, 0.1, a], {a, 0.8, 1.2}]
のように最小点で折れ曲がるキレイな折れ線になるので、
LANG:mathematica optimizeinsegment[t1_, t2_] := FindMinimum[ maxerrorinsegment[t1, t2, a] // Hold, {a, 1} ] // Quiet optimizeinsegment[0, 0.1] Out[]= {0.00619396, {a -> 0.993806}}
のように最適化するのは簡単。
例えば [0, 0.1] の区間での最適値は a = 0.993806 で、そのときの最大誤差は 0.00619396 になる。
これを使って、区間を2つに区切る位置に対する左右区間の誤差の大きさを見てみると、
LANG:mathematica Plot[{ optimizeinsegment[0, t][[1]], optimizeinsegment[t, 1/4][[1]] }, {t, 0.12, 0.13}, GridLines -> Automatic ]
目の子で探した上記のパラメータがほぼ最適値になっていることを確認できる。
より広い範囲で見ると、
LANG:mathematica Plot[{ optimizeinsegment[0, t][[1]], optimizeinsegment[t, 1/4][[1]] }, {t, 0, 1/4} ]
のようになっていて、このグラフは区間をどれほど小さくすればどれだけの精度で近似できるかの指標を与える。
見やすくするために両対数プロットにすると、
LANG:mathematica LogPlot[optimizeinsegment[0, t][[1]], {t, 0.005, 1/4}, GridLines -> Automatic]
横軸が区間の幅、縦軸が誤差である。 これは左端をゼロに取るよう区間を取った際の誤差だけれど、 元の関数の対称性から、実は区間を取る位置には依らず、幅だけで決まる値となる(はず)。
横軸を N 分割の分割数 N にすると、
LANG:mathematica ListLogPlot[Table[{n, optimizeinsegment[0, 1/4/n][[1]]}, {n, 1, 20}], GridLines -> Automatic]
分割数の二乗に反比例して誤差が減ることが分かる。
すなわち精度を1桁よくするのに分割数を約3倍にする必要がある。
数値にして一部抜粋すると、
LANG:mathematica Table[{n, 100 optimizeinsegment[0, 1/4/n][[1]]}, {n, 1, 20}] // Grid
N | error / % |
1 | 3.95661 |
2 | 0.970056 |
3 | 0.429595 |
4 | 0.241345 |
8 | 0.0602639 |
16 | 0.0150615 |
分割数を倍にすると誤差が 1/4 になっていることを確かめられる。
解析的にも解けなくはない†
円弧を直線で近似しているだけなのだから、係数の最適値は閉じた式で書けるはず。
中心角を とすると、
直線の傾きは弧の両端を結んだ直線と平行であるはず。
切片 を持つ直線の方程式は
この直線と原点との距離は、
&math( \frac{\sin t}{\sqrt{(\sin t)^2+(1-\cos t)^2}x_0} );
なので、最大誤差を最小化するには、
となるように を取れば良い。このとき最大誤差は
であり、テーラー展開すると
LANG:mathematica Series[ (Sqrt[2 - 2 Cos[t]] - Sin[t])/ (Sqrt[2 - 2 Cos[t]] + Sin[t]), {t, 0, 8}]
となって、 の小さいところで分割数の二乗に反比例することが確かめられる。また、直線の式は
となる。三次元的な原点を通り、 において上記の直線を通る平面の方程式は、
となる。
例えば全体を1つの平面で近似する のとき、 の係数および最大誤差はそれぞれ、
LANG:mathematica With[{t = Pi/4}, {(Sqrt[2 - 2 Cos[t]] + Sin[t])/(2 Sqrt[2 - 2 Cos[t]]), (Sqrt[2 - 2 Cos[t]] + Sin[t])/(2 Sqrt[2 - 2 Cos[t]]) (1 - Cos[t])/Sin[t], (Sqrt[2 - 2 Cos[t]] - Sin[t])/(Sqrt[2 - 2 Cos[t]] + Sin[t]) }] // N
0.96194, 0.398448, 0.0395661
であり、係数に 128 を掛けると、
123.128, 51.0014
上で求めた値とほぼ一致する。少しずれているのは上での見積もりで誤差の評価に誤差があったためだ。
2つに分割すれば で、
0.990393, 0.197001, 0.00970056
係数に 1024 をかけると、
507.081, 100.865
で、これも上記の値と一致する。
2つ目の区間に対する係数はこれらを だけ回転させて、
LANG:mathematica With[{t = Pi/8}, {{Cos[t], -Sin[t]}, {Sin[t], Cos[t]}} . {507.08103178322705`, 100.86468848532132`} ]
429.883, 287.238
やはり上記の値と一致する。
同様にすると3つに分割した場合 、 誤差は 0.00429595 係数は分母を 512 として、
509.81, 67.1177
475.067, 196.779
407.949, 313.031
を得る。
2つ目、3つ目の区間で係数が上で求めた値と少しずつ違うのは、 上記の近似では区間を完全に3等分できていないためだ。
区間の分割時に近似が入ることを考慮しないと、 最適な係数を求めることはできない。
4分割の係数を求める†
LANG:mathematica In[30]:= With[{t = Pi/4/4}, {(Sqrt[2 - 2 Cos[t]] + Sin[t])/(2 Sqrt[2 - 2 Cos[t]]), (Sqrt[2 - 2 Cos[t]] + Sin[t])/(2 Sqrt[2 - 2 Cos[t]]) (1 - Cos[t])/ Sin[t], (Sqrt[2 - 2 Cos[t]] - Sin[t])/(Sqrt[2 - 2 Cos[t]] + Sin[t])}] // N Out[30]= {0.997592, 0.0982543, 0.00241345} In[32]:= 2048 {0.9975923633360981`, 0.09825427184336248`} Out[32]= {2043.07, 201.225} In[33]:= With[{t = Pi/4/4}, {{Cos[n t], -Sin[n t]}, {Sin[n t], Cos[n t]}}.%32] /. n -> {0, 1, 2, 3} // Grid Out[33]= { {"2043.069160112329`", "1964.555158055597`", "1810.544403010014`", "1586.9554419217554`"}, {"201.22474873520636`", "595.941291959864`", "967.7561455422308`", "1302.3806731753086`"} } In[45]:= Table[{n, n Cos[t]/Sin[t] // Abs[#1 - Round[#1]] &} /. t -> Pi {1, 2, 3}/4/4 // N, {n, 10, 2048}] // SortBy[#1, Function[e, e[[2]] // Total]] & // Take[#1, 10] & Out[45]= { {292., {0.0168683, 0.0496398, 0.0088827}}, {437., {0.0526419, 0.0113268, 0.0167183}}, {145., {0.0357736, 0.0609665, 0.00783559}}, {147., {0.0189053, 0.110606, 0.00104711}}, {729., {0.0695102, 0.038313, 0.025601}}, {584., {0.0337366, 0.0992796, 0.0177654}}, {1760., {0.117506, 0.0158698, 0.0261423}}, {874., {0.105284, 0.0226535, 0.0334366}}, {1941., {0.0659542, 0.0114754, 0.0882147}}, {548., {0.0179583, 0.0109678, 0.139958}} } In[46]:= n Cos[t]/Sin[t] /. {n -> 292, t -> Pi {1, 2, 3}/4/4} // N Out[46]= {1467.98, 704.95, 437.009} In[62]:= Plot[{1, Module[{x, y}, x = Cos[Pi t]; y = Sin[Pi t]; If[292 x < 437 y, 1587 x + 1302 y, If[292 x < 705 y, 1811 x + 968 y, If[292 x < 1468 y, 1965 x + 596 y, 2043 x + 201 y]]]/2048 ]}, {t, 0, 1/4}, GridLines -> Automatic, PlotRange -> {1 - 0.004, 1 + 0.004} ]
これは、次の係数で求めた結果になる。
LANG:c if (292 x < 705 y) { d = 292 x < 437 y ? 1587 x + 1302 y : 1811 x + 968 y ; } else { d = 292 x < 1468 y ? 1965 x + 596 y : 2043 x + 201 y ; } return (d + 1024) / 2048;
ちょっと係数を手で調節することで、
LANG:mathematica Plot[{1, Module[{x, y}, x = Cos[Pi t]; y = Sin[Pi t]; If[73 x < 110 y, 794 x + 651 y, If[73 x < 175 y, 905 x + 484 y, If[73 x < 366 y, 982 x + 299 y, 1022 x + 98 y]]]/1024 ]}, {t, 0, 1/4}, GridLines -> Automatic, PlotRange -> {1 - 0.004, 1 + 0.004} ]
多少係数を小さくできる。
LANG:c if (73 x < 175 y) { d = 73 x < 110 y ? 794 x + 651 y : 905 x + 484 y ; } else { d = 73 x < 366 y ? 982 x + 299 y : 1022 x + 98 y ; } return (d + 512) / 1024;
それにしても†
LANG:c int distance(int x, int y) { x = abs(x); y = abs(y); int max = x>y ? x:y; int min = x>y ? y:x; int d = 2*max<5*min ? 864*max+569*min : 1016*max+190*min; return (d+512)>>10; }
これで「原点からの距離」が 1% の精度で測れるとは驚きだ。
とはいえ†
本当に必要なのは「距離の逆数」の場合が多いんだよなあ。。。
3次元では†
LANG:c int distance3d(int x, int y, int z) { return distance2d(distance2d(x, y), z); }
でいいはず。