/*
 *	clipfix.h
 *
 *	固定小数ライブラリを定義します。
 *
 *	CLiP - Common Library for P/ECE
 *	Copyright (C) 2001-2009 Naoyuki Sawa
 *
 *	* Mon Mar 24 04:30:00 JST 2003 Naoyuki Sawa
 *	- 作成開始。
 *	* Thu Jul  8 06:00:00 JST 2004 Naoyuki Sawa
 *	- fcom()とficom()の引数順序を、(a,b,cond)から(f,cond,b)へ変更しました。
 *	  変更理由は、プログラムの読み易さのためです。
 *	* Mon Jul 12 06:00:00 JST 2004 Naoyuki Sawa
 *	- fild()の定義を((a)*FIXED_SCALE))から((a)<<FRACT_BITS))へ変更しました。
 *	  変更理由は、高速化のためです。
 *	* Thu Jul 22 12:30:00 JST 2004 Naoyuki Sawa
 *	- FRND()とFRND_RANGE()を追加しました。
 *	* Wed Jul 28 04:59:00 JST 2004 Naoyuki Sawa
 *	- FRND()の戻り値をfixedにキャストするのを忘れていたのを修正。
 *	* Wed Sep 15 06:36:00 JST 2004 Naoyuki Sawa
 *	- RND8/16,RND8/16_RANGEのマクロ→関数化に伴い、FRND,FRND_RANGEも関数化しました。
 *	  過去のソースの互換性のため、関数を呼び出すマクロも定義してあります。
 *	* Thu Dec 30 06:52:00 JST 2004 Naoyuki Sawa
 *	- DEG/RAD/GRADマクロを追加しました。
 *	* Sat Jul 08 04:00:13 JST 2006 Naoyuki Sawa
 *	- 逆減算「fsubr()」を追加しました。以下の二つの式は等価です。
 *		fr =  fsub(f1, f2);
 *		fr = fsubr(f2, f1);
 *	  引数の順序を逆にするだけで同じになるので、この関数が役に立つ場面はあまり無いと思います。
 *	- 整数からの逆減算「fisubr()」を追加しました。以下の二つの式は等価です。
 *		fr =   fsub(fild(i1), f2);
 *		fr = fisubr(     f2,  i1);
 *	  引数の順序を逆にすることの他、fild()マクロを一つ省略できるので、少し役に立つと思います。
 *	- 逆除算「fdivr()」を追加しました。以下の二つの式は等価です。
 *		fr =  fdiv(f1, f2);
 *		fr = fdivr(f2, f1);
 *	  引数の順序を逆にするだけで同じになるので、この関数が役に立つ場面はあまり無いと思います。
 *	- 整数からの逆除算「fidivr()」を追加しました。以下の二つの式は等価です。
 *		fr =   fdiv(fild(i1), f2);
 *		fr = fidivr(     f2,  i1);
 *	  引数の順序を逆にすることの他、fild()マクロを一つ省略できるので、少し役に立つと思います。
 *	- ※TODO:ベクトル演算と行列演算にも、fsubr(),fdivr()等と同様の逆演算形式を追加しましょう。
 *	- 固定小数の自己代入演算マクロを定義しました。
 *	  ベクトルと行列の自己代入演算マクロは、今後、必要に応じて追加する予定です。→ 2006/09/14 追加しました。
 *	  ※TODO:自己代入形式「vxform_()」を作成する場合、右辺項へ代入すべきかどうか、要検討です。
 *	- ※TODO:80x87の剰余命令は「FPREM」です。本モジュールにも「fprem()」関数を追加しましょう。 → 2006/07/11 追加しました。
 *	* Tue Jul 11 07:14:35 JST 2006 Naoyuki Sawa
 *	- 剰余演算「fprem()」を追加しました。
 *	  <math.h>のfmod()関数に相当する機能です。
 *	  fprem()の名称は、80x87の剰余命令「FPREM」に由来します。
 *	- 今回は、整数版(fiprem())や逆演算形式(fpremr())は、作成しませんでした。
 *	  作成しないことにした理由は、剰余演算において、それらの形式を利用することは無いだろうと考えたからです。
 *	  また、80x87にも「FPREM」命令の整数版や逆演算形式は存在しないからです。
 *	- FPREMやfmod()の結果の符号は、被序数と同じ符号になります。これは、仕様です。
 *	  一方、C言語の整数「%」演算子の結果の符合は、処理系依存となります。
 *	  GCCやVisual C++の整数「%」演算子は、たまたま、FPREMやfmod()の仕様と同じ動作となっています。
 *	  従って、今回のfprem()の実装には、単純にC言語の整数「%」演算子を利用しました。
 *	  整数「%」演算子の結果が、被序数と同じ符号にならない処理系が存在するのかどうか不明ですが、念のため留意すべきです。
 *	* Thu Sep 14 13:36:26 JST 2006 Naoyuki Sawa
 *	- 浮動小数点数の定数を定義しました。(M_PIなど)
 *	  Borland C++ Compiler 5.5 の <math.h> から拝借したものです。
 *	  Visual C++ 6.0 の <math.h> には定義されていません。
 *	  「fixed d = fld(M_PI * 123.4 / 180.0); // 123.4°を表すfixed値」のように利用します。
 *	- ベクトルと行列の自己代入演算マクロを追加しました。
 *	* Sat Dec 30 17:42:52 JST 2006 Naoyuki Sawa
 *	- M_PI2の定義を追加しました。
 *	  Borland C++ Compiler 5.5はなぜか、２πに相当するシンボルを定義していないのですが、
 *	  使用頻度が高く、シンボル名を混同する危険も少ないので、定義しておくことにしました。
 *	* Sat Feb 10 22:56:55 JST 2007 Naoyuki Sawa
 *	- fceil(),ffloor(),ficeil(),fifloor(),fceil_(),ffloor_()を追加しました。
 *	* Tue Apr 03 00:33:39 JST 2007 Naoyuki Sawa
 *	- POSITIVE_MIN,POSITIVE_MAX,NEGATIVE_MIN,NEGATIVE_MAXを追加しました。
 *	* Tue Apr 03 18:29:20 JST 2007 Naoyuki Sawa
 *	- fcom()マクロとficom()マクロの定義を、より安全なものに変更しました。
 *	* Tue Jan 29 15:21:49 JST 2008 Naoyuki Sawa
 *	- fprem_()を追加しました。
 *	* Thu Apr 02 13:37:51 JST 2009 Naoyuki Sawa
 *	- clipflt1.cモジュールに、xtof()とftox()を追加しました。(clipfixモジュールの変更は有りません)
 *	  xtof(),ftox()関数は、fst(),fld()マクロと同じ働きですが、fst(),fld()よりも効率的な実装としました。
 *	  マクロ展開が必要ない場合は、fst(),fld()よりも、xtof(),ftox()を利用する方が、高速化を期待できます。
 *	  詳細は、clipflt1.cモジュールのコメントを参照してください。
 */
