/*
	from http://http.developer.nvidia.com/GPUGems3/gpugems3_ch37.html
*/


/*
 *	Random nubmers on the GPU
 *
 *	
 *	float RandUniform(unsigned *seeds, unsigned stride);	// float, [0.0 1.0)
 *	unsigned RandUniformui(unsigned *seeds, unsigned stride);	// unsigned,  [0, RAND_MAX]
 *	float RandNormal(unsigned *seeds, unsigned stride);		// float, gaussian mean = 0 std = 1
 *
 *	seeds must point to 4 unsigned values with the specified stride
 */


__device__ __host__ inline unsigned TausStep(unsigned &z, int S1, int S2, int S3, unsigned M)
{
	unsigned b = (((z << S1) ^ z) >> S2);
	return z = (((z & M) << S3) ^ b);
}

__device__ __host__ inline unsigned LCGStep(unsigned &z, unsigned A, unsigned C)
{
	return z = (A*z + C);
}

/* generate a random number, uses an array of 4 unsigned ints */
__device__ __host__ inline float HybridTaus(unsigned *z, unsigned stride)
{
	return 2.3283064365387e-10 * (float)(
		TausStep(z[0], 13, 19, 12, 4294967294UL) ^
		TausStep(z[stride], 2, 25, 4, 4294967288UL) ^
		TausStep(z[2*stride], 3, 11, 17, 4294967280UL) ^
		LCGStep(z[3*stride], 16654525, 1013904223UL)
	);
}

__device__ __host__ inline unsigned HybridTausui(unsigned *z, unsigned stride)
{
	return (
		TausStep(z[0], 13, 19, 12, 4294967294UL) ^
		TausStep(z[stride], 2, 25, 4, 4294967288UL) ^
		TausStep(z[2*stride], 3, 11, 17, 4294967280UL) ^
		LCGStep(z[3*stride], 16654525, 1013904223UL)
	);
}

/*
	Take two random [0,1) floats and produce 2 independent, gaussians with mean=0, var=1
*/
#define PI 3.14159265358979f
//__device__ inline void BoxMuller(float& u1, float& u2){
//    float   r = sqrtf(-2.0f * logf(u1));
//    float phi = 2 * PI * u2;
//    u1 = r * __cosf(phi);
//    u2 = r * __sinf(phi);
//}

__device__ __host__ inline void BoxMuller(float& u1, float& u2){
    float   r = sqrtf(-2.0f * logf(u1));
    float phi = 2 * PI * u2;
    u1 = r * cosf(phi);
    u2 = r * sinf(phi);
}

__device__ __host__ inline float RandUniform(unsigned *z, unsigned stride)
{
	return HybridTaus(z, stride);
}

__device__ __host__ inline unsigned RandUniformui(unsigned *z, unsigned stride)
{
	return HybridTausui(z, stride);
}

__device__ __host__ inline float RandNorm(unsigned *z, unsigned stride)
{
	float u1 = RandUniform(z, stride);
	float u2 = RandUniform(z, stride);
    float r = sqrtf(-2.0f * logf(u1));
    float phi = 2 * PI * u2;
    return r * cosf(phi);
//	return 0.123f;
}


	.file	"tmpxft_0039dc9f_00000000-6_cuda_rand.cudafe1.cpp"
	.text
#APP

#NO_APP
	.type	_ZL26__cudaUnregisterBinaryUtilv, @function
_ZL26__cudaUnregisterBinaryUtilv:
.LFB2018:
	.cfi_startproc
	endbr64
	movq	_ZL20__cudaFatCubinHandle(%rip), %rdi
	jmp	__cudaUnregisterFatBinary@PLT
	.cfi_endproc
.LFE2018:
	.size	_ZL26__cudaUnregisterBinaryUtilv, .-_ZL26__cudaUnregisterBinaryUtilv
	.section	.text.startup,"ax",@progbits
	.type	_ZL24__sti____cudaRegisterAllv, @function
_ZL24__sti____cudaRegisterAllv:
.LFB2041:
	.cfi_startproc
	endbr64
	pushq	%rax
	.cfi_def_cfa_offset 16
	leaq	_ZL15__fatDeviceText(%rip), %rdi
	call	__cudaRegisterFatBinary@PLT
	movq	%rax, %rdi
	movq	%rax, _ZL20__cudaFatCubinHandle(%rip)
	call	__cudaRegisterFatBinaryEnd@PLT
	leaq	_ZL26__cudaUnregisterBinaryUtilv(%rip), %rdi
	popq	%rdx
	.cfi_def_cfa_offset 8
	jmp	atexit@PLT
	.cfi_endproc
.LFE2041:
	.size	_ZL24__sti____cudaRegisterAllv, .-_ZL24__sti____cudaRegisterAllv
	.section	.init_array,"aw"
	.align 8
	.quad	_ZL24__sti____cudaRegisterAllv
	.section	.nvFatBinSegment,"aw"
	.align 8
	.type	_ZL15__fatDeviceText, @object
	.size	_ZL15__fatDeviceText, 24
