/*
 *	framflt2.c
 *
 *	高速単精度浮動小数点演算ライブラリ
 *
 *	CLiP - Common Library for P/ECE
 *	Copyright (C) 2001-2016 Naoyuki Sawa
 *
 *	* Sat Jun 21 22:52:08 JST 2008 Naoyuki Sawa
 *	- 作成開始。
 *	- sqrtf()を作成しました。
 *	* Mon Jun 23 14:59:47 JST 2008 Naoyuki Sawa
 *	- modff()を作成しました。
 *	- fmodf()を作成しました。
 *	* Mon Jun 23 23:18:11 JST 2008 Naoyuki Sawa
 *	- cosf()を作成しました。
 *	- sinf()を作成しました。
 *	* Thu Jun 26 01:31:00 JST 2008 Naoyuki Sawa
 *	- expf()を作成しました。
 *	- ldexpf()を作成しました。
 *	* Sat Jun 28 00:49:08 JST 2008 Naoyuki Sawa
 *	- atanf()を作成しました。
 *	* Sun Jun 29 01:18:14 JST 2008 Naoyuki Sawa
 *	- frexpf()を作成しました。
 *	- logf()を作成しました。
 *	* Mon Jun 30 00:05:52 JST 2008 Naoyuki Sawa
 *	- expf()を高速化しました。
 *	- cosf()を高速化しました。
 *	- sinf()を高速化しました。
 *	* Mon Jun 30 13:54:50 JST 2008 Naoyuki Sawa
 *	- atanf_sub()を高速化しました。
 *	- logf()を高速化しました。
 *	* Mon Jun 30 22:30:41 JST 2008 Naoyuki Sawa
 *	- floorf()を作成しました。
 *	- ceilf()を作成しました。
 *	* Thu Jul 03 23:09:28 JST 2008 Naoyuki Sawa
 *	- atan2f()を作成しました。
 *	* Fri Jul 04 00:09:03 JST 2008 Naoyuki Sawa
 *	- fabsf()を作成しました。
 *	* Sat Jul 12 22:49:09 JST 2008 Naoyuki Sawa
 *	- C99の関数名に変更しました。(float_sqrt→sqrtfなど)
 *	  #defineで対応しようかとも思っていましたが、やっぱり完全に書き換えることにしました。
 *	  #defineでは、アセンブラルーチンから使うときに、元の名前でしか呼び出せないからです。
 *	- 全置換で変換したので、当日以前のコードやコメントも、全部C99の関数名になっています。
 *	  もう、float_*()形式の名前はどこにも残っていないので、忘れてしまって結構です。
 *	* Thu Jul 17 02:03:15 JST 2008 Naoyuki Sawa
 *	- tanf()を作成しました。
 *	* Fri Jul 18 01:31:07 JST 2008 Naoyuki Sawa
 *	- sqrtf()のASM版のバグを修正しました。
 *	* Sat May 02 20:49:01 JST 2009 Naoyuki Sawa
 *	- acosf()を作成しました。
 *	- asinf()を作成しました。
 *	* 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のコメントを参照して下さい。
 *	* Thu Mar 24 21:04:15 JST 2016 Naoyuki Sawa
 *	- signbit(),copysign(),copysignf()を追加しました。
 *	  これらはC99で追加された符号操作用のマクロ,及び,関数で、あまり使う事は無いと思いますが、覚えておくために実装しておく事にしました。
 *	- round(),roundf(),nearbyint(),nearbyintf(),rint(),rintf(),lround(),lroundf(),lrint(),lrintf()を追加しました。
 *	  これらはC99で追加された丸め関数で、特に、lround()が便利だと思います。
 *	  いままでは'(int)((x>=0.0)?(x+0.5):(x-0.5))'のように書いていた処理が、'lround(x)'だけで済むようになりました。
 *	- 宣言は全て、include/math.hに有ります。
 *	  倍精度バージョンの関数、copysign(),lround(),lrint(),round(),nearbyint(),rint(),trunc()の実装は、cliplibc.cに有ります。
 *	  単精度バージョンの関数、copysignf(),lroundf(),lrintf(),roundf(),nearbyintf(),rintf(),truncf()の実装は、framflt2.cに有ります。
 *	  実装を分けた理由は、単精度演算ライブラリが(アプリケーションが使用していないのに)無駄にリンクされて、サイズが増える事を防ぐためです。
 *	* Mon Mar 28 21:40:28 JST 2016 Naoyuki Sawa
 *	- fdim(),fdimf(),fmin(),fminf(),fmax(),fmaxf(),fma(),fmaf()を追加しました。
 *	  これらはC99で追加された関数です。
 *	  fdim()は、何に使うのか良く判らないのですが、覚えておくために実装しておく事にしました。
 *	  fmin(),fmax()は、min(),max()マクロよりも型を明示出来るし副作用が無いので、便利に使用出来ると思います。
 *	  fma()は、積和演算を短く・判り易く書けるので、(使う場面は少ないと思いますが)使い道は有ると思います。
 *	- 宣言は全て、include/math.hに有ります。
 *	  倍精度バージョンの関数、fdim(),fmin(),fmax(),fma()の実装は、cliplibc.cに有ります。
 *	  単精度バージョンの関数、fdimf(),fminf(),fmaxf(),fmaf()の実装は、framflt2.cに有ります。
 *	  実装を分けた理由は、単精度演算ライブラリが(アプリケーションが使用していないのに)無駄にリンクされて、サイズが増える事を防ぐためです。
 *	* Mon Mar 28 22:52:58 JST 2016 Naoyuki Sawa
 *	- fdimf(),fminf(),fmaxf(),fmaf()をアセンブラ化しました。
 *	  fdimf(),fminf(),fmaxf()の共通処理を工夫して、上手くコードサイズを削減出来たと思います。
 *	  尚、__fcmpsと__mulsfが%r4〜%r7を破壊しない事に依存しています。
 *	  もし今後、__fcmpsや__mulsfの実装を変更する場合は、注意して下さい。
 *	* Sun May 29 21:54:56 JST 2016 Naoyuki Sawa
 *	- Visual C++ .NETならば、コンパイラが持っている単精度版の数学関数を使うようにしました。
 *	  詳細は、以下の通りです。
 *	- VC++2015で当ライブラリをコンパイルしてみた所、VC++2015はC99に対応しているので、コンパイラが単精度版の数学関数を持っており、当モジュールで定義している単精度版の数学関数と衝突して、コンパイルエラーになりました。
 *	  VC++6.0の場合は、C++の場合のみインライン版の単精度版の数学関数が定義されるという特殊な構造だったので、コンパイラ側の関数を無効化して常に当モジュールの単精度版の数学関数を使用するようにして対処していたのですが、
 *	  VC++2015の場合は、CでもC++でもインライン版の単精度版の数学関数が定義されるので(c:\Program Files (x86)\Windows Kits\10\Include\10.0.10240.0\ucrt\math.h)、常にコンパイラの単精度版の数学関数を使用する事にしました。
 *	- 上記の対策のために、Visual C++ .NETならば、clipflt.hの一部と、framflt2.cの全部を、条件コンパイルで無効化するようにしました。
 *	  framflt2.cは全部を無効化する代わりにVC++2015のlibclipプロジェクトから外しても良かったのですが、VC++6.0とVC++2015とでlibclipプロジェクトに含めるファイルを同じにしておきたかったので、前者の方法を取る事にしました。
 *	  VC++6.0とVC++2015とでlibclipプロジェクトに含めるファイルを同じにしておきたい理由は、VC++6.0のlibclipプロジェクトを単純に変換して、VC++2015のlibclipプロジェクトを生成したいからです。
 */
#include "clip.h"

//{{2016/05/29追加:Visual C++ .NETならば、コンパイラが持っている単精度版の数学関数を使うようにしました。
#if     !(defined(_MSC_VER) && (_MSC_VER > 1200))	//Visual C++ .NETならば、この範囲のコードを使用しません。↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
//}}2016/05/29追加:Visual C++ .NETならば、コンパイラが持っている単精度版の数学関数を使うようにしました。

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

/* * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
 * - VC++6.0で単精度版の数学関数の定義が衝突する事がある問題を回避しました。
 *   詳細は、clipflt.hの、「ANSIライブラリ数学関数(math.lib)と類似の機能の単精度版」の宣言部分の、2016/02/25のコメントを参照して下さい。
 */
#ifdef  _MSC_VER
float fn_sqrtf(float x)                { return sqrtf(x);          }
float fn_modff(float x, float* intptr) { return modff(x, intptr);  }
float fn_fmodf(float x, float y)       { return fmodf(x, y);       }
float fn_cosf(float x)                 { return cosf(x);           }
float fn_sinf(float x)                 { return sinf(x);           }
float fn_tanf(float x)                 { return tanf(x);           }
float fn_acosf(float x)                { return acosf(x);          }
float fn_asinf(float x)                { return asinf(x);          }
float fn_atanf(float x)                { return atanf(x);          }
float fn_atan2f(float y, float x)      { return atan2f(y, x);      }
float fn_expf(float x)                 { return expf(x);           }
float fn_ldexpf(float x, int exp)      { return ldexpf(x, exp);    }
float fn_frexpf(float x, int* expptr)  { return frexpf(x, expptr); }
float fn_logf(float x)                 { return logf(x);           }
float fn_floorf(float x)               { return floorf(x);         }
float fn_ceilf(float x)                { return ceilf(x);          }
float fn_fabsf(float x)                { return fabsf(x);          }
#endif//_MSC_VER

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

float sqrtf(float x);