#ifndef __CLIP_FIX_H__
#define __CLIP_FIX_H__

/* FIX_NOASMを定義するとCルーチン、定義しなければasmルーチンを使います。 */
#undef FIX_NOASM	/* まず確実に未定義にした上で… */
//#define FIX_NOASM	/* 必要に応じて定義してください（定義すると遅くなります） */

/*****************************************************************************
 *	定数
 *****************************************************************************/

/* Borland C++ Compiler 5.5 <math.h> より */
#define M_E		2.71828182845904523536		// ｅ
#define M_LOG2E		1.44269504088896340736		// log２(ｅ)
#define M_LOG10E	0.434294481903251827651		// log10(ｅ)
#define M_LN2		0.693147180559945309417		// logｅ(２)
#define M_LN10		2.30258509299404568402		// logｅ(10)
#define M_PI		3.14159265358979323846		// π
#define M_PI_2		1.57079632679489661923		// π／２ 
#define M_PI_4		0.785398163397448309616		// π／４
#define M_1_PI		0.318309886183790671538		// １／π
#define M_2_PI		0.636619772367581343076		// ２／π
#define M_1_SQRTPI	0.564189583547756286948		// １／√(π)
#define M_2_SQRTPI	1.12837916709551257390		// ２／√(π)
#define M_SQRT2		1.41421356237309504880		//     √(２)
#define M_SQRT_2	0.707106781186547524401		// １／√(２) = √(１／２)

/* * Sat Dec 30 17:42:52 JST 2006 Naoyuki Sawa
 * - M_PI2の定義を追加しました。
 *   Borland C++ Compiler 5.5はなぜか、２πに相当するシンボルを定義していないのですが、
 *   使用頻度が高く、シンボル名を混同する危険も少ないので、定義しておくことにしました。
 */
#define M_PI2		6.28318530717958647692 		// ２π

/*****************************************************************************
 *	fixed
 *****************************************************************************/

#define FRACT_BITS	10			/* 小数部ビット数 */
#define VALUE_BITS	(32 - FRACT_BITS)	/* 整数部ビット数 */
#define FIXED_SCALE	(1 << FRACT_BITS)	/* 固定小数スケール値 */
#define FRACT_MASK	((1 << FRACT_BITS) - 1)	/* 小数部ビットマスク */

/* 重要: (fchs(POSITIVE_MAX) != NEGATIVE_MAX) です。注意してください!! */
#define POSITIVE_MIN	((fixed)+1)		/* fixed型の正の最小値 */
#define POSITIVE_MAX	((fixed)INT_MAX)	/* fixed型の正の最大値 */
#define NEGATIVE_MIN	((fixed)-1)		/* fixed型の負の最小値 */
#define NEGATIVE_MAX	((fixed)INT_MIN)	/* fixed型の負の最大値 */

/* 固定小数データ型の定義です。
 * * struct _fixed は未定義です。すなわち、fixedの実態は未定義構造体へのポインタです。
 * - fixedをこのように定義した理由は、誤って通常の演算子を適用できないようにするためです。
 *   もしも
 *	(1) typedef int fixed;
 *	(2) typedef void* fixed;
 *   のように定義した場合、次のような誤った演算が通ってしまいます。((1),(2)いずれも)
 *	fixed f;
 *	int n;
 *	...
 *	f += n; // 本当は f = fadd(f, fild(n)) とすべき
 * - 通常の演算子適用を避けるための別の方法として、fixedを構造体として定義する方法もあります。
 *	(3) typedef struct _fixed {
 *		unsigned fract:FRACT_BITS;
 *		int value:VALUE_BITS;
 *	    } fixed;
 *   (3)の方法を採らなかった理由は、実行効率の低下を避けるためです。
 *   gcc33では、構造体引数はレジスタ渡しではなくスタック渡しとなります。
 *   たとえ構造体サイズが32ビット以内であっても、無条件にスタック渡しとなります。
 *   四則演算でスタック渡しが発生すると実行効率が大きく低下するので、この方法は採りませんでした。
 *   また、ポインタに較べて(int)へのキャストが難しく、加減算をマクロ定義するのが難しくなる点も、
 *   (3)の方法を採らなかった理由のひとつです。
 */
