include(`../config.m4')

C		    cycles/limb	  good for cpu?
C AMD K8,K9		n/a
C AMD K10		n/a
C AMD bd1		n/a
C AMD bd2		n/a
C AMD bd3		n/a
C AMD bd4		n/a
C AMD zen		n/a
C AMD bobcat		n/a
C AMD jaguar		n/a
C Intel P4		n/a
C Intel PNR		n/a
C Intel NHM		n/a
C Intel SBR		n/a
C Intel IBR		n/a
C Intel HWL		 0.68		y
C Intel BWL		
C Intel SKL		
C Intel atom		n/a
C Intel SLM		n/a
C VIA nano		 ?

C This code works only for n=64m.

define(`ap',		`%rdi')
define(`n',		`%rsi')

C More commonly used variables are assigned to the first 8 xmm
C registers to maximize VEX2 usage.
define(`x0',	`%ymm0')	define(`x0l',	`%xmm0')
define(`x1',	`%ymm1')	define(`x1l',	`%xmm1')
define(`carry',	`%ymm2')	define(`carryl',`%xmm2')
define(`sum1',	`%ymm3')	define(`sum1l',	`%xmm3')
define(`tmp2',	`%ymm4')
define(`sum2',	`%ymm5')	define(`sum2l',	`%xmm5')
define(`tmp4',	`%ymm6')
define(`sum4',	`%ymm7')	define(`sum4l',	`%xmm7')
define(`tmp8',	`%ymm8')
define(`sum8',	`%ymm9')	define(`sum8l',	`%xmm9')
				define(`sum16l',`%xmm10')

C CSA(carry, a, b, c) sets 2*carry + a = a + b + c.
C The computation is (bitwise)
C	carry = majority(a,b,c) = (a&b) | (a^b) & c = (c&b) | (c^b) & a
C	a = a ^ b ^ c
C Although the functions are symmetric, the implementation in terms
C of dyadic operations uses one input later than the other two.  This
C leads to three variants.  The first uses input a ($2) late:
C	tmp = b & c;		(Where "tmp" is actually "carry')
C	b ^= c;
C	carry = tmp | (a & b);  (using c as a temporary)
C	a ^= b;
C The second uses input c ($4) late:
C	tmp = b & a;
C	b ^= a;
C	carry = tmp | c & b;  (using a as a temporary)
C	a = c ^ b;
C The third uses c as the carry output as well as an input, and
C $1 is a pure temporary:
C	tmp = b & a;
C	b ^= a;
C	a = c ^ b;	a ^ b ^ c, the sum
C	c &= b;
C	c |= tmp;	(b & a) | c & (b ^ a), the carry

define(`CSA1',
	`# ($1, $2) := $2 + $3 + $4
	vpand	$3, $4, $1
	vpxor	$3, $4, $3
	vpand	$2, $3, $4
	vpxor	$2, $3, $2
	vpor	$1, $4, $1'
m4_assert_numargs(4))
define(`CSA2',
	`# ($1, $2) := $3 + $4 + $2
	vpand	$3, $2, $1
	vpxor	$3, $2, $3
	vpand	$4, $3, $2
	vpor	$1, $2, $1
	vpxor	$4, $3, $2'
m4_assert_numargs(4))
define(`CSA3',
	`# ($4, $2) := $2 + $3 + $4 (with $1 as temp)
	vpand	$3, $2, $1
	vpxor	$3, $2, $3
	vpxor	$4, $3, $2
	vpand	$4, $3, $4
	vpor	$4, $1, $4'
m4_assert_numargs(4))

ASM_START()
	TEXT
	.balign 32,,21
PROLOGUE(mpn_popcount2)

	# 4*5 + 2 = 22 bytes.
	vpxor	sum1l, sum1l, sum1l
	vpxor	sum2l, sum2l, sum2l
	vpxor	sum4l, sum4l, sum4l
C These can use VEX2, so are 4 bytes each rather than 5
	vmovdqa	sum1l, sum8l
	vmovdqa	sum2l, sum16l
	xor	%eax, %eax
	.balign 16,,10

# This loop reads 16 256-bit words, totalling 512 bytes,
# per iteration.  The "sum" variables are loop-carried;
# the "tmp" variables are not.
L(top):	vlddqu	(ap), x0
	vlddqu	32(ap), x1
	CSA1(tmp2, sum1, x0, x1)
	vlddqu	64(ap), x0
	vlddqu	96(ap), x1
	add	$256,ap
	CSA1(carry, sum1, x0, x1)
	CSA2(tmp4, sum2, tmp2, carry)

	vlddqu	-128(ap), x0		# Hoist me
	vlddqu	-96(ap), x1
	CSA1(tmp2, sum1, x0, x1)
	vlddqu	-64(ap), x0
	vlddqu	-32(ap), x1
	CSA1(carry, sum1, x0, x1)
	CSA3(x1, sum2, tmp2, carry)
	CSA2(tmp8, sum4, tmp4, carry)

	vlddqu	(ap), x0		# Hoist me
	vlddqu	32(ap), x1
	CSA1(tmp2, sum1, x0, x1)
	vlddqu	64(ap), x0
	vlddqu	96(ap), x1
	add	$256,ap
	CSA1(carry, sum1, x0, x1)

	CSA2(tmp4, sum2, tmp2, carry)
	vlddqu	-128(ap), x0		# Hoist me
	vlddqu	-96(ap), x1
	CSA1(tmp2, sum1, x0, x1)
	vlddqu	-64(ap), x0
	vlddqu	-32(ap), x1
	# This is the point of maximum register pressure.
	CSA1(carry, sum1, x0, x1)
	CSA3(x0, sum2, tmp2, carry)
	CSA3(x0, sum4, tmp4, carry)
	CSA3(x0, sum8, tmp8, carry)
	# And here is minimum pressure; x0, x1, tmp2, tmp4, tmp8 are dead.

	# Now carry is full of weight-16 bits.  Fold and add to sum16l.
	vextracti128 $1, carry, x1l
	CSA2(x0l, sum16l, carryl, x1l)
	# x0l is now weight-32 bits, which we send to the integer ALU.
	vpextrq	$1, x0l, %rdx
	vmovq	x0l, %rcx
	popcnt	%rdx, %rdx
	popcnt	%rcx, %rcx
	add	%rdx, %rax
	add	%rcx, %rax

	sub	$64, n
	ja	L(top)

	# Loop complete (n and ap are dead), now the final summation.
	# All of the sum variables are folded to 128 bits, and summed
	# one last time, then sent to the integer ALU for final popcnt,
	# scaling, and totalling.
	#
	# Folding the 256-bit values is worthwhile, but it's faster
	# to popcnt bits 64..127 on the integer ALU (vextq, popcnt,
	# add) than to fold and sum 64 bits on the SIMD unit.

	vextracti128 $1, sum1, x0l
	# 2-imput sum (carryl, sum1l) = x0l + carryl
	vpand	sum1l, x0l, carryl
	vpxor	sum1l, x0l, sum1l
	# Send to the integer unit
	vpextrq	$1, sum1l, %rdx
	vmovq	sum1l, %rcx
	popcnt	%rdx, n
	popcnt	%rcx, %rcx
	add	%rcx, n
	# Again with the weight-2 sum2
	vextracti128 $1, sum2, x0l
	CSA3(x1l, sum2l, x0l, carryl)
	vpextrq	$1, sum2l, %rdx
	vmovq	sum2l, %rcx
	popcnt	%rdx, %rdx
	popcnt	%rcx, %rcx
	add	%rcx, %rdx
	lea	(n,%rdx,2), n

	# And weight-4 sum4
	vextracti128 $1, sum4, x0l
	CSA3(x1l, sum4l, x0l, carryl)
	vpextrq	$1, sum4l, %rdx
	vmovq	sum4l, %rcx
	popcnt	%rdx, %rdx
	popcnt	%rcx, %rcx
	add	%rcx, %rdx
	lea	(n,%rdx,4), n

	# Weights 8 and 16 are summed separately, into %eax,
	# because LEA is limited to scaling by 8.
	vextracti128 $1, sum8, x0l
	# No more 256-bit math; we could in theory do vzeroupper here.
	CSA3(x1l, sum8l, x0l, carryl)
	vpextrq	$1, sum8l, %rdx
	vmovq	sum8l, %rcx
	popcnt	%rdx, %rdx
	popcnt	%rcx, %rcx
	add	%rcx, %rdx
	lea	(%rdx,%rax,4), %rax

	# (carryl, x0l) := carryl + sum16l
	vpxor	carryl, sum16l, x0l
	vpand	carryl, sum16l, carryl

	# Weight-16 bits
	vpextrq	$1, x0l, %rdx
	vmovq	x0l, %rcx
	popcnt	%rdx, %rdx
	popcnt	%rcx, %rcx
	add	%rcx, %rdx
	lea	(%rax,%rdx,2), %rax

	# And finally, the remaining weight-32 bits
	vpextrq	$1, carryl, %rdx
	vmovq	carryl, %rcx
	vzeroupper
	popcnt	%rdx, %rdx
	popcnt	%rcx, %rcx
	add	%rcx, %rdx
	lea	(%rax,%rdx,4), %rax

	lea	(n,%rax,8), %rax
	ret
EPILOGUE()