/*	
 *	cliprnd.c
 *
 *	乱数
 *
 *	CLiP - Common Library for P/ECE
 *	Copyright (C) 2015 Naoyuki Sawa
 *
 *	* Tue Jan 06 21:26:41 JST 2015 Naoyuki Sawa
 *	- 新規作成。
 *	- 乱数関連の宣言と実装を、clipmisc.*⇒cliprnd.*へ分離しました。
 *	  使用頻度は高いので、無駄なリンクを防ぐためと言うよりも、モジュール分けして分り易くするためです。
 */
#include "clip.h"

/*****************************************************************************
 *	乱数
 *****************************************************************************/

#ifndef PIECE
//int
//rnd8(int* seed)
//{
//	/* 8ビット乱数(『Z80マシン語秘伝の書』より) */
//	*seed = *seed * 5 + 0x3711;
//	return (*seed >> 8) & 0xff;
//}
//↓* Thu Jul 06 11:48:53 JST 2006 Naoyuki Sawa
//↓- Xorshift法による32bit擬似乱数生成関数(rnd32())を新設したことに伴い、8bit擬似乱数も線形合成法→Xorshift法に変更しました。
//↓  線形合成法に問題が有ったわけではありませんが、全体の統一のためと、Xorshift法の方が良い擬似乱数(かな?)と思ったからです。
int
rnd8(int* seed)
{
	/* 32bit擬似乱数値の中央8bitを、8bit擬似乱数値として返します。
	 * - 大きなビット数の擬似乱数値から、小さなビット数の擬似乱数値を取り出す場合、
	 *   最上位や最下位から取り出すよりも、中央部を取り出す方が良いとされています。
	 *
	 *	++++-++++-++++-++++-++++-++++-++++-++++-< rnd32()
	 *	|||| |||| |||| |||| |||| |||| |||| ||||
	 *	1111 1111 1111 1111 0000 0000 0000 0000
	 *	fedc ba98 7654 3210 fedc ba98 7654 3210
	 *	               |||| ||||
	 *	               ++++-++++----------------> rnd8()
	 */
	return (unsigned char)(rnd32(seed) >> ((32 - 8) / 2));
}
//int
//rnd16(int* seed)
//{
//	/* 16ビット乱数(VisualC++ 6.0 \VC98\CRT\SRC\RAND.C を参考にしました) */
//	*seed = *seed * 214013 + 2531011;
//	return (*seed >> 16) & 0xffff;
//}
//↓* Thu Jul 06 11:48:53 JST 2006 Naoyuki Sawa
//↓- Xorshift法による32bit擬似乱数生成関数(rnd32())を新設したことに伴い、16bit擬似乱数も線形合成法→Xorshift法に変更しました。
//↓  線形合成法に問題が有ったわけではありませんが、全体の統一のためと、Xorshift法の方が良い擬似乱数(かな?)と思ったからです。
int
rnd16(int* seed)
{
	/* 32bit擬似乱数値の中央16bitを、16bit擬似乱数値として返します。
	 * - 大きなビット数の擬似乱数値から、小さなビット数の擬似乱数値を取り出す場合、
	 *   最上位や最下位から取り出すよりも、中央部を取り出す方が良いとされています。
	 *
	 *	++++-++++-++++-++++-++++-++++-++++-++++-< rnd32()
	 *	|||| |||| |||| |||| |||| |||| |||| ||||
	 *	1111 1111 1111 1111 0000 0000 0000 0000
	 *	fedc ba98 7654 3210 fedc ba98 7654 3210
	 *	          |||| |||| |||| ||||
	 *	          ++++-++++-++++-++++-----------> rnd16()
	 */
	return (unsigned short)(rnd32(seed) >> ((32 - 16) / 2));
}
/* * Thu Jul 06 11:48:53 JST 2006 Naoyuki Sawa
 * - 32bit擬似乱数生成関数(rnd32())を追加しました。
 *   32bit疑似乱数生成アルゴリズムとして、「Xorshift RNGs」を採用しました。(以下、Xorshift法と呼ぶことにします)
 *   Xorshift法の資料を、「clip\keep\Xorshift RNGs.zip」に保存しました。
 * - 旧いrnd16()関数で利用していた線形合成方法をそのまま使い、16bit右シフトせずにそのまま32bit擬似乱数として返す方法も可能ですが、
 *   その方法では、偶数/奇数が交互に生成されてしまうこと等、良くない擬似乱数であると判断したので、Xorshift法を使うことにしました。
 */