typedef struct _fixed* fixed;

/* 浮動小数or整数→固定小数（グローバル変数の初期化に使うため、インライン関数ではなくマクロとして定義） */
#define  fld(a)		((fixed)(int)((a) * FIXED_SCALE))
/*           整数→固定小数（グローバル変数の初期化に使うため、インライン関数ではなくマクロとして定義） */
#define fild(a)		((fixed)     ((a) << FRACT_BITS))
/* 固定小数→浮動小数 */
static inline float fst(fixed a) { return (float)(int)a / FIXED_SCALE; }
/* 固定小数→整数 */
static inline int fist(fixed a) { return (int)a >> FRACT_BITS; }
/* 符号反転 */
static inline fixed fchs(fixed a) { return (fixed)-(int)a; }
/* 絶対値（fabs()は標準Cライブラリの同名関数とぶつかるので回避） */
static inline fixed ffabs(fixed a) { return (int)a >= 0 ? a : fchs(a); }
/* 加算 */
static inline fixed  fadd(fixed a, fixed b) { return (fixed)((int)a + (int)     b ); }
static inline fixed fiadd(fixed a, int   b) { return (fixed)((int)a + (int)fild(b)); }
/* 減算 */
static inline fixed   fsub(fixed a, fixed b) { return (fixed)((int)a - (int)     b ); }
static inline fixed  fisub(fixed a, int   b) { return (fixed)((int)a - (int)fild(b)); }
static inline fixed  fsubr(fixed a, fixed b) { return (fixed)((int)     b -  (int)a); } /* 2006/07/08 追加 */
static inline fixed fisubr(fixed a, int   b) { return (fixed)((int)fild(b) - (int)a); } /* 2006/07/08 追加 */
/* 比較 */
//static inline int fcom (fixed a, fixed b) { return (int)fsub (a, b); }
//static inline int ficom(fixed a, int   b) { return (int)fisub(a, b); }
//↑これではダメ！aとbの差が大きいときにオーバーフローしてしまい、正しく判定できません。
//#define  fcom(a, cond, b)	((int)(a) cond (int)    (b))
//#define ficom(a, cond, b)	((int)(a) cond (int)fild(b))
//↓* Tue Apr 03 18:29:20 JST 2007 Naoyuki Sawa
//↓- fcom()マクロとficom()マクロの定義を、より安全なものに変更しました。
//↓- 旧い定義では、fcom()マクロやficom()マクロの引数型を間違えても警告が表示されず、誤りを見逃す恐れがありました。
//↓  <例>
//↓	fixed a = ...;
//↓	int b = ...;
//↓	if(fcom(a, >, b)) { ... } ⇒ if((int)a > (int)b) { ... } と展開されてしまい、警告が表示されない。
//↓- 新しい定義では、fcom()マクロやficom()マクロの引数型を間違えると、警告が表示され、誤りを見つけ易くなりました。
//↓  <例>
//↓	fixed a = ...;
//↓	int b = ...;
//↓	if(fcom(a, >, b)) { ... } ⇒ 展開されたコードの、「fixed __b__ = (b);」の部分で警告が表示される。
//↓- VisualC++ではGCC拡張機能が使えないので、旧い定義のままです。
#ifdef PIECE
//#define  fcom(a, cond, b)	({ fixed __a__ = (a); fixed __b__ =     (b); (int)__a__ cond (int)__b__; })
//#define ficom(a, cond, b)	({ fixed __a__ = (a); fixed __b__ = fild(b); (int)__a__ cond (int)__b__; })
//↓* Tue Apr 01 12:10:19 JST 2008 Naoyuki Sawa
//↓- 上のように、GCC拡張構文の中で比較を行うと、冗長なコードが生成されてしまうことに気付きました。
//↓  比較してすぐに分岐するのでなく、比較結果を0/1として一時変数に格納し、その値を見て分岐してしまいます。
//↓  そこでtypecheck()を利用して、GCC拡張構文の外で比較し、かつ、型検査も行う方法に書き換えました。
//↓- ficom()のbを「typecheck(int,b)」で型検査してしまうと、shortやlongとの比較ができなくなってしまいます。
//↓  明示的に型検査しなくても、fild(b)のシフト演算によって型検査できるので、bの型検査は行いません。
#define  fcom(a, cond, b)	(typecheck(fixed, a) && typecheck(fixed, b) && ((int)(a) cond (int)    (b)))
#define ficom(a, cond, b)	(typecheck(fixed, a)                        && ((int)(a) cond (int)fild(b)))
#else /*PIECE*/
#define  fcom(a, cond, b)	((int)(a) cond (int)    (b))
#define ficom(a, cond, b)	((int)(a) cond (int)fild(b))
#endif /*PIECE*/
/* 乗算 */
              fixed  fmul(fixed a, fixed b);
static inline fixed fimul(fixed a, int   b) { return (fixed)((int)a * b); }
/* 除算 */
              fixed   fdiv(fixed a, fixed b);
