/*	
 *	clipmath.c
 *
 *	ユーティリティ：数学関数
 *
 *	CLiP - Common Library for P/ECE
 *	Copyright (C) 2001-2017 Naoyuki Sawa
 *
 *	* Thu Dec 31 09:54:58 JST 2015 Naoyuki Sawa
 *	- 作成開始。
 *	- math.hには含まれない非標準の、数学関数をまとめるためのモジュールです。
 *	  今後、非標準の数学関数を追加する場合は、当モジュールに追加して下さい。
 *	* Sat Mar 19 21:19:21 JST 2016 Naoyuki Sawa
 *	- determinant(),cofactor(),solve()を追加しました。
 *	* Wed May 24 21:26:30 JST 2017 Naoyuki Sawa
 *	- easeIn(),easeOut(),easeInOut()を追加しました。
 *	* Fri Jul 28 21:55:48 JST 2017 Naoyuki Sawa
 *	- lerp_i()を追加しました。
 *	* Sat Jul 29 22:27:07 JST 2017 Naoyuki Sawa
 *	- easeIn_i(),easeOut_i(),easeInOut_i()を追加しました。
 *	* Tue Aug 08 21:44:58 JST 2017 Naoyuki Sawa
 *	- easeOutIn(),easeOutIn_i()を追加しました。
 *	* Tue Aug 08 22:18:34 JST 2017 Naoyuki Sawa
 *	- easeInOut()とeaseOutIn()の処理を共通化して、コードサイズを削減しました。
 *	  easeInOut_i()とeaseOutIn_i()の処理を共通化して、コードサイズを削減しました。
 *	  コードサイズ削減の効果は、前者84バイト,後者12バイトで、合計96バイトも小さくなりました。
 *	* Thu Aug 10 21:41:17 JST 2017 Naoyuki Sawa
 *	- clipmath.cから、整数版の関数の実装を、clipmati.cへ分離しました。
 *	  実際のアプリケーションでは実数版か整数版の片方だけを使用する事が多く、もう片方も無駄にリンクされて容量が増える事を避けるためです。
 *	* Thu Aug 10 22:36:01 JST 2017 Naoyuki Sawa
 *	- easeIn(),easeOut()をアセンブラ化しました。
 *	  コードサイズが、12バイト小さくなりました。
 *	* Fri Aug 11 23:01:44 JST 2017 Naoyuki Sawa
 *	- easeInOut(),easeOutIn()をアセンブラ化しました。
 *	  コードサイズが、64バイトも小さくなりました。
 */
