dnl  AMD64 mpn_modexact_1_odd -- Hensel norm remainder.

dnl  Copyright 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2011, 2012 Free
dnl  Software Foundation, Inc.

dnl  This file is part of the GNU MP Library.

dnl  The GNU MP Library is free software; you can redistribute it and/or modify
dnl  it under the terms of the GNU Lesser General Public License as published
dnl  by the Free Software Foundation; either version 3 of the License, or (at
dnl  your option) any later version.

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 Lesser General Public
dnl  License for more details.

dnl  You should have received a copy of the GNU Lesser General Public License
dnl  along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.

include(`../config.m4')


C	     cycles/limb
C AMD K8,K9	10
C AMD K10	10
C Intel P4	33
C Intel core2	13
C Intel corei	14.5
C Intel atom	35
C VIA nano	 ?


C The dependent chain in the main loop is
C
C                            cycles
C	sub	%rdx, %rax	1
C	imul	%r9, %rax	4
C	mul	%r8		5
C			      ----
C       total		       10
C
C The mov load from src seems to need to be scheduled back before the jz to
C achieve this speed, out-of-order execution apparently can't completely hide
C the latency otherwise.
C
C The l=src[i]-cbit step is rotated back too, since that allows us to avoid it
C for the first iteration (where there's no cbit).
C
C The code alignment used (32-byte) for the loop also seems necessary.  Without
C that the non-PIC case has adc crossing the 0x60 offset, apparently making it
C run at 11 cycles instead of 10.


ABI_SUPPORT(DOS64)
ABI_SUPPORT(STD64)

ASM_START()
	TEXT
	ALIGN(32)
PROLOGUE(mpn_modexact_1_odd)
	FUNC_ENTRY(3)
	mov	$0, R32(%rcx)
IFDOS(`	jmp	L(ent)		')

PROLOGUE(mpn_modexact_1c_odd)
	FUNC_ENTRY(4)
L(ent):
	C rdi	src
	C rsi	size
	C rdx	divisor
	C rcx	carry

	mov	%rdx, %r8		C d
	shr	R32(%rdx)		C d/2

	LEA(	binvert_limb_table, %r9)

	and	$127, R32(%rdx)
	mov	%rcx, %r10		C initial carry

	movzbl	(%r9,%rdx), R32(%rdx)	C inv 8 bits

	mov	(%rdi), %rax		C src[0]
	lea	(%rdi,%rsi,8), %r11	C src end
	mov	%r8, %rdi		C d, made available to imull

	lea	(%rdx,%rdx), R32(%rcx)	C 2*inv
	imul	R32(%rdx), R32(%rdx)	C inv*inv

	neg	%rsi			C -size

	imul	R32(%rdi), R32(%rdx)	C inv*inv*d

	sub	R32(%rdx), R32(%rcx)	C inv = 2*inv - inv*inv*d, 16 bits

	lea	(%rcx,%rcx), R32(%rdx)	C 2*inv
	imul	R32(%rcx), R32(%rcx)	C inv*inv

	imul	R32(%rdi), R32(%rcx)	C inv*inv*d

	sub	R32(%rcx), R32(%rdx)	C inv = 2*inv - inv*inv*d, 32 bits
	xor	R32(%rcx), R32(%rcx)	C initial cbit

	lea	(%rdx,%rdx), %r9	C 2*inv
	imul	%rdx, %rdx		C inv*inv

	imul	%r8, %rdx		C inv*inv*d

	sub	%rdx, %r9		C inv = 2*inv - inv*inv*d, 64 bits
	mov	%r10, %rdx		C initial climb

	ASSERT(e,`	C d*inv == 1 mod 2^64
	mov	%r8, %r10
	imul	%r9, %r10
	cmp	$1, %r10')

	inc	%rsi
	jz	L(one)


	ALIGN(16)
L(top):
	C rax	l = src[i]-cbit
	C rcx	new cbit, 0 or 1
	C rdx	climb, high of last product
	C rsi	counter, limbs, negative
	C rdi
	C r8	divisor
	C r9	inverse
	C r11	src end ptr

	sub	%rdx, %rax		C l = src[i]-cbit - climb

	adc	$0, %rcx		C more cbit
	imul	%r9, %rax		C q = l * inverse

	mul	%r8			C climb = high (q * d)

	mov	(%r11,%rsi,8), %rax	C src[i+1]
	sub	%rcx, %rax		C next l = src[i+1] - cbit
	setc	R8(%rcx)		C new cbit

	inc	%rsi
	jnz	L(top)


L(one):
	sub	%rdx, %rax		C l = src[i]-cbit - climb

	adc	$0, %rcx		C more cbit
	imul	%r9, %rax		C q = l * inverse

	mul	%r8			C climb = high (q * d)

	lea	(%rcx,%rdx), %rax	C climb+cbit
	FUNC_EXIT()
	ret

EPILOGUE(mpn_modexact_1c_odd)
EPILOGUE(mpn_modexact_1_odd)