static inline fixed  fidiv(fixed a, int   b) { return (fixed)((int)a / b); }
static inline fixed  fdivr(fixed a, fixed b) { return fdiv(     b,  a); } /* 2006/07/08 追加 */
static inline fixed fidivr(fixed a, int   b) { return fdiv(fild(b), a); } /* 2006/07/08 追加 */
/* 剰余 */
static inline fixed  fprem(fixed a, fixed b) { return (fixed)((int)a % (int)b); } /* 2006/07/11 追加 */
/* 平方根（負数を指定した場合の結果は未定義です） */
fixed fsqrt(fixed a);
/* 三角関数 */
fixed fsin(fixed a);		//80387以降
fixed fcos(fixed a);		//80387以降
fixed ftan(fixed a);		//80x87ではfptan
fixed fasin(fixed a);		//80x87にはない
fixed facos(fixed a);		//80x87にはない
fixed fatan(fixed a);		//80x87ではfpatan
fixed fatan2(fixed y, fixed x);	//80x87にはない
/* 切り上げ、切り捨て */
static inline fixed  fceil(fixed a) { return (fixed)(((int)a + FRACT_MASK) & ~FRACT_MASK); } /* 2007/02/10追加。80x87には有りません。BREW 2.1に同名のAPIが有るみたいです。 */
static inline fixed ffloor(fixed a) { return (fixed)(((int)a             ) & ~FRACT_MASK); } /* 2007/02/10追加。80x87には有りません。BREW 2.1に同名のAPIが有るみたいです。 */
static inline int  ficeil(fixed a) { return ((int)a + FRACT_MASK) >> FRACT_BITS; } /* 2007/02/10追加 */
static inline int fifloor(fixed a) { return ((int)a             ) >> FRACT_BITS; } /* 2007/02/10追加。ficeil()との対称性のために用意しました。処理内容はfist()と同じです。 */

/* 円周率
 * * 整数倍した時の誤差を防ぐため、π/2を基準に2倍,4倍して定義します。
 * * 三角関数テーブルへのオフセットとしても利用するため、int型とfloat型の両方を定義しました。
 */
#if FIXED_SCALE != 1024
#error 固定小数パラメータが想定値と異なります。
#endif
#define _PI_2	1608	/* 3.14159265358979323846 / 2 * 1024 = 1608.49543863797413809152 */
#define _PI	3216	/* 3.14159265358979323846     * 1024 = 3216.99087727594827618304 */
#define _PI2	6432	/* 3.14159265358979323846 * 2 * 1024 = 6433.98175455189655236608 (本当は6433ですが、_PI_2の4倍に合わせて6432としました) */
#define PI_2	((fixed)_PI_2)
#define PI	((fixed)_PI  )
#define PI2	((fixed)_PI2 )

/* 整数または浮動小数の角度表記を、固定小数のラジアン角に変換します。
 * [in]
 *	a		変換元の角度。
 *			{整数定数,整数変数,浮動小数定数,浮動小数変数}にて指定可能です。
 *			{整数定数,浮動小数定数}を指定した場合は、コンパイル時に計算されます。
 *			{整数変数,浮動小数変数}を指定した場合は、実行時に計算されます。
 * [out]
 *	戻り値		変換後の角度。
 *			固定小数のラジアン角です。
 * [note]
 *	* <例>
 *	  fixed a = DEG(45);   //  45度
 *	  fixed b = RAD(3.14); // 180度
 *	  fixed c = GRAD(200); // 360度
 */
#define  RAD(a)		(fld(a))					/* ラジアン(直角= 1π) */
//#define  DEG(a)	(fidiv(fmul(fld(a), PI_2),  90))		/* デグリー(直角=90度) */
//#define GDEG(a)	(fidiv(fmul(fld(a), PI_2), 100))		/* グラード(直角=100G) */
//↑fmul()は関数呼び出しなので、引数が定数でも実行時計算になってしまいます。
//↓引数が定数の場合は実行時計算を行いたくないので、下のように変形しました。
#define  DEG(a)		((fixed)((int)fld(a) * _PI_2 / (int)fld( 90)))	/* デグリー(直角=90度) */
#define GRAD(a)		((fixed)((int)fld(a) * _PI_2 / (int)fld(100)))	/* グラード(直角=100G) */

///* 0.0以上1.0★未満★の乱数を生成します。 */
//#define FRND(s)		((fixed)(RND16(s) & FRACT_MASK))
///* (a)以上(b)★未満★の乱数を生成します。 */
//#define FRND_RANGE(s, a, b)	fadd(fmul(FRND(s), fsub((b), (a))), (a))
//↓2004/09/15 RND8/16,RND8/16_RANGEのマクロ→関数化に伴い、こちらも関数化しました。
//↓           過去のソースの互換性のため、関数を呼び出すマクロも定義してあります。
/* 0.0以上1.0★未満★の乱数を生成します。 */
fixed frnd(int* seed);
#define FRND(s)			(frnd(&(s)))
/* (a)以上(b)★未満★の乱数を生成します。 */
fixed frnd_range(int* seed, fixed a, fixed b);
#define FRND_RANGE(s, a, b)	(frnd_range(&(s), (a), (b)))

/*****************************************************************************
 *	vector
 *****************************************************************************/

/* ベクトル構造体の定義です。 */
typedef struct _vector {
	fixed x;	/* X成分 */
	fixed y;	/* Y成分 */
	fixed z;	/* Z成分 */
} vector;

/* 定義済みのベクトルです。 */
extern const vector VECTOR_0;	/* ゼロベクトル */
extern const vector VECTOR_X;	/* X軸単位ベクトル */
extern const vector VECTOR_Y;	/* Y軸単位ベクトル */
extern const vector VECTOR_Z;	/* Z軸単位ベクトル */

