平方根を使わずに距離を求める

(133d) 更新


公開メモ

平面上の2点間の距離を平方根を使わずに計算する

こちらの記事で、

http://d.hatena.ne.jp/nowokay/20120604

2点間の距離を平方根を使わず絶対値+最大最小を使って近似的に求める手法を知りました。 浮動小数点演算を積んでいないマイコンやFPGAなどで便利に使えそうです。

そこで、この近似の精度について詳しく調べてみたのがこのページです。

最終的には、

  • どうやら上記のページで与えられているパラメータは最良ではないようです
  • 最適なパラメータはこちら #l9f15e88
  • さらに精度を上げることも可能 #r129c2f4
  • 本当に必要なのは「距離の逆数」の場合が多いんだよなあ。。。

といった感じでした。

リンク先の内容を Mathematica で確かめる

正確な距離関数はこう:

d(x,y)=\sqrt{x^2+y^2}

例えば点 (x_1,y_1) と点 (x_2,y_2) との間の距離を 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]

distance_wo_sqrt1.jpg

最大値、最小値を利用した近似関数のグラフはこう。

\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]

distance_wo_sqrt2.jpg

円錐を多角錐で近似していることになる。

プロットには、上記リンク先で紹介されていた係数 (a=1007/1024, b=441/1024) を使った。

近似値から真値を引いた誤差を表示すると:

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]

distance_wo_sqrt3.jpg

誤差は -2% から +8% くらいまでの間に分布する。

どうやらこのパラメータは近似値が真の距離を大きく下回らないように選ばれているようです? とも思ったのですが、 \theta=n\pi/2 近辺では少々下回っているという謎仕様。。。

リンク先のコメント欄で紹介されている max + (3/8) min という非常に単純化された処理でも 誤差は -2.5% から +7% の間に収まるので、ちょっとこれはいただけない。

実のところ、これらは最良のパラメータではない気がします。

パラメータの最適化を考える

エラーの相対誤差の二乗を最小化してみる。

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]

distance_wo_sqrt4.jpg

横軸が 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]

distance_wo_sqrt5.jpg

誤差は -6% から +2.5% の範囲に分布する。

誤差がゼロの上下に分散したため、上記のパラメータを使うより絶対値が小さくなっている。

ただし、「誤差の二乗の平均値」を最小にしたため、 \theta=n\pi/4 の周辺の狭い範囲では大きなエラーが出るままだ。

最適値の周りでの誤差の変化を見てみる。

LANG:Mathematica
ContourPlot[error[a/1024, b/1024], 
  {a, 964 - 10, 964 + 10}, {b, 420 - 10, 420 + 10}, ImageSize -> Large]

distance_wo_sqrt6.jpg

左下から右上への成分と、左上から右下への成分とに分けられそう。

まず最小値付近を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 の主成分は、

distance_wo_sqrt7.jpg

固有値 0.5 の副次成分は、

distance_wo_sqrt8.jpg

のような形をしている。

主成分のみでほぼ完全に距離関数を復元しているのに対して、
副次成分は \theta=n\pi/2 \theta=n\pi/2+\pi/4 とで異なる符号を持ち、 一方を上げると一方が下がるために平均的な誤差への影響が小さいことが見て取れる。

誤差の最大値を最小化する

誤差の評価を、

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]

distance_wo_sqrt9.jpg

拡大して、

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]

distance_wo_sqrt10.jpg distance_wo_sqrt10b.jpg

最小値を求める(始め 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]

distance_wo_sqrt11.jpg

(983/1024)\mathrm{max}+(407/1024)\mathrm{min} という近似で誤差は \pm 4 % の範囲に収まった。

※ 下で見るように、 (123/128)\mathrm{max}+(51/128)\mathrm{min} でもほぼ 4% に収まる

リンク先の値はどこから来たのか

上記の結果をまとめる。

ab
リンク先1007/1024441/1024
誤差の二乗平均を最小化964/1024420/1024
誤差の最大値を最小化983/1024407/1024

