/*
 *	clipflt1.c
 *
 *	高速単精度浮動小数点演算ライブラリ
 *
 *	CLiP - Common Library for P/ECE
 *	Copyright (C) 2001-2016 Naoyuki Sawa
 *
 *	* Sun Jul 20 13:26:33 JST 2008 Naoyuki Sawa
 *	- 作成開始。
 *	* Sun Jul 20 19:20:18 JST 2008 Naoyuki Sawa
 *	- mat3f_rot()等で、誤って、倍精度版のsin()、cos()を呼び出してしまっていました。
 *	  そのため、mat3f_rot()等を使用すると、大きな性能低下となってしまっていました。
 *	  正しく、単精度版のsinf()、cosf()を呼び出すよう、修正しました。
 *	* Sun Aug 17 22:24:29 JST 2008 Naoyuki Sawa
 *	- vecf→vec3f、matf→mat3f、vecx→vec3x、matx→mat3xに名前変更しました。
 *	  全置換で変換したので、本日以前のコメントなどもすべて名前変更されています。
 *	* Mon Aug 18 00:28:27 JST 2008 Naoyuki Sawa
 *	- vec2f、mat2f、および、それらの関数群をを追加しました。
 *	* Tue Aug 19 23:47:09 JST 2008 Naoyuki Sawa
 *	- vec2f、mat2fの関数群をアセンブラ化し、高速化、省メモリ化しました。
 *	- vec2f_xform,mat2f_xformは、vec3f_xform,mat3f_xformと違って、高速RAM配置候補とはしません。
 *	  vec3f_xform,mat3f_xformは、3D処理のなかでまとまって大量に使用されることが考えられますが、
 *	  vec2f_xform,mat2f_xformは、3D処理ほど速度が重要な場面で使うことは少ないと考えたからです。
 *	* Mon Sep 08 16:31:59 JST 2008 Naoyuki Sawa
 *	- clipflt.cをclipflt1.cとclipflt2.cに分割しました。
 *	- 複雑な関数を使わないプログラムにも、常に全関数がリンクされて実行ファイルが大きくなることを防ぐため、
 *	  使用頻度の高い関数をclipflt1.c、使用頻度の低い関数をclipflt2.cに分けることにしました。
 *	  関数分けの明確な基準は無いので、今後も適宜、関数分けや移動を行なう予定です。
 *	  また、今後も関数が増えたら、必要に応じてclipflt3.c,clipflt4.c,...を追加して行く予定です。
 *	* Thu Apr 02 13:37:51 JST 2009 Naoyuki Sawa
 *	- xtof()とftox()を追加しました。
 *	- xtof(),ftox()は、fst(),fld()と同じ働きですが、以下の二点が異なります。
 *	  - fst(),fld()はマクロ定義されており、コンパイル時の変換も可能で、定数テーブルの初期化に利用できます。
 *	    xtof(),ftox()は関数で実装されており、常に実行時処理となり、定数テーブルの初期化に利用できません。
 *	  - xtof(),ftox()は、固定小数点のバイアスが2^10であることを利用して、効率的な実装としました。
 *	    fst(),fld()は、実行時処理に展開された場合は、上述の特性は考慮されず、非効率的な処理となります。
 *	- xtof()とftox()は、%r4,%r5,%r6,%r7の値を破壊しません。
 *	  xtof()とftox()は独自に定義した関数であり、単精度浮動小数点の基本演算サブルーチンではありませんが、
 *	  framflt1.sモジュールで実装した基本演算サブルーチンと同様のレジスタ規約を前提として、利用可能です。
 *	- xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
 *	* Sun Sep 20 17:20:17 JST 2009 Naoyuki Sawa
 *	- solve_simultaneous_equations()を追加しました。
 *	* Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
 *	- clipflt.hに依存しているモジュールを、libclipに含められるようにするために、clipflt.hをlibclip/extclip.hに含める必要が生じました。
 *	  これまでは、clipflt.hをlibclip/extclip.hに含めると、エラーが出ていたので、エラーが出ないように以下の修正を行いました。
 *	- P/ECE以外では、dief()はdie()と同じにマクロ定義しました。
 *	  P/ECEではdie()がpceFontPrintf()を使っているために、それとは別に汎用版のdief()を用意したのですが、P/ECE以外ではdie()自体が汎用版となっており、同じで問題無いからです。
 *	- clipfix.hが先にインクルードされていなければ、固定小数点に依存する機能を除外するようにしました。
 *	  libclip/extclip.hでclipfix.hをインクルードせずに、clipflt.hだけインクルードした場合のためです。
 *	  もし今後、libclip/extclip.hでclipfix.hをインクルードするように変更したら、このブロックも有効になる予定です。
 *	- cliprend.hが先にインクルードされていなければ、描画に依存する機能を除外するようにしました。
 *	  libclip/extclip.hでcliprend.hをインクルードせずに、clipflt.hだけインクルードした場合のためです。
 *	  もし今後、libclip/extclip.hでcliprend.hをインクルードするように変更したら、このブロックも有効になる予定です。
 *	- P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
 *	  これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
 *	  clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
 *	- VC++6.0で単精度版の数学関数の定義が衝突する事がある問題を回避しました。
 *	  詳細は、clipflt.hの、「ANSIライブラリ数学関数(math.lib)と類似の機能の単精度版」の宣言部分の、2016/02/25のコメントを参照して下さい。
 *	* Tue Mar 08 21:09:23 JST 2016 Naoyuki Sawa
 *	- solve_simultaneous_equations()をマルチスレッド化しました。
 *	  詳細は、clipflt1.c 内の、「マルチスレッド版solve_simultaneous_equations()の検証結果」を参照して下さい。
 */
#include "clip.h"

#ifndef PIECE
#define NOASM	// このシンポルを定義すると、C言語版を使用します。 */
#endif /*PIECE*/

/****************************************************************************
 *	
 ****************************************************************************/

static float matnf_determinant(const float* m/*[n][n]*/, int n);
static float matnf_cofactor(const float* m/*[n][n]*/, int n, int i, int j);

/****************************************************************************
 *	デバッグ表示
 ****************************************************************************/

#ifdef  PIECE	//{{2016/02/25追加:P/ECE以外では、dief()はdie()と同じにマクロ定義しました。P/ECEではdie()がpceFontPrintf()を使っているために、それとは別に汎用版のdief()を用意したのですが、P/ECE以外ではdie()自体が汎用版となっており、同じで問題無いからです。}}

/* * Fri Sep 26 21:35:54 JST 2008 Naoyuki Sawa
 * - 浮動小数点の書式("%f")に対応した、デバッグ表示です。
 *   通常のdie()は、pceFontPrintf()を使用しているため、浮動小数点の書式("%f")に対応していませんでした。
 *   浮動小数点のデバッグ表示を簡単にするため、本関数を追加しました。
 *   浮動小数点の書式("%f")が必要無い場合は、通常のdie()を使う方が安全です。(その理由は、後述のdief()の問題点を参照してください。)
 * - 本関数は、snprintf()を呼び出して、_def_vbuff[]の下端80〜87ラインに、書式付き文字列を展開します。
 *   その後、die()に制御を移して、_def_vbuff[]の0〜79ラインに、展開した文字列を描画します。
 * - 描画が79ラインを超えると、文字列を破壊してしまい、正しくない表示になる危険があります。
 *   また、展開後の文字列に'%'が含まれていると、die()の中のpceFontPrintf()によって再展開され、正しくない表示になる危険があります。
 * - 本関数は、浮動小数点のデバッグ表示という限られた用途なので、省メモリを優先して、これらの問題への対策は採らないことにしました。
 *   長過ぎる文字列をデバッグ表示したり、展開後の文字列に'%'を含まないよう、注意してください。
 */