/* ベクトルの初期化関数は用意していません。
 * GCC拡張機能を利用して、次のように書けます。
 *	vector v;
 *	...
 *	v = (vector){x_value,y_value,z_value};
 *	...
 */

/* ベクトルの各成分を符号反転します。
 * [in]
 *	a		符号反転するベクトル。
 * [out]
 *	戻り値		結果のベクトル。
 */
vector vchs(vector a);

/* ベクトルの各成分を加算します。
 * [in]
 *	a,b		加算するベクトル。
 * [out]
 *	戻り値		結果のベクトル。
 */
vector vadd(vector a, vector b);

/* ベクトルの各成分を減算します。
 * [in]
 *	a,b		減算するベクトル。
 * [out]
 *	戻り値		結果のベクトル。
 */
vector vsub(vector a, vector b);

/* ベクトルの各成分を乗算します。
 * [in]
 *	a,b		乗算するベクトル。
 * [out]
 *	戻り値		結果のベクトル。
 */
vector vmul(vector a, vector b);

/* ベクトルの各成分を同じ値で乗算します。
 * [in]
 *	a		乗算するベクトル。
 *	s		乗算する値。
 * [out]
 *	戻り値		結果のベクトル。
 */
vector vfmul(vector a, fixed b); // 固定小数
vector vimul(vector a, int   b); // 整数

/* ベクトルの各成分を除算します。
 * [in]
 *	a,b		除算するベクトル。
 * [out]
 *	戻り値		結果のベクトル。
 */
vector vdiv(vector a, vector b);

/* ベクトルの各成分を同じ値で除算します。
 * [in]
 *	a		除算するベクトル。
 *	s		除算する値。
 * [out]
 *	戻り値		結果のベクトル。
 */
vector vfdiv(vector a, fixed b); // 固定小数
vector vidiv(vector a, int   b); // 整数

/* ベクトルの内積(dot product)を求めます。
 * [in]
 *	a,b		内積を求めるベクトル。
 * [out]
 *	戻り値		内積の値。
 * [note]
 * 	* 内積の定義: a・b = |a||b|cosθ
 */
fixed vdot(vector a, vector b);

/* ベクトルの外積(cross product)を求めます。
 * [in]
 *	a,b		外積を求めるベクトル。
 * [out]
 *	戻り値		結果のベクトル。
 * [note]
 * 	* 外積の定義: a×b = |a||b|sinθ*c
 * 	  ただし、cはaとbに垂直な単位ベクトルとする。
 */
vector vcross(vector a, vector b);

/* ベクトルを正規化(normalize)します。
 * [in]
 *	a		正規化するベクトル。
 * [out]
 *	戻り値		結果のベクトル。
 * [note]
 *	* Tue Dec 21 06:00:00 JST 2004 Naoyuki Sawa
 *	  ゼロベクトルが指定された場合は、そのままゼロベクトルを返すようにしました。
 *	  (これまでは、ゼロベクトルが指定されると、ゼロ除算エラーになっていました。)
 */
vector vnormal(vector a);

/* ベクトルの長さ(magnitude)を求めます。
 * [in]
 *	a		長さを求めるベクトル。
 * [out]
 * 	戻り値		結果の値。
 */
fixed vmag(vector a);

/* ベクトルの長さの二乗(square magnitude)を求めます。
 * [in]
 *	a		長さの二乗を求めるベクトル。
 * [out]
 * 	戻り値		結果の値。
 */
fixed vmag2(vector a);

/* 反射ベクトルを求めます。
 * [in]
 *	N		反射面の法線ベクトル。
 *	E		反射点から光源へのベクトル。
 *			言い換えると、逆向きの光源ベクトルです。
 * [out]
 *	戻り値		結果のベクトル。
 * [note]
 *	* 反射ベクトルの公式は、次の通りです。
 *		R = 2(E・N)N - E
 *	  考え方:
 *	  ベクトルE,Nが同じ原点を持つとして、EからNへ垂線を下ろした場合、
 *	  その垂線がNと交わるP点は
 *		P = (E・N)N
 *	  です。なぜなら、(E・N)はベクトルE,Nのなす角のcos値だからです。
 *	  Rは、Pを基準にEに点対称な位置、すなわち、
 *	  Eから、E->Pの二倍進んだ位置です。
 *		R = E + 2(P - E)
 *		  = 2P - E
 *		  = 2(E・N)N - E
 *	  以上です。
 *	* ベクトルE,Nは、正規化されていなければいけません。
 *	  正規化されていない場合、間違った結果となります。
 *	  ※実際には、Eだけが正規化されていない場合は、
 *	    単に結果のベクトルの大きさがEに等しくなるだけのように思います。
 *	    (未確認)
 *	* ベクトルEは、"逆向きの"光源ベクトルであることに注意してください。
 *	  この関数は、反射ベクトルを求めると言うよりも、
 *	  Nを基準として、Eに線対称なベクトルを求める、と言った方が正確です。
 *	  Eを光源ベクトルと同じ向きに取った方が使いやすいと思うのですが、
 *	  公式として逆向きの方が一般的なようなので、ならってみました。
 */
vector vreflect(vector N, vector E);

/*****************************************************************************
 *	matrix
 *****************************************************************************/

/* 行列構造体の定義です。 */
typedef struct _matrix {
	/* | a00 a10 a20 |T = 姿勢ベクトルのX成分 */
	/* | a01 a11 a21 |T = 姿勢ベクトルのY成分 */
	/* | a02 a12 a22 |T = 姿勢ベクトルのZ成分 */
	/* | a03 a13 a23 |T = 平行移動成分        */
	fixed a00, a01, a02, a03;
	fixed a10, a11, a12, a13;
	fixed a20, a21, a22, a23;
} matrix;

