dnl  Itanium-2 mpn_modexact_1c_odd -- mpn by 1 exact remainder.

dnl  Contributed to the GNU project by Kevin Ryde.

dnl  Copyright 2003-2005 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 Itanium:      15
C Itanium 2:     8


dnl  Usage: ABI32(`code')
dnl
dnl  Emit the given code only under HAVE_ABI_32.
dnl
define(ABI32,
m4_assert_onearg()
`ifdef(`HAVE_ABI_32',`$1')')


C mp_limb_t mpn_modexact_1c_odd (mp_srcptr src, mp_size_t size,
C                                mp_limb_t divisor, mp_limb_t carry);
C
C The modexact algorithm is usually conceived as a dependent chain
C
C	l = src[i] - c
C	q = low(l * inverse)
C	c = high(q*divisor) + (src[i]<c)
C
C but we can work the src[i]-c into an xma by calculating si=src[i]*inverse
C separately (off the dependent chain) and using
C
C	q = low(c * inverse + si)
C	c = high(q*divisor + c)
C
C This means the dependent chain is simply xma.l followed by xma.hu, for a
C total 8 cycles/limb on itanium-2.
C
C The reason xma.hu works for the new c is that the low of q*divisor is
C src[i]-c (being the whole purpose of the q generated, and it can be
C verified algebraically).  If there was an underflow from src[i]-c, then
C there will be an overflow from (src-c)+c, thereby adding 1 to the new c
C the same as the borrow bit (src[i]<c) gives in the first style shown.
C
C Incidentally, fcmp is not an option for treating src[i]-c, since it
C apparently traps to the kernel for unnormalized operands like those used
C and generated by ldf8 and xma.  On one GNU/Linux system it took about 1200
C cycles.
C
C
C First Limb:
C
C The first limb uses q = (src[0]-c) * inverse shown in the first style.
C This lets us get the first q as soon as the inverse is ready, without
C going through si=s*inverse.  Basically at the start we have c and can use
C it while waiting for the inverse, whereas for the second and subsequent
C limbs it's the other way around, ie. we have the inverse and are waiting
C for c.
C
C At .Lentry the first two instructions in the loop have been done already.
C The load of f11=src[1] at the start (predicated on size>=2), and the
C calculation of q by the initial different scheme.
C
C
C Entry Sequence:
C
C In the entry sequence, the critical path is the calculation of the
C inverse, so this is begun first and optimized.  Apart from that, ar.lc is
C established nice and early so the br.cloop's should predict perfectly.
C And the load for the low limbs src[0] and src[1] can be initiated long
C ahead of where they're needed.
C
C
C Inverse Calculation:
C
C The initial 8-bit inverse is calculated using a table lookup.  If it hits
C L1 (which is likely if we're called several times) then it should take a
C total 4 cycles, otherwise hopefully L2 for 9 cycles.  This is considered
C the best approach, on balance.  It could be done bitwise, but that would
C probably be about 14 cycles (2 per bit beyond the first couple).  Or it
C could be taken from 4 bits to 8 with xmpy doubling as used beyond 8 bits,
C but that would be about 11 cycles.
C
C The table is not the same as binvert_limb_table, instead it's 256 bytes,
C designed to be indexed by the low byte of the divisor.  The divisor is
C always odd, so the relevant data is every second byte in the table.  The
C padding lets us use zxt1 instead of extr.u, the latter would cost an extra
C cycle because it must go down I0, and we're using the first I0 slot to get
C ip.  The extra 128 bytes of padding should be insignificant compared to
C typical ia64 code bloat.
C
C Having the table in .text allows us to use IP-relative addressing,
C avoiding a fetch from ltoff.  .rodata is apparently not suitable for use
C IP-relative, it gets a linker relocation overflow on GNU/Linux.
C
C
C Load Scheduling:
C
C In the main loop, the data loads are scheduled for an L2 hit, which means
C 6 cycles for the data ready to use.  In fact we end up 7 cycles ahead.  In
C any case that scheduling is achieved simply by doing the load (and xmpy.l
C for "si") in the immediately preceding iteration.
C
C The main loop requires size >= 2, and we handle size==1 by an initial
C br.cloop to enter the loop only if size>1.  Since ar.lc is established
C early, this should predict perfectly.
C
C
C Not done:
C
C Consideration was given to using a plain "(src[0]-c) % divisor" for
C size==1, but cycle counting suggests about 50 for the sort of approach
C taken by gcc __umodsi3, versus about 47 for the modexact.  (Both assuming
C L1 hits for their respective fetching.)
C
C Consideration was given to a test for high<divisor and replacing the last
C loop iteration with instead c-=src[size-1] followed by c+=d if underflow.
C Branching on high<divisor wouldn't be good since a mispredict would cost
C more than the loop iteration saved, and the condition is of course data
C dependent.  So the theory would be to shorten the loop count if
C high<divisor, and predicate extra operations at the end.  That would mean
C a gain of 6 when high<divisor, or a cost of 2 if not.
C
C Whether such a tradeoff is a win on average depends on assumptions about
C how many bits in the high and the divisor.  If both are uniformly
C distributed then high<divisor about 50% of the time.  But smallish
C divisors (less chance of high<divisor) might be more likely from
C applications (mpz_divisible_ui, mpz_gcd_ui, etc).  Though biggish divisors
C would be normal internally from say mpn/generic/perfsqr.c.  On balance,
C for the moment, it's felt the gain is not really enough to be worth the
C trouble.
C
C
C Enhancement:
C
C Process two source limbs per iteration using a two-limb inverse and a
C sequence like
C
C	ql  = low (c * il + sil)	quotient low limb
C	qlc = high(c * il + sil)
C	qh1 = low (c * ih + sih)	quotient high, partial
C
C	cl = high (ql * d + c)		carry out of low
C	qh = low (qlc * 1 + qh1)	quotient high limb
C
C	new c = high (qh * d + cl)	carry out of high
C
C This would be 13 cycles/iteration, giving 6.5 cycles/limb.  The two limb
C s*inverse as sih:sil = sh:sl * ih:il would be calculated off the dependent
C chain with 4 multiplies.  The bigger inverse would take extra time to
C calculate, but a one limb iteration to handle an odd size could be done as
C soon as 64-bits of inverse were ready.
C
C Perhaps this could even extend to a 3 limb inverse, which might promise 17
C or 18 cycles for 3 limbs, giving 5.66 or 6.0 cycles/limb.
C

ASM_START()
	.explicit

	.text
	.align	32
.Ltable:
data1	0,0x01, 0,0xAB, 0,0xCD, 0,0xB7, 0,0x39, 0,0xA3, 0,0xC5, 0,0xEF
data1	0,0xF1, 0,0x1B, 0,0x3D, 0,0xA7, 0,0x29, 0,0x13, 0,0x35, 0,0xDF
data1	0,0xE1, 0,0x8B, 0,0xAD, 0,0x97, 0,0x19, 0,0x83, 0,0xA5, 0,0xCF
data1	0,0xD1, 0,0xFB, 0,0x1D, 0,0x87, 0,0x09, 0,0xF3, 0,0x15, 0,0xBF
data1	0,0xC1, 0,0x6B, 0,0x8D, 0,0x77, 0,0xF9, 0,0x63, 0,0x85, 0,0xAF
data1	0,0xB1, 0,0xDB, 0,0xFD, 0,0x67, 0,0xE9, 0,0xD3, 0,0xF5, 0,0x9F
data1	0,0xA1, 0,0x4B, 0,0x6D, 0,0x57, 0,0xD9, 0,0x43, 0,0x65, 0,0x8F
data1	0,0x91, 0,0xBB, 0,0xDD, 0,0x47, 0,0xC9, 0,0xB3, 0,0xD5, 0,0x7F
data1	0,0x81, 0,0x2B, 0,0x4D, 0,0x37, 0,0xB9, 0,0x23, 0,0x45, 0,0x6F
data1	0,0x71, 0,0x9B, 0,0xBD, 0,0x27, 0,0xA9, 0,0x93, 0,0xB5, 0,0x5F
data1	0,0x61, 0,0x0B, 0,0x2D, 0,0x17, 0,0x99, 0,0x03, 0,0x25, 0,0x4F
data1	0,0x51, 0,0x7B, 0,0x9D, 0,0x07, 0,0x89, 0,0x73, 0,0x95, 0,0x3F
data1	0,0x41, 0,0xEB, 0,0x0D, 0,0xF7, 0,0x79, 0,0xE3, 0,0x05, 0,0x2F
data1	0,0x31, 0,0x5B, 0,0x7D, 0,0xE7, 0,0x69, 0,0x53, 0,0x75, 0,0x1F
data1	0,0x21, 0,0xCB, 0,0xED, 0,0xD7, 0,0x59, 0,0xC3, 0,0xE5, 0,0x0F
data1	0,0x11, 0,0x3B, 0,0x5D, 0,0xC7, 0,0x49, 0,0x33, 0,0x55, 0,0xFF


PROLOGUE(mpn_modexact_1c_odd)

	C r32	src
	C r33	size
	C r34	divisor
	C r35	carry

	.prologue
.Lhere:
{ .mmi;	add	r33 = -1, r33		C M0  size-1
	mov	r14 = 2			C M1  2
	mov	r15 = ip		C I0  .Lhere
}{.mmi;	setf.sig f6 = r34		C M2  divisor
	setf.sig f9 = r35		C M3  carry
	zxt1	r3 = r34		C I1  divisor low byte
}	;;

{ .mmi;	add	r3 = .Ltable-.Lhere, r3	C M0  table offset ip and index
	sub	r16 = 0, r34		C M1  -divisor
	.save	ar.lc, r2
	mov	r2 = ar.lc		C I0
}{.mmi;	.body
	setf.sig f13 = r14		C M2  2 in significand
	mov	r17 = -1		C M3  -1
ABI32(`	zxt4	r33 = r33')		C I1  size extend
}	;;

{ .mmi;	add	r3 = r3, r15		C M0  table entry address
ABI32(` addp4	r32 = 0, r32')		C M1  src extend
	mov	ar.lc = r33		C I0  size-1 loop count
}{.mmi;	setf.sig f12 = r16		C M2  -divisor
	setf.sig f8 = r17		C M3  -1
}	;;

{ .mmi;	ld1	r3 = [r3]		C M0  inverse, 8 bits
	ldf8	f10 = [r32], 8		C M1  src[0]
	cmp.ne	p6,p0 = 0, r33		C I0  test size!=1
}	;;

	C Wait for table load.
	C Hope for an L1 hit of 1 cycles to ALU, but could be more.
	setf.sig f7 = r3		C M2  inverse, 8 bits
(p6)	ldf8	f11 = [r32], 8		C M1  src[1], if size!=1
	;;

	C 5 cycles

	C f6	divisor
	C f7	inverse, being calculated
	C f8	-1, will be -inverse
	C f9	carry
	C f10	src[0]
	C f11	src[1]
	C f12	-divisor
	C f13	2
	C f14	scratch

	xmpy.l	f14 = f13, f7		C 2*i
	xmpy.l	f7 = f7, f7		C i*i
	;;
	xma.l	f7 = f7, f12, f14	C i*i*-d + 2*i, inverse 16 bits
	;;

	xmpy.l	f14 = f13, f7		C 2*i
	xmpy.l	f7 = f7, f7		C i*i
	;;
	xma.l	f7 = f7, f12, f14	C i*i*-d + 2*i, inverse 32 bits
	;;

	xmpy.l	f14 = f13, f7		C 2*i
	xmpy.l	f7 = f7, f7		C i*i
	;;

	xma.l	f7 = f7, f12, f14	C i*i*-d + 2*i, inverse 64 bits
	xma.l	f10 = f9, f8, f10	C sc = c * -1 + src[0]
	;;
ASSERT(p6, `
	xmpy.l	f15 = f6, f7 ;;	C divisor*inverse
	getf.sig r31 = f15 ;;
	cmp.eq	p6,p0 = 1, r31	C should == 1
')

	xmpy.l	f10 = f10, f7		C q = sc * inverse
	xmpy.l	f8 = f7, f8		C -inverse = inverse * -1
	br.cloop.sptk.few.clr .Lentry	C main loop, if size > 1
	;;

	C size==1, finish up now
	xma.hu	f9 = f10, f6, f9	C c = high(q * divisor + c)
	mov	ar.lc = r2		C I0
	;;
	getf.sig r8 = f9		C M2  return c
	br.ret.sptk.many b0



.Ltop:
	C r2	saved ar.lc
	C f6	divisor
	C f7	inverse
	C f8	-inverse
	C f9	carry
	C f10	src[i] * inverse
	C f11	scratch src[i+1]

	add	r16 = 160, r32
	ldf8	f11 = [r32], 8		C src[i+1]
	;;
	C 2 cycles

	lfetch	[r16]
	xma.l	f10 = f9, f8, f10	C q = c * -inverse + si
	;;
	C 3 cycles

.Lentry:
	xma.hu	f9 = f10, f6, f9	C c = high(q * divisor + c)
	xmpy.l	f10 = f11, f7		C si = src[i] * inverse
	br.cloop.sptk.few.clr .Ltop
	;;



	xma.l	f10 = f9, f8, f10	C q = c * -inverse + si
	mov	ar.lc = r2		C I0
	;;
	xma.hu	f9 = f10, f6, f9	C c = high(q * divisor + c)
	;;
	getf.sig r8 = f9		C M2  return c
	br.ret.sptk.many b0

EPILOGUE()