;; Register allocation
e	.set	A0	; energy loaded
m0	.set	B0	; input channel mask for channels 0..31
			; (from 2nd input argument)
x	.set	A1	; time/chi2 loaded
k	.set	A1	; condition to dump samples for this channel
m1	.set	B1	; input channel mask for channels 32..63
			; (from 2nd input argument)
hptWR	.set	A2	; condition to fill time histogram and
			; write pointer to histogram contents
hpeWR	.set	B2	; write pointer to energy histogram (or 0)
mask	.set	A3	; output mask for dumping samples
bit	.set	B3	; running bitmask to build the output mask,
			; initially 0x8000 0000
xOut	.set	A4	; pointer to filtered output time/chi2
eIn	.set	B4	; read pointer to energy block
hpt	.set	A5	; base pointer for time histograms
xIn	.set	B5	; read pointer to time/chi2 block
thr_th	.set	A6	; energy cut to fill time histogram
thr4	.set	B6	; energy cut to dump samples
Nbins	.set	A7	; number of bins per channel in energy
			; histogram (16..31) and in time histogram (0..15)
n	.set	B7	; number of channels to skip
xcut	.set	A8	; chi2 cut to dump samples
Offset	.set	B8	; constant containing the biasing exponent
jmpe	.set	A9	; offset to e/t histograms
jmpt	.set	A9	; offset to e/t histograms
N	.set	B9	; current channel being processed, or Sum(n)
binT	.set	A10	; time histogram channel content
binE	.set	B10	; energy histogram channel
hptRD	.set	A11	; read pointer to time histogram
hpeRD	.set	B11	; read pointer to energy histogram
t0	.set	A12	; time extracted from time/chi2
e0	.set	B12	; log2(energy) extracted from e
g	.set	A13	; gain extracted from time/chi2
g2	.set	B13	; copy of g
hpe	.set	B14	; base pointer for energy histograms
SP	.set	B15

	;; The E and t histograms must be in opposite blocks or banks.
	;; The stack and the input data must be in opposite blocks.
_filterHisto:
		MV	B4, B0		; m0
||		MV	A4, B4		; eIn
||		MV	B5, B1		; m1
||		ADDK	256, A4
||		MV	SP, A0

  [!m0]		B	epilog0
||		LMBD	1, m0, n
||		MV	A4, B5		; xIn
||		STW	A13, *SP[0]
||		STW	B3,  *A0[-1]


		LDW	*++xIn[n], x
||		SHL	m0, n, bit
||		ZERO	N
||		STW	B14, *A0[-2]

		LDW	*++eIn[n], e
||		CLR	bit, 31, 31, m0
||		MVK	48, Nbins
|| [m0]		MPY	4, n, N
||		STW	B13, *A0[-3]

		MVKLH	48, Nbins
||		MVKL	_ehist, hpe
||		STW	A12, *SP[-4]
||		STW	B12, *A0[-5]

		XOR	bit, m0, bit
||		ZERO	hptWR
||		MVKL	_thist, hpt
||		MVKH	_ehist, hpe
||		STW	A11, *SP[-6]
||		STW	B11, *A0[-7]

		SHRU	bit, n, bit
||		ZERO	hpeWR
||		ZERO	mask
||		MVKH	_thist, hpt
||		STW	A10, *SP[-8]
||		STW	B10, *A0[-9]

loop0:
		STW	x, *xOut++
|| [m0]		B	loop0
||		LMBD	1, m0, n
||		SHL	x, 18, g:t0
||		CMPGT	x, xcut, k
||		MPY	Nbins, N, jmpt

		LDW	*++xIn[n], x
||		SHL	m0, n, m0
||		SUB	e, Offset, e0
||		SSHL	t0, 2, t0	; set the range, 10 means LSB is bit 0
||		AND	3, g, g
||		MPYHL	Nbins, N, jmpe

		LDW	*++eIn[n], e
||		SSHL	e0, 5, e0
||		MV	g, g2
||		ADD	hpt, jmpt, hptRD
||		SHRU	g:t0, 28, g:t0
|| [!k]		CMPGT	e, thr4, k

   [m0]		ADDAW	N, n, N
