; FILE: gg:src/own/primc/mulmod.ASM          REV: 41 --- (x * y) % z without 64bit math
; History
;

USE_060FPU	EQU	1	; use 68060 FPU for 64bit math

Main
	move.l	#2,d0
	move.l	#3,d1
	move.l	#10,d2
	;bsr	mulmod060

	move.l	#$c05c568a,d0
	move.l	d0,d1
	move.l	#-3,d2
	;bsr	mulmod060

	move.l	#$00193e34,d0
	move.l	#$0025305c,d1
	move.l	#$1a69b118,d2
	bsr	mulmod060

	move.l	#$000037f4,d0
	move.l	#$0003b60c,d1
	move.l	#$02922d18,d2
	bsr	mulmod060
	rts


; optimized 32:32->32 modulo calculation macro
;
MODU32	MACRO	*divisor, dividend(result), tmp
	cmp.l	\2,\1
	bhi.b	.MODU32nodiv\@
	divul.l	\1,\3:\2
	move.l	\3,\2
.MODU32nodiv\@
	ENDM


	IFD	USE_060FPU


	CNOP	0,4
;  IN: d0.l = x (rarely zero)
;      d1.l = y (rarely zero)
;      d2.l = z (*never* zero)
; OUT: d1.l = x * y % z
;NOTE: d0/d6/d7/a4-a5 trashed
mulmod060
	;fmove.l	#%10000000,fpcr

	fmove.l	d0,fp0
	fmul.l	d1,fp0
	fmove.x	fp0,fp1
	fdiv.l	d2,fp0
	fintrz	fp0,fp0
	fmul.l	d2,fp0
	fsub.x	fp0,fp1
	fmove.l	fp1,d1
	illegal


	fmove.l	d0,fp0
	fmove.l	d2,fp2
	fmul.l	d1,fp0
	fmove.x	fp0,fp1
	fdiv.x	fp2,fp0
	fmove.l	fp0,d0
	fmove.l	d0,fp0
	fmul.x	fp2,fp0
	fsub.x	fp0,fp1
	fmove.l	fp1,d1
	illegal

	moveq	#16,d7

	cmp.l	d0,d2
	bhi.b	.nodiv1
	move.l	d0,d6
	divul.l	d2,d0:d6	; x %= z;
.nodiv1
	cmp.l	d1,d2
	bhi.b	.nodiv2
	move.l	d1,d6
	divul.l	d2,d1:d6	; y %= z;
.nodiv2

	move.l	d0,d6
	move.l	d4,a4		; save d4
	move.l	d1,d4
	lsr.l	d7,d6		; hx = x >> 16;
	lsr.l	d7,d4		; hy = y >> 16;

	tst.w	d6
	bne.b	.big
	tst.w	d4
	bne.b	.big

	mulu.w	d0,d1		; l_xy = x * y;
	bra.b	.exit

.big
	; x &= 0x0000ffff;
	; y &= 0x0000ffff;

	move.l	d5,a5		; save d5

	move.w	d0,d5
	mulu.w	d1,d5
	mulu.w	d4,d0		; lxhy = lx * hy;
	mulu.w	d6,d4		; hxhy = hx * hy;
	mulu.w	d1,d6		; lyhx = hx * ly;
	move.l	d5,d1

	lsr.l	d7,d5
	;move.l	d1,d1		; l_xy = (lx * ly) & 0x0000ffff;
	ror.l	d7,d1		; swap l_xy

	add.l	d6,d5		; c = (c >> 16) + lyhx;
	move.l	d5,d6

	and.l	#$0000ffff,d5
	lsr.l	d7,d6		; h_xy = c >> 16;
	add.l	d0,d5		; c = (c & 0x0000ffff) + lxhy;
	move.w	d5,d1		; l_xy |= ((c & 0x0000ffff) << 16);
	lsr.l	d7,d5
	ror.l	d7,d1		; swap l_xy
	add.l	d5,d6
	move.l	a5,d5		; restore d5
	add.l	d4,d6		; h_xy += (c >> 16) + hxhy;
	beq.b	.exit
.adjust
	MODU32	d2,d1,d0
	sub.l	d2,d1		; lxy = (lxy % z) - z;
	subq.l	#1,d6
	bne.b	.adjust
.exit
	move.l	a4,d4		; restore d4

	MODU32	d2,d1,d0	; ret = l_xy % z;
	rts


	ELSE


	CNOP	0,4
;  IN: d0.l = x (rarely zero)
;      d1.l = y (rarely zero)
;      d2.l = z (*never* zero)
; OUT: d1.l = x * y % z
;NOTE: d0/d6/d7/a4-a5 trashed
mulmod060
	moveq	#16,d7

	cmp.l	d0,d2
	bhi.b	.nodiv1
	move.l	d0,d6
	divul.l	d2,d0:d6	; x %= z;
.nodiv1
	cmp.l	d1,d2
	bhi.b	.nodiv2
	move.l	d1,d6
	divul.l	d2,d1:d6	; y %= z;
.nodiv2

	move.l	d0,d6
	move.l	d4,a4		; save d4
	move.l	d1,d4
	lsr.l	d7,d6		; hx = x >> 16;
	lsr.l	d7,d4		; hy = y >> 16;

	tst.w	d6
	bne.b	.big
	tst.w	d4
	bne.b	.big

	mulu.w	d0,d1		; l_xy = x * y;
	bra.b	.exit

.big
	; x &= 0x0000ffff;
	; y &= 0x0000ffff;

	move.l	d5,a5		; save d5

	move.w	d0,d5
	mulu.w	d1,d5
	mulu.w	d4,d0		; lxhy = lx * hy;
	mulu.w	d6,d4		; hxhy = hx * hy;
	mulu.w	d1,d6		; lyhx = hx * ly;
	move.l	d5,d1

	lsr.l	d7,d5
	;move.l	d1,d1		; l_xy = (lx * ly) & 0x0000ffff;
	ror.l	d7,d1		; swap l_xy

	add.l	d6,d5		; c = (c >> 16) + lyhx;
	move.l	d5,d6

	and.l	#$0000ffff,d5
	lsr.l	d7,d6		; h_xy = c >> 16;
	add.l	d0,d5		; c = (c & 0x0000ffff) + lxhy;
	move.w	d5,d1		; l_xy |= ((c & 0x0000ffff) << 16);
	lsr.l	d7,d5
	ror.l	d7,d1		; swap l_xy
	add.l	d5,d6
	move.l	a5,d5		; restore d5
	add.l	d4,d6		; h_xy += (c >> 16) + hxhy;
	beq.b	.exit
.adjust
	MODU32	d2,d1,d0
	sub.l	d2,d1		; lxy = (lxy % z) - z;
	subq.l	#1,d6
	bne.b	.adjust
.exit
	move.l	a4,d4		; restore d4

	MODU32	d2,d1,d0	; ret = l_xy % z;
	rts


	ENDC