/* 定義済みの行列です。 */
extern const matrix MATRIX_I;	/* 単位(identity)行列 */

/* 行列の初期化関数は用意していません。
 * GCC拡張機能を利用して、次のように書けます。
 *	matrix m;
 *	...
 *	m = (matrix){a00_value,a01_value,...,a23_value};
 *	...
 */

/* 転置(transpose)行列を求めます。
 * [in]
 *	a		転置する行列。
 * [out]
 *	戻り値		結果の行列。
 * [note]
 *	* 平行移動成分は除去され、常に(0,0,0)となります。
 */
matrix mxpose(matrix a);

/* 逆行列(inverse)を求めます。
 * [in]
 *	a		逆行列を求める行列。
 * [out]
 *	戻り値		結果の行列。
 * [note]
 *	* 不可逆な行列を指定した場合の結果は、未定義です。
 */
matrix minv(matrix a);

/* 正規直交(normalized orthogonal)行列の逆行列を求めます。
 * [in]
 *	a		逆行列を求める行列。
 * [out]
 *	戻り値		結果の行列。
 * [note]
 *	* minv()よりも、高速です。
 *	* 正規直交でない行列を指定した場合の結果は、未定義です。
 */
matrix minvno(matrix a);

/* ベクトルを変換(transform)します。
 * [in]
 *	a		行列。
 *	b		ベクトル。
 * [out]
 *	戻り値		結果のベクトル。
 */
vector vxform(matrix a, vector b);

/* 行列を変換(transform)します。
 * [in]
 *	a,b		行列。
 *	b00〜b23	行列の成分。
 * [out]
 *	戻り値		結果の行列。
 */
matrix mxform(matrix a, matrix b);

/*  平行移動(translate)を行います。
 * [in]
 *	a		行列。
 *	x,y,z		X,Y,Z軸の平行移動成分。
 * [out]
 *	戻り値		結果の行列。
 */
matrix mxlate(matrix a, fixed x, fixed y, fixed z);

/* 拡大・縮小(scale)を行います。
 * [in]
 *	a		行列。
 *	x,y,z		X,Y,Z軸の拡大・縮小係数。
 * [out]
 *	戻り値		結果の行列。
 */
matrix mscale(matrix a, fixed x, fixed y, fixed z);

/* X,Y,Z軸とも同じ係数で、拡大・縮小を行います。
 * [in]
 *	a		行列。
 *	s		拡大・縮小係数。
 * [out]
 *	戻り値		結果の行列。
 */
matrix mfscale(matrix a, fixed s);

/* 任意軸回転(rotate)を行います。
 * [in]
 *	a		行列。
 *	x,y,z		回転軸のベクトル。
 *	d		回転角度。
 * [out]
 *	戻り値		結果の行列。
 * [note]
 *	* 回転軸のベクトルは、正規化されていなくても構いません。
 */
matrix mrot(matrix a, fixed x, fixed y, fixed z, fixed d);

/* 固定一軸回転(rotate)を行います。
 * [in]
 *	a		行列。
 *	d		回転角度。
 * [out]
 *	戻り値		結果の行列。
 */
matrix mrotx(matrix a, fixed d);
matrix mroty(matrix a, fixed d);
matrix mrotz(matrix a, fixed d);

/* 固定二軸回転(rotate)を行います。
 * [in]
 *	a		行列。
 *	x,y,z		回転角度。
 * [out]
 *	戻り値		結果の行列。
 * [note]
 *	* 三つの回転角度指定の効果は逆順です。
 *	  固定三軸回転の[note]を参照してください。
 */
matrix mrotxyx(matrix a, fixed x, fixed y, fixed z);
matrix mrotxzx(matrix a, fixed x, fixed y, fixed z);
matrix mrotyxy(matrix a, fixed x, fixed y, fixed z);
matrix mrotyzy(matrix a, fixed x, fixed y, fixed z);
matrix mrotzxz(matrix a, fixed x, fixed y, fixed z);
matrix mrotzyz(matrix a, fixed x, fixed y, fixed z);

/* 固定三軸回転(rotate)を行います。
 * [in]
 *	a		行列。
 *	x,y,z		回転角度。
 * [out]
 *	戻り値		結果の行列。
 * [note]
 *	* 行列を転置したので、逆順の効果になりました。
 *	  例えば、
 *		mrotxyz(o, a, x, y, z)
 *	  は、
 *		mrotz(o, a, z);
 *		mroty(o, o, y);
 *		mrotx(o, o, x);
 *	  と同じ結果となります。
 *	  関数の名前を変えた方がいいのかも知れませんが、
 *	  滅多に使わないので、とりあえずそのままにしておきます。
 */
matrix mrotxyz(matrix a, fixed x, fixed y, fixed z);
matrix mrotxzy(matrix a, fixed x, fixed y, fixed z);
matrix mrotyxz(matrix a, fixed x, fixed y, fixed z);
matrix mrotyzx(matrix a, fixed x, fixed y, fixed z);
matrix mrotzxy(matrix a, fixed x, fixed y, fixed z);
matrix mrotzyx(matrix a, fixed x, fixed y, fixed z);

/* 回転軸のベクトルと、回転角度を求めます。
 * [in]
 *	v		回転軸のベクトルを格納するベクトル。
 *	a		調べる行列。
 * [out]
 *	戻り値		回転角度。
 *	*v		回転軸のベクトル。
 */
