dnl  Alpha mpn_modexact_1c_odd -- mpn exact remainder

dnl  Copyright 2003, 2004 Free Software Foundation, Inc.

dnl  This file is part of the GNU MP Library.
dnl
dnl  The GNU MP Library is free software; you can redistribute it and/or modify
dnl  it under the terms of either:
dnl
dnl    * the GNU Lesser General Public License as published by the Free
dnl      Software Foundation; either version 3 of the License, or (at your
dnl      option) any later version.
dnl
dnl  or
dnl
dnl    * the GNU General Public License as published by the Free Software
dnl      Foundation; either version 2 of the License, or (at your option) any
dnl      later version.
dnl
dnl  or both in parallel, as here.
dnl
dnl  The GNU MP Library is distributed in the hope that it will be useful, but
dnl  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
dnl  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
dnl  for more details.
dnl
dnl  You should have received copies of the GNU General Public License and the
dnl  GNU Lesser General Public License along with the GNU MP Library.  If not,
dnl  see https://www.gnu.org/licenses/.

include(`../config.m4')


C      cycles/limb
C EV4:    47
C EV5:    30
C EV6:    15


C mp_limb_t mpn_modexact_1c_odd (mp_srcptr src, mp_size_t size, mp_limb_t d,
C                                mp_limb_t c)
C
C This code follows the "alternate" code in mpn/generic/mode1o.c,
C eliminating cbit+climb from the dependent chain.  This leaves,
C
C        ev4   ev5   ev6
C         1     3     1    subq   y = x - h
C        23    13     7    mulq   q = y * inverse
C        23    14     7    umulh  h = high (q * d)
C        --    --    --
C        47    30    15
C
C In each case, the load latency, loop control, and extra carry bit handling
C hide under the multiply latencies.  Those latencies are long enough that
C we don't need to worry about alignment or pairing to squeeze out
C performance.
C
C For the first limb, some of the loop code is broken out and scheduled back
C since it can be done earlier.
C
C   - The first ldq src[0] is near the start of the routine, for maximum
C     time from memory.
C
C   - The subq y=x-climb can be done without waiting for the inverse.
C
C   - The mulq y*inverse is replicated after the final subq for the inverse,
C     instead of branching to the mulq in the main loop.  On ev4 a branch
C     there would cost cycles, but we can hide them under the mulq latency.
C
C For the last limb, high<divisor is tested and if that's true a subtract
C and addback is done, as per the main mpn/generic/mode1o.c code.  This is a
C data-dependent branch, but we're waiting for umulh so any penalty should
C hide there.  The multiplies saved would be worth the cost anyway.
C
C Enhancements:
C
C For size==1, a plain division (done bitwise say) might be faster than
C calculating an inverse, the latter taking about 130 cycles on ev4 or 70 on
C ev5.  A call to gcc __remqu might be a possibility.

ASM_START()
PROLOGUE(mpn_modexact_1c_odd,gp)

	C r16	src
	C r17	size
	C r18	d
	C r19	c

	LEA(r0, binvert_limb_table)
	srl	r18, 1, r20		C d >> 1

	and	r20, 127, r20		C idx = d>>1 & 0x7F

	addq	r0, r20, r21		C table + idx

ifelse(bwx_available_p,1,
`	ldbu	r20, 0(r21)		C table[idx], inverse 8 bits
',`
	ldq_u	r20, 0(r21)		C table[idx] qword
	extbl	r20, r21, r20		C table[idx], inverse 8 bits
')

	mull	r20, r20, r7		C i*i
	addq	r20, r20, r20		C 2*i

	ldq	r2, 0(r16)		C x = s = src[0]
	lda	r17, -1(r17)		C size--
	clr	r0			C initial cbit=0

	mull	r7, r18, r7		C i*i*d

	subq	r20, r7, r20		C 2*i-i*i*d, inverse 16 bits

	mull	r20, r20, r7		C i*i
	addq	r20, r20, r20		C 2*i

	mull	r7, r18, r7		C i*i*d

	subq	r20, r7, r20		C 2*i-i*i*d, inverse 32 bits

	mulq	r20, r20, r7		C i*i
	addq	r20, r20, r20		C 2*i

	mulq	r7, r18, r7		C i*i*d
	subq	r2, r19, r3		C y = x - climb

	subq	r20, r7, r20		C inv = 2*i-i*i*d, inverse 64 bits

ASSERT(r7, C should have d*inv==1 mod 2^64
`	mulq	r18, r20, r7
	cmpeq	r7, 1, r7')

	mulq	r3, r20, r4		C first q = y * inv

	beq	r17, L(one)		C if size==1
	br	L(entry)


L(top):
	C r0	cbit
	C r16	src, incrementing
	C r17	size, decrementing
	C r18	d
	C r19	climb
	C r20	inv

	ldq	r1, 0(r16)		C s = src[i]
	subq	r1, r0, r2		C x = s - cbit
	cmpult	r1, r0, r0		C new cbit = s < cbit

	subq	r2, r19, r3		C y = x - climb

	mulq	r3, r20, r4		C q = y * inv
L(entry):
	cmpult	r2, r19, r5		C cbit2 = x < climb
	addq	r5, r0, r0		C cbit += cbit2
	lda	r16, 8(r16)		C src++
	lda	r17, -1(r17)		C size--

	umulh	r4, r18, r19		C climb = q * d
	bne	r17, L(top)		C while 2 or more limbs left



	C r0	cbit
	C r18	d
	C r19	climb
	C r20	inv

	ldq	r1, 0(r16)		C s = src[size-1] high limb

	cmpult	r1, r18, r2		C test high<divisor
	bne	r2, L(skip)		C skip if so

	C can't skip a division, repeat loop code

	subq	r1, r0, r2		C x = s - cbit
	cmpult	r1, r0, r0		C new cbit = s < cbit

	subq	r2, r19, r3		C y = x - climb

	mulq	r3, r20, r4		C q = y * inv
L(one):
	cmpult	r2, r19, r5		C cbit2 = x < climb
	addq	r5, r0, r0		C cbit += cbit2

	umulh	r4, r18, r19		C climb = q * d

	addq	r19, r0, r0		C return climb + cbit
	ret	r31, (r26), 1


	ALIGN(8)
L(skip):
	C with high<divisor, the final step can be just (cbit+climb)-s and
	C an addback of d if that underflows

	addq	r19, r0, r19		C c = climb + cbit

	subq	r19, r1, r2		C c - s
	cmpult	r19, r1, r3		C c < s

	addq	r2, r18, r0		C return c-s + divisor

	cmoveq	r3, r2, r0		C return c-s if no underflow
	ret	r31, (r26), 1

EPILOGUE()
ASM_END()