||		SHRU	g2:e0, 28, g2:e0
||		ADD	hpe, jmpe, hpeRD
||		LDW	*++hptRD[t0], binT
||		SADD	1, binT, binT

		LDW	*++hpeRD[e0], binE
||		CLR	m0, 31, 31, m0
||		SADD	1, binE, binE
|| [hptWR]	STW	binT, *hptWR
||		CMPGT	e, thr_th, hptWR


   [hpeWR]	STW	binE, *hpeWR
||		SHRU	bit, n, bit
||		MV	hpeRD, hpeWR
|| [hptWR]	MV	hptRD, hptWR
|| [k]		ADD	bit, mask, mask

epilog0:
		SUB	xIn, N, xIn
||		SUB	eIn, N, eIn
||		LMBD	1, m1, n
|| [!m1]	B	epilog1

		LDW	*++xIn[n], x
||		SHL	m1, n, bit

		LDW	*++eIn[n], e
||		MVK	128, N
||		SADD	1, binT, binT

   [m1]		ADDAW	N, n, N
||		CLR	bit, 31, 31, m1
||		SADD	1, binE, binE
|| [hptWR]	STW	binT, *hptWR

   [hpeWR]	STW	binE, *hpeWR
||		MV	mask, m0
||		XOR	bit, m1, bit
||		ZERO	hptWR

		SHRU	bit, n, bit
||		ZERO	hpeWR
||		ZERO	mask

loop1:
		STW	x, *xOut++
|| [m1]		B	loop1
||		LMBD	1, m1, n
||		SHL	x, 18, g:t0
||		CMPGT	x, xcut, k
||		MPY	Nbins, N, jmpt

		LDW	*++xIn[n], x
||		SHL	m1, n, m1
||		SUB	e, Offset, e0
||		SSHL	t0, 2, t0	; set the histogram range here
||		AND	3, g, g
||		MPYHL	Nbins, N, jmpe

		LDW	*++eIn[n], e
||		SSHL	e0, 5, e0
||		MV	g, g2
||		ADD	hpt, jmpt, hptRD
||		SHRU	g:t0, 28, g:t0
|| [!k]		CMPGT	e, thr4, k

   [m1]		ADDAW	N, n, N
||		SHRU	g2:e0, 28, g2:e0
||		ADD	hpe, jmpe, hpeRD
||		LDW	*++hptRD[t0], binT
||		SADD	1, binT, binT

		LDW	*++hpeRD[e0], binE
||		CLR	m1, 31, 31, m1
||		SADD	1, binE, binE
|| [hptWR]	STW	binT, *hptWR
||		CMPGT	e, thr_th, hptWR


   [hpeWR]	STW	binE, *hpeWR
||		SHRU	bit, n, bit
||		MV	hpeRD, hpeWR
|| [hptWR]	MV	hptRD, hptWR
|| [k]		ADD	bit, mask, mask

epilog1:
		LDW	*SP[-1], B3
||		SUB	xIn, N, xIn	; Optional
||		SUB	eIn, N, eIn	; Optional

		LDW	*SP[-2], B14
||		MV	SP, A0
||		SUB	xIn, xOut, xIn	; Optional, gives # of suppressed bytes

		LDW	*SP[0],  A13
||		LDW	*A0[-3], B13
||		SHR	xIn, 2, xIn	; Optional, gives # of suppressed words

		SADD	1, binE, binE
||		SADD	1, binT, binT
||		LDW	*SP[-4], A12
||		LDW	*A0[-5], B12
||		MVK	32, eIn		; Optional, to compute 32-xIn

   [hpeWR]	STW	binE, *hpeWR
|| [hptWR]	STW	binT, *hptWR
||		SUB	eIn, xIn, xIn	; Optional, the # of selected words

		B	B3
||		LDW	*SP[-6], A11
||		LDW	*A0[-7], B11

		LDW	*SP[-8], A10
||		LDW	*A0[-9], B10
||		MV	mask, A5	; immediately followed by MV m0, A4

*		MVKL	0xdeadbeef, A0
*		MVKH	0xdeadbeef, A0
		MVK	0xe0e, A0
		STW	A0, *xOut
		MV	xIn, A4
*		MV	xOut, A4	; was MV m0, A4

		.global	_filterHisto
		.global	_ehist
		.global	_thist