fixed mrotget(vector* v, matrix a);

/* 姿勢変換(orientate)を行います。
 * [in]
 *	a		行列。
 *	dx,dy,dz	Z軸の姿勢ベクトル。
 *	ux,uy,uz	Y軸の姿勢ベクトル。
 * [out]
 *	戻り値		結果の行列。
 * [note]
 *	* 姿勢ベクトルは、正規化されていなくても構いません。
 */
matrix morient(matrix a,
	fixed dx, fixed dy, fixed dz,
	fixed ux, fixed uy, fixed uz);

/* 視点変換を行います。
 * [in]
 *	a		行列。
 *	ex,ey,ez	視点座標。
 *	cx,cy,cz	注視点座標。
 *	bank		バンク角度。
 * [out]
 *	戻り値		結果の行列。
 */
matrix mlookat(matrix a,
	fixed ex, fixed ey, fixed ez,
	fixed cx, fixed cy, fixed cz,
	fixed bank);

/* 透視変換(frustum)を行います。
 * [in]
 *	a		行列。
 *	width		投影面の幅              (実座標系)
 *	height		投影面の高さ            (実座標系)
 *	front		フロントクリッピング距離(実座標系)
 * [out]
 *	戻り値		結果の行列。
 * [note]
 *	* 完全な透視変換ではなく、Z座標は実座標のまま変化しません。
 *	  なぜなら、このライブラリでは行列を3x4で表しているので、
 *	  Z座標の正規化を考慮した透視変換行列を表せないからです。
 *	  しかし、zバッファではなく、zソートを想定しているので、
 *	  実用上問題はないと思います。
 *	* この行列で座標変換しただけでは、X,Y座標は正規化されません。
 *	  x=x/z,y=y/z を行って、x,yを(-1,1)に投影する必要があります。
 */
matrix mfrustum(matrix a, fixed width, fixed height, fixed front);

/* 透視変換(perspective)を行います。
 * [in]
 *	a		行列。
 *	fovy		垂直視野角度。
 *	aspect		アスペクト比(幅÷高さ)
 *	front		フロントクリッピング距離(実座標系)
 * [out]
 *	戻り値		結果の行列。
 * [note]
 *	* mfrustum()のコメントを参照してください。
 */
matrix mperspect(matrix a, fixed fovy, fixed aspect, fixed front);

/* ビューポート変換(viewport)を行います。
 * [in]
 *	a		行列。
 *	left,top	ビューポート原点座標(スクリーン座標系)
 *	width,height	ビューポートサイズ  (スクリーン座標系)
 *			※引数型は固定小数ではなく整数であることに
 *			  注意してください。
 * [out]
 *	戻り値		結果の行列。
 * [note]
 *	* 透視変換 -> x=x/z,y=y/z -> ビューポート変換
 *	  の順で適用する場合は、mviewport()を使ってください。
 *	* 透視変換 -> ビューポート変換 -> x=x/z,y=y/z
 *	  の順で適用する場合は、mviewportz()を使ってください。
 *	* mviewport()を使った座標変換の手順は、次のとおりです。
 *		1. T = モデル変換 * 視点変換 * 透視変換 による座標変換
 *		2. Zクリッピング
 *		3. x = x / z, y = y / z
 *		4. T = ビューポート変換 による座標変換
 *		5. XYクリッピング
 *	* これに対し、mviewportz()を使った場合は、
 *		1. T = モデル変換 * 視点変換 * 透視変換 * ビューポート変換
 *	  	   による、座標変換
 *		2. Zクリッピング
 *		3. x = x / z, y = y / z
 *		4. XYクリッピング
 *	  となり、手順を一つ減らすことができます。
 */
matrix mviewport (matrix a, int left, int top, int width, int height);
matrix mviewportz(matrix a, int left, int top, int width, int height);

/* 投影変換を行います。
 * [in]
 *	a		行列。
 *	Lx,Ly,Lz	光源ベクトル
 *	Gx,Gy,Gz	投影面上の一点(任意)の座標
 *	Nx,Ny,Nz	投影面の法線ベクトル
 * [out]
 *	戻り値		結果の行列。
 */
matrix mshadow(matrix a,
	fixed Lx, fixed Ly, fixed Lz,
	fixed Gx, fixed Gy, fixed Gz,
	fixed Nx, fixed Ny, fixed Nz);

/*****************************************************************************
 *	自己代入演算
 *****************************************************************************/

/* fixed */
#define   fchs_(a   )	do { (a) =   fchs((a)     ); } while(0)
#define  ffabs_(a   )	do { (a) =  ffabs((a)     ); } while(0)
#define  fceil_(a   )	do { (a) =  fceil((a)     ); } while(0) /* 2007/02/10追加 */
#define ffloor_(a   )	do { (a) = ffloor((a)     ); } while(0) /* 2007/02/10追加 */
#define   fadd_(a, b)	do { (a) =   fadd((a), (b)); } while(0)
#define  fiadd_(a, b)	do { (a) =  fiadd((a), (b)); } while(0)
#define   fsub_(a, b)	do { (a) =   fsub((a), (b)); } while(0)
#define  fisub_(a, b)	do { (a) =  fisub((a), (b)); } while(0)
#define  fsubr_(a, b)	do { (a) =  fsubr((a), (b)); } while(0)
#define fisubr_(a, b)	do { (a) = fisubr((a), (b)); } while(0)
#define   fmul_(a, b)	do { (a) =   fmul((a), (b)); } while(0)
#define  fimul_(a, b)	do { (a) =  fimul((a), (b)); } while(0)
#define   fdiv_(a, b)	do { (a) =   fdiv((a), (b)); } while(0)
#define  fidiv_(a, b)	do { (a) =  fidiv((a), (b)); } while(0)
#define  fdivr_(a, b)	do { (a) =  fdivr((a), (b)); } while(0)
#define fidivr_(a, b)	do { (a) = fidivr((a), (b)); } while(0)
#define  fprem_(a, b)	do { (a) =  fprem((a), (b)); } while(0) /* 2008/01/29追加 */