#include "clip.h"
//*****************************************************************************
//	線形補間
//*****************************************************************************
//Performs a Linear intERPolation.
//[in]
//	x		The first-floating point value.							┬xとyの大小は不問です。(x<y)でも(x=y)でも(x>y)でも構いません。
//	y		The second-floating point value.						┘
//	s		A value that Linearly intERPolates between the x parameter and the y parameter.	─原則的には(0<=s<=1)ですが、(s<0)や(s>1)でも正しく動作します。
//[out]
//	return value	The result of the Linear intERPolation.
//[note]
//	- 単純な線形補間関数です。
//	  標準C関数ではありませんが、多くのゲーム用ライブラリに備わっているので、実装しておく事にしました。
//	- 一般的には'x+(s*(y-x))'で計算しますが、当実装は'x-(s*(x-y))'で計算しました。
//	  'x+(s*(y-x))'より'x-(s*(x-y))'の方が、Asm化した時にレジスタの移動回数が減って効率が良いからです。
#ifndef PIECE
double lerp(double x, double y, double s) {
    //{{式変形して最適化しました。
    //	return (x * (1.0 - s)) + (y * s);	//乗算×2,加減算×2
    //↓式変形して最適化しました。
	return x - (s * (x - y));		//乗算×1,加減算×2
    //}}式変形して最適化しました。
}
#else //PIECE
double lerp(double x, double y, double s);
asm("
		.code
		.align		1
		.global		lerp
lerp:
		pushn		%r1
		;//%r12     := xl
		;//%r13     := xh
		;//%r14     := yl
		;//%r15     := yh
		;//[%sp+ 0] := %r0
		;//[%sp+ 4] := %r1
		;//[%sp+ 8] := retp
		;//[%sp+12] := sl
		;//[%sp+16] := sh
		ld.w		%r0, %r12		;//%r0       := xl
		xcall.d		__subdf3		;//%r11:%r10 :=          (x - y)
		ld.w		%r1, %r13		;//%r1       := xh			*delay*
		xld.w		%r12, [%sp+12]		;//%r12      :=      sl
		xld.w		%r13, [%sp+16]		;//%r13      :=      sh
		ld.w		%r14, %r10		;//%r14      :=          (x - y)l
		xcall.d		__muldf3		;//%r11:%r10 :=     (s * (x - y))
		ld.w		%r15, %r11		;//%r15      :=          (x - y)h	*delay*
		ld.w		%r12, %r0		;//%r12      := xl
		ld.w		%r13, %r1		;//%r13      := xh
		popn		%r1
		ld.w		%r14, %r10		;//%r14      :=     (s * (x - y))l
		xjp.d		__subdf3		;//%r11:%r10 := x - (s * (x - y))
		ld.w		%r15, %r11		;//%r15      :=     (s * (x - y))h	*delay*
");
#endif//PIECE
//*****************************************************************************
//	イージング
//*****************************************************************************
//Performs a EASE-IN interpolation.
//[in]
//	x		The first-floating point value.								┬xとyの大小は不問です。(x<y)でも(x=y)でも(x>y)でも構いません。
//	y		The second-floating point value.							┘
//	s		A value that EASE-IN interpolates between the x parameter and the y parameter.		─(s<0)や(s>1)を指定してもエラーにはなりませんが、恐らく意図と違う結果になります。
//[out]
//	return value	The result of the EASE-IN interpolation.
//[note]
//	- ゆっくり動き出して加速する動作です。
#ifndef PIECE
double easeIn(double x, double y, double s) {
	return lerp(x, y, (s * s)/*調整可*/);
}
#else //PIECE
double easeIn(double x, double y, double s);
asm("
		.code
		.align		1
		.global		easeIn
easeIn:
		xsub		%sp, %sp, 16			;//
		xld.w		[%sp+0], %r12			;//[%sp+0]   := xl
		xld.w		[%sp+4], %r13			;//[%sp+4]   := xh
		xld.w		[%sp+8], %r14			;//[%sp+8]   := yl
		xld.w		[%sp+12], %r15			;//[%sp+12]  := yh
		xld.w		%r12, [%sp+20]			;//%r12      := sl
		xld.w		%r13, [%sp+24]			;//%r13      := sh
		ld.w		%r14, %r12			;//%r14      := sl
		xcall.d		__muldf3			;//%r11:%r10 :=      s * s
		ld.w		%r15, %r13			;//%r15      := sh			*delay*
		xld.w		%r12, [%sp+0]			;//%r12      := xl
		xld.w		%r13, [%sp+4]			;//%r13      := xh
		xld.w		%r14, [%sp+8]			;//%r14      := yl
		xld.w		%r15, [%sp+12]			;//%r15      := yh
		xld.w		[%sp+20], %r10			;//[%sp+20]  :=     (s * s).l				┬スタック渡しされた引数のメモリは破壊しても安全なので、そのメモリを次の関数呼び出しの引数を格納する用途にも使える。
		xld.w		[%sp+24], %r11			;//[%sp+24]  :=     (s * s).h				┘
		xjp.d		lerp				;//return lerp(x, y, s * s)
		xadd		%sp, %sp, 16			;//					*delay*
");
#endif//PIECE
//-----------------------------------------------------------------------------
//Performs a EASE-OUT interpolation.
//[in]
//	x		The first-floating point value.								┬xとyの大小は不問です。(x<y)でも(x=y)でも(x>y)でも構いません。
//	y		The second-floating point value.							┘
//	s		A value that EASE-OUT interpolates between the x parameter and the y parameter.		─(s<0)や(s>1)を指定してもエラーにはなりませんが、恐らく意図と違う結果になります。
//[out]
//	return value	The result of the EASE-OUT interpolation.
//[note]
//	- 高速で動き出して減速する動作です。
#ifndef PIECE
double easeOut(double x, double y, double s) {
	return easeIn(y, x, 1.0 - s);
}
#else //PIECE
double easeOut(double x, double y, double s);
asm("
		.code
		.align		1
		.global		easeOut
easeOut:
		xsub		%sp, %sp, 16
		xld.w		[%sp+0], %r12			;//[%sp+0]   := xl
		xld.w		[%sp+4], %r13			;//[%sp+4]   := xh
		xld.w		[%sp+8], %r14			;//[%sp+8]   := yl
		xld.w		[%sp+12], %r15			;//[%sp+12]  := yh
	;//	ld.w		%r12, 0				;//%r13:%r12 := 1.0
		xld.w		%r13, 0x3FF00000		;//	────────────┐
		xld.w		%r14, [%sp+20]			;//%r14      := sl		│
		xld.w		%r15, [%sp+24]			;//%r15      := sh		│
		xcall.d		__subdf3			;//%r11:%r10 :=        1.0 - s	│
		ld.w		%r12, 0				;//	←───────────┘	*delay*
		xld.w		%r12, [%sp+8]			;//%r12      := yl
		xld.w		%r13, [%sp+12]			;//%r13      := yh
		xld.w		%r14, [%sp+0]			;//%r14      := xl
		xld.w		%r15, [%sp+4]			;//%r15      := xh
		xld.w		[%sp+20], %r10			;//[%sp+20]  :=       (1.0 - s).l			┬スタック渡しされた引数のメモリは破壊しても安全なので、そのメモリを次の関数呼び出しの引数を格納する用途にも使える。
		xld.w		[%sp+24], %r11			;//[%sp+24]  :=       (1.0 - s).h			┘
		xjp.d		easeIn				;//return easeIn(y, x, 1.0 - s)
		xadd		%sp, %sp, 16			;//					*delay*
");
#endif//PIECE
//-----------------------------------------------------------------------------
//Performs a EASE-IN-OUT interpolation.
//[in]
//	x		The first-floating point value.								┬xとyの大小は不問です。(x<y)でも(x=y)でも(x>y)でも構いません。
//	y		The second-floating point value.							┘
//	s		A value that EASE-IN-OUT interpolates between the x parameter and the y parameter.	─(s<0)や(s>1)を指定してもエラーにはなりませんが、恐らく意図と違う結果になります。
//[out]
//	return value	The result of the EASE-IN-OUT interpolation.
//[note]
//	- ゆっくり動き出して加速し、また減速する動作です。
//{{2017/08/08削除:easeInOut()とeaseOutIn()の処理を共通化して、コードサイズを削減しました。
//double easeInOut(double x, double y, double s) {
//	double mid = (x + y) / 2.0;
//	if(s >= 0.5) {
//		x = y;
//		s = 1.0 - s;
//	}
//	return easeIn(x, mid, s * 2.0);										//easeInOut()とeaseOutIn()の違いは、この行だけです。
//}
//}}2017/08/08削除:easeInOut()とeaseOutIn()の処理を共通化して、コードサイズを削減しました。
//-----------------------------------------------------------------------------
//Performs a EASE-OUT-IN interpolation.
//[in]
//	x		The first-floating point value.								┬xとyの大小は不問です。(x<y)でも(x=y)でも(x>y)でも構いません。
//	y		The second-floating point value.							┘
//	s		A value that EASE-OUT-IN interpolates between the x parameter and the y parameter.	─(s<0)や(s>1)を指定してもエラーにはなりませんが、恐らく意図と違う結果になります。
//[out]
//	return value	The result of the EASE-OUT-IN interpolation.
//[note]
//	- 高速で動き出して減速し、また加速する動作です。
//{{2017/08/08削除:easeInOut()とeaseOutIn()の処理を共通化して、コードサイズを削減しました。
//double easeOutIn(double x, double y, double s) {
//	double mid = (x + y) / 2.0;
//	if(s >= 0.5) {
//		x = y;
//		s = 1.0 - s;
//	}
//	return easeOut(x, mid, s * 2.0);									//easeInOut()とeaseOutIn()の違いは、この行だけです。
//}
//}}2017/08/08削除:easeInOut()とeaseOutIn()の処理を共通化して、コードサイズを削減しました。
//-----------------------------------------------------------------------------
//{{2017/08/08追加:easeInOut()とeaseOutIn()の処理を共通化して、コードサイズを削減しました。
#ifndef PIECE
static double easeInOut_easeOutIn(double x, double y, double s, double (*fn)(double x, double y, double s));	//前方参照
double easeInOut(double x, double y, double s) { return easeInOut_easeOutIn(x, y, s, easeIn); }
double easeOutIn(double x, double y, double s) { return easeInOut_easeOutIn(x, y, s, easeOut); }
static double easeInOut_easeOutIn(double x, double y, double s, double (*fn)(double x, double y, double s)) {
#if 0
//□素直な実装。
	double mid = (x + y) / 2.0;
	if(s >= 0.5) {
		x = y;
		s = 1.0 - s;
	}
	return (*fn)(x, mid, s * 2.0);
#else//↑↓どちらのVerでも同じ結果になります。
//□アセンブラ化し易いように式を変形した実装。
	double mid = (x + y) / 2.0;
	   s *= 2.0;
	if(s >= 1.0) {
		s = 2.0 - s;
		x = y;
	}
	return (*fn)(x, mid, s);
#endif
}
#else //PIECE
double easeInOut(double x, double y, double s);
double easeOutIn(double x, double y, double s);
asm("
		.code
		.align		1
		.global		easeInOut
		.global		easeOutIn
easeInOut:
		xld.w		%r9, easeIn
		jp		easeInOut_easeOutIn
easeOutIn:
		xld.w		%r9, easeOut
easeInOut_easeOutIn:
		xsub		%sp, %sp, 28
		;//[%sp+0]  := (fn)
		;//[%sp+4]  := (x.l)
		;//[%sp+8]  := (x.h)
		;//[%sp+12] := (y.l)
		;//[%sp+16] := (y.h)
		;//[%sp+20] := (mid.l)
		;//[%sp+24] := (mid.l)
		;//[%sp+28] := retp
		;//[%sp+32] := s.l
		;//[%sp+36] := s.h
		xld.w		[%sp+0], %r9			;//[%sp+0]   := fn
		xld.w		[%sp+4], %r12			;//[%sp+4]   := x.l
		xld.w		[%sp+8], %r13			;//[%sp+8]   := x.h
		xld.w		[%sp+12], %r14			;//[%sp+12]  := y.l
		xld.w		[%sp+16], %r15			;//[%sp+16]  := y.h
		xcall		__adddf3			;//%r11:%r10 :=       (x + y)							[mod] %r4,%r5,%r6,%r7,%r12,%r13,%r14,%r15
		ld.w		%r12, %r10			;//%r12      :=       (x + y).l
		ld.w		%r13, %r11			;//%r13      :=       (x + y).h
	;//	ld.w		%r14, 0				;//%r14      :=                 (2.0).l		─┐
	;//{{最適化:3命令⇒2命令										　│
	;//	xld.w		%r15, 0x40000000		;//%r15      :=                 (2.0).h		　│
	;//↓最適化:3命令⇒2命令										　│
		ld.w		%r15, 1				;//%r15      :=                 (2.0).h		　│
		rr		%r15, 2				;//						　│
	;//}}最適化:3命令⇒2命令										　│
		xcall.d		__divdf3			;//%r11:%r10 := mid = (x + y) /  2.0		　│				[mod] %r4,%r5,%r12,%r13,%r14,%r15
		ld.w		%r14, 0				;//%r14      :=                 (2.0).l		←┘		*delay*
		xld.w		[%sp+20], %r10			;//[%sp+20]  := mid.l
		xld.w		[%sp+24], %r11			;//[%sp+24]  := mid.h
		xld.w		%r12, [%sp+32]			;//%r12      := s.l
		xld.w		%r13, [%sp+36]			;//%r13      := s.h
	;//	ld.w		%r14, 0				;//%r14      :=     (2.0).l			─┐
	;//{{最適化:3命令⇒2命令										　│
	;//	xld.w		%r15, 0x40000000		;//%r15      :=     (2.0).h			　│
	;//↓最適化:3命令⇒2命令										　│
		ld.w		%r15, 1				;//%r15      :=     (2.0).h			　│
		rr		%r15, 2				;//						　│
	;//}}最適化:3命令⇒2命令										　│
		xcall.d		__muldf3			;//%r11:%r10 := s *= 2.0			　│─┐─┐			[mod] %r4,%r5,%r12,%r13,%r14,%r15
		ld.w		%r14, 0				;//%r14      :=     (2.0).l			←┘　│　│	*delay*
		ld.w		%r12, %r10			;//%r12      := s.l				　　←┤　│
		ld.w		%r13, %r11			;//%r13      := s.h				　　←┘　│
	;//	ld.w		%r14, 0				;//%r14      :=     (1.0).l			─┐　│　│
		xld.w		%r15, 0x3FF00000		;//%r15      :=     (1.0).h			　│　│　│
		xcall.d		__fcmpd				;//%psr      := s -  1.0			　│　│　│			[mod] %r12,%r13,%r14,%r15	__fcmpdが%r10,%r11を破壊しない事に依存しています。
		ld.w		%r14, 0				;//%r14      :=     (1.0).l			←┘　│　│	*delay*
		xld.w		%r12, [%sp+4]			;//%r12      := x.l				　　　│　│
		xld.w		%r13, [%sp+8]			;//%r13      := x.h				　　　│　│
		jrlt		easeInOut_easeOutIn_L10		;//if(          s >= 1.0) {			　　　│　│
		ld.w		%r12, 0				;//  %r12      :=    (2.0).l			　　　│　│
	;//{{最適化:3命令⇒2命令										　　　│　│
	;//	xld.w		%r13, 0x40000000		;//  %r13      :=    (2.0).h			　　　│　│
	;//↓最適化:3命令⇒2命令										　　　│　│
		ld.w		%r13, 1				;//  %r13      :=    (2.0).h			　　　│　│
		rr		%r13, 2				;//						　　　│　│
	;//}}最適化:3命令⇒2命令										　　　│　│
		ld.w		%r14, %r10			;//  %r14      :=           s.l			　　←┤　│
	;//	ld.w		%r15, %r11			;//  %r15      :=           s.h			─┐←┘　│
		xcall.d		__subdf3			;//  %r11:%r10 := s = 2.0 - s			　│─┐　│			[mod] %r4,%r5,%r6,%r7,%r12,%r13,%r14,%r15
		ld.w		%r15, %r11			;//  %r15      :=           s.h			←┘　│　│	*delay*
		xld.w		%r12, [%sp+12]			;//  %r12      := x.l = y.l			　　　│　│
		xld.w		%r13, [%sp+16]			;//  %r13      := x.h = y.h			　　　│　│
easeInOut_easeOutIn_L10:					;//}						　　　│　│
		xld.w		%r14, [%sp+20]			;//%r14      :=      mid.l			　　　│　│
		xld.w		%r15, [%sp+24]			;//%r15      :=      mid.h			　　　│　│
		xld.w		[%sp+32], %r10			;//[%sp+32]  :=           s.l			　　←┤←┤			┬スタック渡しされた引数のメモリは破壊しても安全なので、そのメモリを次の関数呼び出しの引数を格納する用途にも使える。
		xld.w		[%sp+36], %r11			;//[%sp+36]  :=           s.h			　　←┘←┘			┘
		;//[%sp+0]  := fn	←──①ret.dによってfnが読み出され…
		;//[%sp+4]  := 不要	┐←─②…%spはここを指す。
		;//[%sp+8]  := 不要	├──③xadd %sp,%sp,24によってこの範囲を破棄する。
		;//[%sp+12] := 不要	│
		;//[%sp+16] := 不要	│
		;//[%sp+20] := 不要	│
		;//[%sp+24] := 不要	┘
		;//[%sp+28] := retp	─┐
		;//[%sp+32] := s.l	┬│─④分岐先のfnにて引数として参照される。
		;//[%sp+36] := s.h	┘└─⑤分岐先のfnからの戻り先となる。
		ret.d						;//return (*fn)(  x, mid, s)							┬この二行はちょっと高等な(^^)テクニックを使っています。詳細は、この上のコメント(①～⑤)を参照して下さい。
		xadd		%sp, %sp, 24			;//								*delay*		┘
");
#endif//PIECE
//}}2017/08/08追加:easeInOut()とeaseOutIn()の処理を共通化して、コードサイズを削減しました。
//*****************************************************************************
//	連立方程式
//*****************************************************************************
//m[n行][n列]の行列式(determinant)を求めます。
//[note]
//	- 既存のmatnf_determinant()を、単制度⇒倍制度に変更した関数です。
//	  使い方やアルゴリズムは、matnf_determinant()と同じです。
//	  詳細は、clipflt1.cのmatnf_determinant()の実装を参照して下さい。
//	- 関数名を短くしても衝突しなさそうだったので、思い切って短くしました。
//	  もし今後関数名が衝突した場合は、変更して下さい。
double determinant(const double* m/*[n][n]*/, int n) {
	int i;
	double d;
	if(n < 1) { DIE(); }
	//余因子展開を用いて、行列式(determinant)を求めます。
	d = 0.0;
	for(i = 0; i < n; i++) {
		d += (m[(n * i) + 0/*j*/] * cofactor(m, n, i, 0/*j*/));
	}
	return d;
}
//-----------------------------------------------------------------------------
//m[n行][n列]からi行j列を除いた、余因子(cofactor)を求めます。
//[note]
//	- 既存のmatnf_cofactor()を、単制度⇒倍制度に変更した関数です。
//	  使い方やアルゴリズムは、matnf_cofactor()と同じです。
//	  詳細は、clipflt1.cのmatnf_cofactor()の実装を参照して下さい。
//	- 関数名を短くしても衝突しなさそうだったので、思い切って短くしました。
//	  もし今後関数名が衝突した場合は、変更して下さい。
double cofactor(const double* m/*[n][n]*/, int n, int i, int j) {
	int k, l;
	double d;
	double* s/*[n-1][n-1]*/;
	double* 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列になる場合
#if     !(defined(_WIN32)&&defined(_MT))
//Win32シングルスレッドの場合,P/ECEの場合,その他の環境の場合に有効なブロック
	s = malloc(sizeof(double) * ((n - 1) * (n - 1)));
	if(!s) { DIE(); }
#else //!(defined(_WIN32)&&defined(_MT))
//Win32マルチスレッドの場合に有効なブロック
	s = alloca(sizeof(double) * ((n - 1) * (n - 1)));
#endif//!(defined(_WIN32)&&defined(_MT))
	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 = determinant(s, (n - 1));
	//余因子(cofactor)を求めます。
	if((i ^ j) & 1) { d = -d; }
	//小行列を解放します。
#if     !(defined(_WIN32)&&defined(_MT))
//Win32シングルスレッドの場合,P/ECEの場合,その他の環境の場合に有効なブロック
	free(s);
#endif//!(defined(_WIN32)&&defined(_MT))
	return d;
}
//-----------------------------------------------------------------------------
//n次連立方程式の解を求めます。
//[note]
//	- 既存のsolve_simultaneous_equations()を、単制度⇒倍制度に変更した関数です。
//	  使い方やアルゴリズムは、solve_simultaneous_equations()と同じです。
//	  詳細は、clipflt1.cのsolve_simultaneous_equations()の実装を参照して下さい。
//	- ただし、solve_simultaneous_equations()と較べて、一つだけ変更した点が有ります。
//	  solve_simultaneous_equations()は戻り値無しでしたが、当関数は処理の途中で求めた行列式の値を返すようにしました。
//	  連立方程式の解が有効かを行列式の値を見て判断するために、別途determinant()を呼び出すよりも効率が良いからです。
//	  solve_simultaneous_equations()では、有効な解が無い連立方程式は指定されない前提で戻り値無しにしていましたが、
//	  実際には、有効な解が無い連立方程式も指定される可能性が有ると考えて上記の戻り値の仕様変更を行う事にしました。
//	  あくまでも処理の途中で求めた行列式の値を返すようにしただけの変更であり、処理内容自体は一切変更していません。
//	  また、当関数が返した値を参照するかどうかも、呼び出し側次第です。
//	  有効な解が無い連立方程式を指定しない前提ならば、戻り値を無視しても構いません。
//	- 関数名を短くしても衝突しなさそうだったので、思い切って短くしました。
//	  もし今後関数名が衝突した場合は、変更して下さい。
//	  ちなみにMATLABにも、連立方程式を解く関数としてsolve()という名前の関数があるそうです。
//	  用途は似ていますが、引数仕様等は互換ではありません。たまたま名前が一致しただけです。
#if     !(defined(_WIN32)&&defined(_MT))
//Win32シングルスレッドの場合,P/ECEの場合,その他の環境の場合に有効なブロック
double solve(const double* A/*[n][n]*/, const double* b/*[n]*/, double* x/*[n]*/, int n) {
	int i, j;
	double d;
	//行列式(determinant)を求めます。
	d = determinant(A, n);
	//クラメルの公式を用いて、Aの逆行列を求めながら、bを掛けて、xを求めます。
	for(i = 0; i < n; i++) {
		x[i] = 0.0;
		for(j = 0; j < n; j++) {
			x[i] += (cofactor(A, n, j, i) / d) * b[j];
			//            順序注意!!~~~~
		}
	}
	//行列式を返します。
	return d;
}
#else //!(defined(_WIN32)&&defined(_MT))
//Win32マルチスレッドの場合に有効なブロック
struct _ST_solve {
	const double*		A/*[n][n]*/;	//┬solve()の引数と同じです。solve_start_routine()からは参照専用です。書き換えてはいけません。
	const double*		b/*[n]*/;	//│
	double*			x/*[n]*/;	//│
	int			n;		//┘
	double			d;		//─solve()が値を格納します。solve_start_routine()からは参照専用です。書き換えてはいけません。
	volatile int		i;		//┬i＝0→n。solve_start_routine()が読み書きします。iを読み書きする時はcsをロックして下さい。
	CRITICAL_SECTION 	cs;		//┘
};
static unsigned __stdcall solve_start_routine(void* _arg) {
	struct _ST_solve* arg = _arg;
	const double* const A = arg->A;
	const double* const b = arg->b;
	double*       const x = arg->x;
	int           const n = arg->n;
	double        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.0;
		for(j = 0; j < n; j++) {
			x[i] += (cofactor(A, n, j, i) / d) * b[j];
			//            順序注意!!~~~~
		}
	}
	return 0;
}
double solve(const double* A/*[n][n]*/, const double* b/*[n]*/, double* x/*[n]*/, int n) {
	struct _ST_solve arg;
	HANDLE hThread[16/*調整可*/];	//当関数が作成するスレッドの数はこの配列サイズが上限となる。2016/03/19現在は16スレッドを上限とした。もし今後16コア以上のCPUが一般的になった場合はこの配列サイズを増やせ。
	int nThread, iThread;
	SYSTEM_INFO si;
	//構造体を初期化する。
	arg.A = A;
	arg.b = b;
	arg.x = x;
	arg.n = n;
	arg.d = 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_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);	//忘れないで!!
	//行列式を返します。
	return arg.d;
}
#endif//!(defined(_WIN32)&&defined(_MT))