上のプロットにリンク先の点を赤で記入すると以下のようになる。

distance_wo_sqrt40.png  distance_wo_sqrt41.png

リンク先の値はどうやって決めたものなのか・・・ちょっとよく分かりません。

もうちょっとちゃんと考えてみる

もともと、どうして min と max で距離が求まるか、というと、、

d(x,y) &=\sqrt{x^2+y^2}\\ &=\sqrt{|x|^2+|y|^2}\\ &=\sqrt{|y|^2+|x|^2}\\

という関係から、

\begin{cases} x'=\mathrm{max}(|x|,|y|)\\ y'=\mathrm{min}(|x|,|y|)\\ \end{cases}

と置けば、

d(x,y)=d(x',y')

が成り立つところがキモなのだった。

この (x',y') は、 x'\ge 0,y'\ge 0 かつ、 x'\ge y' より、 第一象限の 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]
}]

distance_wo_sqrt32.svg

図中の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]

distance_wo_sqrt33.jpg distance_wo_sqrt34.jpg

正確な距離関数は本来円錐状なのだけれど、 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
]

distance_wo_sqrt35.jpg

本当はたわんでいる黄色のグラフをビシッと伸ばしてしまい青の三角形で近似する。

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}
]

distance_wo_sqrt36.svg

正しいグラフは r=1 の 1/8 円弧になるけれど、 近似グラフは直線になる。

昔から粗い近似として使われていたとリンク先のコメントで書かれていた \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}
]

distance_wo_sqrt37.svg

ずれてると言えばずれてるけど、まあそこそあってる。

こんな近似式をどこから導いたんだろう、と疑問だったのだけれど、 こうしてちゃんと考えればまずまず納得できる。(でもやっぱり賢いなあ)

近似精度の評価

近似の精度を評価することを考えると、

\begin{cases} x=r\cos\theta\\ y=r\sin\theta \end{cases}

と置けば、

d(x,y)=\sqrt{x^2+y^2}=|r|

\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]\\

となることから、相対的な近似精度を見るには r=1 に固定して \theta の関数として扱えば良い。

そこで、

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]

distance_wo_sqrt12.svg

となる。(横軸は \pi を単位とした方位角)

今度は正しいグラフが直線に、近似グラフが曲線になっている。 (さっきまでとグラフの色が変わっていてわかりにくいので気をつけて)

どうしてこんな関数になるかを理解するため、 \mathrm{max}(|\cos\theta|,|\sin\theta|) \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}]

distance_wo_sqrt13.svg

青が max、黄色が min。sin と cos の大きい方、小さい方がそれぞれ max, min になる。

これらの線形結合で最も 1 に近い関数を作ったのが上記の結果になる。

max も min も \theta=0\sim \pi/4 の範囲が反転しながら繰り返し現れるから、 誤差を最小化するには [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}]

distance_wo_sqrt14.svg

このプロットから、このパラメータに対する誤差の最大値は約 4% であることが分かる。

さらに誤差を小さくするには、例えば min の2倍の周期の三角波を

\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}]

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]
]

distance_wo_sqrt16.svg

とすると、最大誤差を 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
]

distance_wo_sqrt19.jpg distance_wo_sqrt18.jpg

左が正確な距離関数、右が近似関数。

相対誤差を評価すると、

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}
]

distance_wo_sqrt17.svg

横軸は \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/

以下のコードを試したところ、 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;
}

ビット幅などいろいろ工夫の余地は残る

リンク先のコメント欄にある、

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% 程度に収まる。

\mathrm{max}+(3/8)\mathrm{min}
distance_wo_sqrt20.svg


分母を 32 にすれば -5% ~ 6% の範囲にできる。

\mathrm{max}+(11/32)\mathrm{min}
distance_wo_sqrt21.svg


max にも係数を付ければ分母が 128 の段階で誤差をほぼ 4% 以内に収められる。

