平方根を使わずに距離を求める の変更点
更新- 追加された行はこの色です。
- 削除された行はこの色です。
- プログラミング/平方根を使わずに距離を求める へ行く。
- プログラミング/平方根を使わずに距離を求める の差分を削除
[[公開メモ]]
#contents
&katex();
* 平面上の２点間の距離を平方根を使わずに計算する [#v08053c6]
こちらの記事で、
http://d.hatena.ne.jp/nowokay/20120604
２点間の距離を平方根を使わず絶対値＋最大最小を使って近似的に求める手法を知りました。
浮動小数点演算を積んでいないマイコンやFPGAなどで便利に使えそうです。
そこで、この近似の精度について詳しく調べてみたのがこのページです。
最終的には、
- どうやら上記のページで与えられているパラメータは最良ではないようです
- 最適なパラメータはこちら [[#l9f15e88]]
- さらに精度を上げることも可能 [[#r129c2f4]]
- 本当に必要なのは「距離の逆数」の場合が多いんだよなあ。。。
といった感じでした。
* リンク先の内容を Mathematica で確かめる [#y6b0a391]
正確な距離関数はこう：
&math(d(x,y)=\sqrt{x^2+y^2});
例えば点 &math((x_1,y_1)); と点 &math((x_2,y_2)); との間の距離を
&math(d(x_2-x_1,y_2-y_1)); として求められる。
この関数をグラフにすると、頂点を下に向けた円錐が現れる。
LANG:Mathematica
distance[x_, y_] := Sqrt[x^2 + y^2]
Plot3D[distance[x, y], {x, -1, 1}, {y, -1, 1}, ImageSize -> Large]
&attachref(distance_wo_sqrt1.jpg,,50%);
最大値、最小値を利用した近似関数のグラフはこう。
&math(\tilde d(x,y;a,b)=a\,\mathrm{max}(|x|,|y|)+b\,\mathrm{min}(|x|,|y|));
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]
&attachref(distance_wo_sqrt2.jpg,,50%);
円錐を多角錐で近似していることになる。
プロットには、上記リンク先で紹介されていた係数 &math((a=1007/1024, b=441/1024)); を使った。
近似値から真値を引いた誤差を表示すると：
&math(e(x,y;a,b)=\tilde d(x,y;a,b)-d(x,y));
LANG:Mathematica
Plot3D[distance2[x, y, 1007/1024, 441/1024] - distance[x, y],
{x, -1, 1}, {y, -1, 1}, ImageSize -> Large, PlotPoints -> 80]
&attachref(distance_wo_sqrt3.jpg,,50%);
誤差は -2% から +8% くらいまでの間に分布する。
どうやらこのパラメータは近似値が真の距離を大きく下回らないように選ばれているようです？
とも思ったのですが、&math(\theta=n\pi/2); 近辺では少々下回っているという謎仕様。。。
リンク先のコメント欄で紹介されている max + (3/8) min という非常に単純化された処理でも
誤差は -2.5% から +7% の間に収まるので、ちょっとこれはいただけない。
実のところ、これらは最良のパラメータではない気がします。
* パラメータの最適化を考える [#p8c95460]
エラーの相対誤差の二乗を最小化してみる。
&math(e_1(a,b)=\int_{10^{-5}}^1 dx\,\int_{10^{-5}}^1 dy\,|e(x,y;a,b)|^2/d(x,y)^2);
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]
&attachref(distance_wo_sqrt4.jpg,,50%);
横軸が 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]
&attachref(distance_wo_sqrt5.jpg,,50%);
誤差は -6% から +2.5% の範囲に分布する。
誤差がゼロの上下に分散したため、上記のパラメータを使うより絶対値が小さくなっている。
ただし、「誤差の二乗の平均値」を最小にしたため、&math(\theta=n\pi/4); の周辺の狭い範囲では大きなエラーが出るままだ。
最適値の周りでの誤差の変化を見てみる。
LANG:Mathematica
ContourPlot[error[a/1024, b/1024],
{a, 964 - 10, 964 + 10}, {b, 420 - 10, 420 + 10}, ImageSize -> Large]
&attachref(distance_wo_sqrt6.jpg,,50%);
左下から右上への成分と、左上から右下への成分とに分けられそう。
まず最小値付近を２次関数でフィッティングして、
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 の主成分は、
&attachref(distance_wo_sqrt7.jpg,,50%);
固有値 0.5 の副次成分は、
&attachref(distance_wo_sqrt8.jpg,,50%);
のような形をしている。
主成分のみでほぼ完全に距離関数を復元しているのに対して、~
副次成分は &math(\theta=n\pi/2); と &math(\theta=n\pi/2+\pi/4); とで異なる符号を持ち、
一方を上げると一方が下がるために平均的な誤差への影響が小さいことが見て取れる。
* 誤差の最大値を最小化する [#jf7830a7]
誤差の評価を、
&math(e_2(a,b)=\mathrm{max}\Big(\ \big|\tilde d(x,y;a,b)-d(x,y)\big|\ \Big|\ \sqrt{x^2+y^2}=1\ \Big));
のようにして、誤差の最大値を最小化することを考える。
パラメータ 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]
&attachref(distance_wo_sqrt9.jpg,,50%);
拡大して、
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]
&attachref(distance_wo_sqrt10.jpg,,50%);
&attachref(distance_wo_sqrt10b.jpg,,50%);
最小値を求める（始め &math(a=1000, b=400); から探したところうまく見つからなかった）。
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]
&attachref(distance_wo_sqrt11.jpg,,50%);
&math((983/1024)\mathrm{max}+(407/1024)\mathrm{min});
という近似で誤差は &math(\pm 4);% の範囲に収まった。
※ 下で見るように、&math((123/128)\mathrm{max}+(51/128)\mathrm{min}); でもほぼ 4% に収まる
* リンク先の値はどこから来たのか [#l9f15e88]
上記の結果をまとめる。
| | a | b|コメント|
| リンク先 | 1007/1024| 441/1024|あんまり良くないような？|
|誤差の二乗平均を最小化|964/1024| 420/1024||
|誤差の最大値を最小化|983/1024| 407/1024|普通はこれが最良|
上のプロットにリンク先の点を赤で記入すると以下のようになる。
&attachref(distance_wo_sqrt40.png,,50%);
&attachref(distance_wo_sqrt41.png,,50%);
リンク先の値はどうやって決めたものなのか・・・ちょっとよく分かりません。
* もうちょっとちゃんと考えてみる [#jf9e0004]
もともと、どうして 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}
);
と置けば、
&math(d(x,y)=d(x',y'));
が成り立つところがキモなのだった。
この &math((x',y')); は、&math(x'\ge 0,y'\ge 0); かつ、&math(x'\ge y'); より、
第一象限の &math(y=x); 以下の部分に存在する点となる。
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]
}]
&attachref(distance_wo_sqrt32.svg);
図中の８個の点はすべて、網掛け領域内の (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]
&attachref(distance_wo_sqrt33.jpg,,66%);
&attachref(distance_wo_sqrt34.jpg,,66%);
正確な距離関数は本来円錐状なのだけれど、
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
]
&attachref(distance_wo_sqrt35.jpg,,66%);
本当はたわんでいる黄色のグラフをビシッと伸ばしてしまい青の三角形で近似する。
&math(z=1); との交線をプロットすると、
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}
]
&attachref(distance_wo_sqrt36.svg);
正しいグラフは &math(r=1); の 1/8 円弧になるけれど、
近似グラフは直線になる。
昔から粗い近似として使われていたとリンク先のコメントで書かれていた
&math(\mathrm{max}+(3/8)\mathrm{min}); であればこうなる。
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}
]
&attachref(distance_wo_sqrt37.svg);
ずれてると言えばずれてるけど、まあそこそあってる。
こんな近似式をどこから導いたんだろう、と疑問だったのだけれど、
こうしてちゃんと考えればまずまず納得できる。（でもやっぱり賢いなあ）
* 近似精度の評価 [#ec7a92df]
近似の精度を評価することを考えると、
&math(
\begin{cases}
x=r\cos\theta\\
y=r\sin\theta
\end{cases}
);
と置けば、
&math(d(x,y)=\sqrt{x^2+y^2}=|r|);
&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]\\
);
となることから、相対的な近似精度を見るには &math(r=1); に固定して &math(\theta); の関数として扱えば良い。
そこで、
&math(d^*(\theta;a,b)=a\,\mathrm{max}(|\cos\theta|,|\sin\theta|)+b\,\mathrm{min}(|\cos\theta|,|\sin\theta|));
を、上記で求めた「誤差の最大値を最小化するパラメータ」に対してプロットすると、
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]
&attachref(distance_wo_sqrt12.svg);
となる。（横軸は &math(\pi); を単位とした方位角）
今度は正しいグラフが直線に、近似グラフが曲線になっている。
（さっきまでとグラフの色が変わっていてわかりにくいので気をつけて）
どうしてこんな関数になるかを理解するため、
&math(\mathrm{max}(|\cos\theta|,|\sin\theta|)); と
&math(\mathrm{min}(|\cos\theta|,|\sin\theta|)); とを別々にプロットしてみると、
LANG:mathematica
Plot[{Max[Abs[Cos[t Pi]], Abs[Sin[t Pi]]], Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]]}, {t, 0, 2}]
&attachref(distance_wo_sqrt13.svg);
青が max、黄色が min。sin と cos の大きい方、小さい方がそれぞれ max, min になる。
これらの線形結合で最も 1 に近い関数を作ったのが上記の結果になる。
max も min も &math(\theta=0\sim \pi/4); の範囲が反転しながら繰り返し現れるから、
誤差を最小化するには &math([0,\pi/4]); の範囲だけ考えれば良い。
これは、上で 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}]
&attachref(distance_wo_sqrt14.svg);
このプロットから、このパラメータに対する誤差の最大値は約 4% であることが分かる。
さらに誤差を小さくするには、例えば min の２倍の周期の三角波を
&math(\big|\mathrm{min}(x,y)-(49/128) \tilde d(x,y;a,b)\big|);
のようにして作って、
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}]
&attachref(distance_wo_sqrt15.svg);
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]
]
&attachref(distance_wo_sqrt16.svg);
とすると、最大誤差を 1% 程度に抑えられる。
ただ、計算量が増えてしまうのでもう少し何とかならないか？
* 計算量を減らす [#qbba8f2b]
計算量を増加させずに精度を上げるには、
リンク先と同様に 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
]
&attachref(distance_wo_sqrt19.jpg,,50%);
&attachref(distance_wo_sqrt18.jpg,,50%,ogp);
左が正確な距離関数、右が近似関数。
相対誤差を評価すると、
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}
]
&attachref(distance_wo_sqrt17.svg);
横軸は &math(\pi); を単位とする方位角、縦軸は相対誤差です。
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; // 四捨五入
}
こちらで
https://code.sololearn.com/
以下のコードを試したところ、&math(1000\le x,y\le 10000); の範囲で
相対誤差の最大値は 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;
}
* ビット幅などいろいろ工夫の余地は残る [#r129c2f4]
リンク先のコメント欄にある、
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% 程度に収まる。
&math(\mathrm{max}+(3/8)\mathrm{min}); ~
&attachref(distance_wo_sqrt20.svg);
~
分母を 32 にすれば -5% ～ 6% の範囲にできる。
&math(\mathrm{max}+(11/32)\mathrm{min}); ~
&attachref(distance_wo_sqrt21.svg);
~
max にも係数を付ければ分母が 128 の段階で誤差をほぼ 4% 以内に収められる。
&math((123/128)\mathrm{max}+(51/128)\mathrm{min});~
&attachref(distance_wo_sqrt22.svg);
区間を分けない限り、これがほぼ最良なのは、この周辺での誤差分布をプロットすれば分かる。
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}}
&attachref(distance_wo_sqrt25.png);
最小値は (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]
&attachref(distance_wo_sqrt26.gif);
パラメータを右上に動かすと中央部分、左や左下へ動かすと両端部分のエラーが増加する。
~
区間を２つに分ける場合、次の式なら分母を 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 回
&attachref(distance_wo_sqrt23.svg);
~
区間を３分割すれば 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;
&attachref(distance_wo_sqrt24.svg);
* 区分法による最適化の方針 [#p4f2d6ff]
区間を分けた場合、ベストな結果を与えるのは区間の両端の値が等しい時なので、
ここから a と b との関係が決まる。
つまり、区間 &math([t_1, t_2]); での最良値を与える a, b の間には、
&math(d^*(t_1;a,b)=d^*(t_2;a,b));
が成り立つので、ここから &math(b = A a); のような線形関係を導ける。
例えば &math([0,0.1]); に対してなら、
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}]
&attachref(distance_wo_sqrt27.svg);
のように最小点で折れ曲がるキレイな折れ線になるので、
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 になる。
これを使って、区間を２つに区切る位置に対する左右区間の誤差の大きさを見てみると、
LANG:mathematica
Plot[{
optimizeinsegment[0, t][[1]],
optimizeinsegment[t, 1/4][[1]]
}, {t, 0.12, 0.13}, GridLines -> Automatic
]
&attachref(distance_wo_sqrt28.svg);
目の子で探した上記のパラメータがほぼ最適値になっていることを確認できる。
より広い範囲で見ると、
LANG:mathematica
Plot[{
optimizeinsegment[0, t][[1]],
optimizeinsegment[t, 1/4][[1]]
}, {t, 0, 1/4}
]
&attachref(distance_wo_sqrt29.svg);
のようになっていて、このグラフは区間をどれほど小さくすればどれだけの精度で近似できるかの指標を与える。
見やすくするために両対数プロットにすると、
LANG:mathematica
LogPlot[optimizeinsegment[0, t][[1]], {t, 0.005, 1/4},
GridLines -> Automatic]
&attachref(distance_wo_sqrt30.svg);
横軸が区間の幅、縦軸が誤差である。
これは左端をゼロに取るよう区間を取った際の誤差だけれど、
元の関数の対称性から、実は区間を取る位置には依らず、幅だけで決まる値となる（はず）。
横軸を N 分割の分割数 N にすると、
LANG:mathematica
ListLogPlot[Table[{n, optimizeinsegment[0, 1/4/n][[1]]}, {n, 1, 20}], GridLines -> Automatic]
&attachref(distance_wo_sqrt31.svg);
分割数の二乗に反比例して誤差が減ることが分かる。
すなわち精度を１桁よくするのに分割数を約３倍にする必要がある。
数値にして一部抜粋すると、
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 になっていることを確かめられる。
* 解析的にも解けなくはない [#f4332555]
円弧を直線で近似しているだけなのだから、係数の最適値は閉じた式で書けるはず。
&ref(distance_wo_sqrt36.svg);
中心角を &math(t); とすると、
直線の傾きは弧の両端を結んだ直線と平行であるはず。
&math(\frac{-\sin t}{1-\cos t});
&math(x); 切片 &math(x_0); を持つ直線の方程式は
&math(y=\frac{-\sin t}{1-\cos t}(x-x_0));
この直線と原点との距離は、
&math(
\frac{\sin t}{\sqrt{(\sin t)^2+(1-\cos t)^2}x_0}
);
なので、最大誤差を最小化するには、
&math(x_0-1 = 1 - \frac{\sin t}{\sqrt{(\sin t)^2+(1-\cos t)^2}}x_0);
となるように &math(x_0); を取れば良い。このとき最大誤差は
&math(\frac{\sqrt{2-2\cos t}-\sin t}{\sqrt{2-2\cos t}+\sin t});
であり、テーラー展開すると
LANG:mathematica
Series[
(Sqrt[2 - 2 Cos[t]] - Sin[t])/
(Sqrt[2 - 2 Cos[t]] + Sin[t]), {t, 0, 8}]
&math(\frac{t^2}{16}+\frac{t^4}{384}+\frac{17 t^6}{184320}+\frac{31 t^8}{10321920}+O\left(t^9\right));
となって、&math(t); の小さいところで分割数の二乗に反比例することが確かめられる。また、直線の式は
&math(y=\frac{-\sin t}{1-\cos t}\left(x-\frac{2\sqrt{2-2\cos t}}{\sqrt{2-2\cos t}+\sin t}\right));
となる。三次元的な原点を通り、&math(z=1); において上記の直線を通る平面の方程式は、
&math(z=\frac{\sqrt{2-2\cos t}+\sin t}{2\sqrt{2-2\cos t}}\left(x+\frac{1-\cos t}{\sin t}y\right));
となる。
例えば全体を１つの平面で近似する &math(t=\pi/4); のとき、&math(x,y); の係数および最大誤差はそれぞれ、
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
上で求めた値とほぼ一致する。少しずれているのは上での見積もりで誤差の評価に誤差があったためだ。
２つに分割すれば &math(t=\pi/8); で、
0.990393, 0.197001, 0.00970056
係数に 512 をかけると、
507.081, 100.865
で、これも上記の値と一致する。
２つ目の区間に対する係数はこれらを &math(\pi/8); だけ回転させて、
LANG:mathematica
With[{t = Pi/8},
{{Cos[t], -Sin[t]},
{Sin[t], Cos[t]}} . {507.08103178322705`, 100.86468848532132`}
]
429.883, 287.238
やはり上記の値と一致する。
同様にすると３つに分割した場合 &math(t=\pi/12);、
誤差は 0.00429595 係数は分母を 512 として、
509.81, 67.1177~
475.067, 196.779~
407.949, 313.031~
を得る。
２つ目、３つ目の区間で係数が上で求めた値と少しずつ違うのは、
上記の近似では区間を完全に３等分できていないためだ。
区間の分割時に近似が入ることを考慮しないと、
最適な係数を求めることはできない。
** ４分割の係数を求める [#t1bf0874]
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}
]
&attachref(distance_wo_sqrt38.svg);
これは、次の係数で求めた結果になる。
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}
]
&attachref(distance_wo_sqrt39.svg);
多少係数を小さくできる。
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;
* それにしても [#u4233dc8]
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% の精度で測れるとは驚きだ。
* ３次元では [#f64e777f]
LANG:c
int distance3d(int x, int y, int z)
{
return distance2d(distance2d(x, y), z);
}
でいいはず。
* とはいえ [#tc1de7de]
本当に必要なのは「距離の逆数」の場合が多いんだよなあ。。。
ニュートン法を使うやり方と、その初期値の求め方はこちらに説明されていた（英語）：~
https://en.wikipedia.org/wiki/Fast_inverse_square_root
* 質問・コメント [#uabcd372]
#article_kcaptcha
**係数の値について [#s8b5ea68]
>[[tkzn_io]] (&timetag(2023-11-05T15:12:32+09:00, 2023-11-06 (月) 00:12:32);)~
~
983/1024,407/1024という係数の値について、~
私も同様となりました。また、正4n角形での係数についても考えてみました。~
https://qiita.com/tkzn_io/items/ed84ebed79779560de00~
//
- コメントをありがとうございます。正$4n$角形に対応して [[#f4332555]] で $t=2\pi/4n=\pi/2n$ としてみたところ tkzn_io さんと同じく $\tan^2(\pi/8n)$ という誤差が得られました。１つ目の区間の係数は $\displaystyle z=\cos^2(\pi/8n)x+\frac12\big(\sin(\pi/4n)+\tan(\pi/4n)\big)y$ となるようでした。 -- [[武内（管理人）]] &new{2023-11-06 (月) 05:28:50};
- コメントをありがとうございます。正$4n$角形に対応して [[#f4332555]] で $t=2\pi/4n=\pi/2n$ としてみたところ tkzn_io さんと同じく $\tan^2(\pi/8n)$ という誤差を得ることができました。１つ目の区間の係数は $\displaystyle z=\cos^2(\pi/8n)x+\frac12\big(\sin(\pi/4n)+\tan(\pi/4n)\big)y$ となるようでした。 -- [[武内（管理人）]] &new{2023-11-06 (月) 05:28:50};
#comment_kcaptcha
Counter: 18509 (from 2010/06/03),
today: 7,
yesterday: 4