/* vector */
#define    vchs_(a   )	do { (a) =    vchs((a)     ); } while(0)
#define    vadd_(a, b)	do { (a) =    vadd((a), (b)); } while(0)
#define    vsub_(a, b)	do { (a) =    vsub((a), (b)); } while(0)
#define    vmul_(a, b)	do { (a) =    vmul((a), (b)); } while(0)
#define   vfmul_(a, b)	do { (a) =   vfmul((a), (b)); } while(0)
#define   vimul_(a, b)	do { (a) =   vimul((a), (b)); } while(0)
#define    vdiv_(a, b)	do { (a) =    vdiv((a), (b)); } while(0)
#define   vfdiv_(a, b)	do { (a) =   vfdiv((a), (b)); } while(0)
#define   vidiv_(a, b)	do { (a) =   vidiv((a), (b)); } while(0)
#define vnormal_(a   )	do { (a) = vnormal((a)     ); } while(0)

/* matrix */
#define     mxpose_(a)						do { (a) =     mxpose((a)); } while(0)
#define       minv_(a)						do { (a) =       minv((a)); } while(0)
#define     minvno_(a)						do { (a) =     minvno((a)); } while(0)
#define     mxform_(a, b)					do { (a) =     mxform((a), (b)); } while(0)
#define     mxlate_(a, x, y, z)					do { (a) =     mxlate((a), (x), (y), (z)); } while(0)
#define     mscale_(a, x, y, z)					do { (a) =     mscale((a), (x), (y), (z)); } while(0)
#define    mfscale_(a, s)					do { (a) =    mfscale((a), (s)); } while(0)
#define       mrot_(a, x, y, z, d)				do { (a) =       mrot((a), (x), (y), (z), (d)); } while(0)
#define      mrotx_(a, d)					do { (a) =      mrotx((a), (d)); } while(0)
#define      mroty_(a, d)					do { (a) =      mroty((a), (d)); } while(0)
#define      mrotz_(a, d)					do { (a) =      mrotz((a), (d)); } while(0)
#define    mrotxyx_(a, x, y, z)					do { (a) =    mrotxyx((a), (x), (y), (z)); } while(0)
#define    mrotxzx_(a, x, y, z)					do { (a) =    mrotxzx((a), (x), (y), (z)); } while(0)
#define    mrotyxy_(a, x, y, z)					do { (a) =    mrotyxy((a), (x), (y), (z)); } while(0)
#define    mrotyzy_(a, x, y, z)					do { (a) =    mrotyzy((a), (x), (y), (z)); } while(0)
#define    mrotzxz_(a, x, y, z)					do { (a) =    mrotzxz((a), (x), (y), (z)); } while(0)
#define    mrotzyz_(a, x, y, z)					do { (a) =    mrotzyz((a), (x), (y), (z)); } while(0)
#define    mrotxyz_(a, x, y, z)					do { (a) =    mrotxyz((a), (x), (y), (z)); } while(0)
#define    mrotxzy_(a, x, y, z)					do { (a) =    mrotxzy((a), (x), (y), (z)); } while(0)
#define    mrotyxz_(a, x, y, z)					do { (a) =    mrotyxz((a), (x), (y), (z)); } while(0)
#define    mrotyzx_(a, x, y, z)					do { (a) =    mrotyzx((a), (x), (y), (z)); } while(0)
#define    mrotzxy_(a, x, y, z)					do { (a) =    mrotzxy((a), (x), (y), (z)); } while(0)
#define    mrotzyx_(a, x, y, z)					do { (a) =    mrotzyx((a), (x), (y), (z)); } while(0)
#define    morient_(a, dx, dy, dz, ux, uy, uz)			do { (a) =    morient((a), (dx), (dy), (dz), (ux), (uy), (uz)); } while(0)
#define    mlookat_(a, ex, ey, ez, cx, cy, cz, bank)		do { (a) =    mlookat((a), (ex), (ey), (ez), (cx), (cy), (cz), (bank)); } while(0)
#define   mfrustum_(a, width, height, front)			do { (a) =   mfrustum((a), (width), (height), (front)); } while(0)
#define  mperspect_(a, fovy, aspect, front)			do { (a) =  mperspect((a), (fovy), (aspect), (front)); } while(0)
#define  mviewport_(a, left, top, width, height)		do { (a) =  mviewport((a), (left), (top), (width), (height)); } while(0)
#define mviewportz_(a, left, top, width, height)		do { (a) = mviewportz((a), (left), (top), (width), (height)); } while(0)
#define    mshadow_(a, Lx, Ly, Lz, Gx, Gy, Gz, Nx, Ny, Nz)	do { (a) =    mshadow((a), (Lx), (Ly), (Lz), (Gx), (Gy), (Gz), (Nx), (Ny), (Nz)); } while(0)

#endif /*__CLIP_FIX_H__*/