_ZL15__fatDeviceText:
	.long	1180844977
	.long	1
	.quad	fatbinData
	.quad	0
	.local	_ZL20__cudaFatCubinHandle
	.comm	_ZL20__cudaFatCubinHandle,8,8
	.ident	"GCC: (Ubuntu 13.3.0-6ubuntu2~24.04) 13.3.0"
	.section	.note.GNU-stack,"",@progbits
	.section	.note.gnu.property,"a"
	.align 8
	.long	1f - 0f
	.long	4f - 1f
	.long	5
0:
	.string	"GNU"
1:
	.align 8
	.long	0xc0000002
	.long	3f - 2f
2:
	.long	0x3
3:
	.align 8
4:


#include <hip/hip_runtime.h>
/*
	from http://http.developer.nvidia.com/GPUGems3/gpugems3_ch37.html
*/


/*
 *	Random nubmers on the GPU
 *
 *	
 *	float RandUniform(unsigned *seeds, unsigned stride);	// float, [0.0 1.0)
 *	unsigned RandUniformui(unsigned *seeds, unsigned stride);	// unsigned,  [0, RAND_MAX]
 *	float RandNormal(unsigned *seeds, unsigned stride);		// float, gaussian mean = 0 std = 1
 *
 *	seeds must point to 4 unsigned values with the specified stride
 */


__device__ __host__ inline unsigned TausStep(unsigned &z, int S1, int S2, int S3, unsigned M)
{
	unsigned b = (((z << S1) ^ z) >> S2);
	return z = (((z & M) << S3) ^ b);
}

__device__ __host__ inline unsigned LCGStep(unsigned &z, unsigned A, unsigned C)
{
	return z = (A*z + C);
}

/* generate a random number, uses an array of 4 unsigned ints */
__device__ __host__ inline float HybridTaus(unsigned *z, unsigned stride)
{
	return 2.3283064365387e-10 * (float)(
		TausStep(z[0], 13, 19, 12, 4294967294UL) ^
		TausStep(z[stride], 2, 25, 4, 4294967288UL) ^
		TausStep(z[2*stride], 3, 11, 17, 4294967280UL) ^
		LCGStep(z[3*stride], 16654525, 1013904223UL)
	);
}

__device__ __host__ inline unsigned HybridTausui(unsigned *z, unsigned stride)
{
	return (
		TausStep(z[0], 13, 19, 12, 4294967294UL) ^
		TausStep(z[stride], 2, 25, 4, 4294967288UL) ^
		TausStep(z[2*stride], 3, 11, 17, 4294967280UL) ^
		LCGStep(z[3*stride], 16654525, 1013904223UL)
	);
}

/*
	Take two random [0,1) floats and produce 2 independent, gaussians with mean=0, var=1
*/
#define PI 3.14159265358979f
//__device__ inline void BoxMuller(float& u1, float& u2){
//    float   r = sqrtf(-2.0f * logf(u1));
//    float phi = 2 * PI * u2;
//    u1 = r * __cosf(phi);
//    u2 = r * __sinf(phi);
//}

__device__ __host__ inline void BoxMuller(float& u1, float& u2){
    float   r = sqrtf(-2.0f * logf(u1));
    float phi = 2 * PI * u2;
    u1 = r * cosf(phi);
    u2 = r * sinf(phi);
}

__device__ __host__ inline float RandUniform(unsigned *z, unsigned stride)
{
	return HybridTaus(z, stride);
}

__device__ __host__ inline unsigned RandUniformui(unsigned *z, unsigned stride)
{
	return HybridTausui(z, stride);
}

__device__ __host__ inline float RandNorm(unsigned *z, unsigned stride)
{
	float u1 = RandUniform(z, stride);
	float u2 = RandUniform(z, stride);
    float r = sqrtf(-2.0f * logf(u1));
    float phi = 2 * PI * u2;
    return r * cosf(phi);
//	return 0.123f;
}


	.text
	.file	"cuda_rand.hip"
	.type	__hip_cuid_,@object             # @__hip_cuid_
	.bss
	.globl	__hip_cuid_
__hip_cuid_:
	.byte	0                               # 0x0
	.size	__hip_cuid_, 1

	.ident	"AMD clang version 18.0.0git (https://github.com/RadeonOpenCompute/llvm-project roc-6.3.2 25012 e5bf7e55c91490b07c49d8960fa7983d864936c4)"
	.section	".note.GNU-stack","",@progbits
	.addrsig
	.addrsig_sym __hip_cuid_

	.text
	.p2alignl 7, 3214868480
	.fill 96, 4, 3214868480
	.type	__hip_cuid_,@object             ; @__hip_cuid_
	.section	.bss,"aw",@nobits
	.globl	__hip_cuid_
__hip_cuid_:
	.byte	0                               ; 0x0
	.size	__hip_cuid_, 1

	.ident	"AMD clang version 18.0.0git (https://github.com/RadeonOpenCompute/llvm-project roc-6.3.2 25012 e5bf7e55c91490b07c49d8960fa7983d864936c4)"
	.section	".note.GNU-stack","",@progbits
	.addrsig
	.addrsig_sym __hip_cuid_
	.amdgcn_target "amdgcn-amd-amdhsa--gfx1100"
	.amdgpu_metadata
---
amdhsa.kernels:  []
amdhsa.target:   amdgcn-amd-amdhsa--gfx1100
amdhsa.version:
  - 1
  - 2
...

	.end_amdgpu_metadata