#ifdef NOASM
void dief(const char* fmt, ...) {
	va_list ap;
	va_start(ap, fmt);
	vsnprintf(&_def_vbuff[DISP_X * (DISP_Y - 8)], DISP_X * 8, fmt, ap);
	va_end(ap);
	die("%s", &_def_vbuff[DISP_X * (DISP_Y - 8)]);
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		dief
dief:
		;// [%sp+0]  := retp (noreturn)
		;// [%sp+4]  := fmt
		;// [%sp+8]  := ...
		popn		%r0
		xld.w		%r0, _def_vbuff + 10240			;// %r0     := &_def_vbuff[DISP_X * (DISP_Y - 8)]
		;// %r0      := &_def_vbuff[DISP_X * (DISP_Y - 8)]
		;// [%sp+0]  := fmt
		;// [%sp+4]  := ...
		ld.w		%r12, %r0				;// %r12    := &_def_vbuff[DISP_X * (DISP_Y - 8)]
		xld.w		%r13, 1024				;// %r13    :=                                     DISP_X * 8
		xcall		snprintf				;//   snprintf(&_def_vbuff[DISP_X * (DISP_Y - 8)], DISP_X * 8, fmt, ...)
		;// %r0      := &_def_vbuff[DISP_X * (DISP_Y - 8)]
		pushn		%r0					;// [%sp+0] := &_def_vbuff[DISP_X * (DISP_Y - 8)]
		xcall		die					;//        die(&_def_vbuff[DISP_X * (DISP_Y - 8)])
		;// [%sp+0]  := retp (noreturn)
		;// [%sp+4]  := &_def_vbuff[DISP_X * (DISP_Y - 8)]
");
#endif /*NOASM*/

#endif//PIECE	//{{2016/02/25追加:P/ECE以外では、dief()はdie()と同じにマクロ定義しました。P/ECEではdie()がpceFontPrintf()を使っているために、それとは別に汎用版のdief()を用意したのですが、P/ECE以外ではdie()自体が汎用版となっており、同じで問題無いからです。}}

/****************************************************************************
 *	浮動小数点乱数
 ****************************************************************************/

/* 0.0f以上、1.0f未満の乱数 */
#ifdef NOASM
float rndf(int* seed) {
	int m = rnd32(seed);
	m <<= 8;
	if(!m) {
		return 0.0f;
	} else {
		int e = 127;
		int c;
		do {
			c = m & (1 << 31);
			m <<= 1;
			e--;
		} while(!c);
		m = ((unsigned)m >> 9) | (e << 23);
		return *(float*)&m;
	}
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		rndf
rndf:
		xcall		rnd32			;// %r10 := m = ???????? ???????? ???????? ????????
		sll		%r10, 8			;// %r10 := m = ???????? ???????? ???????? 00000000
		jreq		rndf_RET		;//             |
		;//					;//             +----------------------------------------------+
		xld.w		%r11, 127		;// %r11 := e = bias                                           |
		 add		%r10, %r10		;// %r10 := m = ???????? ???????? ???????0 00000000, %psr(C) <-+
		jruge.d		-1			;//             (一回で抜けても、%r10[8:0]は確実に0になります。)
		sub		%r11, 1			;// %r11 := e--								*delay*
		;//
		or		%r10, %r11		;// %r10 := mmmmmmmm mmmmmmmm mmmmmmm0 eeeeeeee
		rr		%r10, 8			;// %r10 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmm0
rndf_RET:
		ret.d
		rr		%r10, 1			;// %r10 := 0eeeeeee emmmmmmm mmmmmmmm mmmmmmmm				*delay*
");
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

/* a以上、b未満の乱数 */
#ifdef NOASM
float rndf_range(int* seed, float a, float b) {
	return rndf(seed) * (b - a) + a;
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		rndf_range
rndf_range:
		pushn		%r1
		ld.w		%r0, %r13		;// %r0  :=                   a
		xcall.d		rndf			;// %r10 := rndf(seed)
		ld.w		%r1, %r14		;// %r1  :=               b			*delay*
		;//
		ld.w		%r12, %r1		;// %r12 :=               b
		ld.w		%r1, %r10		;// %r1  := rndf(seed)
		xcall.d		__subsf3		;// %r10 :=              (b - a)
		ld.w		%r13, %r0		;// %r13 :=                   a			*delay*
		;//
		ld.w		%r12, %r1		;// %r12 := rndf(seed)
		xcall.d		__mulsf3		;// %r10 := rndf(seed) * (b - a)
		ld.w		%r13, %r10		;// %r13 :=              (b - a)		*delay*
		;//
		ld.w		%r13, %r0		;// %r13 :=                        a
		popn		%r1
		xjp.d		__addsf3		;// %r10 := rndf(seed) * (b - a) + a
		ld.w		%r12, %r10		;// %r12 := rndf(seed) * (b - a)		*delay*
");
#endif /*NOASM*/

/****************************************************************************
 *	固定小数点{値,ベクトル,行列}⇔浮動小数点{値,ベクトル,行列}
 ****************************************************************************/

#ifdef  __CLIP_FIX_H__	//{{2016/02/25追加:clipfix.hが先にインクルードされていなければ、固定小数点に依存する機能を除外するようにしました。libclip/extclip.hでclipfix.hをインクルードせずに、clipflt.hだけインクルードした場合のためです。もし今後、libclip/extclip.hでclipfix.hをインクルードするように変更したら、このブロックも有効になる予定です。}}

/* 固定小数値と浮動小数値の相互変換を行ないます。
 * [in]
 *	xtof:	x		変換する固定小数値。
 *	ftox:	f		変換する浮動小数値。
 * [out]
 *	xtof:	戻り値		変換した浮動小数値。
 *	ftox:	戻り値		変換した固定小数値。
 * [note]
 *	* xtof(),ftox()は、fst(),fld()と同じ働きですが、以下の二点が異なります。
 *	  - fst(),fld()はマクロ定義されており、コンパイル時の変換も可能で、定数テーブルの初期化に利用できます。
 *	    xtof(),ftox()は関数で実装されており、常に実行時処理となり、定数テーブルの初期化に利用できません。
 *	  - xtof(),ftox()は、固定小数点のバイアスが2^10であることを利用して、効率的な実装としました。
 *	    fst(),fld()は、実行時処理に展開された場合は、上述の特性は考慮されず、非効率的な処理となります。
 *	* xtof()とftox()は、%r4,%r5,%r6,%r7の値を破壊しません。
 *	  xtof()とftox()は独自に定義した関数であり、単精度浮動小数点の基本演算サブルーチンではありませんが、
 *	  framflt1.sモジュールで実装した基本演算サブルーチンと同様のレジスタ規約を前提として、利用可能です。
 */
float xtof(fixed x);	/* %r4,%r5,%r6,%r7の値を破壊しません */
fixed ftox(float f);	/* %r4,%r5,%r6,%r7の値を破壊しません */
asm("
		.code
		.align		1
;//-----------------------------------------------------;
		.global		xtof
xtof:
		xcall		__floatsisf		;// %r10 := seeeeeee emmmmmmm mmmmmmmm mmmmmmmm
		swap		%r12, %r10		;// %r12 := emmmmmmm mmmmmmmm mmmmmmmm seeeeeee
		sla		%r10, 1			;// %r10 := eeeeeeem mmmmmmmm mmmmmmmm mmmmmmm0
		jreq		2			;// if(|f|)
		 sub		%r12, 5			;//   exp  -= 10 (アンダーフローは有り得ません)
		ret.d
		swap		%r10, %r12		;// %r10 := seeeeeee emmmmmmm mmmmmmmm mmmmmmmm		*delay*
;//-----------------------------------------------------;
		.global		ftox
ftox:
		swap		%r10, %r12		;// %r10 := emmmmmmm mmmmmmmm mmmmmmmm seeeeeee
		sla		%r12, 1			;// %r12 := eeeeeeem mmmmmmmm mmmmmmmm mmmmmmm0
		jreq		2			;// if(|f|)
		 add		%r10, 5			;//   exp  += 10 (オーバーフローは考慮しません)
		xjp.d		__fixsfsi		;// %r10 := (int)f
		swap		%r12, %r10		;// %r12 := seeeeeee emmmmmmm mmmmmmmm mmmmmmmm		*delay*
");

#ifdef NOASM
void vec3f_conv(vec3f* vec1, const vec3x* vec2) {
	      float* dst = &vec1->x;
	const fixed* src = &vec2->x;
	int cnt = 3;
	do {
//{{2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
//		*dst++ = fst(*src++);
//↓2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
		*dst++ = xtof(*src++);
//}}2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
	} while(--cnt);
}
void mat3f_conv(mat3f* mat1, const mat3x* mat2) {
	      float* dst = &mat1->a00;
	const fixed* src = &mat2->a00;
	int cnt = 12;
	do {
//{{2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
//		*dst++ = fst(*src++);
//↓2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
		*dst++ = xtof(*src++);
//}}2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
	} while(--cnt);
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		vec3f_conv
		.global		mat3f_conv
vec3f_conv:
		jp.d		vec3f_mat3f_conv
		ld.w		%r4, 3			;// %r4  := cnt = 3					*delay*
		;//-------------------------------------;
mat3f_conv:
		ld.w		%r4, 12			;// %r4  := cnt = 12
		;//-------------------------------------;
vec3f_mat3f_conv:
		ld.w		%r5, %r12		;// %r5  := (float*)dst
		ld.w		%r6, %r13		;// %r6  := (fixed*)src
;//{{2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
;//		xld.w		%r7, 0x3a800000		;// %r7  := (1.0f / 1024.0f)
;//vec3f_mat3f_conv_loop:
;//		ld.w		%r12, [%r6]+		;// %r12 := x = (fixed)*src++
;//		xcall		__floatsisf		;// %r10 := f = (float)x
;//		;//
;//		ld.w		%r12, %r10		;// %r12 := f
;//		xcall.d		__mulsf3		;// %r10 := f *= (1.0f / 1024.0f)
;//		ld.w		%r13, %r7		;// %r13 :=      (1.0f / 1024.0f)			*delay*
;//		;//
;//		ld.w		[%r5]+, %r10		;// *dst++ = f
;//↓2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
vec3f_mat3f_conv_loop:
		ld.w		%r12, [%r6]+		;// %r12 := x = (fixed)*src++
		xcall		xtof			;// %r10 := f = xtof(x)
		ld.w		[%r5]+, %r10		;// *dst++ = f
;//}}2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
		sub		%r4, 1			;// %r4  := cnt--
		jrne		vec3f_mat3f_conv_loop
		ret
");
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

#ifdef NOASM
void mat3x_conv(mat3x* mat1, const mat3f* mat2) {
	      fixed* dst = &mat1->a00;
	const float* src = &mat2->a00;
	int cnt = 12;
	do {
//{{2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
//		*dst++ = fld(*src++);
//↓2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
		*dst++ = ftox(*src++);
//}}2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
	} while(--cnt);
}
void vec3x_conv(vec3x* vec1, const vec3f* vec2) {
	      fixed* dst = &vec1->x;
	const float* src = &vec2->x;
	int cnt = 3;
	do {
//{{2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
//		*dst++ = fld(*src++);
//↓2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
		*dst++ = ftox(*src++);
//}}2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
	} while(--cnt);
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		vec3x_conv
		.global		mat3x_conv
vec3x_conv:
		jp.d		vec3x_mat3x_conv
		ld.w		%r4, 3			;// %r4  := cnt = 3					*delay*
		;//-------------------------------------;
mat3x_conv:
		ld.w		%r4, 12			;// %r4  := cnt = 12
		;//-------------------------------------;
vec3x_mat3x_conv:
		ld.w		%r5, %r12		;// %r5  := (fixed*)dst
		ld.w		%r6, %r13		;// %r6  := (float*)src
;//{{2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
;//		xld.w		%r7, 0x44800000		;// %r7  := 1024.0f
;//vec3x_mat3x_conv_loop:
;//		ld.w		%r12, [%r6]+		;// %r12 := f = (float)*src++
;//		xcall.d		__mulsf3		;// %r10 := f *= 1024.0f
;//		ld.w		%r13, %r7		;// %r13 :=      1024.0f				*delay*
;//		;//
;//		xcall.d		__fixsfsi		;// %r10 := x = (fixed)f
;//		ld.w		%r12, %r10		;// %r12 := x						*delay*
;//		;//
;//		ld.w		[%r5]+, %r10		;// *dst++ = x
;//↓2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
vec3x_mat3x_conv_loop:
		ld.w		%r12, [%r6]+		;// %r12 := f = (float)*src++
		xcall		ftox			;// %r10 := x = ftox(f)
		ld.w		[%r5]+, %r10		;// *dst++ = x
;//}}2009/04/02:xtof()とftox()の追加に伴い、vec3f_conv(),vec3x_conv(),mat3f_conv(),mat3x_conv()が、xtof(),ftox()を利用するように変更しました。
		sub		%r4, 1			;// %r4  := cnt--
		jrne		vec3x_mat3x_conv_loop
		ret
");
#endif /*NOASM*/

#endif//__CLIP_FIX_H__	//{{2016/02/25追加:clipfix.hが先にインクルードされていなければ、固定小数点に依存する機能を除外するようにしました。libclip/extclip.hでclipfix.hをインクルードせずに、clipflt.hだけインクルードした場合のためです。もし今後、libclip/extclip.hでclipfix.hをインクルードするように変更したら、このブロックも有効になる予定です。}}

/****************************************************************************
 *	浮動小数点ベクトル(3D)
 ****************************************************************************/

const vec3f vec3f_0 = { 0.0f, 0.0f, 0.0f };
const vec3f vec3f_1 = { 1.0f, 1.0f, 1.0f };
const vec3f vec3f_x = { 1.0f, 0.0f, 0.0f };
const vec3f vec3f_y = { 0.0f, 1.0f, 0.0f };
const vec3f vec3f_z = { 0.0f, 0.0f, 1.0f };

/*--------------------------------------------------------------------------*/

#ifdef NOASM
void vec3f_add(vec3f* vec1, const vec3f* vec2) {
	vec1->x += vec2->x;
	vec1->y += vec2->y;
	vec1->z += vec2->z;
}
void vec3f_sub(vec3f* vec1, const vec3f* vec2) {
	vec1->x -= vec2->x;
	vec1->y -= vec2->y;
	vec1->z -= vec2->z;
}
void vec3f_mul(vec3f* vec1, const vec3f* vec2) {
	vec1->x *= vec2->x;
	vec1->y *= vec2->y;
	vec1->z *= vec2->z;
}
void vec3f_div(vec3f* vec1, const vec3f* vec2) {
	vec1->x /= vec2->x;
	vec1->y /= vec2->y;
	vec1->z /= vec2->z;
}
void vec3f_neg(vec3f* vec) {
	vec->x = -vec->x;
	vec->y = -vec->y;
	vec->z = -vec->z;
}
#else /*NOASM*/
//{{2008/08/19:vec2fとvec3fの処理を統合するため削除しました。
//asm("
//		.code
//		.align		1
//		.global		vec3f_add
//		.global		vec3f_sub
//		.global		vec3f_mul
//		.global		vec3f_div
//		.global		vec3f_neg
//vec3f_neg:
//		xld.w		%r4, __negsf2
//		jp.d		vec3f_op
//		ld.w		%r13, %r12		;// vec2はダミーの読み込みになります。			*delay*
//		;//-------------------------------------;
//vec3f_div:
//		xld.w		%r4, __divsf3
//		jp		vec3f_op
//		;//-------------------------------------;
//vec3f_mul:
//		xld.w		%r4, __mulsf3
//		jp		vec3f_op
//		;//-------------------------------------;
//vec3f_sub:
//		xld.w		%r4, __subsf3
//		jp		vec3f_op
//		;//-------------------------------------;
//vec3f_add:
//		xld.w		%r4, __addsf3
//		;//-------------------------------------;
//vec3f_op:
//		ld.w		%r5, %r12		;// %r5  := vec1
//		ld.w		%r6, %r13		;// %r6  := vec2
//		ld.w		%r7, 3			;// %r7  := cnt
//vec3f_op_loop:
//		ld.w		%r12, [%r5]		;// %r12 := vec1->x
//		ld.w		%r13, [%r6]+		;// %r13 :=           vec2->x
//		call		%r4			;// %r10 := vec1->x + vec2->x
//		ld.w		[%r5]+, %r10
//		sub		%r7, 1			;// %r7  := cnt--
//		jrne		vec3f_op_loop
//		ret
//");
//}}2008/08/19:vec2fとvec3fの処理を統合するため削除しました。
//↓* Tue Aug 19 23:47:09 JST 2008 Naoyuki Sawa
//↓- vec2fのアセンブラ化に併せて、vec2fとvec3fの処理を統合しました。
asm("
		.code
		.align		1
		.global		vec2f_add
		.global		vec2f_sub
		.global		vec2f_mul
		.global		vec2f_div
		.global		vec2f_neg
		.global		vec3f_add
		.global		vec3f_sub
		.global		vec3f_mul
		.global		vec3f_div
		.global		vec3f_neg
;//-----------------------------------------------------;
vec2f_neg:
		ld.w		%r13, %r12		;// vec2はダミーの読み込みになります。
		xld.w		%r4, __negsf2
		jp.d		vec2f_vec3f_op
		ld.w		%r5, 2			;// %r5  := cnt						*delay*
		;//-------------------------------------;
vec2f_div:
		xld.w		%r4, __divsf3
		jp.d		vec2f_vec3f_op
		ld.w		%r5, 2			;// %r5  := cnt						*delay*
		;//-------------------------------------;
vec2f_mul:
		xld.w		%r4, __mulsf3
		jp.d		vec2f_vec3f_op
		ld.w		%r5, 2			;// %r5  := cnt						*delay*
		;//-------------------------------------;
vec2f_sub:
		xld.w		%r4, __subsf3
		jp.d		vec2f_vec3f_op
		ld.w		%r5, 2			;// %r5  := cnt						*delay*
		;//-------------------------------------;
vec2f_add:
		xld.w		%r4, __addsf3
		jp.d		vec2f_vec3f_op
		ld.w		%r5, 2			;// %r5  := cnt						*delay*
;//-----------------------------------------------------;
vec3f_neg:
		ld.w		%r13, %r12		;// vec2はダミーの読み込みになります。
		xld.w		%r4, __negsf2
		jp.d		vec2f_vec3f_op
		ld.w		%r5, 3			;// %r5  := cnt						*delay*
		;//-------------------------------------;
vec3f_div:
		xld.w		%r4, __divsf3
		jp.d		vec2f_vec3f_op
		ld.w		%r5, 3			;// %r5  := cnt						*delay*
		;//-------------------------------------;
vec3f_mul:
		xld.w		%r4, __mulsf3
		jp.d		vec2f_vec3f_op
		ld.w		%r5, 3			;// %r5  := cnt						*delay*
		;//-------------------------------------;
vec3f_sub:
		xld.w		%r4, __subsf3
		jp.d		vec2f_vec3f_op
		ld.w		%r5, 3			;// %r5  := cnt						*delay*
		;//-------------------------------------;
vec3f_add:
		xld.w		%r4, __addsf3
		ld.w		%r5, 3			;// %r5  := cnt
;//-----------------------------------------------------;
vec2f_vec3f_op:
		ld.w		%r6, %r12		;// %r6  := vec1
		ld.w		%r7, %r13		;// %r7  := vec2
vec2f_vec3f_op_loop:
		ld.w		%r12, [%r6]		;// %r12 := vec1->x
		ld.w		%r13, [%r7]+		;// %r13 :=           vec2->x
		call		%r4			;// %r10 := vec1->x + vec2->x
		ld.w		[%r6]+, %r10
		sub		%r5, 1			;// %r5  := cnt--
		jrne		vec2f_vec3f_op_loop
		ret
");
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

#ifdef NOASM
void vec3f_normal(vec3f* vec) {
	float mag = vec3f_mag(vec);
	vec->x /= mag;
	vec->y /= mag;
	vec->z /= mag;
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		vec3f_normal
vec3f_normal:
		xsub		%sp, %sp, 12
		xld.w		[%sp+0], %r12		;// [%sp+0] := vec
		xcall		vec3f_mag		;// %r10    := mag
		xld.w		%r12, [%sp+0]		;// %r12    := vec
		xld.w		[%sp+0], %r10
		xld.w		[%sp+4], %r10
		xld.w		[%sp+8], %r10
		ld.w		%r13, %sp		;// %r13    :=       [mag,mag,mag]
		xcall		vec3f_div		;// %r10    := vec / [mag,mag,mag]
		xadd		%sp, %sp, 12
		ret
");
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

#ifdef NOASM
float vec3f_mag(const vec3f* vec) {
	return sqrtf(vec3f_mag2(vec));
}
float vec3f_mag2(const vec3f* vec) {
	return vec3f_dot(vec, vec);
}
float vec3f_dot(const vec3f* vec1, const vec3f* vec2) {
	return (vec1->x * vec2->x) + (vec1->y * vec2->y) + (vec1->z * vec2->z);
}
#else /*NOASM*/
//{{2008/08/19:vec2fとvec3fの処理を統合するため削除しました。
//asm("
//		.code
//		.align		1
//		.global		vec3f_mag
//		.global		vec3f_mag2
//		.global		vec3f_dot
//vec3f_mag:
//		call		vec3f_mag2
//		xjp.d		sqrtf
//		ld.w		%r12, %r10		;//							*delay*
//		;//-------------------------------------;
//vec3f_mag2:
//		ld.w		%r13, %r12
//		;//-------------------------------------;
//vec3f_dot:
//		ld.w		%r4, %r12		;// %r4  := vec1
//		ld.w		%r5, %r13		;// %r5  := vec2
//		ld.w		%r6, 3			;// %r6  := cnt
//		ld.w		%r7, 0			;// %r7  := sum
//vec3f_dot_loop:
//		ld.w		%r12, [%r4]+		;// %r12 := vec1->x
//		ld.w		%r13, [%r5]+		;// %r13 :=           vec2->x
//		xcall		__mulsf3		;// %r10 := vec1->x * vec2->x
//		ld.w		%r12, %r10		;// %r13 := vec1->x * vec2->x
//		xcall.d		__addsf3		;// %r10 := vec1->x * vec2->x + sum
//		ld.w		%r13, %r7		;// %r12 :=                     sum			*delay*
//		sub		%r6, 1			;// %r6  := cnt--
//		jrne.d		vec3f_dot_loop
//		ld.w		%r7, %r10		;// %r7  := sum						*delay*
//		ret.d
//		ld.w		%r10, %r7		;// %r10 := sum						*delay*
//");
//}}2008/08/19:vec2fとvec3fの処理を統合するため削除しました。
//↓* Tue Aug 19 23:47:09 JST 2008 Naoyuki Sawa
//↓- vec2fのアセンブラ化に併せて、vec2fとvec3fの処理を統合しました。
asm("
		.code
		.align		1
		.global		vec2f_mag
		.global		vec2f_mag2
		.global		vec2f_dot
		.global		vec3f_mag
		.global		vec3f_mag2
		.global		vec3f_dot
;//-----------------------------------------------------;
vec2f_mag:
		call		vec2f_mag2
		xjp.d		sqrtf
		ld.w		%r12, %r10		;//							*delay*
		;//-------------------------------------;
vec2f_mag2:
		ld.w		%r13, %r12
		;//-------------------------------------;
vec2f_dot:
		jp.d		vec2f_vec3f_dot
		ld.w		%r4, 2			;// %r4  := cnt						*delay*
;//-----------------------------------------------------;
vec3f_mag:
		call		vec3f_mag2
		xjp.d		sqrtf
		ld.w		%r12, %r10		;//							*delay*
		;//-------------------------------------;
vec3f_mag2:
		ld.w		%r13, %r12
		;//-------------------------------------;
vec3f_dot:
		ld.w		%r4, 3			;// %r4  := cnt
;//-----------------------------------------------------;
vec2f_vec3f_dot:
		ld.w		%r5, %r12		;// %r5  := vec1
		ld.w		%r6, %r13		;// %r6  := vec2
		ld.w		%r7, 0			;// %r7  := sum
vec2f_vec3f_dot_loop:
		ld.w		%r12, [%r5]+		;// %r12 := vec1->x
		ld.w		%r13, [%r6]+		;// %r13 :=           vec2->x
		xcall		__mulsf3		;// %r10 := vec1->x * vec2->x
		ld.w		%r12, %r10		;// %r13 := vec1->x * vec2->x
		xcall.d		__addsf3		;// %r10 := vec1->x * vec2->x + sum
		ld.w		%r13, %r7		;// %r12 :=                     sum			*delay*
		sub		%r4, 1			;// %r4  := cnt--
		jrne.d		vec2f_vec3f_dot_loop
		ld.w		%r7, %r10		;// %r7  := sum						*delay*
		ret.d
		ld.w		%r10, %r7		;// %r10 := sum						*delay*
");
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

#ifdef NOASM
vec3f vec3f_cross(const vec3f* vec1, const vec3f* vec2) {
	vec3f tmp;
	tmp.x = (vec1->y * vec2->z) - (vec1->z * vec2->y);
	tmp.y = (vec1->z * vec2->x) - (vec1->x * vec2->z);
	tmp.z = (vec1->x * vec2->y) - (vec1->y * vec2->x);
	return tmp;
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		vec3f_cross
vec3f_cross:
		pushn		%r3
		ld.w		%r0, %r12		;// %r0  := 結果を格納する構造体のアドレス
		ld.w		%r1, [%r13]+		;// %r1  := vec1->x
		ld.w		%r2, [%r13]+		;// %r2  := vec1->y
		ld.w		%r3, [%r13]		;// %r3  := vec1->z
		ld.w		%r4, [%r14]+		;// %r4  := vec2->x
		ld.w		%r5, [%r14]+		;// %r5  := vec2->y
		ld.w		%r6, [%r14]		;// %r6  := vec2->z
		;//
		ld.w		%r12, %r2		;// %r12 := vec1->y
		xcall.d		__mulsf3		;// %r10 := vec1->y * vec2->z
		ld.w		%r13, %r6		;// %r13 :=           vec2->z				*delay*
		ld.w		%r7, %r10		;// %r7  := vec1->y * vec2->z
		ld.w		%r12, %r3		;// %r12 :=                     vec1->z
		xcall.d		__mulsf3		;// %r10 :=                     vec1->z * vec2->y
		ld.w		%r13, %r5		;// %r13 :=                               vec2->y	*delay*
		ld.w		%r12, %r7		;// %r12 := vec1->y * vec2->z
		xcall.d		__subsf3		;// %r10 := vec1->y * vec2->z - vec1->z * vec2->y
		ld.w		%r13, %r10		;// %r13 :=                     vec1->z * vec2->y	*delay*
		ld.w		[%r0]+, %r10
		;//
		ld.w		%r12, %r3		;// %r12 := vec1->z
		xcall.d		__mulsf3		;// %r10 := vec1->z * vec2->x
		ld.w		%r13, %r4		;// %r13 :=           vec2->x				*delay*
		ld.w		%r7, %r10		;// %r7  := vec1->z * vec2->x
		ld.w		%r12, %r1		;// %r12 :=                     vec1->x
		xcall.d		__mulsf3		;// %r10 :=                     vec1->x * vec2->z
		ld.w		%r13, %r6		;// %r13 :=                               vec2->z	*delay*
		ld.w		%r12, %r7		;// %r12 := vec1->z * vec2->x
		xcall.d		__subsf3		;// %r10 := vec1->z * vec2->x - vec1->x * vec2->z
		ld.w		%r13, %r10		;// %r13 :=                     vec1->x * vec2->z	*delay*
		ld.w		[%r0]+, %r10
		;//
		ld.w		%r12, %r1		;// %r12 := vec1->x
		xcall.d		__mulsf3		;// %r10 := vec1->x * vec2->y
		ld.w		%r13, %r5		;// %r13 :=           vec2->y				*delay*
		ld.w		%r7, %r10		;// %r7  := vec1->x * vec2->y
		ld.w		%r12, %r2		;// %r12 :=                     vec1->y
		xcall.d		__mulsf3		;// %r10 :=                     vec1->y * vec2->x
		ld.w		%r13, %r4		;// %r13 :=                               vec2->x	*delay*
		ld.w		%r12, %r7		;// %r12 := vec1->x * vec2->y
		xcall.d		__subsf3		;// %r10 := vec1->x * vec2->y - vec1->y * vec2->x
		ld.w		%r13, %r10		;// %r13 :=                     vec1->y * vec2->x	*delay*
		ld.w		[%r0]+, %r10
		;//
		ld.w		%r10, %r0		;// %r10 := 結果を格納する構造体のアドレス + 12
		popn		%r3
		ret.d
		sub		%r10, 12		;// %r10 := 結果を格納する構造体のアドレス		*delay*
");
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

/* 二点間の距離を求める */
#ifdef NOASM
float vec3f_dist(const vec3f* vec1, const vec3f* vec2) {
	return sqrtf(vec3f_dist2(vec1, vec2));
}
/* 二点間の距離の二乗を求める */
float vec3f_dist2(const vec3f* vec1, const vec3f* vec2) {
	vec3f tmp = *vec1;
	vec3f_sub(&tmp, vec2);
	return vec3f_mag2(&tmp);
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		vec3f_dist
		.global		vec3f_dist2
vec3f_dist:
		call		vec3f_dist2		;// %r10    := vec3f_dist2(vec1, vec2)
		xjp.d		sqrtf
		ld.w		%r12, %r10		;//										*delay*
vec3f_dist2:
		xsub		%sp, %sp, 12
		;//
		ld.w		%r4, [%r12]+		;// %r4     := vec1->x
		ld.w		%r5, [%r12]+		;// %r5     := vec1->y
		ld.w		%r6, [%r12]		;// %r6     := vec1->z
		xld.w		[%sp+0], %r4		;// [%sp+0] := tmp.x = vec1->x
		xld.w		[%sp+4], %r5		;// [%sp+4] := tmp.y = vec1->y
		xld.w		[%sp+8], %r6		;// [%sp+8] := tmp.z = vec1->z
		;//
		ld.w		%r12, %sp		;// %r12    :=           &tmp
		xcall		vec3f_sub		;//            vec3f_sub(&tmp, vec2)
		;//
		ld.w		%r12, %sp		;// %r12    :=            &tmp
		xcall		vec3f_mag2		;// %r10    := vec3f_mag2(&tmp)
		;//
		xadd		%sp, %sp, 12
		ret
");
#endif /*NOASM*/

/****************************************************************************
 *	浮動小数点行列(3D)
 ****************************************************************************/

const mat3f mat3f_0 = {
	0.0f, 0.0f, 0.0f, 0.0f,
	0.0f, 0.0f, 0.0f, 0.0f,
	0.0f, 0.0f, 0.0f, 0.0f,
};
const mat3f mat3f_1 = {
	1.0f, 0.0f, 0.0f, 0.0f,
	0.0f, 1.0f, 0.0f, 0.0f,
	0.0f, 0.0f, 1.0f, 0.0f,
};

/*--------------------------------------------------------------------------*/

void mat3f_xlate(mat3f* mat, float x, float y, float z) {
	mat3f tmp = {
		1.0f, 0.0f, 0.0f, x,
		0.0f, 1.0f, 0.0f, y,
		0.0f, 0.0f, 1.0f, z,
	};
	mat3f_xform(mat, &tmp);
}

/*--------------------------------------------------------------------------*/

void mat3f_scale(mat3f* mat, float x, float y, float z) {
	mat3f tmp = {
		   x, 0.0f, 0.0f, 0.0f,
		0.0f,    y, 0.0f, 0.0f,
		0.0f, 0.0f,    z, 0.0f,
	};
	mat3f_xform(mat, &tmp);
}

/*--------------------------------------------------------------------------*/

void mat3f_rot(mat3f* mat, const vec3f* axis, float ang) {
	float s = sinf(ang);
	float c = cosf(ang);
	float ic = 1.0f - c;
	vec3f _axis = *axis;
	vec3f_normal(&_axis);
	{
		mat3f tmp = {
			(_axis.x * _axis.x * ic) + (          c), (_axis.x * _axis.y * ic) - (_axis.z * s), (_axis.x * _axis.z * ic) + (_axis.y * s), 0.0f,
			(_axis.y * _axis.x * ic) + (_axis.z * s), (_axis.y * _axis.y * ic) + (          c), (_axis.y * _axis.z * ic) - (_axis.x * s), 0.0f,
			(_axis.z * _axis.x * ic) - (_axis.y * s), (_axis.z * _axis.y * ic) + (_axis.x * s), (_axis.z * _axis.z * ic) + (          c), 0.0f,
		};
		mat3f_xform(mat, &tmp);
	}
}

/*--------------------------------------------------------------------------*/

void mat3f_rotx(mat3f* mat, float ang) {
	float s = sinf(ang);
	float c = cosf(ang);
	mat3f tmp = {
		1.0f, 0.0f, 0.0f, 0.0f,
		0.0f,    c,   -s, 0.0f,
		0.0f,    s,    c, 0.0f,
	};
	mat3f_xform(mat, &tmp);
}

/*--------------------------------------------------------------------------*/

void mat3f_roty(mat3f* mat, float ang) {
	float s = sinf(ang);
	float c = cosf(ang);
	mat3f tmp = {
		   c, 0.0f,    s, 0.0f,
		0.0f, 1.0f, 0.0f, 0.0f,
		  -s, 0.0f,    c, 0.0f,
	};
	mat3f_xform(mat, &tmp);
}

/*--------------------------------------------------------------------------*/

void mat3f_rotz(mat3f* mat, float ang) {
	float s = sinf(ang);
	float c = cosf(ang);
	mat3f tmp = {
		   c,   -s, 0.0f, 0.0f,
		   s,    c, 0.0f, 0.0f,
		0.0f, 0.0f, 1.0f, 0.0f,
	};
	mat3f_xform(mat, &tmp);
}

/*--------------------------------------------------------------------------*/

void mat3f_invno(mat3f* mat) {
	mat3f tmp = {
		mat->a00, mat->a10, mat->a20, -((mat->a00 * mat->a03) + (mat->a10 * mat->a13) + (mat->a20 * mat->a23)),
		mat->a01, mat->a11, mat->a21, -((mat->a01 * mat->a03) + (mat->a11 * mat->a13) + (mat->a21 * mat->a23)),
		mat->a02, mat->a12, mat->a22, -((mat->a02 * mat->a03) + (mat->a12 * mat->a13) + (mat->a22 * mat->a23)),
	};
	*mat = tmp;
}

void mat3f_inv(mat3f* mat) { /* 逆行列を求めます。正規直交行列でなくても構いません。 */
	int i;
	int j;
	float d;
	float m[4 * 4];

	/* 一時的な、正方行列を作成します。 */
	memcpy(m, mat, sizeof(mat3f));
	m[(4 * 3) + 0] = 0.0f;
	m[(4 * 3) + 1] = 0.0f;
	m[(4 * 3) + 2] = 0.0f;
	m[(4 * 3) + 3] = 1.0f;

	/* クラメルの公式を用いて、逆行列を求めます。 */
	d = matnf_determinant(m, 4);
	for(i = 0; i < 3; i++) {
		for(j = 0; j < 4; j++) {
			((float(*)[4])mat)[i][j] = matnf_cofactor(m, 4, j, i) / d;
		}	//                                    順序注意!!~~~~
	}
}

/*--------------------------------------------------------------------------*/

void mat3f_orient(mat3f* mat, const vec3f* forward, const vec3f* up) {
	vec3f _forward;
	vec3f _right;
	vec3f _up;
	_forward = *forward;
	vec3f_normal(&_forward);
	_right = vec3f_cross(up, &_forward);
	vec3f_normal(&_right);
	_up = vec3f_cross(&_forward, &_right);
	{
		mat3f tmp = {
			_right.x, _up.x, _forward.x, 0.0f,
			_right.y, _up.y, _forward.y, 0.0f,
			_right.z, _up.z, _forward.z, 0.0f,
		};
		mat3f_xform(mat, &tmp);
	}
}

/*--------------------------------------------------------------------------*/

void mat3f_lookat(mat3f* mat, const vec3f* eye, const vec3f* center, float ang) {
	vec3f forward;
	mat3f tmp = mat3f_1;
	mat3f_xlate(&tmp, eye->x, eye->y, eye->z);
	forward = *center;
	vec3f_sub(&forward, eye);
	mat3f_orient(&tmp, &forward, &vec3f_y);
	mat3f_rotz(&tmp, ang);
	mat3f_invno(&tmp);
	mat3f_xform(mat, &tmp);
}

/*--------------------------------------------------------------------------*/

void mat3f_frustum(mat3f* mat, float width, float height, float front) {
	mat3f tmp = {
		2.0f * front / width,                  0.0f, 0.0f, 0.0f,
		                0.0f, 2.0f * front / height, 0.0f, 0.0f,
		                0.0f,                  0.0f, 1.0f, 0.0f,
	};
	mat3f_xform(mat, &tmp);
}

/*--------------------------------------------------------------------------*/

void mat3f_perspect(mat3f* mat, float fovy, float aspect, float front) {
	float height = front * tanf(fovy / 2.0f) * 2.0f;
	float width = height * aspect;
	mat3f_frustum(mat, width, height, front);
}

/*--------------------------------------------------------------------------*/

void mat3f_viewportz(mat3f* mat, int left, int top, int width, int height) {
	mat3f tmp = {
		width / 2.0f,           0.0f, left + width / 2.0f, 0.0f,
		        0.0f, -height / 2.0f, top + height / 2.0f, 0.0f,
		        0.0f,           0.0f,                1.0f, 0.0f,
	};
	mat3f_xform(mat, &tmp);
}

/****************************************************************************
 *	浮動小数点ベクトル(2D)
 ****************************************************************************/

const vec2f vec2f_0 = { 0.0f, 0.0f };
const vec2f vec2f_1 = { 1.0f, 1.0f };
const vec2f vec2f_x = { 1.0f, 0.0f };
const vec2f vec2f_y = { 0.0f, 1.0f };

/*--------------------------------------------------------------------------*/

#ifdef NOASM
void vec2f_xform(const mat2f* mat, vec2f* vec) {
	vec2f tmp = *vec;
	vec->x = (mat->a00 * tmp.x) + (mat->a01 * tmp.y) + mat->a02;
	vec->y = (mat->a10 * tmp.x) + (mat->a11 * tmp.y) + mat->a12;
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		vec2f_xform
vec2f_xform:
		pushn		%r1
		;//
		ld.w		%r0, %r12		;// %r0  := mat
		ld.w		%r1, %r13		;// %r1  := vec
		ld.w		%r4, 2			;// %r4  := cnt
		ld.w		%r5, [%r13]+		;// %r5  := tmp.x
		ld.w		%r6, [%r13]		;// %r6  := tmp.y
vec2f_xform_loop:
		ld.w		%r12, [%r0]+		;// %r12 :=        mat->a00        , %r0  := &mat->a01
		xcall.d		__mulsf3		;// %r10 := sum  = mat->a00 * tmp.x
		ld.w		%r13, %r5		;// %r13 :=                   tmp.x			*delay*
		ld.w		%r7, %r10		;// %r7  := sum
		;//
		ld.w		%r12, [%r0]+		;// %r12 :=        mat->a01        , %r0  := &mat->a02
		xcall.d		__mulsf3		;// %r10 :=        mat->a01 * tmp.y
		ld.w		%r13, %r6		;// %r13 :=                   tmp.y			*delay*
		ld.w		%r12, %r10		;// %r12 :=        mat->a01 * tmp.y
		xcall.d		__addsf3		;// %r10 := sum += mat->a01 * tmp.y
		ld.w		%r13, %r7		;// %r13 := sum						*delay*
		;//
		ld.w		%r12, [%r0]+		;// %r12 :=        mat->a02        , %r0  := &mat->a10
		xcall.d		__addsf3		;// %r10 := sum += mat->a02
		ld.w		%r13, %r10		;// %r13 := sum						*delay*
		;//
		ld.w		[%r1]+, %r10		;// vec->x = sum
		sub		%r4, 1			;// %r4  := cnt--
		jrne		vec2f_xform_loop
		;//
		popn		%r1
		ret
");
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

#ifdef NOASM
void vec2f_add(vec2f* vec1, const vec2f* vec2) {
	vec1->x += vec2->x;
	vec1->y += vec2->y;
}
void vec2f_sub(vec2f* vec1, const vec2f* vec2) {
	vec1->x -= vec2->x;
	vec1->y -= vec2->y;
}
void vec2f_mul(vec2f* vec1, const vec2f* vec2) {
	vec1->x *= vec2->x;
	vec1->y *= vec2->y;
}
void vec2f_div(vec2f* vec1, const vec2f* vec2) {
	vec1->x /= vec2->x;
	vec1->y /= vec2->y;
}
void vec2f_neg(vec2f* vec) {
	vec->x = -vec->x;
	vec->y = -vec->y;
}
#else /*NOASM*/
// * Tue Aug 19 23:47:09 JST 2008 Naoyuki Sawa
// - vec2fのアセンブラ化に併せて、vec2fとvec3fの処理を統合しました。
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

#ifdef NOASM
void vec2f_normal(vec2f* vec) {
	float mag = vec2f_mag(vec);
	vec->x /= mag;
	vec->y /= mag;
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		vec2f_normal
vec2f_normal:
		xsub		%sp, %sp, 8
		xld.w		[%sp+0], %r12		;// [%sp+0] := vec
		xcall		vec2f_mag		;// %r10    := mag
		xld.w		%r12, [%sp+0]		;// %r12    := vec
		xld.w		[%sp+0], %r10
		xld.w		[%sp+4], %r10
		ld.w		%r13, %sp		;// %r13    :=       [mag,mag]
		xcall		vec2f_div		;// %r10    := vec / [mag,mag]
		xadd		%sp, %sp, 8
		ret
");
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

#ifdef NOASM
float vec2f_mag(const vec2f* vec) {
	return sqrtf(vec2f_mag2(vec));
}
float vec2f_mag2(const vec2f* vec) {
	return vec2f_dot(vec, vec);
}
float vec2f_dot(const vec2f* vec1, const vec2f* vec2) {
	return (vec1->x * vec2->x) + (vec1->y * vec2->y);
}
#else /*NOASM*/
// * Tue Aug 19 23:47:09 JST 2008 Naoyuki Sawa
// - vec2fのアセンブラ化に併せて、vec2fとvec3fの処理を統合しました。
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

#ifdef NOASM
float vec2f_cross(const vec2f* vec1, const vec2f* vec2) {
	return (vec1->x * vec2->y) - (vec1->y * vec2->x);
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		vec2f_cross
vec2f_cross:
		ld.w		%r4, [%r12]+		;// %r4  := vec1->x
		ld.w		%r5, [%r12]		;// %r5  := vec1->y
		ld.w		%r6, [%r13]+		;// %r6  := vec2->x
		;//
		ld.w		%r13, [%r13]		;// %r13 :=            vec2->y
		xcall.d		__mulsf3		;// %r10 := (vec1->x * vec2->y)
		ld.w		%r12, %r4		;// %r12 :=  vec1->x					*delay*
		ld.w		%r4, %r10		;// %r4  := (vec1->x * vec2->y)
		;//
		ld.w		%r12, %r5		;// %r12 :=  vec1->y
		xcall.d		__mulsf3		;// %r10 := (vec1->y * vec2->x)
		ld.w		%r13, %r6		;// %r13 :=            vec2->x				*delay*
		;//
		ld.w		%r12, %r4		;// %r12 := (vec1->x * vec2->y)
		xjp.d		__subsf3		;// %r10 := (vec1->x * vec2->y) - (vec1->y * vec2->x)
		ld.w		%r13, %r10		;// %r13 :=                       (vec1->y * vec2->x)	*delay*
");
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

/* 二点間の距離を求める */
#ifdef NOASM
float vec2f_dist(const vec2f* vec1, const vec2f* vec2) {
	return sqrtf(vec2f_dist2(vec1, vec2));
}
/* 二点間の距離の二乗を求める */
float vec2f_dist2(const vec2f* vec1, const vec2f* vec2) {
	vec2f tmp = *vec1;
	vec2f_sub(&tmp, vec2);
	return vec2f_mag2(&tmp);
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		vec2f_dist
		.global		vec2f_dist2
vec2f_dist:
		call		vec2f_dist2		;// %r10    := vec2f_dist2(vec1, vec2)
		xjp.d		sqrtf
		ld.w		%r12, %r10		;//										*delay*
vec2f_dist2:
		xsub		%sp, %sp, 8
		;//
		ld.w		%r4, [%r12]+		;// %r4     := vec1->x
		ld.w		%r5, [%r12]		;// %r5     := vec1->y
		xld.w		[%sp+0], %r4		;// [%sp+0] := tmp.x = vec1->x
		xld.w		[%sp+4], %r5		;// [%sp+4] := tmp.y = vec1->y
		;//
		ld.w		%r12, %sp		;// %r12    :=           &tmp
		xcall		vec2f_sub		;//            vec2f_sub(&tmp, vec2)
		;//
		ld.w		%r12, %sp		;// %r12    :=            &tmp
		xcall		vec2f_mag2		;// %r10    := vec2f_mag2(&tmp)
		;//
		xadd		%sp, %sp, 8
		ret
");
#endif /*NOASM*/

/****************************************************************************
 *	浮動小数点行列(2D)
 ****************************************************************************/

const mat2f mat2f_0 = {
	0.0f, 0.0f, 0.0f,
	0.0f, 0.0f, 0.0f,
};
const mat2f mat2f_1 = {
	1.0f, 0.0f, 0.0f,
	0.0f, 1.0f, 0.0f,
};

/*--------------------------------------------------------------------------*/

#ifdef NOASM
void mat2f_xform(mat2f* mat1, const mat2f* mat2) {
	mat2f tmp = *mat1;
	mat1->a00 = (tmp.a00 * mat2->a00) + (tmp.a01 * mat2->a10);
	mat1->a01 = (tmp.a00 * mat2->a01) + (tmp.a01 * mat2->a11);
	mat1->a02 = (tmp.a00 * mat2->a02) + (tmp.a01 * mat2->a12) + tmp.a02;
	mat1->a10 = (tmp.a10 * mat2->a00) + (tmp.a11 * mat2->a10);
	mat1->a11 = (tmp.a10 * mat2->a01) + (tmp.a11 * mat2->a11);
	mat1->a12 = (tmp.a10 * mat2->a02) + (tmp.a11 * mat2->a12) + tmp.a12;
}
#else /*NOASM*/
asm("
		.code
		.align		1
		.global		mat2f_xform
mat2f_xform:
		pushn		%r2
		xsub		%sp, %sp, 24		;// %sp  := &tmp
		ld.w		%r0, %r12		;// %r0  := mat1
		ld.w		%r1, %r13		;// %r1  := mat2
		ld.w		%r12, %sp		;// %r12 := &tmp
		xld.w		%r14, 24		;// %r14 :=             sizeof(mat2f)
		xcall.d		memcpy			;//  memcpy(&tmp, mat1, sizeof(mat2f))
		ld.w		%r13, %r0		;// %r13 :=       mat1					*delay*
		;//
		xld.w		%r2, __mulsf3		;// %r2  := __mulsf3
		xld.w		%r4, __addsf3		;// %r4  := __addsf3
		ld.w		%r5, 2			;// %r5  := cnt1
		;//-------------------------------------;
mat2f_xform_lp1:
		ld.w		%r6, 3			;// %r6  := cnt2
		;//-------------------------------------;
mat2f_xform_lp2:
		xld.w		%r12, [%sp+0]		;// %r12 :=        tmp.a00
		ld.w		%r13, [%r1]		;// %r13 :=                  mat2->a00
		call.d		%r2			;// %r10 := sum  = tmp.a00 * mat2->a00 -+
		add		%r1, 12			;// %r1  :=                 &mat2->a10  |		*delay*
		;//					;//                                     |
		xld.w		%r12, [%sp+4]		;// %r12 :=        tmp.a01              |
		ld.w		%r13, [%r1]		;// %r13 :=                  mat2->a10  |
		call.d		%r2			;// %r10 :=        tmp.a01 * mat2->a10 -|-+
		ld.w		%r7, %r10		;// %r7  := sum <-----------------------+ |		*delay*
		;//					;//                                       |
		ld.w		%r12, %r10		;// %r12 :=        tmp.a01 * mat2->a10 <--+
		call.d		%r4			;// %r10 := sum += tmp.a01 * mat2->a10 -----+
		ld.w		%r13, %r7		;// %r13 := sum                             |		*delay*
		;//					;//                                         |
		sub		%r6, 1			;// %r6  := cnt2--                          |
		jreq		4			;// if(cnt2)                                |
		 ld.w		[%r0]+, %r10		;//   mat1->a00 = sum, %r0  := &mat1->a01 <-+
		 jp.d		mat2f_xform_lp2		;//                                         |
		 sub		%r1, 8			;//                    %r1  := &mat2->a01   |		*delay*
		;//-------------------------------------;// else                                    |
		xld.w		%r12, [%sp+8]		;//   %r12 :=        tmp.a02                |
		call.d		%r4			;//   %r10 := sum += tmp.a02 ---------------|-+
		ld.w		%r13, %r10		;//   %r13 := sum <-------------------------+ |		*delay*
		ld.w		[%r0]+, %r10		;//   mat1->a02 = sum, %r0  := &mat1->a10 <---+
		;//
		xadd		%sp, %sp, 12		;//                    %sp  :=   &tmp.a10
		sub		%r5, 1			;//   %r5  := cnt1--
		jrne.d		mat2f_xform_lp1
		sub		%r1, 20			;//                    %r1  := &mat2->a00		*delay*
		;//-------------------------------------;
		;// %spは2回+12して既に元に戻っています。
		popn		%r2
		ret
");
#endif /*NOASM*/

/*--------------------------------------------------------------------------*/

void mat2f_xlate(mat2f* mat, float x, float y) {
	mat2f tmp = {
		1.0f, 0.0f, x,
		0.0f, 1.0f, y,
	};
	mat2f_xform(mat, &tmp);
}
void mat2f_scale(mat2f* mat, float x, float y) {
	mat2f tmp = {
		   x, 0.0f, 0.0f,
		0.0f,    y, 0.0f,
	};
	mat2f_xform(mat, &tmp);
}
void mat2f_rot(mat2f* mat, float ang) {
	float s = sinf(ang);
	float c = cosf(ang);
	mat2f tmp = {
		c, -s, 0.0f,
		s,  c, 0.0f,
	};
	mat2f_xform(mat, &tmp);
}

/*--------------------------------------------------------------------------*/

void mat2f_invno(mat2f* mat) {
	mat2f tmp = {
		mat->a00, mat->a10, -((mat->a00 * mat->a02) + (mat->a10 * mat->a12)),
		mat->a01, mat->a11, -((mat->a01 * mat->a02) + (mat->a11 * mat->a12)),
	};
	*mat = tmp;
}

void mat2f_inv(mat2f* mat) { /* 逆行列を求めます。正規直交行列でなくても構いません。 */
	int i;
	int j;
	float d;
	float m[3 * 3];

	/* 一時的な、正方行列を作成します。 */
	memcpy(m, mat, sizeof(mat2f));
	m[(3 * 2) + 0] = 0.0f;
	m[(3 * 2) + 1] = 0.0f;
	m[(3 * 2) + 2] = 1.0f;

	/* クラメルの公式を用いて、逆行列を求めます。 */
	d = matnf_determinant(m, 3);
	for(i = 0; i < 2; i++) {
		for(j = 0; j < 3; j++) {
			((float(*)[3])mat)[i][j] = matnf_cofactor(m, 3, j, i) / d;
		}	//                                    順序注意!!~~~~
	}
}

/****************************************************************************
 *	逆行列を求めるために利用するサブルーチン
 ****************************************************************************/

//m[n行][n列]の行列式(determinant)を求めます。
static float matnf_determinant(const float* m/*[n][n]*/, int n) {
	int i;
	float d;
	if(n < 1) { DIE(); }
	//余因子展開を用いて、行列式(determinant)を求めます。
	d = 0.0f;
	for(i = 0; i < n; i++) {
		d += (m[(n * i) + 0/*j*/] * matnf_cofactor(m, n, i, 0/*j*/));
	}
	return d;
}
//m[n行][n列]からi行j列を除いた、余因子(cofactor)を求めます。
static float matnf_cofactor(const float* m/*[n][n]*/, int n, int i, int j) {
	int k, l;
	float d;
	float* s/*[n-1][n-1]*/;
	float* t/*[n-1][n-1]*/;
	if((n < 1) ||
	   (i < 0) || (i > (n - 1)) ||
	   (j < 0) || (j > (n - 1))) { DIE(); }
	//m[n行][n列]からi行j列を除いた、小行列(submatrix)を作成します。
	if(n == 1) { return 1; } //小行列が0行0列になる場合
//{{2016/03/08変更:solve_simultaneous_equations()をマルチスレッド化しました。
//	s = malloc(sizeof(float) * ((n - 1) * (n - 1)));
//	if(!s) { DIE(); }
//↓2016/03/08変更:solve_simultaneous_equations()をマルチスレッド化しました。
//│ * Tue Mar 08 21:16:31 JST 2016 Naoyuki Sawa
//│ - solve_simultaneous_equations()をマルチスレッド化しました。
//│ - 上記の変更のために、Win32環境でマルチスレッドを使う場合は、matnf_cofactor()内で、malloc()ではなくalloca()を使うように修正しました。
//│   matnf_cofactor()が再帰で何度も呼び出され、マルチスレッドがmalloc()の排他制御で直列化され、マルチスレッドの効果が得られないからです。
//│   こうしないと、マルチスレッド化した方が却って遅くなってしまいます。
//│   詳細は、後述の「マルチスレッド版solve_simultaneous_equations()の検証結果」を参照して下さい。
//│ - 尚、実際には、Win32環境でマルチスレッドを使う場合だけでなく、Win32環境でシングルスレッドの場合も、alloca()を使う方が早いのですが、
//│   特にマルチスレッドの場合はalloca()を使う事が必須である事を強調するために、シングルスレッドの場合はmalloc()のままにする事にしました。
//│ - P/ECE環境では、これまで通り、malloc()を使います。
//│   P/ECE環境の独自実装のalloca() (clipalca.c参照)は、互換性のためだけに実装したものであり、malloc()よりも効率が悪いからです。
#if     !(defined(_WIN32)&&defined(_MT))
//Win32シングルスレッドの場合,P/ECEの場合,その他の環境の場合に有効なブロック
	s = malloc(sizeof(float) * ((n - 1) * (n - 1)));
	if(!s) { DIE(); }
#else //!(defined(_WIN32)&&defined(_MT))
//Win32マルチスレッドの場合に有効なブロック
	s = alloca(sizeof(float) * ((n - 1) * (n - 1)));
#endif//!(defined(_WIN32)&&defined(_MT))
//}}2016/03/08変更:solve_simultaneous_equations()をマルチスレッド化しました。
	t = s;
	for(k = 0; k < n; k++) {
		for(l = 0; l < n; l++) {
			if((k != i) && (l != j)) {
				*t++ = *m;
			}
			m++;
		}
	}
	//小行列式(minor determinant)を求めます。
	d = matnf_determinant(s, (n - 1));
	//余因子(cofactor)を求めます。
	if((i ^ j) & 1) { d = -d; }
	//小行列を解放します。
//{{2016/03/08変更:solve_simultaneous_equations()をマルチスレッド化しました。
//	free(s);
//↓2016/03/08変更:solve_simultaneous_equations()をマルチスレッド化しました。
#if     !(defined(_WIN32)&&defined(_MT))
//Win32シングルスレッドの場合,P/ECEの場合,その他の環境の場合に有効なブロック
	free(s);
#endif//!(defined(_WIN32)&&defined(_MT))
//}}2016/03/08変更:solve_simultaneous_equations()をマルチスレッド化しました。
	return d;
}

/****************************************************************************
 *	連立方程式
 ****************************************************************************/

//{{2016/03/08変更:solve_simultaneous_equations()をマルチスレッド化しました。
//void solve_simultaneous_equations(const float* A/*[n][n]*/, const float* b/*[n]*/, float* x/*[n]*/, int n) {
//	int i;
//	int j;
//	float d;
//	//クラメルの公式を用いて、Aの逆行列を求めながら、bを掛けて、xを求めます。
//	d = matnf_determinant(A, n);
//	for(i = 0; i < n; i++) {
//		x[i] = 0.0f;
//		for(j = 0; j < n; j++) {
//			x[i] += (matnf_cofactor(A, n, j, i) / d) * b[j];
//			//                  順序注意!!~~~~
//		}
//	}
//}
//↓2016/03/08変更:solve_simultaneous_equations()をマルチスレッド化しました。
#if     !(defined(_WIN32)&&defined(_MT))
//Win32シングルスレッドの場合,P/ECEの場合,その他の環境の場合に有効なブロック
void solve_simultaneous_equations(const float* A/*[n][n]*/, const float* b/*[n]*/, float* x/*[n]*/, int n) {
	int i, j;
	float d;
	//行列式(determinant)を求めます。
	d = matnf_determinant(A, n);
	//クラメルの公式を用いて、Aの逆行列を求めながら、bを掛けて、xを求めます。
	for(i = 0; i < n; i++) {
		x[i] = 0.0f;
		for(j = 0; j < n; j++) {
			x[i] += (matnf_cofactor(A, n, j, i) / d) * b[j];
			//                  順序注意!!~~~~
		}
	}
}
#else //!(defined(_WIN32)&&defined(_MT))
//Win32マルチスレッドの場合に有効なブロック
struct _ST_solve_simultaneous_equations {
	const float*		A/*[n][n]*/;	//┬solve_simultaneous_equations()の引数と同じです。solve_simultaneous_equations_start_routine()からは参照専用です。書き換えてはいけません。
	const float*		b/*[n]*/;	//│
	float*			x/*[n]*/;	//│
	int			n;		//┘
	float			d;		//─solve_simultaneous_equations()が値を格納します。solve_simultaneous_equations_start_routine()からは参照専用です。書き換えてはいけません。
	volatile int		i;		//┬i＝0→n。solve_simultaneous_equations_start_routine()が読み書きします。iを読み書きする時はcsをロックして下さい。
	CRITICAL_SECTION 	cs;		//┘
};
static unsigned __stdcall solve_simultaneous_equations_start_routine(void* _arg) {
	struct _ST_solve_simultaneous_equations* arg = _arg;
	const float* const A = arg->A;
	const float* const b = arg->b;
	float*       const x = arg->x;
	int          const n = arg->n;
	float        const d = arg->d;
	for(;;) {
		int i, j;
	    EnterCriticalSection(&arg->cs);	//{{arg->iの読み書き:ここから
		i = arg->i;
		if(i != n) { arg->i = (i + 1); }
	    LeaveCriticalSection(&arg->cs);	//}}arg->iの読み書き:ここまで
		if(i == n) { break; }
		x[i] = 0.0f;
		for(j = 0; j < n; j++) {
			x[i] += (matnf_cofactor(A, n, j, i) / d) * b[j];
			//                  順序注意!!~~~~
		}
	}
	return 0;
}
void solve_simultaneous_equations(const float* A/*[n][n]*/, const float* b/*[n]*/, float* x/*[n]*/, int n) {
	struct _ST_solve_simultaneous_equations arg;
	HANDLE hThread[16/*調整可*/];	//当関数が作成するスレッドの数はこの配列サイズが上限となる。2016/03/08現在は16スレッドを上限とした。もし今後16コア以上のCPUが一般的になった場合はこの配列サイズを増やせ。
	int nThread, iThread;
	SYSTEM_INFO si;
	//構造体を初期化する。
	arg.A = A;
	arg.b = b;
	arg.x = x;
	arg.n = n;
	arg.d = matnf_determinant(A, n);	//行列式(determinant)を求めます。
	arg.i = 0;
	//クリティカルセクションを作成する。
	InitializeCriticalSection(&arg.cs);
	//スレッド数を決めるために、CPUの論理コア数を取得する。
	GetSystemInfo(&si);
	nThread = si.dwNumberOfProcessors;
	//CPUの論理コア数が、hThread[]配列のサイズよりも多い場合は、切り詰める。
	if(nThread > ARRAY_SIZE(hThread)) { nThread = ARRAY_SIZE(hThread); }
	//CPUの論理コア数が、配列の次元数よりも多い場合は、切り詰める。
	// - 配列の次元数より多いスレッドを生成しても、動作上は問題無いが、何もしないスレッドが生じて無駄である。
	if(nThread > n) { nThread = n; }
	//スレッドを作成する。
	for(iThread = 0; iThread < nThread; iThread++) {
		hThread[iThread] = (HANDLE)_beginthreadex(NULL,
			0,
			solve_simultaneous_equations_start_routine,
			&arg,
			0,
			NULL);
		if(!hThread[iThread]) { DIE(); }	//バグ,又は,リソース不足
	}
	//スレッドの完了を待つ。
	iThread = WaitForMultipleObjects(nThread, hThread, TRUE, INFINITE);
	if((unsigned)(iThread - WAIT_OBJECT_0) >= (unsigned)nThread) { DIE(); }	//バグ
	//スレッドを開放する。
	for(iThread = 0; iThread < nThread; iThread++) {
		if(!CloseHandle(hThread[iThread])) { DIE(); }	//バグ
	}
	//クリティカルセクションを削除する。
	DeleteCriticalSection(&arg.cs);	//忘れないで!!
}
#endif//!(defined(_WIN32)&&defined(_MT))
//}}2016/03/08変更:solve_simultaneous_equations()をマルチスレッド化しました。

#if 0
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
■マルチスレッド版solve_simultaneous_equations()の検証結果
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
□テストプログラム
│void test() {
│  static const float A[10*10]={1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
│                               1, 1, 1,-1, 1, 1, 1, 1, 1, 1,
│                               1, 1,-1, 1, 2, 3, 4, 4, 4, 4,
│                               1,-1, 1, 1, 1, 1, 2, 2, 2, 2,
│                               1, 1, 9, 1, 3, 1, 1, 1, 1, 1,
│                               1, 1, 8, 1, 1, 1, 4, 4, 4, 4,
│                               1, 1, 7, 1, 1, 1, 1, 5, 5, 5,
│                               1, 1, 7, 1, 1, 1, 1, 1, 1, 3,
│                               1, 1, 7, 1, 1, 1, 1, 1, 1, 1,
│                               1, 1, 7, 1, 1, 1, 2, 3, 1, 1};
│  static const float b[   10]={0, 4,-4, 2, 2, 2, 3, 5, 6, 7};
│               float x[   10];
│  int i,t=GetTickCount();
│  solve_simultaneous_equations(A,b,x,ARRAY_SIZE(x));
│  t=GetTickCount()-t;
│  for(i=0;i<ARRAY_SIZE(x);i++){printf("x[%d]=% f\n",i,x[i]);}
│  printf("time=%d\n",t);
│}
□結果
│x[0]= 4.500000
│x[1]=-1.833333
│x[2]= 1.000000
│x[3]=-2.000000
│x[4]=-3.000000
│x[5]= 3.000000
│x[6]=-0.916667
│x[7]= 0.958333
│x[8]=-1.208333
│x[9]=-0.500000
│上記の結果が正しいことを、「高精度計算サイト - 行列式(n次元)」(http://keisan.casio.jp/exec/system/1279265553)でも計算して確認した。
│以下の、どのテスト環境・どのスレッド構成でも、上記の結果になる事を確認した。
───────────────────────────────────────
□テスト環境
│OS  Windows 10 Home
│CPU Core i7-4500U
│RAM 8 GB
□時間
│△ 1スレッド、malloc()を使用	time=13719	┬malloc()を使うとスレッドが増えるほど遅くなっている。当関数はmalloc()の使用頻度が高いので、malloc()で直列化されてマルチスレッドの効果が得られず、オーバーヘッドだけが増えているのだろう。
│・ 2スレッド、malloc()を使用	time=22578	│尚、現在のビルド環境ではBoehm GCも併用しており、上記のテストプログラムで使用したメモリアロケータは、VC++6.0標準のLIBCMTのmalloc()ではなく、Boehm GCのGC_MALLOC()である。
│・ 4スレッド、malloc()を使用	time=38766	┘Boehm GCがGC_MALLOC()がLIBCMTのmalloc()よりも遅い可能性は有るが、LIBCMTのmalloc()でも何らかの排他制御を行っているはずであり、LIBCMTのmalloc()を使ったとしても今回のテスト結果と同じ傾向になると思う。
│・ 1スレッド、alloca()を使用	time= 2781	┐
│・ 2スレッド、alloca()を使用	time= 1750	├alloca()を使うとスレッドが増えるほど早くなっている。理論上は、'4スレッド、alloca()を使用'が、当CPUでの最適な構成であるはず。ただし、2⇒4の差が少ない事から考えて、HyperThreadingの効果は小さいようだ。
│○ 4スレッド、alloca()を使用	time= 1625	┘
│・10スレッド、alloca()を使用	time= 1437	─試しに、行列式の次数である10スレッドに分割してみると、何故か4スレッド時よりも早くなった。CPUの論理コア数以上に分割しても早くなるはずがないのだが… TurboBoostでも掛かっているのだろうか?
│2016/03/08変更前のプログラム構成では、'△'の動作でした。
│2016/03/08変更後のプログラム構成では、'○'の動作になりました。
───────────────────────────────────────
□テスト環境
│OS  Windows XP Home
│CPU Pentium 4 530
│RAM 1.5 GB
□時間
│△ 1スレッド、malloc()を使用	time=65641	┐┐
│・ 2スレッド、malloc()を使用	time=66343	┘├Pentium 4では、当関数のマルチスレッド化の効果は得られず、逆に遅くなってしまいました。FP演算器が一個しか無い(?)ので、当関数のようなFP演算を多用する処理では、FP演算器で競合してしまうのだと思います。
│・ 1スレッド、alloca()を使用	time=12984	┐│冷却の悪い環境でテストしたので、もしかしたらTM2(ThermalMonitoring2)が働いて、普通よりも遅い結果になっているかも知れません。しかし、スレッド構成ごとの相対的な比較としては、このとおりだと思います。
│○ 2スレッド、alloca()を使用	time=13297	┘┘
│2016/03/08変更前のプログラム構成では、'△'の動作でした。
│2016/03/08変更後のプログラム構成では、'○'の動作になりました。
#endif

/****************************************************************************
 *	
 ****************************************************************************/