#ifndef PIECE
// * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
// - P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
//   これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
//   clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
float sqrtf(float x) {
	return sqrt(x);
}
#else //PIECE
//{{2008/06/21:高速化のため、削除しました。
////---------------------------------------------------------------------------
////	/* アルゴリズム確認用です。asm版も同じアルゴリズムです。 */
////	float sqrtf(float _x) {
////		FLT_INT x;
////		int exp;
////		unsigned long long man1;
////		unsigned long long man2;
////		unsigned long long bit;
////		unsigned long long tmp;
////	
////		x.f = _x;
////	
////		exp = x.i.e;
////		if(!exp) {
////			return 0;
////		}
////	
////		man1 = x.i.m | (1 << 23);
////		man1 <<= 23;
////	
////		exp -= 127;
////		if(exp & 1) {
////			man1 <<= 1;
////		}
////		exp >>= 1;
////		exp += 127;
////	
////		/* man1 min = 0x0000400000000000 */
////		/*      max = 0x0000ffffff000000 */
////	
////		man2 = (1 << 23);
////		bit = man2;
////		while(bit >>= 1) {
////			tmp = man2 * man2;
////			if(tmp < man1) {
////				man2 += bit;
////			} else if(tmp > man1) {
////				man2 -= bit;
////			} else {
////				break;
////			}
////		}
////	
////		/* man2 min = 0x00800000 */
////		/*      max = 0x00ffffff */
////	
////		x.i.e = exp;
////		x.i.m = man2;
////	
////		return x.f;
////	}
////---------------------------------------------------------------------------
//asm("
//		.code
//		.align		1
//		.global		sqrtf
//sqrtf:
//		sll		%r12, 1				; %r12 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmm0
//		rl		%r12, 8				; %r12 := mmmmmmmm mmmmmmmm mmmmmmm0 eeeeeeee
//		xand		%r10, %r12, 255			; %r10 := 00000000 00000000 00000000 eeeeeeee = exp
//		jreq		sqrtf_ret			; if(!exp) return 0 (rr %r10,1はダミーになる)
//		;//
//		xor		%r12, %r10			; %r12 := mmmmmmmm mmmmmmmm mmmmmmm0 00000000
//		or		%r12, 1				; %r12 := mmmmmmmm mmmmmmmm mmmmmmm0 00000001
//		xrl		%r12, 14			; %r12 := mmmmmmmm m0000000 01mmmmmm mmmmmmmm = man1
//		;//
//		xsub		%r10, %r10, 127			; %r10 := exp -= bias
//		xand		%r15, %r10, 1			; if(exp & 1)
//		rl		%r12, %r15			;   %r12 := mmmmmmmm 00000000 1mmmmmmm mmmmmmmm = man1
//		sra		%r10, 1				; %r10 := exp >>= 1
//		xadd		%r10, %r10, 127			; %r10 := exp += bias
//		;//
//		xand		%r13, %r12, 0xffff		; %r13 := HI(man1)
//		xor		%r12, %r13			; %r12 := LO(man1)
//		;//
//		xld.w		%r11, 0x800000			; %r11 := 00000000 10000000 00000000 00000000 = man2
//		ld.w		%r14, %r11			; %r14 := 00000000 10000000 00000000 00000000 = bit
//sqrtf_loop:
//		srl		%r14, 1				; while(bit >>= 1)
//		jreq		sqrtf_exit
//		mltu.w		%r11, %r11			;   {%alr,%ahr} := man2 ^ 2
//		ld.w		%r15, %ahr			;   %r15 := HI(man2 ^ 2)
//		cmp		%r15, %r13
//		jrult		sqrtf_ult			;   if(man2 ^ 2 < man1) man2 += bit, continue
//		jrugt		sqrtf_ugt			;   if(man2 ^ 2 > man1) man2 -= bit, continue
//		ld.w		%r15, %alr			;   %r15 := LO(man2 ^ 2)
//		cmp		%r15, %r12
//		jrult		sqrtf_ult			;   if(man2 ^ 2 < man1) man2 += bit, continue
//		jrugt		sqrtf_ugt			;   if(man2 ^ 2 > man1) man2 -= bit, continue
//		;//						;   if(man2 ^ 2 == man1) break
//sqrtf_exit:
//		xsll		%r11, 9				; %r11 := mmmmmmmm mmmmmmmm mmmmmmm0 00000000
//		or		%r10, %r11			; %r10 := mmmmmmmm mmmmmmmm mmmmmmm0 eeeeeeee
//		rr		%r10, 8				; %r10 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmm0
//sqrtf_ret:
//		ret.d
//		rr		%r10, 1				; %r10 := 0eeeeeee emmmmmmm mmmmmmmm mmmmmmmm			*delay*
//sqrtf_ult:
//		jp.d		sqrtf_loop
//		add		%r11, %r14			;   %r11 := man2 += bit						*delay*
//sqrtf_ugt:
//		jp.d		sqrtf_loop
//		sub		%r11, %r14			;   %r11 := man2 -= bit						*delay*
//");
//}}2008/06/21:高速化のため、削除しました。
//↓* Sat Jun 21 23:12:05 JST 2008 Naoyuki Sawa
//↓- 高速化しました。
//↓  ・比較対象の2乗数を、64bitの上位32bit内にずらして、下位32bitを0としました。
//↓    上位32bit、下位32bitの順で比較する際、最初の上位32bit比較で大小が確定する可能性が高くなります。
//↓  ・比較対象の2乗数の下位32bitは0なので、.GT.か.EQ.のみで、.LT.の場合を考慮する必要が無くなります。
//↓  ・以上の工夫により、約10%速くなり、コードサイズもわずかに小さくなりました。
//---------------------------------------------------------------------------
//	/* アルゴリズム確認用です。asm版も同じアルゴリズムです。 */
//	float sqrtf(float _x) {
//		FLT_INT x;
//		int exp;
//		unsigned long long man1;
//		unsigned long long man2;
//		unsigned long long bit;
//		unsigned long long tmp;
//		int i;
//	
//		x.f = _x;
//	
//		exp = x.i.e;
//		if(!exp) {
//			return 0;
//		}
//	
//		man1 = x.i.m | (1 << 23);
//		man1 <<= 33;
//	
//		exp -= 127;
//		if(exp & 1) {
//			man1 <<= 1;
//		}
//		exp >>= 1;
//		exp += 127;
//	
//		/* man1 min = 0x0100000000000000 */
//		/*      max = 0x03fffffc00000000 */
//	
//		man2 = (1 << 28);
//		bit = man2;
//		i = 23;
//		do {
//			bit >>= 1;
//			tmp = man2 * man2;
//			if(tmp < man1) {
//				man2 += bit;
//			} else if(tmp > man1) {
//				man2 -= bit;
//			} else {
//				break;
//			}
//		} while(--i);
//	
//		/* man2 min = 0x10000000 */
//		/*      max = 0x1fffffe0 */
//	
//		x.i.e = exp;
//		x.i.m = man2 >> 5;
//	
//		return x.f;
//	}
//---------------------------------------------------------------------------
//{{2008/07/18:バグ修正のため、削除しました。
//asm("
//		.code
//		.align		1
//		.global		sqrtf
//sqrtf:
//		sll		%r12, 1				; %r12 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmm0
//		rl		%r12, 8				; %r12 := mmmmmmmm mmmmmmmm mmmmmmm0 eeeeeeee
//		xand		%r10, %r12, 255			; %r10 := 00000000 00000000 00000000 eeeeeeee = exp
//		jreq		sqrtf_ret			; if(!exp) return 0 (rr %r10,1はダミーになる)
//		;//
//		xor		%r12, %r10			; %r12 := mmmmmmmm mmmmmmmm mmmmmmm0 00000000
//		or		%r12, 1				; %r12 := mmmmmmmm mmmmmmmm mmmmmmm0 00000001
//		rr		%r12, 8				; %r12 := 00000001 mmmmmmmm mmmmmmmm mmmmmmm0 = HI(man1)
//		;//
//		xsub		%r10, %r10, 127			; %r10 := exp -= bias
//		xand		%r15, %r10, 1			; if(exp & 1)
//		rl		%r12, %r15			;   %r12 := 0000001m mmmmmmmm mmmmmmmm mmmmmm00 = HI(man1)
//		sra		%r10, 1				; %r10 := exp >>= 1
//		xadd		%r10, %r10, 127			; %r10 := exp += bias
//		;//
//		ld.w		%r11, 1				; %r11 := 00010000 00000000 00000000 00000000 = man2
//		rr		%r11, 4
//		ld.w		%r13, %r11			; %r13 := 00001000 00000000 00000000 00000000 = bit
//		rr		%r13, 1
//		ld.w		%r14, 23			; %r14 := i = 23
//sqrtf_loop:							; do {
//		mltu.w		%r11, %r11			;   {%alr,%ahr} := man2 ^ 2
//		ld.w		%r15, %ahr			;   %r15 := HI(man2 ^ 2)
//		cmp		%r15, %r12													←{{2008/07/18:この比較命令で反映したフラグを…}}
//		jrult.d		sqrtf_ult			;   if(man2 ^ 2 < man1) man2 += bit
//		add		%r11, %r13			;   %r11 := man2 + bit						*delay*		←{{2008/07/18:条件分岐の前に破壊してしまうバグでした}}
//		jrugt.d		sqrtf_ugt			;   if(man2 ^ 2 > man1) man2 -= bit						←{{2008/07/18:本当はここまで維持しなければなりません}}
//		sub		%r11, %r13			;   %r11 := man2						*delay*
//		ld.w		%r15, %alr			;   %r15 := LO(man2 ^ 2)
//		cmp		%r15, 0
//		jreq		sqrtf_exit			;   if(man2 ^ 2 == man1) break
//sqrtf_ugt:
//		sub		%r11, %r13			;   %r11 := man2 - bit
//sqrtf_ult:
//		sub		%r14, 1				; } while(--i)
//		jrne.d		sqrtf_loop
//		rr		%r13, 1				;								*delay*
//sqrtf_exit:
//		sll		%r11, 4				; %r11 := mmmmmmmm mmmmmmmm mmmmmmm0 00000000
//		or		%r10, %r11			; %r10 := mmmmmmmm mmmmmmmm mmmmmmm0 eeeeeeee
//		rr		%r10, 8				; %r10 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmm0
//sqrtf_ret:
//		ret.d
//		rr		%r10, 1				; %r10 := 0eeeeeee emmmmmmm mmmmmmmm mmmmmmmm			*delay*
//");
//}}2008/07/18:バグ修正のため、削除しました。
// * Fri Jul 18 01:31:07 JST 2008 Naoyuki Sawa
// - 上のASM版にはバグがあり、「sqrtf(1.0f)->2.0f」や「sqrtf(4.0f)->4.0f」になってしまっていました。
//   バグの原因は、比較命令で反映したフラグを、条件分岐の前に破壊してしまっていたことでした。
//   これはASM版特有のバグであり、アルゴリズム確認用のC言語版は正しい結果が得られます。
// - 上述のバグを修正しました。
//   遅延分岐が使えなくなりましたが、ほとんど性能低下は無く、コードはさらに1命令小さくなっています。
asm("
		.code
		.align		1
		.global		sqrtf
sqrtf:
		sll		%r12, 1				; %r12 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmm0
		rl		%r12, 8				; %r12 := mmmmmmmm mmmmmmmm mmmmmmm0 eeeeeeee
		xand		%r10, %r12, 255			; %r10 := 00000000 00000000 00000000 eeeeeeee = exp
		jreq		sqrtf_ret			; if(!exp) return 0 (rr %r10,1はダミーになる)
		;//
		xor		%r12, %r10			; %r12 := mmmmmmmm mmmmmmmm mmmmmmm0 00000000
		or		%r12, 1				; %r12 := mmmmmmmm mmmmmmmm mmmmmmm0 00000001
		rr		%r12, 8				; %r12 := 00000001 mmmmmmmm mmmmmmmm mmmmmmm0 = HI(man1)
		;//
		xsub		%r10, %r10, 127			; %r10 := exp -= bias
		xand		%r15, %r10, 1			; if(exp & 1)
		rl		%r12, %r15			;   %r12 := 0000001m mmmmmmmm mmmmmmmm mmmmmm00 = HI(man1)
		sra		%r10, 1				; %r10 := exp >>= 1
		xadd		%r10, %r10, 127			; %r10 := exp += bias
		;//
		ld.w		%r11, 1				; %r11 := 00010000 00000000 00000000 00000000 = man2
		rr		%r11, 4
		ld.w		%r13, %r11			; %r13 := 00010000 00000000 00000000 00000000 = bit * 2
		ld.w		%r14, 23			; %r14 := i = 23
sqrtf_loop:							; do {
		mltu.w		%r11, %r11			;   {%alr,%ahr} := man2 ^ 2
		ld.w		%r15, %ahr			;   %r15 := HI(man2 ^ 2)
		cmp		%r15, %r12
		jrult		sqrtf_ult			;   if(man2 ^ 2 < man1) man2 += bit
		jrugt		sqrtf_ugt			;   if(man2 ^ 2 > man1) man2 -= bit
		ld.w		%r15, %alr			;   %r15 := LO(man2 ^ 2)
		cmp		%r15, 0
		jreq		sqrtf_exit			;   if(man2 ^ 2 == man1) break
sqrtf_ugt:
		sub		%r11, %r13			;   %r11 := man2 - bit * 2 (下で相殺して(man2 - bit)になります)
sqrtf_ult:
		rr		%r13, 1				;   %r13 := bit * 1
		sub		%r14, 1				; } while(--i)
		jrne.d		sqrtf_loop
		add		%r11, %r13			; %r11 := man2 + bit						*delay*
sqrtf_exit:
		sll		%r11, 4				; %r11 := mmmmmmmm mmmmmmmm mmmmmmm0 00000000
		or		%r10, %r11			; %r10 := mmmmmmmm mmmmmmmm mmmmmmm0 eeeeeeee
		rr		%r10, 8				; %r10 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmm0
sqrtf_ret:
		ret.d
		rr		%r10, 1				; %r10 := 0eeeeeee emmmmmmm mmmmmmmm mmmmmmmm			*delay*
");
#endif//PIECE

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

float modff(float x, float* intptr);

#ifndef PIECE
// * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
// - P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
//   これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
//   clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
float modff(float x, float* intptr) {
	double intval;
	double retval = modf(x, &intval);
	*intptr = intval;
	return retval;
}
#else //PIECE
asm("
		.code
		.align		1
		.global		modff
modff:
		;// %r12 := xval
		;// %r13 := intptr
		ld.w		%r10, %r12				; %r10 := x = xval
;//		ld.w		%r11, 0					; %r11 := y = 0 ----------------------------------+
		;// %r10 := x = xval					;                                                 |
		;// %r11 := y = 0					;                                                 |
		;// %r12 := xval					;                                                 |
		;// %r13 := intptr					;                                                 |
		rl		%r12, 1					; %r12 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmms     |
		xand		%r4, %r12, 1				; %r4  := 00000000 00000000 00000000 0000000s = s |
		rl		%r12, 8					; %r12 := mmmmmmmm mmmmmmmm mmmmmmms eeeeeeee     |
		xand		%r5, %r12, 255				; %r5  := 00000000 00000000 00000000 eeeeeeee = e |
		xcmp		%r5, 127				; if(e < 127) return [x, 0]                       |
		jrlt.d		modff_ret				;                                                 |
		ld.w		%r11, 0					; <-----------------------------------------------+		*delay*
		;// %r4  := s (0 or 1)
		;// %r5  := e (min = 127, max = 255)
		;// %r10 := x = xval
		;// %r11 := y = 0
		;// %r12 := mmmmmmmm mmmmmmmm mmmmmmms eeeeeeee
		;// %r13 := intptr
		ld.w		%r11, %r10				; %r11 := y = xval
;//		ld.w		%r10, 0					; %r10 := x = 0 -----------+
		;// %r4  := s (0 or 1)					;                          |
		;// %r5  := e (min = 127, max = 255)			;                          |
		;// %r10 := x = 0					;                          |
		;// %r11 := y = xval					;                          |
		;// %r12 := mmmmmmmm mmmmmmmm mmmmmmms eeeeeeee		;                          |
		;// %r13 := intptr					;                          |
		xld.w		%r6, 150				; %r6  := i = 150 - e      |
		sub		%r6, %r5				;                          |
		jrle.d		modff_ret				; if(i <= 0) return [0, x] |
		ld.w		%r10, 0					; <------------------------+					*delay*
		;// %r4  := s (0 or 1)
		;// %r5  := e (min = 127, max = 149)
		;// %r6  := i (min =   1, max =  23)
		;// %r10 := x = 0
		;// %r11 := y = xval
		;// %r12 := mmmmmmmm mmmmmmmm mmmmmmms eeeeeeee
		;// %r13 := intptr
		ld.w		%r7, -1					; %r7  := 11111111 11111111 11111111 11111111 = mask
		;//{{VARSHIFT
		cmp		%r6, 8					; %r6  := mask <<= i
		jrle		4
		sll		%r7, 8
		jp.d		-3
		sub		%r6, 8					;								*delay*
		sll		%r7, %r6
		;//}}VARSHIFT
		not		%r6, %r7				; %r6  := ~mask << 9
		xsll		%r6, 9
		;// %r4  := s (0 or 1)
		;// %r5  := e (min = 127, max = 149)
		;// %r6  := 00000000 00000000 00000010 00000000 (i =  1)
		;//         00000000 00000000 00000110 00000000 (i =  2)
		;//         00000000 00000000 00001110 00000000 (i =  3)
		;//                          .
		;//                          .
		;//                          .
		;//         00111111 11111111 11111110 00000000 (i = 21)
		;//         01111111 11111111 11111110 00000000 (i = 22)
		;//         11111111 11111111 11111110 00000000 (i = 23)
		;// %r7  := 11111111 11111111 11111111 11111110 (i =  1)
		;//         11111111 11111111 11111111 11111100 (i =  2)
		;//         11111111 11111111 11111111 11111000 (i =  3)
		;//                          .
		;//                          .
		;//                          .
		;//         11111111 11100000 00000000 00000000 (i = 21)
		;//         11111111 11000000 00000000 00000000 (i = 22)
		;//         11111111 10000000 00000000 00000000 (i = 23)
		;// %r10 := x = 0
		;// %r11 := y
		;// %r12 := mmmmmmmm mmmmmmmm mmmmmmms eeeeeeee
		;// %r13 := intptr
		and		%r12, %r6				; %r12 := mmmmmmmm mmmmmmmm mmmmmmms eeeeeeee & (~mask << 9)
		jreq.d		modff_ret				; if(!m) return [0, y]
		and		%r11, %r7				; %r11 := y = x & mask						*delay*
		;// %r4  := s (0 or 1)
		;// %r5  := e (min = 127, max = 149)
		;// %r10 := x = 0
		;// %r11 := y
		;// %r12 := mmmmmmmm mmmmmmmm mmmmmmm0 00000000
		;// %r13 := intptr
		;//{{VARSCAN
		scan1		%r6, %r12
		sll		%r12, %r6
		jrult.d		-2
		sub		%r5, %r6				;								*delay*
		;//}}VARSCAN
		;// %r4  := s (0 or 1)
		;// %r5  := e (min = 127 - 22, max = 149 - 0)
		;// %r10 := x = 0
		;// %r11 := y
		;// %r12 := 1mmmmmmm mmmmmmmm mmmmmmm0 00000000
		;// %r13 := intptr
		sub		%r5, 1					; %r5  := e-- (min = 127 - 23, max = 149 - 1)
		sll		%r12, 1					; %r12 := mmmmmmmm mmmmmmmm mmmmmm00 00000000
		or		%r12, %r5				; %r12 := mmmmmmmm mmmmmmmm mmmmmm00 eeeeeeee
		rr		%r12, 8					; %r12 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmm00
		or		%r12, %r4				; %r12 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmm0s
		rr		%r12, 1					; %r12 := seeeeeee emmmmmmm mmmmmmmm mmmmmmm0
		ld.w		%r10, %r12				; %r10 := seeeeeee emmmmmmm mmmmmmmm mmmmmmm0
modff_ret:
		;// %r10 := x
		;// %r11 := y
		;// %r13 := intptr
		ld.w		[%r13], %r11				; *intptr = y
		ret
");
#endif//PIECE

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

float fmodf(float x, float y);

#ifndef PIECE
// * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
// - P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
//   これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
//   clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
float fmodf(float x, float y) {
	return fmod(x, y);
}
#else //PIECE
asm("
		.code
		.align		1
		.global		fmodf
fmodf:
		;// %r12    := x
		;// %r13    := y
		;//-----------------------------------------------------;
		ld.w		%r10, %r12				; %r10 :=  x
		ld.w		%r11, %r13				; %r11 :=  y
		sll		%r10, 1					; %r10 := |x| << 1
		sll		%r11, 1					; %r11 := |y| << 1
		sub		%r10, %r11				; if(|x| >  |y|) goto ugt
		jrugt		fmodf_ugt
		;//-----------------------------------------------------;
		jreq		fmodf_ret				; if(|x| == |y|) return 0 (%r10は既に0です)
		;//-----------------------------------------------------;
		jp.d		fmodf_ret				; if(|x| <  |y|) return x
		ld.w		%r10, %r12				; %r10 := x							*delay*
fmodf_ugt:	;//-----------------------------------------------------;
		xsub		%sp, %sp, 8
		xld.w		[%sp+4], %r13				; [%sp+4] := y
		;// %r12    := x
		;// %r13    := y
		;// [%sp+4] := y
		xcall		__divsf3				; %r10 :=      x / y
		;// %r10    := x / y
		;// [%sp+4] := y
		ld.w		%r13, %sp				; %r13 :=             &dummy
		xcall.d		modff					; %r10 := modf(x / y, &dummy)
		ld.w		%r12, %r10				; %r12 :=      x / y						*delay*
		;// %r10    := modf(x / y, &dummy)
		;// [%sp+0] := dummy
		;// [%sp+4] := y
		xld.w		%r12, [%sp+4]				; %r12 :=                       y
		xcall.d		__mulsf3				; %r10 := modf(x / y, &dummy) * y
		ld.w		%r13, %r10				; %r13 := modf(x / y, &dummy)					*delay*
		;// %r10    := modf(x / y, &dummy) * y
		xadd		%sp, %sp, 8
fmodf_ret:	;//-----------------------------------------------------;
		ret
");
#endif//PIECE

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

float cosf(float x);
float sinf(float x);
float tanf(float x);

#ifndef PIECE
// * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
// - P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
//   これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
//   clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
float cosf(float x) {
	return cos(x);
}
float sinf(float x) {
	return sin(x);
}
float tanf(float x) {
	return tan(x);
}
#else //PIECE
//{{2008/06/30:高速化のため、削除しました。
////	float cosf(float x) {
////		return cosf_sinf(x, 0/*cos*/);
////	}
//asm("
//		.code
//		.align		1
//		.global		cosf
//cosf:
//		xjp.d		cosf_sinf
//		ld.w		%r13, 0			; *delay*
//");
////	float sinf(float x) {
////		return cosf_sinf(x, 1/*sin*/);
////	}
//asm("
//		.code
//		.align		1
//		.global		sinf
//sinf:
//		xjp.d		cosf_sinf
//		ld.w		%r13, 1			; *delay*
//");
//static float cosf_sinf(float x, int cos_sin/*0=cos,1=sin*/) __attribute__((unused));
//static float cosf_sinf(float x, int cos_sin/*0=cos,1=sin*/) {
//	float aa;
//	float a;
//	float b;
//	float c;
//	int d;
//
//	/* |x|が大きいほど収束に時間がかかるので、-2π〜+2πに切り詰める。 */
//	x = fmodf(x, 6.283185307179586476925286766559f);
//
//	a = 0;
//	c = 1;
//	if(!cos_sin) {
//		/* cos(x) = x^0/0! - x^2/2! + x^4/4! - x^6/6! + ... */
//		b = 1;
//		d = 0;
//	} else {
//		/* sin(x) = x^1/1! - x^3/3! + x^5/5! - x^7/7! + ... */
//		b = x;
//		d = 1;
//	}
//	x *= x;
//
//	for(;;) {
//		aa = a;
//		a += b / c;
//		if(a == aa) {
//			break;
//		}
//		b *= x;
//		c *= (float)-((d + 1) * (d + 2));
//		d += 2;
//	}
//
//	return a;
//}
//}}2008/06/30:高速化のため、削除しました。
//↓* Mon Jun 30 00:17:28 JST 2008 Naoyuki Sawa
//↓- アセンブラ化し、高速化しました。
asm("
		.code
		.align		1
		.global		cosf
		.global		sinf
		.global		tanf
tanf:		;//-------------------------------------;
		pushn		%r0
		call.d		sinf			; %r10 := sinf(x)
		ld.w		%r0, %r12		; %r0  :=                x			*delay*
		ld.w		%r12, %r0		; %r12 :=                x
		call.d		cosf			; %r10 :=           cosf(x)
		ld.w		%r0, %r10		; %r0  := sinf(x)				*delay*
		ld.w		%r12, %r0		; %r12 := sinf(x)
		popn		%r0
		xjp.d		__divsf3		; %r10 := sinf(x) / cosf(x)
		ld.w		%r13, %r10		; %r13 :=           cosf(x)			*delay*
cosf:		;//-------------------------------------;
		jp.d		cosf_sinf
		ld.w		%r4, 0			; %r4  := d = cos_sin = 0			*delay*
sinf:		;//-------------------------------------;
		ld.w		%r4, 1			; %r4  := d = cos_sin = 1
cosf_sinf:	;//-------------------------------------;
		pushn		%r0
		;// %r4  := d
		;// %r12 := x
		xld.w		%r13,0x40c90fdb		; %r10 := x = fmodf(x, 6.283185307179586476925286766559f)
		xcall.d		fmodf
		ld.w		%r0, %r4		; %r0  := d					*delay*
		;// %r0  := d
		;// %r10 := x
		ld.w		%r4, 0			; %r4  := a = 0.0f
		xld.w		%r5, 0x3f800000		; %r5  := c = 1.0f
		cmp		%r0, 0			; if(d)
		jreq.d		3
		 ld.w		%r6, %r5		; %r6  := b = 1.0f				*delay*
		 ld.w		%r6, %r10		;   %r6  := b = x
		;// %r0  := d
		;// %r4  := a
		;// %r5  := c
		;// %r6  := b
		;// %r10 := x
		ld.w		%r12, %r10		; %r10 := x *= x
		xcall.d		__mulsf3
		ld.w		%r13, %r10		;						*delay*
		ld.w		%r7, %r10		; %r7  := x
		;// %r0  := d
		;// %r4  := a
		;// %r5  := c
		;// %r6  := b
		;// %r7  := x
cosf_sinf_loop:
		;// %r0  := d
		;// %r4  := a
		;// %r5  := c
		;// %r6  := b
		;// %r7  := x
		ld.w		%r12, %r6		; %r10 := b / c
		xcall.d		__divsf3
		ld.w		%r13, %r5		;						*delay*
		;// %r0  := d
		;// %r4  := a
		;// %r5  := c
		;// %r6  := b
		;// %r7  := x
		;// %r10 := b / c
		ld.w		%r12, %r4		; %r4  := a -> aa
		xcall.d		__addsf3		; %r10 := a += b / c
		ld.w		%r13, %r10		;						*delay*
		;// %r0  := d
		;// %r4  := aa
		;// %r5  := c
		;// %r6  := b
		;// %r7  := x
		;// %r10 := a
		cmp		%r10, %r4		; if(a == aa) break
		jreq		cosf_sinf_ret
		;//
		ld.w		%r4, %r10		; %r4  := a
		;// %r0  := d
		;// %r4  := a
		;// %r5  := c
		;// %r6  := b
		;// %r7  := x
		ld.w		%r12, %r6		; %r10 := b *= x
		xcall.d		__mulsf3
		ld.w		%r13, %r7		;						*delay*
		ld.w		%r6, %r10		; %r6  := b
		;//
		not		%r12, %r0		; %r12 := ~d =   - (d + 1)
		add		%r0, 2			; %r0  :=  d =               (d + 2)
		mlt.h		%r12, %r0		; %alr :=        -((d + 1) * (d + 2))
		xcall.d		__floatsisf		; %r10 := (float)-((d + 1) * (d + 2))
		ld.w		%r12, %alr		;						*delay*
		;//
		ld.w		%r12, %r5		; %r10 := c *= (float)-((d + 1) * (d + 2))
		xcall.d		__mulsf3
		ld.w		%r13, %r10		;						*delay*
		;//
		jp.d		cosf_sinf_loop
		ld.w		%r5, %r10		; %r5  := c					*delay*
cosf_sinf_ret:
		;// %r10 := a
		popn		%r0
		ret
");
#endif//PIECE

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

float acosf(float x);
float asinf(float x);
float atanf(float x);
float atan2f(float y, float x);

#ifndef PIECE
// * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
// - P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
//   これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
//   clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
float acosf(float x) {
	return acos(x);
}
float asinf(float x) {
	return asin(x);
}
float atanf(float x) {
	return atan(x);
}
float atan2f(float y, float x) {
	return atan2(y, x);
}
#else //PIECE
//---------------------------------------------------------------------------
//	/* アルゴリズム確認用です。asm版も同じアルゴリズムです。 */
//	float acosf(float x) {
//		return atan2f(sqrtf(1.0f - x * x), x);
//	}
//	float asinf(float x) {
//		return atan2f(x, sqrtf(1.0f - x * x));
//	}
//---------------------------------------------------------------------------
asm("
		.code
		.align		1
		.global		acosf
		.global		asinf
acosf:		;//-------------------------------------;
		pushn		%r0
		call.d		acosf_asinf		;// %r10 := atan2f(   sqrtf(1.0f - x * x), x)
		ld.w		%r0, %r12		;// %r0  :=                                x	*delay*
		ld.w		%r13, %r0		;// %r13 :=                                x
		popn		%r0
		xjp.d		atan2f			;// %r10 := atan2f(   sqrtf(1.0f - x * x), x)
		ld.w		%r12, %r10		;// %r12 :=           sqrtf(1.0f - x * x)	*delay*
asinf:		;//-------------------------------------;
		pushn		%r0
		call.d		acosf_asinf		;// %r10 := atan2f(x, sqrtf(1.0f - x * x))
		ld.w		%r0, %r12		;// %r0  :=        x                     	*delay*
		ld.w		%r12, %r0		;// %r13 :=        x                     
		popn		%r0
		xjp.d		atan2f			;// %r10 := atan2f(x, sqrtf(1.0f - x * x))
		ld.w		%r13, %r10		;// %r12 :=           sqrtf(1.0f - x * x)	*delay*
acosf_asinf:	;//-------------------------------------;
		xcall.d		__mulsf3		;// %r10 :=                        x * x
		ld.w		%r13, %r12		;// %r13 :=                            x	*delay*
		xld.w		%r12, 0x3f800000	;// %r12 :=                 1.0f
		xcall.d		__subsf3		;// %r10 :=                 1.0f - x * x
		ld.w		%r13, %r10		;// %r13 :=                        x * x	*delay*
		xjp.d		sqrtf			;// %r10 :=           sqrtf(1.0f - x * x)
		ld.w		%r12, %r10		;// %r12 :=                 1.0f - x * x	*delay*
");
/* * Thu Jun 26 23:29:28 JST 2008 Naoyuki Sawa
 * - arctan x をテイラー展開で解くと、|x| が 1.0 に近いとき、極端に収束が遅くなります。(数千回のループになります)
 *   この問題を解決するために、以下の資料を参考にさせていただきました。
 *	-------------------------------------------------------------------------------------------------------
 *	「MSN相談箱 QNo.3628429 三角関数、逆三角関数の算出方法がわかりません」の回答から一部抜粋
 *	  URL: http://questionbox.jp.msn.com/qa3628429.html (2008/06/26。keepフォルダに保存しておきました)
 *	＞arctan x は tan π/8 = √2-1 を使って, tan (θ-π/8) = (tan θ - (√2-1)) / (1 + (√2-1) tan θ) から
 *	＞arctan x = arctan [(x - (√2-1)) / (1 + (√2-1)x)] + π/8 とする方法はあったような気がします.
 *	＞つまり, |x| < 0.42 なら普通に arctan の展開級数を使い, 0.42 ≦ |x| ≦ 1 ならこの式で小さくする, と.
 *	-------------------------------------------------------------------------------------------------------
 *   上記資料では、|x| < 0.42 で条件分けしていますが、実際には、条件分けせず常にこの方法を使って大丈夫みたいです。
 *   この方法が行っていることは、要するに、「tan(0) <= x <= tan(π/4) の x に対して arctan x を求めようとすると、
 *   x ≒ tan(π/4) の時に極端に収束が遅くなる。そこで、x を tan(-π/8) <= x' <= tan(-π/8) の範囲の x' にずらして
 *   arctan x' を求め、求めた角度 -π/8 <= θ' <= π/8 を、0 <= θ <= π/8 にずらし戻して返す。」ということです。
 *   (上記資料の中の、「tan (θ-π/8) = (tan θ - (√2-1)) / (1 + (√2-1) tan θ)」は、三角関数の加法定理です。)
 */
//---------------------------------------------------------------------------
//	/* アルゴリズム確認用です。asm版も同じアルゴリズムです。 */
//	#define	M_PI_2		1.57079632679489661923f		/* π/2 */
//	float atanf(float x) {
//		int flg = 0;
//		if(x < 0.0f) {
//			flg |= 1;
//			x = -x;
//		}
//		if(x > 1.0f) {
//			flg |= 2;
//			x = 1.0f / x;
//		}
//		x = atanf_sub(x);
//		if(flg & 2) {
//			x = M_PI_2 - x;
//		}
//		if(flg & 1) {
//			x = -x;
//		}
//		return x;
//	}
//---------------------------------------------------------------------------
asm("
		.code
		.align		1
		.global		atanf
atanf:
		pushn		%r0
		ld.w		%r10, %r12		; %r10 := seeeeeee emmmmmmm mmmmmmmm mmmmmmmm = x
		;//
		rl		%r10, 1			; %r10 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmms
		ld.w		%r0, %r10		; %r0  := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmms
		srl		%r10, 1			; %r10 := 0eeeeeee emmmmmmm mmmmmmmm mmmmmmmm = x = |x|
		;//
		xld.w		%r12, 0x3f800000	; %r12 := 1.0f
		cmp		%r10, %r12		; if(x > 1.0f)
		jrle.d		atanf_L1
		 and		%r0, 1			; %r0  := 00000000 00000000 00000000 0000000s		*delay*
		 or		%r0, 2			;   %r0  := 00000000 00000000 00000000 000000ts
		 xcall.d	__divsf3		;   %r10 := x = 1.0f / x
		 ld.w		%r13, %r10		;							*delay*
atanf_L1:
		;//
		xcall.d		atanf_sub		; %r10 := x = atanf_sub(x)
		ld.w		%r12, %r10		;							*delay*
		;//
		rr		%r0, 2			; %r0  := ts000000 00000000 00000000 00000000
		add		%r0, %r0		; %r0  := s0000000 00000000 00000000 00000000, C = t
		jruge		atanf_L2		; if(t)
		 xld.w		%r12, 0x3fc90fdb	;   %r10 := x = π/2 - x
		 xcall.d	__subsf3
		 ld.w		%r13, %r10		;							*delay*
atanf_L2:
		;//
		xor		%r10, %r0		; if(s) x = -x
		popn		%r0
		ret
");
//---------------------------------------------------------------------------
//	/* アルゴリズム確認用です。asm版も同じアルゴリズムです。 */
//	#define M_PI	3.1415926535897932384626433832795f
//	#define M_PI_2	1.5707963267948966192313216916398f
//	float atan2f(float y, float x) {
//		int revx;
//		int revy;
//		int xchg;
//		float r;
//	
//		revx = x < 0.0f;
//		if(revx) {
//			x = -x;
//		}
//		revy = y < 0.0f;
//		if(revy) {
//			y = -y;
//		}
//		xchg = y > x;
//		if(xchg) {
//			r = x;
//			    x = y;
//			        y = r;
//		}
//	
//		if(x == 0.0f) {
//			return 0.0f;
//		}
//		r = atanf_sub(y / x);
//	
//		if(xchg) {
//			r = M_PI_2 - r;
//		}
//		if(revx) {
//			r = M_PI - r;
//		}
//		if(revy) {
//			r = -r;
//		}
//	
//		return r;
//	}
//---------------------------------------------------------------------------
asm("
		.code
		.align		1
		.global		atan2f
atan2f:
		pushn		%r0
		;//
		add		%r13, %r13		; %r13 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmm0 = (|x| << 1), C = s
		adc		%r0, %r0		; %r0  := ???????? ???????? ???????? ???????x (x = revx)
		srl		%r13, 1			; %r13 := 0eeeeeee emmmmmmm mmmmmmmm mmmmmmmm = |x|
		;//
		add		%r12, %r12		; %r12 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmm0 = (|y| << 1), C = s
		adc		%r0, %r0		; %r0  := ???????? ???????? ???????? ??????xy (y = revy)
		srl		%r12, 1			; %r12 := 0eeeeeee emmmmmmm mmmmmmmm mmmmmmmm = |y|
		;//
		cmp		%r12, %r13		; if(|y| > |x|) {
		jrule.d		6
		 and		%r0, 3			; %r13 := 00000000 00000000 00000000 000000xy				*delay*
		 or		%r0, 4			;   %r13 := 00000000 00000000 00000000 00000zxy (z = xchg)
		 xor		%r12, %r13
		 xor		%r13, %r12		;   %r13 := |y|
		 xor		%r12, %r13		;   %r12 := |x|
		;---------------------------------------;
;//{{自作の__divsf3は、0/0=0を返すように作ってあるので、この判定を省略して継続しても同じ結果が得られます。
;//		cmp		%r13, 0			; if(|x| == 0.0f) return 0.0f
;//		jreq.d		atan2f_ret
;//		ld.w		%r10, 0			; 分岐しないときはダミー処理。						*delay*
;//}}自作の__divsf3は、0/0=0を返すように作ってあるので、この判定を省略して継続しても同じ結果が得られます。
		;//
		xcall.d		__divsf3		; %r10 := tan(θ) = |y| / |x| (0.0f <= tan(θ) <= 1.0f)
		rr		%r0, 3			; %r0  := zxy00000 00000000 00000000 00000000				*delay*
		;//
		xcall.d		atanf_sub		; %r10 := θ (0.0f <= θ <= π/4)
		ld.w		%r12, %r10		;									*delay*
		;---------------------------------------;
		add		%r0, %r0		; %r0  := xy000000 00000000 00000000 00000000, C = xchg
		jruge		atan2f_L1		; if(xchg)
		 xld.w		%r12, 0x3fc90fdb	;   %r10 := θ = π/2 - θ (0.0f <= θ <= π/2)
		 xcall.d	__subsf3
		 ld.w		%r13, %r10		;									*delay*
atan2f_L1:
		;//
		add		%r0, %r0		; %r0  := y0000000 00000000 00000000 00000000, C = revx
		jruge		atan2f_L2		; if(revx)
		 xld.w		%r12, 0x40490fdb	;   %r10 := θ = π - θ (0.0f <= θ <= π)
		 xcall.d	__subsf3
		 ld.w		%r13, %r10		;									*delay*
atan2f_L2:
		;//
		xor		%r10, %r0		; if(revy) θ = -θ (-π <= θ <= π)
;//atan2f_ret:
		popn		%r0
		ret
");
//{{2008/06/30:高速化のため、削除しました。
//#define	M_SQRT2_1	0.414213562373095048802f	/* √2-1 */
//#define	M_PI_8		0.392699081698724154808f	/* π/8 */
//static float atanf_sub(float x) __attribute__((unused));
//static float atanf_sub(float x) {	/* 0.0 <= x <= 1.0 */
//	float aa;
//	float a;
//	float b;
//	int c;
//
//	/* 0.0 <= x <= 1.0 */
//	x = (x - M_SQRT2_1) / (x * M_SQRT2_1 + 1);	/* x=tan(0〜π/4) → x'=tan(-π/8〜π/8) */
//	/* |x| <= √2-1 */
//
//	a = x;
//	b = x;
//	c = 1;
//	x = -(x * x);
//
//	for(;;) {
//		aa = a;
//		b *= x;
//		c += 2;
//		a += b / c;
//		if(a == aa) {
//			break;
//		}
//	}
//
//	/* |a| <= π/8 */
//	return a + M_PI_8;	/* θ'=-π/8〜π/8 → θ=0〜π/4 */
//	/* 0.0 <= a <= π/4 */
//}
//}}2008/06/30:高速化のため、削除しました。
//↓* Mon Jun 30 13:54:50 JST 2008 Naoyuki Sawa
//↓- アセンブラ化し、高速化しました。
asm("
		.code
		.align		1
atanf_sub:
		xld.w		%r4, 0x3ed413cd			; %r4  := √2-1
		;// %r4  := √2-1
		ld.w		%r5, %r12			; %r5  := x
		xcall.d		__subsf3			; %r10 := x - (√2-1)
		ld.w		%r13, %r4			;						*delay*
		;// %r4  := √2-1
		;// %r5  := x
		;// %r10 := x - (√2-1)
		ld.w		%r12, %r5			; %r10 := x * (√2-1)
		ld.w		%r5, %r10			; %r5  := x - (√2-1)
		xcall.d		__mulsf3
		ld.w		%r13, %r4			;						*delay*
		;// %r5  := x - (√2-1)
		;// %r10 := x * (√2-1)
		xld.w		%r13, 0x3f800000		; %r10 := x * (√2-1) + 1
		xcall.d		__addsf3
		ld.w		%r12, %r10			;						*delay*
		;// %r5  := x - (√2-1)
		;// %r10 := x * (√2-1) + 1
		ld.w		%r12, %r5			; %r10 := x = (x - (√2-1)) / (x * (√2-1) + 1)
		xcall.d		__divsf3
		ld.w		%r13, %r10			;						*delay*
		;// %r10 := x
		ld.w		%r4, %r10			; %r4  := a = x
		ld.w		%r5, %r10			; %r5  := b = x
		ld.w		%r6, 1				; %r6  := c = 1
		ld.w		%r12, %r10			; %r7  := x = -(x * x)
		xcall.d		__mulsf3
		ld.w		%r13, %r10			;						*delay*
		xxor		%r7, %r10, 0x80000000
		;// %r4  := a
		;// %r5  := b
		;// %r6  := c
		;// %r7  := x
atanf_sub_loop:
		ld.w		%r12, %r5			; %r5  := b *= x
		xcall.d		__mulsf3
		ld.w		%r13, %r7			;						*delay*
		ld.w		%r5, %r10
		;// %r4  := a
		;// %r5  := b
		;// %r6  := c
		;// %r7  := x
		add		%r6, 2				; %r6  := c += 2
		xcall.d		__floatsisf			; %r10 := (float)c
		ld.w		%r12, %r6			;						*delay*
		;// %r4  := a
		;// %r5  := b
		;// %r6  := c
		;// %r7  := x
		;// %r10 := (float)c
		ld.w		%r12, %r5			; %r10 := b / (float)c
		xcall.d		__divsf3
		ld.w		%r13, %r10			;						*delay*
		;// %r4  := a
		;// %r5  := b
		;// %r6  := c
		;// %r7  := x
		;// %r10 := b / (float)c
		ld.w		%r12, %r4			; %r4  := a -> aa
		xcall.d		__addsf3			; %r10 := a += b / (float)c
		ld.w		%r13, %r10			;						*delay*
		;// %r4  := aa
		;// %r5  := b
		;// %r6  := c
		;// %r7  := x
		;// %r10 := a
		cmp		%r10, %r4			; if(a == aa) break
		jrne.d		atanf_sub_loop
		ld.w		%r4, %r10			; %r4  := a					*delay*
		;// %r4  := a
		xld.w		%r13, 0x3ec90fdb		; %r10 := a + π/8
		xjp.d		__addsf3
		ld.w		%r12, %r4			;						*delay*
");
#endif//PIECE

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

float expf(float x);

#ifndef PIECE
// * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
// - P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
//   これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
//   clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
float expf(float x) {
	return exp(x);
}
#else //PIECE
//---------------------------------------------------------------------------
//	/* アルゴリズム確認用です。asm版も同じアルゴリズムです。 */
//	float expf(float x) {
//		float aa;
//		float a;
//		float b;
//		int c;
//		int d;
//		int y;
//	
//		/* |x|が大きいままテイラー展開で解くと、収束に時間がかかるので、
//		 *	exp(x) = e^x * 2^y (|x|=ln(2)以下の小数、y=整数)
//		 * と変換しておいてから、e^xの部分だけをテイラー展開で解きます。
//		 */
//		x *= 1.4426950408889634073599247e0f;	/* x /= ln(2) */
//		//{{選択
//		// 後者のほうが、|x|の範囲を小さくできて収束が早いのですが、
//		// 実際には、前者は11回、後者は9回、程度の違いしか無いので、
//		// コードサイズの低減を優先して、単純な前者を選択しました。
//		y = (int)x;				/* こうすれば  0.0<=x<=1.0 になります。 */
//		//y = (int)floorf(x + 0.5f);		/* こうすれば -0.5<=x<=0.5 になります。 */
//		//}}選択
//		x -= y;
//		x *= 6.93147180559945309417e-1f;	/* x *= ln(2) */
//	
//		/* exp(x) = x^0/0! + x^1/1! + x^2/2! + x^3/3! + ... */
//		a = 1.0f; /* = x^0/0! */
//		b = 1.0f;
//		c = 1;
//		d = 1;
//		for(;;) {
//			aa = a;
//			b *= x;
//			c *= d;
//			a += b / c;
//			if(a == aa) {
//				break;
//			}
//			d++;
//		}
//	
//		return ldexpf(a, y); /* e^x * 2^y */
//	}
//---------------------------------------------------------------------------
//{{2008/06/30:高速化のため、削除しました。
//asm("
//		.code
//		.align		1
//		.global		expf
//expf:
//		pushn		%r3
//		xsub		%sp, %sp, 12
//		;//
//		xld.w		%r13, 0x3fb8aa3b	; %r10    := x *= 1.4426950408889634073599247e0f
//		xcall		__mulsf3
//		xld.w		[%sp+0], %r10		; [%sp+0] := x
//		;//
//		xcall.d		__fixsfsi		; %r10    := y = (int)x
//		ld.w		%r12, %r10		;							*delay*
//		xld.w		[%sp+4], %r10		; [%sp+4] := y
//		;//
//		xcall.d		__floatsisf		; %r10    :=      (float)y
//		ld.w		%r12, %r10		;							*delay*
//		;//
//		xld.w		%r12, [%sp+0]		; %r10    := x -= (float)y
//		xcall.d		__subsf3
//		ld.w		%r13, %r10		;							*delay*
//		;//
//		xld.w		%r13, 0x3f317218	; %r10    := x *= 6.93147180559945309417e-1f
//		xcall.d		__mulsf3
//		ld.w		%r12, %r10		;							*delay*
//		xld.w		[%sp+0], %r10		; [%sp+0] := x
//		;//
//		xld.w		%r0, 0x3f800000		; %r0     := a = 1.0f
//		ld.w		%r1, %r0		; %r1     := b = 1.0f
//		ld.w		%r2, 1			; %r2     := c = 1
//		ld.w		%r3, 1			; %r3     := d = 1
//expf_loop:
//		xld.w		[%sp+8], %r0		; [%sp+8] := aa = a
//		;//
//		ld.w		%r13, [%sp+0]		; %r10    := b *= x
//		xcall.d		__mulsf3
//		ld.w		%r12, %r1		;							*delay*
//		ld.w		%r1, %r10		; %r1     := b
//		;//
//		mlt.w		%r2, %r3		; %r2     := c *= d
//		ld.w		%r2, %alr
//		;//
//		xcall.d		__floatsisf		; %r10    :=          (float)c
//		ld.w		%r12, %r2		;							*delay*
//		;//
//		ld.w		%r12, %r1		; %r10    :=      b / (float)c
//		xcall.d		__divsf3
//		ld.w		%r13, %r10		;							*delay*
//		;//
//		ld.w		%r12, %r0		; %r10    := a += b / (float)c
//		xcall.d		__addsf3
//		ld.w		%r13, %r10		;							*delay*
//		ld.w		%r0, %r10		; %r0     := a
//		;//
//		xld.w		%r12, [%sp+8]		; if(a != aa) goto loop
//		cmp		%r0, %r12		;							!INTERLOCK!
//		jrne.d		expf_loop
//		add		%r3, 1			; %r3     := d++					*delay*
//		;//
//		xld.w		%r13, [%sp+4]		; %r10    := ldexpf(a, y)
//		xcall.d		ldexpf
//		ld.w		%r12, %r0		;							*delay*
//		;//
//		xadd		%sp, %sp, 12
//		popn		%r3
//		ret
//");
//}}2008/06/30:高速化のため、削除しました。
//↓ * Mon Jun 30 00:08:49 JST 2008 Naoyuki Sawa
//↓ - 単精度浮動小数点の基本演算サブルーチンが、%r4,%r5,%r6,%r7の値を変更しないことを仕様としたことに追従し、
//↓   レジスタ使用を最適化して、高速化しました。
asm("
		.code
		.align		1
		.global		expf
expf:
		pushn		%r1
		;// %r12 := x
		xld.w		%r13, 0x3fb8aa3b	; %r10 := x *= 1.4426950408889634073599247e0f
		xcall		__mulsf3
		ld.w		%r4, %r10		; %r4  := x
		;// %r4  := x
		xcall.d		__fixsfsi		; %r10 := y = (int)x
		ld.w		%r12, %r4		;							*delay*
		ld.w		%r5, %r10		; %r5  := y
		;// %r4  := x
		;// %r5  := y
		xcall.d		__floatsisf		; %r10 :=      (float)y
		ld.w		%r12, %r5		;							*delay*
		;// %r4  := x
		;// %r5  := y
		;// %r10 := (float)y
		ld.w		%r12, %r4		; %r10 := x -= (float)y
		xcall.d		__subsf3
		ld.w		%r13, %r10		;							*delay*
		;// %r5  := y
		;// %r10 := x
		xld.w		%r13, 0x3f317218	; %r10 := x *= 6.93147180559945309417e-1f
		xcall.d		__mulsf3
		ld.w		%r12, %r10		;							*delay*
		ld.w		%r4, %r10		; %r4  := x
		;// %r4  := x
		;// %r5  := y
		xld.w		%r6, 0x3f800000		; %r6  := a = 1.0f
		ld.w		%r7, %r6		; %r7  := b = 1.0f
		ld.w		%r0, 1			; %r0  := c = 1
		ld.w		%r1, 1			; %r1  := d = 1
		;// %r0  := c
		;// %r1  := d
		;// %r4  := x
		;// %r5  := y
		;// %r6  := a
		;// %r7  := b
expf_loop:
		;// %r0  := c
		;// %r1  := d
		;// %r4  := x
		;// %r5  := y
		;// %r6  := a
		;// %r7  := b
		ld.w		%r13, %r4		; %r10 := b *= x
		xcall.d		__mulsf3
		ld.w		%r12, %r7		;							*delay*
		ld.w		%r7, %r10		; %r7  := b
		;//
		mlt.w		%r0, %r1		; %r0  := c *= d	(mlt.hでは桁不足になるのでダメ)
		ld.w		%r0, %alr
		;//
		xcall.d		__floatsisf		; %r10 :=          (float)c
		ld.w		%r12, %r0		;							*delay*
		;// %r0  := c
		;// %r1  := d
		;// %r4  := x
		;// %r5  := y
		;// %r6  := a
		;// %r7  := b
		;// %r10 := (float)c
		ld.w		%r12, %r7		; %r10 :=      b / (float)c
		xcall.d		__divsf3
		ld.w		%r13, %r10		;							*delay*
		;// %r0  := c
		;// %r1  := d
		;// %r4  := x
		;// %r5  := y
		;// %r6  := a
		;// %r7  := b
		;// %r10 := b / (float)c
		ld.w		%r12, %r6		; %r6  := a -> aa
		xcall.d		__addsf3		; %r10 := a += b / (float)c
		ld.w		%r13, %r10		;							*delay*
		;// %r0  := c
		;// %r1  := d
		;// %r4  := x
		;// %r5  := y
		;// %r6  := aa
		;// %r7  := b
		;// %r10 := a
		cmp		%r10, %r6		; if(a != aa) goto loop
		ld.w		%r6, %r10		; %r6  := a
		jrne.d		expf_loop
		add		%r1, 1			; %r1  := d++						*delay*
		;// %r5  := y
		;// %r6  := a
		popn		%r1
		;//
		ld.w		%r12, %r6		; %r10 := ldexpf(a, y)
		xjp.d		ldexpf
		ld.w		%r13, %r5		;							*delay*
");
#endif//PIECE

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

float ldexpf(float x, int exp);

#ifndef PIECE
// * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
// - P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
//   これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
//   clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
float ldexpf(float x, int exp) {
	return ldexp(x, exp);
}
#else //PIECE
//---------------------------------------------------------------------------
//	/* アルゴリズム確認用です。asm版も同じアルゴリズムです。 */
//	float ldexpf(float _x, int exp) {
//		FLT_INT x;
//		int e;
//	
//		x.f = _x;
//	
//		e = x.i.e;
//		if(!e) {
//			return 0.0f;
//		}
//		e += exp;
//		if(e <= 0) {
//			return 0.0f;
//		}
//		x.i.e = e;
//	
//		return x.f;
//	}
//---------------------------------------------------------------------------
asm("
		.code
		.align		1
		.global		ldexpf
ldexpf:
		xrl		%r12, 9			; %r12 := mmmmmmmm mmmmmmmm mmmmmmms eeeeeeee
		xand		%r10, %r12, 255		; %r10 := 00000000 00000000 00000000 eeeeeeee = e
		jreq		ldexpf_ret0		; if(!e) return 0.0f
		xor		%r12, %r10		; %r12 := mmmmmmmm mmmmmmmm mmmmmmms 00000000
		add		%r10, %r13		; %r10 := e += exp
		jrle		ldexpf_ret0		; if(e <= 0) return 0.0f
		or		%r10, %r12		; %r12 := mmmmmmmm mmmmmmmm mmmmmmms eeeeeeee
		rr		%r10, 8			; %r10 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmms
		ret.d
		rr		%r10, 1			; %r10 := seeeeeee emmmmmmm mmmmmmmm mmmmmmmm		*delay*
ldexpf_ret0:
		ret.d
		ld.w		%r10, 0			;							*delay*
");
#endif//PIECE

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

float frexpf(float x, int* expptr);

#ifndef PIECE
// * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
// - P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
//   これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
//   clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
float frexpf(float x, int* expptr) {
	return frexp(x, expptr);
}
#else //PIECE
//---------------------------------------------------------------------------
//	/* アルゴリズム確認用です。asm版も同じアルゴリズムです。 */
//	float frexpf(float _x, int* expptr) {
//		FLT_INT x;
//		int e;
//	
//		x.f = _x;
//	
//		e = x.i.e;		/* e = bias + exp */
//		e--;			/* e = bias + exp - 1 */
//		if(e <= 0) {
//			*expptr = 0;
//			return 0.0f;
//		}
//		e -= (127 - 2);		/* e = exp + 1 */
//		*expptr = e;
//		x.i.e -= e;		/* x.i.e = bias - 1 */
//	
//		return x.f;
//	}
//---------------------------------------------------------------------------
asm("
		.code
		.align		1
		.global		frexpf
frexpf:
		xrl		%r12, 9			; %r12 := mmmmmmmm mmmmmmmm mmmmmmms eeeeeeee (e = bias + exp)
		ld.ub		%r10, %r12		; %r10 := e = bias + exp
		sub		%r10, 1			; %r10 := e = bias + exp - 1
		jrle		frexpf_ret0		; if(e <= 0) return 0
		xsub		%r10, %r10, 125		; %r10 := e = exp + 1
		ld.w		[%r13], %r10		; [%r13] := *expptr = exp + 1
		sub		%r12, %r10		; %r12 := mmmmmmmm mmmmmmmm mmmmmmms eeeeeeee (e = bias - 1)
		xrr		%r12, 9			; %r12 := seeeeeee emmmmmmm mmmmmmmm mmmmmmmm
		ret.d
		ld.w		%r10, %r12		; %r10 := seeeeeee emmmmmmm mmmmmmmm mmmmmmmm			*delay*
frexpf_ret0:
		ld.w		[%r13], %r8		; [%r13] := *expptr = 0
		ret.d
		ld.w		%r10, 0			; %r10 := return 0						*delay*
");
#endif//PIECE

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

float logf(float x);

#ifndef PIECE
// * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
// - P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
//   これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
//   clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
float logf(float x) {
	return log(x);
}
#else //PIECE
//{{2008/06/30:高速化のため、削除しました。
//float
//logf(float x)
//{
//	float aa;
//	float a;
//	float b;
//	int c;
//	int y;
//
//	/* ln(x) = ln(x * 2^y) = ln(x) + ln(2^y) = ln(x) + ln(2) * y */
//	x = frexpf(x, &y);
//	if(x < 0.5f) {
//		return 0.0f; /* 元のxが0.0以下だった場合は解けません */
//	}
//	/* 0.5 <= x < 1.0 */
//
//	/* log(x)の基本的なテイラー展開は、
//	 *
//	 *	log(1+x) =   x^1/1 - x^2/2 + x^3/3 - ...	(|x| < 1)
//	 *
//	 * または
//	 *
//	 *	log(1-x) = -(x^1/1 + x^2/2 + x^3/3 + ...)	(|x| < 1)
//	 *
//	 * ですが、収束が遅く、たとえばx=0.5のとき、120回ぐらいのループになります。
//	 * そこで、以下のような展開式を使って、収束を早める方法があるのだそうです。
//	 *
//	 *	log(x)/2 = 1/1 * ((x-1)/(x+1))^1 +
//	 *	           1/3 * ((x-1)/(x+1))^3 +
//	 *	           1/5 * ((x-1)/(x+1))^5 + ...		(1/2√2 <= x < √2)
//	 *
//	 * 以下のコードは、展開式を使った方法で実装しました。
//	 * 遅くても、25回ぐらいのループで収束するみたいです。
//	 */
//
//	if(x < 0.70710678118654752440084436210485f) {	/* sqrt(2)/2 */
//		x *= 2.0f;
//		y--;
//	}
//	/* sqrt(2)/2 <= x < sqrt(2) */
//
//	a = 0.69314718055994530941723212145818f * y;
//	x = (x - 1.0f) / (x + 1.0f);
//	b = x * 2.0f;
//	x *= x;
//	c = 1;
//
//	for(;;) {
//		aa = a;
//		a += b / c;
//		if(a == aa) {
//			break;
//		}
//		b *= x;
//		c += 2;
//	}
//
//	return a;
//}
//}}2008/06/30:高速化のため、削除しました。
//↓ * Mon Jun 30 13:54:50 JST 2008 Naoyuki Sawa
//↓- アセンブラ化し、高速化しました。
asm("
		.code
		.align		1
		.global		logf
logf:
		xsub		%sp, %sp, 4
		;// %r12 := x
		ld.w		%r13, %sp			; %r10 := x = frexpf(x, &y)
		xcall		frexpf
		ld.w		%r4, %r10			; %r4  := x
		;// %r4   := x
		;// [%sp] := y
		xld.w		%r13, 0x3f000000		; if(x < 0.5f) return 0
		xcall.d		__fcmps
		ld.w		%r12, %r4			;					*delay*
		jrlt		logf_ret0
		;//
		xld.w		%r13, 0x3f3504f3		; if(x < 1/2√2)
		xcall.d		__fcmps
		ld.w		%r12, %r4			;					*delay*
		xld.w		%r12, [%sp]			; %r12 := y
		jrge		logf_L1
		 xadd		%r4, %r4, 0x00800000		;   %r4  := x *= 2.0f
		 sub		%r12, 1				;   %r12 := y--
logf_L1:
		;// %r4  := x
		;// %r12 := y
		xcall		__floatsisf			; %r10 := (float)y
		;// %r4  := x
		;// %r10 := (float)y
		xld.w		%r12, 0x3f317218		; %r10 := a = ln(2) * (float)y
		xcall.d		__mulsf3
		ld.w		%r13, %r10			;					*delay*
		ld.w		%r5, %r10			; %r5  := a
		;// %r4  := x
		;// %r5  := a
		xld.w		%r6, 0x3f800000			; %r6  :=          1.0f
		ld.w		%r12, %r4			; %r10 :=     (x - 1.0f)
		xcall.d		__subsf3
		ld.w		%r13, %r6			;					*delay*
		;// %r4  := x
		;// %r5  := a
		;// %r6  := 1.0f
		;// %r10 := (x - 1.0f)
		ld.w		%r12, %r4			; %r10 :=                  (x + 1.0f)
		ld.w		%r4, %r10			; %r4  :=     (x - 1.0f)
		xcall.d		__addsf3
		ld.w		%r13, %r6			;					*delay*
		;// %r4  := (x - 1.0f)
		;// %r5  := a
		;// %r10 := (x + 1.0f)
		ld.w		%r12, %r4			; %r10 := x = (x - 1.0f) / (x + 1.0f)
		xcall.d		__divsf3
		ld.w		%r13, %r10			;					*delay*
		;// %r5  := a
		;// %r10 := x
		xadd		%r6, %r10, 0x00800000		; %r6  := b = x * 2.0f
		;// %r5  := a
		;// %r6  := b
		;// %r10 := x
		ld.w		%r12, %r10			; %r10 := x *= x
		xcall.d		__mulsf3
		ld.w		%r13, %r10			;					*delay*
		ld.w		%r4, %r10			; %r4  := x
		;// %r4  := x
		;// %r5  := a
		;// %r6  := b
		ld.w		%r7, 1				; %r7  := c = 1
		;// %r4  := x
		;// %r5  := a
		;// %r6  := b
		;// %r7  := c
logf_loop:
		xcall.d		__floatsisf			; %r10 := (float)c
		ld.w		%r12, %r7			;					*delay*
		;// %r4  := x
		;// %r5  := a
		;// %r6  := b
		;// %r7  := c
		;// %r10 := (float)c
		ld.w		%r12, %r6			; %r10 := b / (float)c
		xcall.d		__divsf3
		ld.w		%r13, %r10			;					*delay*
		;// %r4  := x
		;// %r5  := a
		;// %r6  := b
		;// %r7  := c
		;// %r10 := b / (float)c
		ld.w		%r12, %r5			; %r5  := a -> aa
		xcall.d		__addsf3			; %r10 := a += b / (float)c
		ld.w		%r13, %r10			;					*delay*
		;// %r4  := x
		;// %r5  := aa
		;// %r6  := b
		;// %r7  := c
		;// %r10 := a
		cmp		%r10, %r5			; if(a == aa) break
		jreq		logf_ret
		ld.w		%r5, %r10			; %r5  := a
		;// %r4  := x
		;// %r5  := a
		;// %r6  := b
		;// %r7  := c
		ld.w		%r12, %r6			; %r10 := b *= x
		xcall.d		__mulsf3
		ld.w		%r13, %r4			;					*delay*
		ld.w		%r6, %r10			; %r6  := b
		;//
		jp.d		logf_loop
		add		%r7, 2				; %r7  := c += 2			*delay*
		;//
logf_ret0:
		ld.w		%r10, 0				; %r10 := a = 0
logf_ret:
		;// %r10 := a
		xadd		%sp, %sp, 4
		ret
");
#endif//PIECE

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

float floorf(float x);
float ceilf(float x);

#ifndef PIECE
// * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
// - P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
//   これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
//   clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
float floorf(float x) {
	return floor(x);
}
float ceilf(float x) {
	return ceil(x);
}
#else //PIECE
//---------------------------------------------------------------------------
//	/* アルゴリズム確認用です。asm版も同じアルゴリズムです。 */
//	float floorf(float x) {
//		float y;
//	
//		x = modff(x, &y);
//		if(x < 0.0f) {
//			y -= 1.0f;
//		}
//	
//		return y;
//	}
//---------------------------------------------------------------------------
asm("
		.code
		.align		1
		.global		floorf
floorf:
		xsub		%sp, %sp, 4
		ld.w		%r13, %sp		; %r10 := x = modff(x, &y)
		xcall		modff
		cmp		%r10, 0			; %psr := x < 0.0f
		xld.w		%r10, [%sp]		; %r10 := y
		xadd		%sp, %sp, 4
		;//
		jrlt		2			; if(x >= 0.0f) return y
		 ret
		xld.w		%r13, 0x3f800000	; return y - 1.0f
		xjp.d		__subsf3
		ld.w		%r12, %r10		;				*delay*
");
//---------------------------------------------------------------------------
//	/* アルゴリズム確認用です。asm版も同じアルゴリズムです。 */
//	float ceilf(float x) {
//		return -floorf(-x);
//	}
//---------------------------------------------------------------------------
asm("
		.align		1
		.global		ceilf
ceilf:
		rl		%r12, 1			; %r10 := floorf(-x)
		xor		%r12, 1
		xcall.d		floorf
		rr		%r12, 1			;				*delay*
		;//
		rl		%r10, 1			; return -floorf(-x)
		xor		%r10, 1
		ret.d
		rr		%r10, 1			;				*delay*
");
#endif//PIECE

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

float fabsf(float x);

#ifndef PIECE
// * Thu Feb 25 21:08:23 JST 2016 Naoyuki Sawa
// - P/ECE以外では、単精度版の数学関数は、倍精度版の数学関数を呼び出すようにしました。
//   これまでは、当モジュール内の単精度版の数学関数を、C言語版が無く直接アセンブラ版で作成していたので、P/ECE以外に含められませんでした。
//   clipflt.hに依存しているモジュールを、libclipに含められるようにするためだけの目的で、動作すれば良いだけなので速度は考慮していません。
float fabsf(float x) {
	return fabs(x);
}
#else //PIECE
asm("
		.code
		.align		1
		.global		fabsf
fabsf:
		sll		%r12, 1			; %r12 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmm0
		srl		%r12, 1			; %r12 := 0eeeeeee emmmmmmm mmmmmmmm mmmmmmmm
		ret.d
		ld.w		%r10, %r12		; %r10 := 0eeeeeee emmmmmmm mmmmmmmm mmmmmmmm	*delay*
");
#endif//PIECE

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

#ifndef PIECE
float  copysignf(float x, float y) {
	int* px = (int*)&x;					//*px := seeeeeee emmmmmmm mmmmmmmm mmmmmmmm
	int* py = (int*)&y;					//*py := SEEEEEEE EMMMMMMM MMMMMMMM MMMMMMMM
	*px &= 0x7FFFFFFF;					//*px := 0eeeeeee emmmmmmm mmmmmmmm mmmmmmmm
	*py &= 0x80000000;					//*py := S0000000 00000000 00000000 00000000
	*px |= *py;						//*px := Seeeeeee emmmmmmm mmmmmmmm mmmmmmmm
	return x;						//return Seeeeeee emmmmmmm mmmmmmmm mmmmmmmm
}
long     lroundf(float x) { return (long)(x + ((x >= 0.0f) ? 0.5f : -0.5f)); }
long      lrintf(float x) { return lroundf(x); }		//本当はfesetround()の設定に従うべきだが、簡単に実装するためにlroundf()と同じ動作にした。実用上はこれでも支障無いと思うし、実際の所当関数を使う事が無い気もする。
float     roundf(float x) { return lroundf(x); }
float nearbyintf(float x) { return  roundf(x); }		//本当はfesetround()の設定に従うべきだが、簡単に実装するためにroundf()と同じ動作にした。実用上はこれでも支障無いと思うし、実際の所当関数を使う事が無い気もする。
float      rintf(float x) { return  roundf(x); }		//本当はfesetround()の設定に従うべきだが、簡単に実装するためにroundf()と同じ動作にした。実用上はこれでも支障無いと思うし、実際の所当関数を使う事が無い気もする。
float     truncf(float x) { return (long)x; }
#else //PIECE
float  copysignf(float x, float y);
long     lroundf(float x);
long      lrintf(float x);
float     roundf(float x);
float nearbyintf(float x);
float      rintf(float x);
float     truncf(float x);
asm("
		.align		1
		.global		copysignf
		.global		lroundf
		.global		lrintf
		.global		roundf
		.global		nearbyintf
		.global		rintf
		.global		truncf
copysignf:
		add		%r13, %r13			;//%r13 := EEEEEEEE MMMMMMMM MMMMMMMM MMMMMMM0, %psr(C) := S
		adc		%r12, %r12			;//%r12 := eeeeeeee mmmmmmmm mmmmmmmm mmmmmmmS
		rr		%r12, 1				;//%r12 := Seeeeeee emmmmmmm mmmmmmmm mmmmmmmm
		ret.d						;//return  Seeeeeee eeeemmmm mmmmmmmm mmmmmmmm
		ld.w		%r10, %r12			;//								*delay*
		;//---------------------------------------------;//
lroundf:
lrintf:
		xand		%r13, %r12, 0x80000000		;//%r13 :=      ((x >= 0.0f) ? 0.5f : -0.5f)
		xoor		%r13, %r13, 0x3F000000		;//
		xcall		__addsf3			;//%r10 := x += ((x >= 0.0f) ? 0.5f : -0.5f)
		xjp.d		__fixsfsi			;//return (long)x
		ld.w		%r12, %r10			;//								*delay*
		;//---------------------------------------------;//
roundf:
nearbyintf:
rintf:
		call		lroundf				;//%r10 :=       lroundf(x)
		xjp.d		__floatsisf			;//return (float)lroundf(x)
		ld.w		%r12, %r10			;//								*delay*
		;//---------------------------------------------;//
truncf:
		xcall		__fixsfsi			;//%r10 :=       (int)x
		xjp.d		__floatsisf			;//return (float)(int)x
		ld.w		%r12, %r10			;//								*delay*
");
#endif//PIECE

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

#ifndef PIECE
float fdimf(float x, float y) {
	return fmaxf(x - y, 0.0f);
}
float fminf(float x, float y) {
	return (x <= y) ? x : y;
}
float fmaxf(float x, float y) {
	return (x <= y) ? y : x;
}
float fmaf(float x, float y, float z) {
	return (x * y) + z;
}
#else //PIECE
// * Mon Mar 28 22:52:58 JST 2016 Naoyuki Sawa
// - fdimf(),fminf(),fmaxf(),fmaf()をアセンブラ化しました。
//   fdimf(),fminf(),fmaxf()の共通処理を工夫して、上手くコードサイズを削減出来たと思います。
//   尚、__fcmpsと__mulsfが%r4〜%r7を破壊しない事に依存しています。
//   もし今後、__fcmpsや__mulsfの実装を変更する場合は、注意して下さい。
float fdimf(float x, float y);
float fminf(float x, float y);
float fmaxf(float x, float y);
float fmaf(float x, float y, float z);
asm("
		.code
		.align		1
		.global		fdimf
		.global		fminf
		.global		fmaxf
		.global		fmaf
fdimf:
		xcall		__subsf3		;//%r10 :=      x - y
		ld.w		%r12, %r10		;//%r12 :=      x - y
		ld.w		%r13, 0			;//%r13 :=             0
		;//* FALLTHRU *//			;//return fmaxf(x - y, 0)
fmaxf:
		ld.w		%r4, %r13						;//%r4  :=         y
		jp.d		4							;//			─┐
		 ld.w		%r5, %r12						;//%r5  :=   x		　│	*delay*
		;//* FALLTHRU *//							;//			　│
fminf:											;//			　│
		 ld.w		%r4, %r12		;//%r4  :=   x			;//			　│			;//┐
		 ld.w		%r5, %r13		;//%r5  :=         y		;//			　│			;//│
		xcall		__fcmps			;//%psr :=  (x <=> y)		;//%psr :=  (x <=> y)	←┘			;//│
		jrle.d		3			;//       if(x <=  y) {		;//       if(x <=  y) {				;//├__fcmpsが%r4〜%r7を破壊しない事に依存しています。
		 ld.w		%r10, %r4		;//%r10 :=   x        } else {	;//%r10 :=         y  } else {	*delay*		;//│
		 ld.w		%r10, %r5		;//%r11 :=         y  }		;//%r11 :=   x        }				;//│
		ret					;//										;//┘
		;//-------------------------------------;//
fmaf:
		xcall.d		__mulsf3		;//%r10 := (x * y)								;//┐
		ld.w		%r4, %r14		;//%r4  :=           z						*delay*		;//│
		ld.w		%r12, %r10		;//%r12 := (x * y)								;//├__mulsfが%r4〜%r7を破壊しない事に依存しています。
		xjp.d		__addsf3		;//return  (x * y) + z								;//│
		ld.w		%r13, %r4		;//%r13 :=           z						*delay*		;//┘
");
#endif//PIECE

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

//{{2016/05/29追加:Visual C++ .NETならば、コンパイラが持っている単精度版の数学関数を使うようにしました。
#endif//!(defined(_MSC_VER) && (_MSC_VER > 1200))	//Visual C++ .NETならば、この範囲のコードを使用しません。↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑
//}}2016/05/29追加:Visual C++ .NETならば、コンパイラが持っている単精度版の数学関数を使うようにしました。
