`circle.c` のバグ

『C言語による最新アルゴリズム事典』p. 65 の circle.c では小さな円がいびつになってしまいます。これについて読者のかたから非常に緻密な分析をご報告いただきましたので，ほぼそのまま掲載させていただきます。

なお，同じことがその次の楕円 ellipse.c にも当てはまります。

1. 現象

半径の小さな（半径5以下の）円を描くと，いびつな円になってしまう。

        ■■■■■■■
        ■□□□□□■
        ■□□□□□■
        ■□□□□□■
        ■□□□□□■
        ■□□□□□■
        ■■■■■■■
                半径３

        □■■■■■■■□
        ■□□□□□□□■
        ■□□□□□□□■
        ■□□□□□□□■
        ■□□□□□□□■
        ■□□□□□□□■
        ■□□□□□□□■
        ■□□□□□□□■
        □■■■■■■■□
                　半径４

        □□■■■■■■■□□
        □■□□□□□□□■□
        ■□□□□□□□□□■
        ■□□□□□□□□□■
        ■□□□□□□□□□■
        ■□□□□□□□□□■
        ■□□□□□□□□□■
        ■□□□□□□□□□■
        ■□□□□□□□□□■
        □■□□□□□□□■□
        □□■■■■■■■□□
                        半径５

2. 修正方法

関数 gr_circle() 内の次の if 文を修正する。

if ((r -= (y++ << 1) - 1) < 0)          /* 誤 */
if ((r -= (y++ << 1) + 1) <= 0)         /* 正 */

3. 理由

円の公式を示す。

r² = x² + y² (1)

このロジックは，円を書くために(1)式を変形させて使用する。

x² = r² - y² (2)

y に対応する x を(2)式から求め，整数値に丸めて点 (x,y) を書く。これを y=0 から y<=x の間 y++ しながらループすると 1/8 円が書ける。これを８方向に写像すれば円となる。

細かい説明は省くが，x座標が変化するのは，以下の式が成り立つときである。 y² = (2n + 1)r - n² - n - 0.25 n=0,1,2,3... (3)

このロジックでは，次にx座標が変化する y²の値(4)を求めておき，後述する仕掛けで r=(4)-y² となるようにしてある。このため，r<0 となったときに x（と次のx座標が変化する y²の値）を変更することで，円が書けるようになる。ここで，(3)式の最後に -0.25 という項がある。これは，r=0 となった時点で既にx座標を変化させなければならないことを示している。従って，ソース中の条件式は <= でなければならない。

次に，rを求める仕掛けを説明する。 y=0,1,2,3... のとき、y²=0,1,4,9...である。 y² のn番目(n=0,1,2...)の項を f(n) で表すと， f(n+1)-f(n)=2n+1 となる。 y の値は y++ により随時更新されているから， r=(4)-y² を求めるためには r-=(2y+1) を繰り返せばよい。従って，ソース中の条件式では，1引くのではなく1足さなければならない。

4. 修正結果

修正後の実行結果を示す。

        □□■■■□□
        □■□□□■□
        ■□□□□□■
        ■□□□□□■
        ■□□□□□■
        □■□□□■□
        □□■■■□□
                半径３

        □□□■■■□□□
        □■■□□□■■□
        □■□□□□□■□
        ■□□□□□□□■
        ■□□□□□□□■
        ■□□□□□□□■
        □■□□□□□■□
        □■■□□□■■□
        □□□■■■□□□
                　半径４

        □□□■■■■■□□□
        □□■□□□□□■□□
        □■□□□□□□□■□
        ■□□□□□□□□□■
        ■□□□□□□□□□■
        ■□□□□□□□□□■
        ■□□□□□□□□□■
        ■□□□□□□□□□■
        □■□□□□□□□■□
        □□■□□□□□■□□
        □□□■■■■■□□□
                        半径５

5. 追記

これも同じ読者の方からご指摘いただいたのですが，ソースコード中に

r += (x-- - 1) << 1;

とあるのは

r += --x << 1;

と書く方がスマートでした。

リンクはご自由にどうぞ。

松阪大学奥村晴彦 okumura@matsusaka-u.ac.jp