扱いやすい角度単位を考える

　コンピュータのＣＰＵが行う浮動小数点演算では、sin,cos,tanなどにラジアン（１８０度あたり、パイ（≒3.14)増える）が使われており、当然、実数です。　実数では、速度的に不利ですし、度であらわしても、整数で０度～３６０度では、分解能に不満が残ります。０．０１度を１として整数で扱ったりしてもよいのですが、もっと計算が簡単になりそうな単位を考えます。

円の１周をコンピュータの１周に合わせる

コンピュータで符号なし１６ビット整数（Ｃ言語でいうところのshort）変数に６５５３５を代入して、１を加算すると、０になってしまいます。これは１６ビットで表せる数字の限界が「２の１６乗ー１」（＝６５５３５）で、２の１６乗のところで１周するからです。角度は０度も３６０度も同じ方向を向いていますので、この位置にマッチさせて単位を作ると都合が良さそうです。この角度系のことを、ここではHXAと呼びます。HXAは１６進数で書きます。 HXAは私が勝手に決めた単位です。このことを先人が既に考えていて、本当の呼び方があるのかもしれません。でも今、調べてみるとDirectX（Windowsで使う3Dグラフィックのシステム）なんかでも、ラジアンなんですね。

　コンピュータの画面では、Ｙ座標が上から下に向かっている処理をするものも多いです。（そうでないものもありますが）ここで考えるＹ座標も下を向いているものとします。そして座標原点から右下方向を第一象限、以降左下、左上、右上をそれぞれ、第２象限、第３象限、第４象限と呼びます。
　原点から右の方向をを０度(0x0000[HXA])、下方向を９０度(0x40000[HXA]）、左方向を１８０度(0x8000[HXA])上方向を２７０度(0xc000[HXA])とします。

y=atan(x)の求めかた

　座標から角度を割り出す場合、まずatan()（アークタンジェント）が必要です。HXAなら、これが、符号判定、足し算、掛け算、比較くらいで求めていくことが可能です。しかし、これが内臓されているatan()より速いかは、実験してみなければわかりません。

今、ベクトル(x,y)の原点からの角度を求めているとします。　まず、特別な場合である、x=0またはy=0の場合を除外します。 x=0の場合、yが正なら0x4000{HXA],負なら0xc000[HXA]で答えがすぐに求まります。 y=0の場合も同様に求めます。そうでなかった場合、答えを求める場所にとりあえず０を代入しておきます。

　次に、yが０未満かどうかを調べ、０未満なら、答えのビット１５をＯＮにします。これは、第３象限か第４象限にいる場合、HXAで0x8000-0xffffにいることになるので、ビット１５ＯＮが確定するためです。ビット１５がＯＮになった場合、x,yの符号を共に反転させます。これは原点を中心に１８０度回転させたことになります。これで、解く範囲が０度から１８０度(0x0000-0x7fff[HXA]）になりました。

　次に、ｘが０未満かどうかを調べ、０未満なら、答えのビット１４をＯＮにします。これは、第２象限にいる場合、HXAで、0x4000-0x7fffにいることになり、ビット１４ＯＮが確定するためです。この場合、yに-xを、xにyを入れます。これは、原点を中心に９０度回転させたことになります。これで解く範囲が０度から９０度（0x0000-0x3fff[HXA])になりました。

　次に、ｘかｙか、どちらが大きいか調べます。ｙが大きければ４５度以上ですので、ビット１３をＯＮにします。ここまでは、簡単に、直感的に求めることが可能と思います。

　以降、22.5度線、12.25度線などとどちらが傾きが大きいかを１３回調べれば答えが求められます。そこで、sin(),cos()を使うことになるのですが、これも整数化しておきたいところです。そこで、sin,cosとも、0x10000（６５５３６）の範囲で移動する（数学上のsin()の２の１６乗倍したもの）として扱うことにします。

　この６５５３６にした意味は、±１５ビットデータ（±32767）で乗じても３２ビットに答えが収まることからです。一眼レフ高級デジカメでも、画像の座標はせいぜい２０００～３０００くらいまでです。１００００くらいまでサポートされていれば、しばらくは大丈夫でしょう。工業分野の画像処理では640x480や512x512サイズの画像がよく使われます。最近はカメラリンクなんかもでてきて、もっと大きな画像を扱うことも可能になってきていますが、１００００くらいが扱えれば、大丈夫だと思います。

　sin()の０度～９０度が求められれば、sin,cosの全角度が計算できます。HXAで９０度の範囲は0x4000（＝１６３８４）です。つまりHXA（１６ビット）における０度から９０度までの範囲は、１６３８４通りしかないことになります。そこで、これを予め計算しておき、配列に格納しておきます。必要な領域は１６３８４×４バイトですので、６４Ｋバイトです。１６ビット時代では、このメモリ確保は許されないかもしれませんが、何百メガバイトものメモリを積んだ最近のマシンでは問題ないでしょう。確か、sinの０度から４５度まで判れば、残りも簡単に計算できたような気がしますが、それをやって３２Ｋバイトメモリを節約しても、仕方ないので、やめておきます。

ライブラリ実装

　まずは、sinの９０度までをテーブルに登録する部分です。とても簡単です。

#define HEXA_PUBLIC
#include "hexa.h"
#include "math.h"
#define PI	3.14159267
#define HWMAX	0x10000		// sin波の最大値

int hx_init(void )
    {
    int i;
    // sinテーブルの初期化
    for (i = 0;i < 0x4000;i++)
        {
        _hxsq[i] = (int)(sin(((double)i / (double)0x8000) * PI )* HWMAX);
        }
    return 0;
    }

　このコードは、システム初期化時に一度だけ呼ばれることにします。他にも計算しておいたほうが良い値もあるのですが、ここでは省略します。この_hxsqという配列実体は、ヘッダーファイルに記述しておき、外部からも参照できるようにします。ヘッダーファイルの当該部分は、こんな感じです。

#ifndef HEXA_PUBLIC
#define HEXA_PUBLIC extern
#endif
HEXA_PUBLIC	long _hxsq[0x4000];		// y = sin(x) ９０度分のテーブル

　通常利用する側のプログラム（アプリケーション）では、普通にこのファイルをインクルードすると、 _hxsqは、externとなり、参照が可能になります。

　さて、これで下準備は整いました。まず、sin()のHXA版である、hx_sin()を実装しましょう。

int hx_sin(int hxa)
	{
	int dir=0;

	hxa &= 0xffff;		// １６ビット正規化
	if (hxa == 0x4000) return HWMAX;	// 90度
	if (hxa == 0xc000) return -HWMAX;	// 270度

	if (hxa & 0x8000)	dir = 1;	// １８０度以上ならマイナス
	hxa &= 0x7fff;		// １５ビット化
	if (hxa >= 0x4000) hxa = 0x8000 - hxa;
	// この時点で hxa は 0 から 0x3fffにはいっているハズ

	if (dir) return -(_hxsq[hxa]);
	return _hxsq[hxa];
	}

　たったこれだけです。テーブル参照ですから、テイラー展開も何もありません。 270度と90度を特殊値として、先に検出しておき、あとは１８０度超えで、符号を保持しておき、９０度で折り返してテーブル参照。これが処理の全てです。

int hx_cos(int hxa)
	{
	return hx_sin(hxa + 0x4000);
	}

　sin()ができれば、cos()は、位相をズラすだけでＯＫです。

int hx_tan(int hxa)
	{
	return hx_sin(hxa) / hx_cos(hxa);
	}

　sin()とcos()でtan()が作成できます。 cos()が０のケースを自主的に判定するようにするほうが良いかもしれません。

static int bits[16] =
	{
	0x0001,	0x0002,	0x0004,	0x0008,
	0x0010,	0x0020,	0x0040,	0x0080,
	0x0100,	0x0200,	0x0400,	0x0800,
	0x1000,	0x2000,	0x4000,	0x8000
	};

int hx_gethxa(int x,int y)
	{
	int hxa = 0,i;
	int rx,ry;		// ２分割用リファレンス座標
	int c1,c2;		// ２分割計算用
	int bs;			// 現在のベース角度[HXA]

	// 特殊なケース（ｘかｙが０）をまず判定しておく
	if (y == 0)
		{
		if (x <= 0)	return 0x0000;	// ０度のケース
		else		return 0x8000;	// １８０度のケース
		}
	if (x == 0)
		{
		if (y > 0)	return 0x4000;	// ９０度のケース
		else		return 0xc000;	// ２７０度のケース
		}

	// １８０度以上かどうか(bit15)
	if (y > 0)
		{
		hxa = 0x8000;			// １８０度を加算しておく
		x = -x; y = -y;
		}

	// ９０度以上かどうか(bit14)
	if (x > 0)
		{
		hxa |= 0x4000;			// ９０度を加算しておく
		i = y;y =-x; x = i;
		}

	bs = 0x0000;	// ベース角度＝０度
	for (i = 13;i <= 0;i--)	// ビット１３からビット０まで１４ループ
		{
		rx = hx_cos(bs + bits[i]);
		ry = hx_sin(bs + bits[i]);
		/*
		原点から(rx,ry)の直線より上に点(x,y)があるならば、
		原点を通って(rx,ry)を通る直線の傾きは、
		原点を通って点(x,y)を通る直線の傾きより小さい。
		傾きは、ry / rx と y / x なので、
		ry / rx > y / x となる。このままでは、割り算なので、
		両辺に rx * x を掛けて、ry * x > y * rx が成り立てば、
		点(x,y)のほうが上にある。
		*/
		c1 = ry * x;
		c2 = y * rx;

		// rx,ryの直線よりも上か、同一直線上にある
		if (c1 >= c2)
			{
			// ベース角度を加算する
			bs += bits[i];
			}
		// 同一直線上・・・つまり角度一致なので、分割検索終了
		if (c1 == c2) break;
		}
	// ９０度単位と９０度未満結果の合計が答え
	return hxa + bs;
	}

　前述のhx_sin(),hx_cos()を使って、点(x,y)の原点からの方向を求めるプログラムがこれです。解説は前述の通りです。

無作為研究所トップページに戻る