; RXS-M-XS 32bit->32bit Permuted Congruential Generator
; Written by Harry Sintonen <sintonen@iki.fi> 2024-05-30
; (C) 2024 Harry Sintonen. This work is openly licensed via CC BY 4.0.
; https://creativecommons.org/licenses/by/4.0/

OPT020	EQU	0		; Enable to optimize for 68020+
RANDMAX	EQU	1		; Enable to get limited to RAND_MAX ($7fffffff)

	; PCG constants
MUL1	EQU	747796405
MUL2	EQU	277803737
ADDV	EQU	2891336453

	ifne	OPT020
	MC68020
	endc

	ifd	TEST
main:	moveq	#1,d0
	bsr	pcg_srand
	bsr	pcg_rand
	bsr	pcg_rand
	bsr	pcg_rand
	bsr	pcg_rand
	rts
	endc

	XDEF	pcg_rand

; OUT: d0.l = random value
pcg_rand:
	lea	pcg_state(pc),a0
	; fall through

	XDEF	pcg_rand_r

;  IN: a0.l = state pointer
; OUT: d0.l = random value
pcg_rand_r:
	move.l	(a0),d0
	ifne	OPT020
	move.l	d0,d1		; get old state
	mulu.l	#MUL1,d0
	add.l	#ADDV,d0
	else
	; AB * CD = (A*C << 32) + ((B*C + A*D) << 16) + B*D
	move.w	d0,d1
	mulu.w	#MUL1&$ffff,d1
	move.l	d1,a1		; a1 = B*D
	move.l	d0,d1
	swap	d0
	mulu.w	#MUL1&$ffff,d0	; d0 = B*C
	mulu.w	#MUL1>>16,d1	; d1 = A*D
	add.l	d1,d0
	swap	d0
	clr.w	d0
	add.l	a1,d0
	add.l	#ADDV,d0
	move.l	(a0),d1		; get old state
	endc
	move.l	d0,(a0)		; store new state

	move.l	d1,a1
	move.l	d1,d0
	lsr.l	#4,d1
	rol.l	#4,d0
	and.w	#$f,d0
	lsr.l	d0,d1
	move.l	a1,d0
	eor.l	d1,d0		; (state >> ((state >> 28) + 4)) ^ state

	ifne	OPT020
	mulu.l	#MUL2,d0
	else
	; AB * CD = (A*C << 32) + ((B*C + A*D) << 16) + B*D
	move.w	d0,d1
	mulu.w	#MUL2&$ffff,d1
	move.l	d1,a1		; a1 = B*D
	move.l	d0,d1
	swap	d0
	mulu.w	#MUL2&$ffff,d0	; d0 = B*C
	mulu.w	#MUL2>>16,d1	; d1 = A*D
	add.l	d1,d0
	swap	d0
	clr.w	d0
	add.l	a1,d0
	endc
	move.l	d0,d1
	swap	d1
	lsr.w	#22-16,d1
	and.w	#(1<<(32-22))-1,d1
	eor.w	d1,d0		; v ^ (v >> 22)

	ifne	RANDMAX
	bclr	#31,d0
	endc
	rts


	XDEF	pcg_srand

; IN: d0.l = seed
pcg_srand:
	lea	pcg_state(pc),a0
	; fall through

	XDEF	pcg_srand_r

; IN: d0.l = seed
;     a0.l = state pointer
pcg_srand_r:
	add.l	#ADDV,d0
	ifne	OPT020
	mulu.l	#MUL1,d0
	else
	; AB * CD = (A*C << 32) + ((B*C + A*D) << 16) + B*D
	move.w	d0,d1
	mulu.w	#MUL1&$ffff,d1
	move.l	d1,a1		; a1 = B*D
	move.l	d0,d1
	swap	d0
	mulu.w	#MUL1&$ffff,d0	; d0 = B*C
	mulu.w	#MUL1>>16,d1	; d1 = A*D
	add.l	d1,d0
	swap	d0
	clr.w	d0
	add.l	a1,d0
	endc
	add.l	#ADDV,d0
	move.l	d0,(a0)
	rts

pcg_state:
	dc.l	1940201539	; PCG state init with seed 1