(123/128)\mathrm{max}+(51/128)\mathrm{min}
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}}

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]

distance_wo_sqrt26.gif

パラメータを右上に動かすと中央部分、左や左下へ動かすと両端部分のエラーが増加する。


区間を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 回

distance_wo_sqrt23.svg


区間を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;

distance_wo_sqrt24.svg

区分法による最適化の方針

区間を分けた場合、ベストな結果を与えるのは区間の両端の値が等しい時なので、 ここから a と b との関係が決まる。

つまり、区間 [t_1, t_2] での最良値を与える a, b の間には、

d^*(t_1;a,b)=d^*(t_2;a,b)

が成り立つので、ここから b = A a のような線形関係を導ける。

例えば [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}]

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 になる。

これを使って、区間を2つに区切る位置に対する左右区間の誤差の大きさを見てみると、

LANG:mathematica
Plot[{
    optimizeinsegment[0, t][[1]], 
    optimizeinsegment[t, 1/4][[1]]
  }, {t, 0.12, 0.13}, GridLines -> Automatic
]

distance_wo_sqrt28.svg

目の子で探した上記のパラメータがほぼ最適値になっていることを確認できる。

より広い範囲で見ると、

LANG:mathematica
Plot[{
    optimizeinsegment[0, t][[1]], 
    optimizeinsegment[t, 1/4][[1]]
  }, {t, 0, 1/4}
]

distance_wo_sqrt29.svg

のようになっていて、このグラフは区間をどれほど小さくすればどれだけの精度で近似できるかの指標を与える。

見やすくするために両対数プロットにすると、

LANG:mathematica
LogPlot[optimizeinsegment[0, t][[1]], {t, 0.005, 1/4}, 
  GridLines -> Automatic]

distance_wo_sqrt30.svg

横軸が区間の幅、縦軸が誤差である。 これは左端をゼロに取るよう区間を取った際の誤差だけれど、 元の関数の対称性から、実は区間を取る位置には依らず、幅だけで決まる値となる(はず)。

横軸を N 分割の分割数 N にすると、

LANG:mathematica
ListLogPlot[Table[{n, optimizeinsegment[0, 1/4/n][[1]]}, {n, 1, 20}], GridLines -> Automatic]

distance_wo_sqrt31.svg

分割数の二乗に反比例して誤差が減ることが分かる。

すなわち精度を1桁よくするのに分割数を約3倍にする必要がある。

数値にして一部抜粋すると、

LANG:mathematica
Table[{n, 100 optimizeinsegment[0, 1/4/n][[1]]}, {n, 1, 20}] // Grid
Nerror / %
13.95661
20.970056
30.429595
40.241345
80.0602639
160.0150615

分割数を倍にすると誤差が 1/4 になっていることを確かめられる。

解析的にも解けなくはない

円弧を直線で近似しているだけなのだから、係数の最適値は閉じた式で書けるはず。

distance_wo_sqrt36.svg

中心角を t とすると、

直線の傾きは弧の両端を結んだ直線と平行であるはず。

  \frac{-\sin t}{1-\cos t}

x 切片 x_0 を持つ直線の方程式は

  y=\frac{-\sin t}{1-\cos t}(x-x_0)

この直線と原点との距離は、

  \frac{\sin t}{\sqrt{(\sin t)^2+(1-\cos t)^2}x_0}

なので、最大誤差を最小化するには、

  x_0-1 = 1 - \frac{\sin t}{\sqrt{(\sin t)^2+(1-\cos t)^2}}x_0

となるように x_0 を取れば良い。このとき最大誤差は

  \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}]

\frac{t^2}{16}+\frac{t^4}{384}+\frac{17 t^6}{184320}+\frac{31 t^8}{10321920}+O\left(t^9\right)

となって、 t の小さいところで分割数の二乗に反比例することが確かめられる。また、直線の式は

  y=\frac{-\sin t}{1-\cos t}\left(x-\frac{2\sqrt{2-2\cos t}}{\sqrt{2-2\cos t}+\sin t}\right)

