1 dnl AMD K7 mpn_mod_1 -- mpn by limb remainder.
3 dnl K7: 17.0 cycles/limb.
6 dnl Copyright (C) 1999, 2000 Free Software Foundation, Inc.
8 dnl This file is part of the GNU MP Library.
10 dnl The GNU MP Library is free software; you can redistribute it and/or
11 dnl modify it under the terms of the GNU Lesser General Public License as
12 dnl published by the Free Software Foundation; either version 2.1 of the
13 dnl License, or (at your option) any later version.
15 dnl The GNU MP Library is distributed in the hope that it will be useful,
16 dnl but WITHOUT ANY WARRANTY; without even the implied warranty of
17 dnl MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 dnl Lesser General Public License for more details.
20 dnl You should have received a copy of the GNU Lesser General Public
21 dnl License along with the GNU MP Library; see the file COPYING.LIB. If
22 dnl not, write to the Free Software Foundation, Inc., 59 Temple Place -
23 dnl Suite 330, Boston, MA 02111-1307, USA.
26 include(`../config.m4')
29 C mp_limb_t mpn_mod_1 (mp_srcptr src, mp_size_t size, mp_limb_t divisor);
30 C mp_limb_t mpn_mod_1c (mp_srcptr src, mp_size_t size, mp_limb_t divisor,
33 C The code here is the same as mpn_divrem_1, but with the quotient
34 C discarded. See mpn/x86/k7/mmx/divrem_1.c for some comments.
37 dnl MUL_THRESHOLD is the size at which the multiply by inverse method is
38 dnl used, rather than plain "divl"s. Minimum value 2.
40 dnl The inverse takes about 50 cycles to calculate, but after that the
41 dnl multiply is 17 c/l versus division at 41 c/l.
43 dnl Using mul or div is about the same speed at 3 limbs, so the threshold
44 dnl is set to 4 to get the smaller div code used at 3.
46 deflit(MUL_THRESHOLD, 4)
49 defframe(PARAM_CARRY, 16)
50 defframe(PARAM_DIVISOR,12)
51 defframe(PARAM_SIZE, 8)
52 defframe(PARAM_SRC, 4)
54 defframe(SAVE_EBX, -4)
55 defframe(SAVE_ESI, -8)
56 defframe(SAVE_EDI, -12)
57 defframe(SAVE_EBP, -16)
59 defframe(VAR_NORM, -20)
60 defframe(VAR_INVERSE, -24)
61 defframe(VAR_SRC_STOP,-28)
63 deflit(STACK_SPACE, 28)
70 movl PARAM_CARRY, %edx
72 subl $STACK_SPACE, %esp
73 deflit(`FRAME',STACK_SPACE)
76 movl PARAM_DIVISOR, %ebp
80 jmp LF(mpn_mod_1,start_1c)
90 movl $0, %edx C initial carry (if can't skip a div)
91 subl $STACK_SPACE, %esp
92 deflit(`FRAME',STACK_SPACE)
98 movl PARAM_DIVISOR, %ebp
103 movl -4(%esi,%ecx,4), %eax C src high limb
105 cmpl %ebp, %eax C carry flag if high<divisor
107 cmovc( %eax, %edx) C src high limb as initial carry
108 sbbl $0, %ecx C size-1 to skip one div
122 cmpl $MUL_THRESHOLD, %ecx
123 jae L(mul_by_inverse)
127 C With a MUL_THRESHOLD of 4, this "loop" only ever does 1 to 3 iterations,
128 C but it's already fast and compact, and there's nothing to gain by
131 C Using PARAM_DIVISOR in the divl is a couple of cycles faster than %ebp.
138 C eax scratch (quotient)
140 C ecx counter, limbs, decrementing
141 C edx scratch (remainder)
146 movl -4(%esi,%ecx,4), %eax
157 addl $STACK_SPACE, %esp
165 C -----------------------------------------------------------------------------
176 bsrl %ebp, %eax C 31-l
181 movl %ebx, VAR_SRC_STOP
184 movl %ecx, %ebx C size
187 movl %edx, %edi C carry
195 shll %cl, %ebp C d normalized
201 subl %ebp, %edx C (b-d)-1 so edx:eax = b*(b-d)-1
203 divl %ebp C floor (b*(b-d)-1) / d
207 movl %eax, VAR_INVERSE
208 leal -12(%esi,%ebx,4), %eax C &src[size-3]
210 movl 8(%eax), %esi C src high limb
211 movl 4(%eax), %edx C src second highest limb
213 shldl( %cl, %esi, %edi) C n2 = carry,high << l
215 shldl( %cl, %edx, %esi) C n10 = high,second << l
217 movl %eax, %ecx C &src[size-3]
220 ifelse(MUL_THRESHOLD,2,`
222 je L(inverse_two_left)
226 C The dependent chain here is the same as in mpn_divrem_1, but a few
227 C instructions are saved by not needing to store the quotient limbs.
228 C Unfortunately this doesn't get the code down to the theoretical 16 c/l.
230 C There's four dummy instructions in the loop, all of which are necessary
231 C for the claimed 17 c/l. It's a 1 to 3 cycle slowdown if any are removed,
232 C or changed from load to store or vice versa. They're not completely
233 C random, since they correspond to what mpn_divrem_1 has, but there's no
234 C obvious reason why they're necessary. Presumably they induce something
235 C good in the out of order execution, perhaps through some load/store
236 C ordering and/or decoding effects.
238 C The q1==0xFFFFFFFF case is handled here the same as in mpn_divrem_1. On
239 C on special data that comes out as q1==0xFFFFFFFF always, the loop runs at
245 C ebx scratch (nadj, q1)
246 C ecx src pointer, decrementing
252 C mm0 scratch (src qword)
253 C mm7 rshift for normalization
255 cmpl $0x80000000, %esi C n1 as 0=c, 1=nc
257 movl PARAM_SIZE, %ebx C dummy
259 leal (%ebp,%esi), %ebx
260 cmovc( %esi, %ebx) C nadj = n10 + (-n1 & d), ignoring overflow
261 sbbl $-1, %eax C n2+n1
263 mull VAR_INVERSE C m*(n2+n1)
265 movq (%ecx), %mm0 C next src limb and the one below it
268 movl %ecx, PARAM_SIZE C dummy
272 addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
273 leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
278 adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
285 leal 0(%ecx), %ecx C dummy
292 movl VAR_SRC_STOP, %eax
296 sbbl %edx, %edi C n - (q1+1)*d
297 movl %esi, %edi C remainder -> n2
298 leal (%ebp,%esi), %edx
302 cmovc( %edx, %edi) C n - q1*d if underflow from using q1+1
307 L(inverse_loop_done):
310 C -----------------------------------------------------------------------------
314 C ebx scratch (nadj, q1)
321 C mm0 scratch (src dword)
324 cmpl $0x80000000, %esi C n1 as 0=c, 1=nc
327 leal (%ebp,%esi), %ebx
328 cmovc( %esi, %ebx) C nadj = n10 + (-n1 & d), ignoring overflow
329 sbbl $-1, %eax C n2+n1
331 mull VAR_INVERSE C m*(n2+n1)
333 movd 4(%ecx), %mm0 C src low limb
339 addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
340 leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
343 adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
359 sbbl %edx, %edi C n - (q1+1)*d
360 movl %esi, %edi C remainder -> n2
361 leal (%ebp,%esi), %edx
365 cmovc( %edx, %edi) C n - q1*d if underflow from using q1+1
371 C ebx scratch (nadj, q1)
378 C mm0 src limb, shifted
381 cmpl $0x80000000, %esi C n1 as 0=c, 1=nc
384 leal (%ebp,%esi), %ebx
385 cmovc( %esi, %ebx) C nadj = n10 + (-n1 & d), ignoring overflow
386 sbbl $-1, %eax C n2+n1
388 mull VAR_INVERSE C m*(n2+n1)
390 movl VAR_NORM, %ecx C for final denorm
396 addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
397 leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
402 adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
416 movl %esi, %eax C remainder
419 sbbl %edx, %edi C n - (q1+1)*d
420 leal (%ebp,%eax), %edx
423 cmovc( %edx, %eax) C n - q1*d if underflow from using q1+1
426 shrl %cl, %eax C denorm remainder
427 addl $STACK_SPACE, %esp
433 C -----------------------------------------------------------------------------
435 C Special case for q1=0xFFFFFFFF, giving q=0xFFFFFFFF meaning the low dword
436 C of q*d is simply -d and the remainder n-q*d = n10+d
447 movl VAR_SRC_STOP, %edx
448 leal (%ebp,%esi), %edi C n-q*d remainder -> next n2
451 movd %mm0, %esi C next n10
455 jmp L(inverse_loop_done)