int
rnd32(int* seed)
{
#define MAGIC_NUMBER 2463534242u /* 0以外ならば何でも可 */

	/* 前回の乱数値を取り出します。
	 * - 以下のシフト演算は、算術シフトでなく論理シフトで行う必要があるので、符号無し値として取り出します。
	 *   関数インターフェイス(seed引数)は、利用側の使い易さを考えて、符号付き値へのポインタとしてあります。
	 */
	unsigned y = *seed;

	/* 前回の乱数値が0だったならば、0以外の適当な値に置き換えます。
	 * - Xorshift法の乱数列は、M系列乱数に性格が似ていて、0を生成できません。
	 *   32ビット乱数列を生成する場合、周期は「2^32」でなく「(2^32)-1」です。
	 *   すなわち、0x00000000〜0xFFFFFFFFの範囲の乱数列を生成するのではなく、0x00000001〜0xFFFFFFFFの範囲の乱数列を生成します。
	 *   0を含めるには、乱数列の適当な位置に、明示的に0を挟む必要があります。
	 * - MAGIC_NUMBERの値に意味は無く、0以外の値ならば何でも構いません。
	 *   Xorshift法の原典のサンプルコードが、たまたま初期シード値として「2463534242」を使っていたので、それにちなんだだけです。
	 *   初期シード値も、0以外ならば何でも良く、Xorshift法の原典資料のサンプルコードにおいても、「2463534242」に意味はありません。
	 * - 本関数は、Xorshift法の原典のサンプルコードに手を加えて、0にも対応しているので、初期シード値として0を指定することも可能です。
	 *   MAGIC_NUMBERにだけは、0を指定することは不可です。
	 */
	if(y == 0) y = MAGIC_NUMBER;

	/* Xorshift法を使って、前回の乱数値から、今回の乱数値を生成します。
	 * - シフト数[13,17,5]の組み合わせは、Xorshift法の原典のサンプルコードを参考にしました。
	 *   Xorshift法の原典には、これ以外にも、様々なシフト数の組み合わせが載っています。
	 */
	y ^= (y << 13);
	y ^= (y >> 17);
	y ^= (y <<  5);

	/* 今回の乱数値が特定の値になったら、その前に、明示的に0を挟みます。
	 * - Xorshift法は0を生成できないので、乱数列の適当な位置に、明示的に0を挟む必要があります。
	 *   上述のコメントを参照してください。
	 */
	if(y == MAGIC_NUMBER) y = 0;

	/* 今回の乱数値を、次回の乱数シードとして格納します。
	 * - 関数の戻り値である乱数値と格納する乱数シードは同じ値なので、seed引数を無くすことも可能ですが、
	 *   たとえば、以下のような呼び出し方ができると便利なので、あえてseed引数を設けてあります。
	 *
	 *	retval = rnd32(&seed) % 100; // 0〜99の乱数値を生成する
	 *
	 *   また、今後、疑似乱数生成に利用するアルゴリズムを変更し、戻り値と乱数シードが同じでなくなった場合にも対応できます。
	 *   実際に今回、疑似乱数生成アルゴリズムを、線形合成法→Xorshift法に変更したわけですが、関数インターフェイスに変化はありません。
	 * - rnd32()の戻り値は(int)なので、0x00000000〜0xFFFFFFFFではなく、-0x80000000〜0x7FFFFFFFとなります。
	 *   rnd8()の戻り値が0x00〜0xFFであることや、rnd16()の戻り値が0x0000〜0xFFFFであることに較べると、一貫性に欠けます。
	 *   それでもあえてこの仕様とした理由は、rnd8()もrnd16()も戻り値が(int)であり、rnd32()だけ戻り値の型を変えるのが嫌だったからです。
	 *   0x00000000〜0xFFFFFFFFの乱数値が必要ならば、呼び出し側にて戻り値を(unsigned)にキャストしてください。
	 *   明示的な(unsigned)へのキャストを要求する方が、見た目で判らない(unsigned)値が式に与える影響の伝播よりも、安全だと考えました。
	 */
	return (*seed = y);

#undef MAGIC_NUMBER
}
int
rnd8_range(int* seed, int a, int b)
{
	return rnd8(seed) % (b - a) + a;
}
int
rnd16_range(int* seed, int a, int b)
{
	return rnd16(seed) % (b - a) + a;
}
int
rnd32_range(int* seed, int a, int b)
{
	return (unsigned)rnd32(seed) % (unsigned)(b - a) + a;
	//     ~~~~~~~~~~←─必要!!─→~~~~~~~~~~
}
#else /*PIECE*/
asm("
		.code
		.align		1
		.global		rnd8
		.global		rnd16
		.global		rnd32
		.global		rnd8_range
		.global		rnd16_range
		.global		rnd32_range
		;//-------------------------------------;//
rnd8:
		call		rnd16			;// %r10 := y = rnd32(seed) >> 12
		sra		%r10, 4			;//           = rnd16(seed) >> 4
		ret.d
		ld.ub		%r10, %r10		;// return (unsigned char)y		*delay*
		;//-------------------------------------;//
rnd16:
		call		rnd32			;// %r10 := y = rnd32(seed) >> 8
		sra		%r10, 8
		ret.d
		ld.uh		%r10, %r10		;// return (unsigned short)y		*delay*
		;//-------------------------------------;//
rnd32:
		ld.w		%r10, [%r12]		;// %r10 := y = *seed
		xld.w		%r11, 2463534242	;// %r11 := MAGIC_NUMBER		*anti-interlock*
		cmp		%r10, 0			;// if(y == 0) {
		jrne		2			;//    y  = MAGIC_NUMBER
		 ld.w		%r10, %r11		;// }
		;//
		ld.w		%r9, %r10		;// %r9  :=       y
		xsll		%r9, 13			;// %r9  :=      (y << 13)
		xor		%r10, %r9		;// %r10 := y ^= (y << 13)
		;//
		ld.w		%r9, %r10		;// %r9  :=       y
		xsrl		%r9, 17			;// %r9  :=      (y >> 17)
		xor		%r10, %r9		;// %r10 := y ^= (y >> 17)
		;//
		ld.w		%r9, %r10		;// %r9  :=       y
		xsll		%r9, 5			;// %r9  :=      (y <<  5)
		xor		%r10, %r9		;// %r10 := y ^= (y <<  5)
		;//
		cmp		%r10, %r11		;// if(y == MAGIC_NUMBER) {
		jrne		2			;//    y  = 0
		 ld.w		%r10, 0			;// }
		;//
		ld.w		[%r12], %r10		;// *seed = y
		ret					;// return  y
		;//-------------------------------------;//
rnd8_range:
		ld.w		%r15, %r13		;// %r15 :=                                             a
		ld.w		%r13, %r14		;// %r13 :=                                    b
		call.d		rnd8			;// %r10 :=            rnd8(seed)				アセンブラ版の rnd8は %r9,%r10,%r11,%psr 以外を破壊しません。
		sub		%r13, %r15		;// %r13 :=                                   (b - a)		*delay*
		;//
		xcall.d		__modsi3		;// %r10 :=            rnd8(seed) %           (b - a)		CLiP版の__modsi3は %r9,%r10,%alr,%ahr,%psr 以外を破壊しません。(clipidiv.s 参照)
		ld.w		%r12, %r10		;// %r12 :=            rnd8(seed)				*delay*
		;//
		ret.d
		add		%r10, %r15		;// %r10 :=            rnd8(seed) %           (b - a) + a	*delay*
		;//-------------------------------------;//
rnd16_range:
		ld.w		%r15, %r13		;// %r15 :=                                             a
		ld.w		%r13, %r14		;// %r13 :=                                    b
		call.d		rnd16			;// %r10 :=           rnd16(seed)				アセンブラ版のrnd16は %r9,%r10,%r11,%psr 以外を破壊しません。
		sub		%r13, %r15		;// %r13 :=                                   (b - a)		*delay*
		;//
		xcall.d		__modsi3		;// %r10 :=           rnd16(seed) %           (b - a)		CLiP版の__modsi3は %r9,%r10,%alr,%ahr,%psr 以外を破壊しません。(clipidiv.s 参照)
		ld.w		%r12, %r10		;// %r12 :=           rnd16(seed)				*delay*
		;//
		ret.d
		add		%r10, %r15		;// %r10 :=           rnd16(seed) %           (b - a) + a	*delay*
		;//-------------------------------------;//
rnd32_range:
		ld.w		%r15, %r13		;// %r15 :=                                             a
		ld.w		%r13, %r14		;// %r13 :=                                    b
		call.d		rnd32			;// %r10 :=           rnd32(seed)				アセンブラ版のrnd32は %r9,%r10,%r11,%psr 以外を破壊しません。
		sub		%r13, %r15		;// %r13 :=                                   (b - a)		*delay*
		;//
		xcall.d		__umodsi3		;// %r10 := (unsigned)rnd32(seed) % (unsigned)(b - a)		CLiP版の__modsi3は %r9,%r10,%alr,%ahr,%psr 以外を破壊しません。(clipidiv.s 参照)
		ld.w		%r12, %r10		;// %r12 :=           rnd32(seed)				*delay*
		;//
		ret.d
		add		%r10, %r15		;// %r10 := (unsigned)rnd32(seed) % (unsigned)(b - a) + a	*delay*
");
#endif/*PIECE*/

/*****************************************************************************
 *	ユーティリティ：シャフル
 *****************************************************************************/

//{{2010/10/21:高速化・省メモリ化に伴い削除
//void
//shuffle(void* base, int num, int width, int* seed)
//{
//	int i;
//	int c;
//	int i_dst;
//	int i_src;
//	unsigned char* p_dst;
//	unsigned char* p_src;
//
//	///* 0要素未満の配列は不正です。また、
//	// * 現在の実装は16bit乱数を利用しているため、65536要素超過の配列には未対応です。
//	// */
//	//if((num < 0) || (num > (1<<16))) {
//	//	DIE();
//	//}
//	//↓2006/07/06 numの範囲検査を削除しました。
//	// * Thu Jul 06 11:48:53 JST 2006 Naoyuki Sawa
//	// - rnd32()の新設に伴い、16bit乱数→32bit乱数を利用するよう変更しました。
//	//   この変更により、65536要素以上の配列をシャッフルできなかった制限が無くなりました。
//	//   (num < 0)の制限は引き続き存在しますが、自明なので、検査を省略することにしました。
//
//	/* 先頭要素から最終要素まで順番に... */
//	p_dst = base;
//	for(i_dst = 0; i_dst < num; i_dst++) {
//
//		/* 交換相手の要素を決定します。 */
//		//i_src = RND16_RANGE(*seed, i_dst, num);
//		//↓* Thu Jul 06 11:48:53 JST 2006 Naoyuki Sawa
//		//↓- rnd32()の新設に伴い、16bit乱数→32bit乱数を利用するよう変更しました。
//		//↓  この変更により、65536要素以上の配列をシャッフルできなかった制限が無くなりました。
//		i_src = RND32_RANGE(*seed, i_dst, num);
//		p_src = (unsigned char*)base + width * i_src;
//
//		/* 1バイトづつ交換します。 */
//		for(i = 0; i < width; i++) {
//			c = *p_dst;
//			    *p_dst = *p_src;
//			             *p_src = c;
//			p_dst++;
//			p_src++;
//		}
//	}
//}
//}}2010/10/21:高速化・省メモリ化に伴い削除
//↓ * Thu Oct 21 03:01:12 JST 2010 Naoyuki Sawa
//↓ - shuffle()関数を、高速化・省メモリ化しました。
//↓   高速化の効果は、要素数32768個のint型配列をシャフルした場合、約10%の高速化です。
//↓   省メモリ化の効果は、4バイト小さくなっただけです。
void
shuffle(void* base, int num, int width, int* seed)
{
	int i;
	int c;
	unsigned char* p_dst;
	unsigned char* p_src;

	/* 交換先要素は、配列の先頭要素から、配列の(最終要素-1)まで。 */
	p_dst = base;
	while(num > 1) {
		/* 交換元要素は、交換先要素から、配列の最終要素まで。 */
		p_src = &p_dst[((unsigned)rnd32(seed) % num) * width];
		//              ~~~~~~~~~~必要!
		/* 交換先要素と、交換元要素を、一文字づつ交換します。 */
		for(i = 0; i < width; i++) {
			c = *p_dst;
			    *p_dst = *p_src;
			             *p_src = c;
			p_dst++;
			p_src++;
		}
		num--;
	}
}

/*****************************************************************************
 *	ユーティリティ：範囲乱数
 *****************************************************************************/

#ifndef PIECE
unsigned minimum_mask_greater_than_or_equal_to(unsigned x) {
	unsigned y = 0;
	while(y < x) {
		y = (y << 1) | 1;
	}
	return y;
}
unsigned random_value_less_than_or_equal_to(unsigned x, unsigned (*generator)(void* param), void* param) {
	unsigned mask = minimum_mask_greater_than_or_equal_to(x);
	unsigned y;
	do {
		y = (*generator)(param) & mask;
	} while(y > x);
	return  y;
}
#else /*PIECE*/
asm("
		.align		1
		.global		minimum_mask_greater_than_or_equal_to
		.global		random_value_less_than_or_equal_to
		;//-----------------------------------------------------;//
minimum_mask_greater_than_or_equal_to:
		;//							;// %r12  := x     ┌────────────────┐
		ld.w		%r10, 0					;// %r10  := y = 0 │  ┌────────┐          │
		;//							;// do {           │  ↓                │          │
		cmp		%r10, %r12				;//     %psr[C] := c = (y < x) ? 1 : 0   │          │
		jrult.d		-1					;//                              ┌───┴───┐  │
		adc		%r10, %r10				;//     %psr[C] := y[31], %r10 := (y[30:0]<<1) | c ←┘		*delay*
		;//							;// } while(c == 1)  └──────┐
		adc		%r10, %r8				;// %r10 := (y[30:0]<<1) | y[31] ←┘
		ret.d							;// return   y[31:0]
		rr		%r10, 1					;// %r10 :=  y[31:0]						*delay*
		;//-----------------------------------------------------;//
random_value_less_than_or_equal_to:
		pushn		%r3
		ld.w		%r0, %r12				;// %r0  := x
		ld.w		%r1, %r13				;// %r1  := generator
		call.d		minimum_mask_greater_than_or_equal_to	;// %r10 := mask = minimum_mask_greater_than_or_equal_to(x)
		ld.w		%r2, %r14				;// %r2  := param						*delay*
		ld.w		%r3, %r10				;// %r3  := mask
		;//							;// do {
		call.d		%r1					;//   %r10 := y = (*generator)(param)
		 ld.w		%r12, %r2				;//   %r12 :=                  param				*delay*
		 and		%r10, %r3				;//   %r10 := y &= mask
		 cmp		%r10, %r0				;// } while(y > x)
		jrugt		-4
		popn		%r3
		ret
		;//-----------------------------------------------------;//
");
#endif/*PIECE*/