となる。三次元的な原点を通り、 z=1 において上記の直線を通る平面の方程式は、

  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)

となる。

例えば全体を1つの平面で近似する t=\pi/4 のとき、 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

上で求めた値とほぼ一致する。少しずれているのは上での見積もりで誤差の評価に誤差があったためだ。

2つに分割すれば t=\pi/8 で、

0.990393, 0.197001, 0.00970056

係数に 1024 をかけると、

507.081, 100.865

で、これも上記の値と一致する。

2つ目の区間に対する係数はこれらを \pi/8 だけ回転させて、

LANG:mathematica
With[{t = Pi/8},
 {{Cos[t], -Sin[t]},
  {Sin[t], Cos[t]}} . {507.08103178322705`, 100.86468848532132`}
]

429.883, 287.238

やはり上記の値と一致する。

同様にすると3つに分割した場合 t=\pi/12 、 誤差は 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}
]

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}
]

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;

それにしても

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);
}

でいいはず。

とはいえ

本当に必要なのは「距離の逆数」の場合が多いんだよなあ。。。

質問・コメント





添付ファイル: filedistance_wo_sqrt41.png 85件 [詳細] filedistance_wo_sqrt40.png 92件 [詳細] filedistance_wo_sqrt39.svg 164件 [詳細] filedistance_wo_sqrt38.svg 162件 [詳細] filedistance_wo_sqrt37.svg 216件 [詳細] filedistance_wo_sqrt36.svg 228件 [詳細] filedistance_wo_sqrt35.jpg 222件 [詳細] filedistance_wo_sqrt34.jpg 264件 [詳細] filedistance_wo_sqrt33.jpg 245件 [詳細] filedistance_wo_sqrt32.svg 217件 [詳細] filedistance_wo_sqrt31.svg 207件 [詳細] filedistance_wo_sqrt30.svg 182件 [詳細] filedistance_wo_sqrt29.svg 190件 [詳細] filedistance_wo_sqrt28.svg 178件 [詳細] filedistance_wo_sqrt27.svg 219件 [詳細] filedistance_wo_sqrt26.gif 227件 [詳細] filedistance_wo_sqrt25.png 227件 [詳細] filedistance_wo_sqrt24.svg 207件 [詳細] filedistance_wo_sqrt23.svg 205件 [詳細] filedistance_wo_sqrt10b.jpg 243件 [詳細] filedistance_wo_sqrt22.svg 189件 [詳細] filedistance_wo_sqrt21.svg 204件 [詳細] filedistance_wo_sqrt20.svg 201件 [詳細] filedistance_wo_sqrt19.jpg 228件 [詳細] filedistance_wo_sqrt18.jpg 256件 [詳細] filedistance_wo_sqrt17.svg 211件 [詳細] filedistance_wo_sqrt16.svg 267件 [詳細] filedistance_wo_sqrt15.svg 226件 [詳細] filedistance_wo_sqrt14.svg 231件 [詳細] filedistance_wo_sqrt13.svg 225件 [詳細] filedistance_wo_sqrt12.svg 240件 [詳細] filedistance_wo_sqrt11.jpg 247件 [詳細] filedistance_wo_sqrt10.jpg 234件 [詳細] filedistance_wo_sqrt9.jpg 229件 [詳細] filedistance_wo_sqrt8.jpg 239件 [詳細] filedistance_wo_sqrt7.jpg 233件 [詳細] filedistance_wo_sqrt6.jpg 209件 [詳細] filedistance_wo_sqrt5.jpg 252件 [詳細] filedistance_wo_sqrt4.jpg 248件 [詳細] filedistance_wo_sqrt3.jpg 236件 [詳細] filedistance_wo_sqrt2.jpg 225件 [詳細] filedistance_wo_sqrt1.jpg 236件 [詳細]

Counter: 4907 (from 2010/06/03), today: 10, yesterday: 0