/*	
 *	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バイトも小さくなりました。
 */
#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]
//	- ゆっくり動き出して加速する動作です。
double easeIn(double x, double y, double s) {
	return lerp(x, y, (s * s)/*調整可*/);
}
//-----------------------------------------------------------------------------
//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]
//	- 高速で動き出して減速する動作です。
double easeOut(double x, double y, double s) {
	return easeIn(y, x, 1.0 - s);
}
//-----------------------------------------------------------------------------
//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()の処理を共通化して、コードサイズを削減しました。
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)) {
	double mid = (x + y) / 2.0;
	if(s >= 0.5) {
		x = y;
		s = 1.0 - s;
	}
	return (*fn)(x, mid, s * 2.0);
}
//}}2017/08/08追加:easeInOut()とeaseOutIn()の処理を共通化して、コードサイズを削減しました。
//*****************************************************************************
//	線形補間(整数版)
//*****************************************************************************
//x～yを、n/dで補間します。
//[in]
//	x		開始値。		INT_MIN～INT_MAX
//	y		終了値。		INT_MIN～INT_MAX
//	n		補間係数の分子。	INT_MIN～INT_MAX
//	d		補間係数の分母。	1～INT_MAX,-1～INT_MIN	0のみ不可
//[out]
//	戻り値		x+(y-x)*(n/d)		結果がINT_MIN～INT_MAXの範囲に収まらない場合は、戻り値は不定となる。<例>{x=0,y=1<<30,n=5,d=1}を指定した場合など。
//[note]
//	- 当関数は、他の補完処理からも利用する基本処理になると思ったので、汎用性を重視して、引数の値の範囲の制限を設けないように作成しました。
//	  内部的には演算途中で64ビット演算を行っているので、(結果以外は)演算途中での桁あふれを心配する必要は有りません。
//	  実数版の'lerp()'と、ほぼ同じ感覚で使えると思います。
//	  詳細は、以下の通りです。
//	- xとyの大小は不問です。
//	  (x<y)でも(x=y)でも(x>y)でも構いません。
//	- xとyはそれぞれ、INT_MIN～INT_MAXの任意の値です。
//	  (|x-y|>INT_MAX)でも構いません。				一般的には(|x-y|≦INT_MAX)だとは思いますが、汎用性のために(|x-y|>INT_MAX)も対応しました。この対応のために、やや処理が多くなっています。
//	- nは一般的には、0以上d以下(d<0の場合は0以下d以上)の値です。
//	  ただし必須ではなく、上記の範囲外の値でも構いません。
//	- dは一般的には、1～INT_MAXの値です。
//	  ただし必須ではなく、-1～INT_MINの値でも構いません。		一般的には(1≦d≦INT_MAX)だとは思いますが、汎用性のために(-1≧d≧INT_MIN)も対応しました。この対応のために、やや処理が多くなっています。
#ifndef PIECE
int lerp_i(int x, int y, int n, int d) {
	if(d < 0) {
		n = -n;
		d = -d;
	}
	if(y < x) {
		y -= x;
		y = -y;
	    //{{共通処理
		if(n < 0) {
			n = -n;
			y = (uint64_t)(unsigned)y * (uint64_t)(unsigned)n / (uint64_t)(unsigned)d;
			y = -y;	//    ~~~~~~~~~~              ~~~~~~~~~~              ~~~~~~~~~~	元々の変数y,n,dは符号付きの値ですが、既に符号無しの値に補正済みなので、bit31をbit63～32へ符号拡張しないように、まず(unsigned)へのキャストが必要です。
		} else {
			y = (uint64_t)(unsigned)y * (uint64_t)(unsigned)n / (uint64_t)(unsigned)d;
		}		//    ~~~~~~~~~~              ~~~~~~~~~~              ~~~~~~~~~~	元々の変数y,n,dは符号付きの値ですが、既に符号無しの値に補正済みなので、bit31をbit63～32へ符号拡張しないように、まず(unsigned)へのキャストが必要です。
	    //}}共通処理
		y = -y;
	} else {
		y -= x;
	    //{{共通処理
		if(n < 0) {
			n = -n;
			y = (uint64_t)(unsigned)y * (uint64_t)(unsigned)n / (uint64_t)(unsigned)d;
			y = -y;	//    ~~~~~~~~~~              ~~~~~~~~~~              ~~~~~~~~~~	元々の変数y,n,dは符号付きの値ですが、既に符号無しの値に補正済みなので、bit31をbit63～32へ符号拡張しないように、まず(unsigned)へのキャストが必要です。
		} else {
			y = (uint64_t)(unsigned)y * (uint64_t)(unsigned)n / (uint64_t)(unsigned)d;
		}		//    ~~~~~~~~~~              ~~~~~~~~~~              ~~~~~~~~~~	元々の変数y,n,dは符号付きの値ですが、既に符号無しの値に補正済みなので、bit31をbit63～32へ符号拡張しないように、まず(unsigned)へのキャストが必要です。
	    //}}共通処理
	}
	return x + y;
}
#else //PIECE
//□アセンブラ版のlerp_i()が保存するレジスタについて:
//│アセンブラ版のlerp_i()は、%r9,%r10,%r11,%r13,%r14,%r15,%alr,%ahr,%psr だけを破壊します。
//│呼び出し規約によって保存が義務付けられている %r0～%r3,及び,%r8 に加えて、%r4～%r7,及び,%r12 も保存します。	←──────────────────────┬注目!!
//│呼び出し側の関数は、上記のレジスタが保存される事を前提として、上記のレジスタの退避を省略しても構いません。	←──────────────────────┘
//□アセンブラ版の性能について(2017/07/28現在)
//│C言語版よりもアセンブラ版の方が、64バイト小さく、約10倍も高速です。	←──────────────────────────────────────────┬注目!!
//│C言語版でframdi.oを高速RAMに置いた場合と比べても、約3倍も高速です。	←──────────────────────────────────────────┘
//│C言語版とアセンブラ版の速度が大きく違う理由は、C言語版では、必要無い部分までフル64ビット乗除算を行ってしまっていますが、
//│アセンブラ版では、必要最低限の32ビット×32ビット⇒64ビット乗算，及び，64ビット÷32ビット⇒32ビット除算に最適化したので、
//│C言語版と違って、ソフトウェアによる64ビット乗除算ルーチン(framdi.o)は呼び出さず、機械語命令だけで済ませられたからです。
//│速度計測の結果の詳細は、以下の通りです。
//│【速度計測に使用したテストプログラム】
//│#define N 10000
//│int data[N][4];
//│int app_main(int argc, char* argv[]) {
//│  int i, j, t;
//│  srand48(0);
//│  for(i = 0; i < N; i++) {
//│    double d;
//│         data[i][0] = mrand48();
//│         data[i][1] = mrand48();
//│         data[i][2] = mrand48();
//│    do { data[i][3] = mrand48(); } while(!data[i][3]);
//│    d = lerp(data[i][0], data[i][1], (double)data[i][2] / (double)data[i][3]);
//│    if((d < INT_MIN) || (d > INT_MAX)) { i--; }
//│  }
//│  t = pceTimerGetCount();
//│  for(j = 0; j < 10; j++) {
//│    for(i = 0; i < N; i++) {
//│      lerp_i(data[i][0], data[i][1], data[i][2], data[i][3]);
//│    }
//│  }
//│  t = pceTimerGetCount() - t;
//│  die("%d", t);
//│  return 0;
//│}
//│【速度計測の結果】TURBOモード、スタックは高速RAMに配置して実行しました。
//│・C言語版(framdi.oをSRAMに配置)     7555ミリ秒
//│・C言語版(framdi.oを高速RAMに配置)  2272ミリ秒
//│・アセンブラ版(framdi.oは関係無し)   767ミリ秒
int lerp_i(int x, int y, int n, int d);
asm("
		.code
		.align		1
		.global		lerp_i
lerp_i:
		;//%r12 := x
		;//%r13 := y
		;//%r14 := n
		;//%r15 := d
		cmp		%r15, 0			;//if(d < 0) {
		jrge.d		lerp_i_L10		;//  
		ld.w		%r9, 8-1		;//%r9       := i = 8-1		*delay*	─┐
		not		%r14, %r14		;//  %r14    := n = -n			　│
		add		%r14, 1			;//  					　│
		not		%r15, %r15		;//  %r15    := d = -d			　│
		add		%r15, 1			;//  					　│
lerp_i_L10:						;//}					　│
		sub		%r13, %r12		;//%r13      := y -= x			　│
		jrge.d		lerp_i_L20		;//if(変更前のy < x) {			　│
		ld.w		%r11, %psr		;//%psr保存①			*delay*	　│
		not		%r13, %r13		;//  %r13    := y = -y			　│
		add		%r13, 1			;//  					　│
lerp_i_L20:						;//}					　│
		cmp		%r14, -1		;//if(変更前のn <= -1) {		　│		(n>-1)ならば%psr(Z)=0で分岐する。(n<=-1)ならば%psr(Z)=1の可能性も有るが、この下の処理で%psr(Z)=0になる。
		jrgt.d		lerp_i_L30		;//  					　│
		ld.w		%r10, %psr		;//%psr保存②			*delay*	　│
		sub		%r14, 1			;//  %r14    := n = -n			　│		この処理で必ず%psr(Z)=0になる。-1して0になるのは(n=1)だった場合のみだが、(n=1)ならばここへは来ないから。
		not		%r14, %r14		;//					　│
lerp_i_L30:						;//}					　│	─┐	分岐した場合も分岐しなかった場合も、ここで必ず%psr(Z)=0になっている。この後のループの初回の%psr(Z)判定まで変化しないよう注意せよ。
		mltu.w		%r13, %r14		;//%ahr:%alr := (y * n)			　│	　│
		ld.w		%r14, %ahr		;//%r14      := (y * n)[63:32]		　│	　│
		div0u		%r15			;//%ahr      :=    0			　│	　│	%psr(DS,N)=0
		ld.w		%ahr, %r14		;//%ahr:%alr := (y * n)			　│	　│
	;//	ld.w		%r9, 8-1		;//%r9       := i = 8-1			←┘	　│
lerp_i_L40:						;//do {						　│
		div1		%r15			;//  1						　│	┐
		div1		%r15			;//  2						　│	├div1を実行する時は'%psr(DS,N)=0'を維持している必要がある事に注意せよ。
		div1		%r15			;//  3						　│	│
		div1		%r15			;//  4						　│	┘
		jrne.d		lerp_i_L40		;//} while(i != 0)				←┤	初回は%psr(Z)=0であるはずなので必ずループする。二回目以降はこの下で変化した%psr(Z)を参照してループする。
		sub		%r9, 1			;//  %r9     := i--		*delay*		─┘
		ld.w		%psr, %r10		;//%psr復元②
		jrgt.d		lerp_i_L50		;//if(変更前のn <= -1) {
		ld.w		%r10, %alr		;//%r10      := y = (y * n) / d	*delay*
		not		%r10, %r10		;//  %r10    := y = -y
		add		%r10, 1			;//  
lerp_i_L50:						;//}
		ld.w		%psr, %r11		;//%psr復元①
		jrge		lerp_i_L60		;//if(変更前のn < 0) {
		not		%r10, %r10		;//  %r10    := y = -y
		add		%r10, 1			;//  
lerp_i_L60:						;//}
		ret.d					;//return       x + y
		add		%r10, %r12		;//%r10      := x + y		*delay*
");
#endif//PIECE
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#if 0
//□lerp_i()のテストスイート
int app_main(int argc, char* argv[]) {
	int x, y, n, d, a;
	double b;
	for(;;) {
		     x = mrand48();
		     y = mrand48();
		     n = mrand48();
		do { d = mrand48(); } while(!d);
		a = lerp_i(x, y,         n,          d);
		b = lerp(  x, y, (double)n / (double)d);
	#ifndef PIECE
		       printf(                     "% 11d % 11d % 2d\n", a, lround(b), a - lround(b));
	#else //PIECE
		surface_clear(&surface, 0);
		render_printf(&render, 0, 0, 2, 3, "% 11d % 11d % 2d\n", a, lround(b), a - lround(b));
		update();
		if(pcePadGetDirect() & PAD_SELECT) { break; }
	#endif//PIECE
		if((b >= INT_MIN) && (b <= INT_MAX)) {	//結果がINT_MIN～INT_MAXの範囲に収まらない場合を除けば…
			if(abs(a - lround(b)) > 1) {	//小数部の切り捨て誤差の範囲で正しい結果が得られるはず。
				die("x = % 11d\n"
				    "y = % 11d\n"
				    "n = % 11d\n"
				    "d = % 11d\n"
				    "a = % 11d\n"
				    "b = % 11d\n", x, y, n, d, a, lround(b));
			}
		}
	}
	return 0;
}
#endif
//*****************************************************************************
//	イージング(整数版)
//*****************************************************************************
//ゆっくり動き出して加速する動作です。
//[in]
//	x		開始値。		INT_MIN～INT_MAX
//	y		終了値。		INT_MIN～INT_MAX
//	n		補間係数の分子。	0～d
//	d		補間係数の分母。	1～SHRT_MAX
//[out]
//	戻り値		補間結果。
//[note]
//	- lerp_i()は基本処理なので汎用性を重視して引数の値の範囲の制限を設けないようにしましたが、当関数(easeIn_i())は末端処理だと思うので汎用性は重視しませんでした。
//	  そのため、lerp_i()よりも引数d(ひいては引数nも)の値の範囲の制限が厳しいです。
//	  一般的な使い方では充分だと思いますが、引数の値が5桁(=?万)以上になるケースでは、注意して下さい。
#ifndef PIECE
int easeIn_i(int x, int y, int n, int d) {
	n *= n;	//調整可
	d *= d;	//調整可
	return lerp_i(x, y, n, d);
}
#else //PIECE
int easeIn_i(int x, int y, int n, int d);
asm("
		.code
		.align		1
		.global		easeIn_i
easeIn_i:
		mlt.h		%r14, %r14		;//%r14 := n *= n
		ld.w		%r14, %alr		;//
		mlt.h		%r15, %r15		;//%r15 := d *= d
		xjp.d		lerp_i			;//return lerp_i(x, y, n, d)
		ld.w		%r15, %alr		;//				*delay*
");
#endif//PIECE
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#if 0
//□easeIn_i()のテストスイート
int app_main(int argc, char* argv[]) {
	int x, y, n, d, a;
	double b;
	for(;;) {
		     x = mrand48();
		     y = mrand48();
		do { d = lrand48() % (SHRT_MAX + 1); } while(!d);
		     n = lrand48() % d;
		a = easeIn_i(x, y,         n,          d);
		b = easeIn(  x, y, (double)n / (double)d);
	#ifndef PIECE
		       printf(                     "% 11d % 11d % 2d\n", a, lround(b), a - lround(b));
	#else //PIECE
		surface_clear(&surface, 0);
		render_printf(&render, 0, 0, 2, 3, "% 11d % 11d % 2d\n", a, lround(b), a - lround(b));
		update();
		if(pcePadGetDirect() & PAD_SELECT) { break; }
	#endif//PIECE
		if((b >= INT_MIN) && (b <= INT_MAX)) {	//結果がINT_MIN～INT_MAXの範囲に収まらない場合を除けば…
			if(abs(a - lround(b)) > 1) {	//小数部の切り捨て誤差の範囲で正しい結果が得られるはず。
				die("x = % 11d\n"
				    "y = % 11d\n"
				    "n = % 11d\n"
				    "d = % 11d\n"
				    "a = % 11d\n"
				    "b = % 11d\n", x, y, n, d, a, lround(b));
			}
		}
	}
	return 0;
}
#endif
//-----------------------------------------------------------------------------
//高速で動き出して減速する動作です。
//[in]
//	x		開始値。		INT_MIN～INT_MAX
//	y		終了値。		INT_MIN～INT_MAX
//	n		補間係数の分子。	0～d
//	d		補間係数の分母。	1～SHRT_MAX
//[out]
//	戻り値		補間結果。
//[note]
//	- lerp_i()は基本処理なので汎用性を重視して引数の値の範囲の制限を設けないようにしましたが、当関数(easeOut_i())は末端処理だと思うので汎用性は重視しませんでした。
//	  そのため、lerp_i()よりも引数d(ひいては引数nも)の値の範囲の制限が厳しいです。
//	  一般的な使い方では充分だと思いますが、引数の値が5桁(=?万)以上になるケースでは、注意して下さい。
#ifndef PIECE
int easeOut_i(int x, int y, int n, int d) {
	n = -n;
	n += d;
	return easeIn_i(y, x, n, d);
}
#else //PIECE
int easeIn_i(int x, int y, int n, int d);
asm("
		.code
		.align		1
		.global		easeOut_i
easeOut_i:
		not		%r14, %r14		;//%r14 := n = -n
		add		%r14, 1			;//
		add		%r14, %r15		;//%r14 := n += d
		ld.w		%r9, %r12		;//%r9  :=            x
		ld.w		%r12, %r13		;//%r12 :=         y
		xjp.d		easeIn_i		;//return easeIn_i(y, x, n, d)
		ld.w		%r13, %r9		;//%r13 :=            x		*delay*
");
#endif//PIECE
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#if 0
//□easeOut_i()のテストスイート
int app_main(int argc, char* argv[]) {
	int x, y, n, d, a;
	double b;
	for(;;) {
		     x = mrand48();
		     y = mrand48();
		do { d = lrand48() % (SHRT_MAX + 1); } while(!d);
		     n = lrand48() % d;
		a = easeOut_i(x, y,         n,          d);
		b = easeOut(  x, y, (double)n / (double)d);
	#ifndef PIECE
		       printf(                     "% 11d % 11d % 2d\n", a, lround(b), a - lround(b));
	#else //PIECE
		surface_clear(&surface, 0);
		render_printf(&render, 0, 0, 2, 3, "% 11d % 11d % 2d\n", a, lround(b), a - lround(b));
		update();
		if(pcePadGetDirect() & PAD_SELECT) { break; }
	#endif//PIECE
		if((b >= INT_MIN) && (b <= INT_MAX)) {	//結果がINT_MIN～INT_MAXの範囲に収まらない場合を除けば…
			if(abs(a - lround(b)) > 1) {	//小数部の切り捨て誤差の範囲で正しい結果が得られるはず。
				die("x = % 11d\n"
				    "y = % 11d\n"
				    "n = % 11d\n"
				    "d = % 11d\n"
				    "a = % 11d\n"
				    "b = % 11d\n", x, y, n, d, a, lround(b));
			}
		}
	}
	return 0;
}
#endif
//-----------------------------------------------------------------------------
//ゆっくり動き出して加速し、また減速する動作です。
//[in]
//	x		開始値。		INT_MIN～INT_MAX	┬処理の都合上、演算途中で、(y + x)がINT_MIN～INT_MAXの範囲を超えた場合も、結果が不正になります。
//	y		終了値。		INT_MIN～INT_MAX	┘
//	n		補間係数の分子。	0～d
//	d		補間係数の分母。	1～SHRT_MAX
//[out]
//	戻り値		補間結果。
//[note]
//	- lerp_i()は基本処理なので汎用性を重視して引数の値の範囲の制限を設けないようにしましたが、当関数(easeInOut_i())は末端処理だと思うので汎用性は重視しませんでした。
//	  そのため、lerp_i()よりも引数d(ひいては引数nも)の値の範囲の制限が厳しいです。
//	  一般的な使い方では充分だと思いますが、引数の値が5桁(=?万)以上になるケースでは、注意して下さい。
//{{2017/08/08削除:easeInOut_i()とeaseOutIn_i()の処理を共通化して、コードサイズを削減しました。
//#ifndef PIECE
//int easeInOut_i(int x, int y, int n, int d) {
//	y += x;								//処理の都合上、演算途中で、(y + x)がINT_MIN～INT_MAXの範囲を超えた場合も、結果が不正になります。
//	n <<= 1;
//	if(n >= d) {
//		x = -x;
//		x += y;
//		n = -n;
//		n += d << 1;
//	}
//	y >>= 1;
//	return easeIn_i(x, y, n, d);										//easeInOut_i()とeaseOutIn_i()の違いは、この行だけです。
//}
//#else //PIECE
//int easeInOut_i(int x, int y, int n, int d);
//asm("
//		.code
//		.align		1
//		.global		easeInOut_i
//easeInOut_i:
//		sla		%r14, 1			;//%r14 := n <<= 1
//		cmp		%r14, %r15		;//if(n >= d) {
//		jrlt.d		easeInOut_i_L10		;//  
//		add		%r13, %r12		;//%r13 := y += x		*delay*
//		not		%r12, %r12		;//  %r12 := x = -x
//		add		%r12, 1			;//  
//		add		%r12, %r13		;//  %r12 := x += y
//		not		%r14, %r14		;//  %r14 := n = -n
//		add		%r14, 1			;//  
//		add		%r14, %r15		;//  %r14 := n += d << 1
//		add		%r14, %r15		;//  
//easeInOut_i_L10:					;//}
//		xjp.d		easeIn_i		;//return easeIn_i(x, y, n, d)				//easeInOut_i()とeaseOutIn_i()の違いは、この行だけです。
//		sra		%r13, 1			;//%r13 := y >>= 1		*delay*
//");
//#endif//PIECE
//}}2017/08/08削除:easeInOut_i()とeaseOutIn_i()の処理を共通化して、コードサイズを削減しました。
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#if 0
//□easeInOut_i()のテストスイート
int app_main(int argc, char* argv[]) {
	int x, y, n, d, a;
	double b;
	for(;;) {
		     x = mrand48();
		     y = mrand48();
		do { d = lrand48() % (SHRT_MAX + 1); } while(!d);
		     n = lrand48() % d;
		if((((double)x + (double)y) < (double)INT_MIN) || (((double)x + (double)y) > (double)INT_MAX)) { continue; }
		a = easeInOut_i(x, y,         n,          d);
		b = easeInOut(  x, y, (double)n / (double)d);
	#ifndef PIECE
		       printf(                     "% 11d % 11d % 2d\n", a, lround(b), a - lround(b));
	#else //PIECE
		surface_clear(&surface, 0);
		render_printf(&render, 0, 0, 2, 3, "% 11d % 11d % 2d\n", a, lround(b), a - lround(b));
		update();
		if(pcePadGetDirect() & PAD_SELECT) { break; }
	#endif//PIECE
		if((b >= INT_MIN) && (b <= INT_MAX)) {	//結果がINT_MIN～INT_MAXの範囲に収まらない場合を除けば…
			if(abs(a - lround(b)) > 1) {	//小数部の切り捨て誤差の範囲で正しい結果が得られるはず。
				die("x = % 11d\n"
				    "y = % 11d\n"
				    "n = % 11d\n"
				    "d = % 11d\n"
				    "a = % 11d\n"
				    "b = % 11d\n", x, y, n, d, a, lround(b));
			}
		}
	}
	return 0;
}
#endif
//-----------------------------------------------------------------------------
//高速で動き出して減速し、また加速する動作です。
//[in]
//	x		開始値。		INT_MIN～INT_MAX	┬処理の都合上、演算途中で、(y + x)がINT_MIN～INT_MAXの範囲を超えた場合も、結果が不正になります。
//	y		終了値。		INT_MIN～INT_MAX	┘
//	n		補間係数の分子。	0～d
//	d		補間係数の分母。	1～SHRT_MAX
//[out]
//	戻り値		補間結果。
//[note]
//	- lerp_i()は基本処理なので汎用性を重視して引数の値の範囲の制限を設けないようにしましたが、当関数(easeOutIn_i())は末端処理だと思うので汎用性は重視しませんでした。
//	  そのため、lerp_i()よりも引数d(ひいては引数nも)の値の範囲の制限が厳しいです。
//	  一般的な使い方では充分だと思いますが、引数の値が5桁(=?万)以上になるケースでは、注意して下さい。
//{{2017/08/08削除:easeInOut_i()とeaseOutIn_i()の処理を共通化して、コードサイズを削減しました。
//#ifndef PIECE
//int easeOutIn_i(int x, int y, int n, int d) {
//	y += x;								//処理の都合上、演算途中で、(y + x)がINT_MIN～INT_MAXの範囲を超えた場合も、結果が不正になります。
//	n <<= 1;
//	if(n >= d) {
//		x = -x;
//		x += y;
//		n = -n;
//		n += d << 1;
//	}
//	y >>= 1;
//	return easeOut_i(x, y, n, d);										//easeInOut_i()とeaseOutIn_i()の違いは、この行だけです。
//}
//#else //PIECE
//int easeOutIn_i(int x, int y, int n, int d);
//asm("
//		.code
//		.align		1
//		.global		easeOutIn_i
//easeOutIn_i:
//		sla		%r14, 1			;//%r14 := n <<= 1
//		cmp		%r14, %r15		;//if(n >= d) {
//		jrlt.d		easeOutIn_i_L10		;//  
//		add		%r13, %r12		;//%r13 := y += x		*delay*
//		not		%r12, %r12		;//  %r12 := x = -x
//		add		%r12, 1			;//  
//		add		%r12, %r13		;//  %r12 := x += y
//		not		%r14, %r14		;//  %r14 := n = -n
//		add		%r14, 1			;//  
//		add		%r14, %r15		;//  %r14 := n += d << 1
//		add		%r14, %r15		;//  
//easeOutIn_i_L10:					;//}
//		xjp.d		easeOut_i		;//return easeIn_i(x, y, n, d)				//easeInOut_i()とeaseOutIn_i()の違いは、この行だけです。
//		sra		%r13, 1			;//%r13 := y >>= 1		*delay*
//");
//#endif//PIECE
//}}2017/08/08削除:easeInOut_i()とeaseOutIn_i()の処理を共通化して、コードサイズを削減しました。
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#if 0
//□easeOutIn_i()のテストスイート
int app_main(int argc, char* argv[]) {
	int x, y, n, d, a;
	double b;
	for(;;) {
		     x = mrand48();
		     y = mrand48();
		do { d = lrand48() % (SHRT_MAX + 1); } while(!d);
		     n = lrand48() % d;
		if((((double)x + (double)y) < (double)INT_MIN) || (((double)x + (double)y) > (double)INT_MAX)) { continue; }
		a = easeOutIn_i(x, y,         n,          d);
		b = easeOutIn(  x, y, (double)n / (double)d);
	#ifndef PIECE
		       printf(                     "% 11d % 11d % 2d\n", a, lround(b), a - lround(b));
	#else //PIECE
		surface_clear(&surface, 0);
		render_printf(&render, 0, 0, 2, 3, "% 11d % 11d % 2d\n", a, lround(b), a - lround(b));
		update();
		if(pcePadGetDirect() & PAD_SELECT) { break; }
	#endif//PIECE
		if((b >= INT_MIN) && (b <= INT_MAX)) {	//結果がINT_MIN～INT_MAXの範囲に収まらない場合を除けば…
			if(abs(a - lround(b)) > 1) {	//小数部の切り捨て誤差の範囲で正しい結果が得られるはず。
				die("x = % 11d\n"
				    "y = % 11d\n"
				    "n = % 11d\n"
				    "d = % 11d\n"
				    "a = % 11d\n"
				    "b = % 11d\n", x, y, n, d, a, lround(b));
			}
		}
	}
	return 0;
}
#endif
//-----------------------------------------------------------------------------
//{{2017/08/08追加:easeInOut_i()とeaseOutIn_i()の処理を共通化して、コードサイズを削減しました。
#ifndef PIECE
static int easeInOut_easeOutIn_i(int x, int y, int n, int d, int (*fn)(int x, int y, int n, int d));		//前方参照
int easeInOut_i(int x, int y, int n, int d) { return easeInOut_easeOutIn_i(x, y, n, d, easeIn_i); }
int easeOutIn_i(int x, int y, int n, int d) { return easeInOut_easeOutIn_i(x, y, n, d, easeOut_i); }
static int easeInOut_easeOutIn_i(int x, int y, int n, int d, int (*fn)(int x, int y, int n, int d)) {
	y += x;								//処理の都合上、演算途中で、(y + x)がINT_MIN～INT_MAXの範囲を超えた場合も、結果が不正になります。
	n <<= 1;
	if(n >= d) {
		x = -x;
		x += y;
		n = -n;
		n += d << 1;
	}
	y >>= 1;
	return (*fn)(x, y, n, d);
}
#else //PIECE
int easeInOut_i(int x, int y, int n, int d);
int easeOutIn_i(int x, int y, int n, int d);
asm("
		.code
		.align		1
		.global		easeInOut_i
		.global		easeOutIn_i
easeInOut_i:
		xld.w		%r9, easeIn_i			;//%r9  := fn = easeIn_i
		jp		easeInOut_easeOutIn_i
easeOutIn_i:
		xld.w		%r9, easeOut_i			;//%r9  := fn = easeOut_i
easeInOut_easeOutIn_i:
		sla		%r14, 1				;//%r14 := n <<= 1
		cmp		%r14, %r15			;//if(n >= d) {
		jrlt.d		easeInOut_easeOutIn_i_L10	;//  
		add		%r13, %r12			;//%r13 := y += x		*delay*
		not		%r12, %r12			;//  %r12 := x = -x
		add		%r12, 1				;//  
		add		%r12, %r13			;//  %r12 := x += y
		not		%r14, %r14			;//  %r14 := n = -n
		add		%r14, 1				;//  
		add		%r14, %r15			;//  %r14 := n += d << 1
		add		%r14, %r15			;//  
easeInOut_easeOutIn_i_L10:					;//}
		sra		%r13, 1				;//%r13 := y >>= 1
		jp		%r9				;//return (*fn)(x, y, n, d)	(jp.d %rbは不可)
");
#endif//PIECE
//}}2017/08/08追加:easeInOut_i()とeaseOutIn_i()の処理を共通化して、コードサイズを削減しました。
//*****************************************************************************
//	連立方程式
//*****************************************************************************
//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))
