dnl Intel P5 mpn_mod_1 -- mpn by limb remainder.
dnl Copyright 1999, 2000, 2002 Free Software Foundation, Inc.
dnl
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
dnl modify it under the terms of the GNU Lesser General Public License as
dnl published by the Free Software Foundation; either version 2.1 of the
dnl License, or (at your option) any later version.
dnl
dnl The GNU MP Library is distributed in the hope that it will be useful,
dnl but WITHOUT ANY WARRANTY; without even the implied warranty of
dnl MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
dnl Lesser General Public License for more details.
dnl
dnl You should have received a copy of the GNU Lesser General Public
dnl License along with the GNU MP Library; see the file COPYING.LIB. If
dnl not, write to the Free Software Foundation, Inc., 59 Temple Place -
dnl Suite 330, Boston, MA 02111-1307, USA.
include(`../config.m4')
C P5: 28.0 cycles/limb
C mp_limb_t mpn_mod_1 (mp_srcptr src, mp_size_t size, mp_limb_t divisor);
C mp_limb_t mpn_mod_1c (mp_srcptr src, mp_size_t size, mp_limb_t divisor,
C mp_limb_t carry);
C mp_limb_t mpn_preinv_mod_1 (mp_srcptr src, mp_size_t size, mp_limb_t divisor,
C mp_limb_t inverse);
C
C This code is not unlike mpn/x86/p6/mod_1.asm, it does the same sort of
C multiply by inverse without on-the-fly shifts. See that code for some
C general comments.
C
C Alternatives:
C
C P5 shldl is 4 cycles, so shifting on the fly would be at least 5 cycles
C slower, probably more depending what it did to register usage. Using MMX
C on P55 would be better, but still at least 4 or 5 instructions and so 2 or
C 3 cycles.
dnl These thresholds are the sizes where the multiply by inverse method is
dnl used, rather than plain "divl"s. Minimum value 2.
dnl
dnl MUL_NORM_THRESHOLD is for an already normalized divisor (high bit set),
dnl MUL_UNNORM_THRESHOLD for an unnormalized divisor.
dnl
dnl With the divl loop at 44 c/l and the inverse at 28 c/l with about 70
dnl cycles to setup, the threshold should be about ceil(70/16)==5, which is
dnl what happens in practice.
dnl
dnl An unnormalized divisor gets an extra 40 cycles at the end for the
dnl final (r*2^n)%(d*2^n) and shift. This increases the threshold by about
dnl 40/16=3.
dnl
dnl PIC adds between 4 and 7 cycles (not sure why it varies), but this
dnl doesn't change the thresholds.
dnl
dnl The entry sequence code that chooses between MUL_NORM_THRESHOLD and
dnl MUL_UNNORM_THRESHOLD is a bit horrible, but it adds only 2 cycles
dnl (branch free) and ensures the choice between div or mul is optimal.
deflit(MUL_NORM_THRESHOLD, ifdef(`PIC',5,5))
deflit(MUL_UNNORM_THRESHOLD, ifdef(`PIC',8,8))
deflit(MUL_NORM_DELTA, eval(MUL_NORM_THRESHOLD - MUL_UNNORM_THRESHOLD))
defframe(PARAM_INVERSE, 16) dnl mpn_preinv_mod_1
defframe(PARAM_CARRY, 16) dnl mpn_mod_1c
defframe(PARAM_DIVISOR, 12)
defframe(PARAM_SIZE, 8)
defframe(PARAM_SRC, 4)
dnl re-using parameter space
define(VAR_NORM, `PARAM_DIVISOR')
define(VAR_INVERSE, `PARAM_SIZE')
TEXT
ALIGN(8)
PROLOGUE(mpn_preinv_mod_1)
deflit(`FRAME',0)
pushl %ebp FRAME_pushl()
pushl %esi FRAME_pushl()
movl PARAM_SRC, %esi
movl PARAM_SIZE, %edx
pushl %edi FRAME_pushl()
pushl %ebx FRAME_pushl()
movl PARAM_DIVISOR, %ebp
movl PARAM_INVERSE, %eax
movl -4(%esi,%edx,4), %edi C src high limb
leal -8(%esi,%edx,4), %esi C &src[size-2]
movl $0, VAR_NORM
decl %edx
jnz L(start_preinv)
subl %ebp, %edi C src-divisor
popl %ebx
sbbl %ecx, %ecx C -1 if underflow
movl %edi, %eax C src-divisor
andl %ebp, %ecx C d if underflow
popl %edi
addl %ecx, %eax C remainder, with possible addback
popl %esi
popl %ebp
ret
EPILOGUE()
ALIGN(8)
PROLOGUE(mpn_mod_1c)
deflit(`FRAME',0)
movl PARAM_DIVISOR, %eax
movl PARAM_SIZE, %ecx
sarl $31, %eax C d highbit
movl PARAM_CARRY, %edx
orl %ecx, %ecx
jz L(done_edx) C result==carry if size==0
andl $MUL_NORM_DELTA, %eax
pushl %ebp FRAME_pushl()
addl $MUL_UNNORM_THRESHOLD, %eax C norm or unnorm thresh
pushl %esi FRAME_pushl()
movl PARAM_SRC, %esi
movl PARAM_DIVISOR, %ebp
cmpl %eax, %ecx
jb L(divide_top)
movl %edx, %eax C carry as pretend src high limb
leal 1(%ecx), %edx C size+1
cmpl $0x1000000, %ebp
jmp L(mul_by_inverse_1c)
EPILOGUE()
ALIGN(8)
PROLOGUE(mpn_mod_1)
deflit(`FRAME',0)
movl PARAM_SIZE, %ecx
pushl %ebp FRAME_pushl()
orl %ecx, %ecx
jz L(done_zero)
movl PARAM_SRC, %eax
movl PARAM_DIVISOR, %ebp
sarl $31, %ebp C -1 if divisor normalized
movl -4(%eax,%ecx,4), %eax C src high limb
movl PARAM_DIVISOR, %edx
pushl %esi FRAME_pushl()
andl $MUL_NORM_DELTA, %ebp
cmpl %edx, %eax C carry flag if high<divisor
sbbl %edx, %edx C -1 if high<divisor
addl $MUL_UNNORM_THRESHOLD, %ebp C norm or unnorm thresh
addl %edx, %ecx C size-1 if high<divisor
jz L(done_eax)
cmpl %ebp, %ecx
movl PARAM_DIVISOR, %ebp
movl PARAM_SRC, %esi
jae L(mul_by_inverse)
andl %eax, %edx C high as initial carry if high<divisor
L(divide_top):
C eax scratch (quotient)
C ebx
C ecx counter, limbs, decrementing
C edx scratch (remainder)
C esi src
C edi
C ebp divisor
movl -4(%esi,%ecx,4), %eax
divl %ebp
decl %ecx
jnz L(divide_top)
popl %esi
popl %ebp
L(done_edx):
movl %edx, %eax
ret
L(done_zero):
xorl %eax, %eax
popl %ebp
ret
C -----------------------------------------------------------------------------
C
C The divisor is normalized using the same code as the pentium
C count_leading_zeros in longlong.h. Going through the GOT for PIC costs a
C couple of cycles, but is more or less unavoidable.
ALIGN(8)
L(mul_by_inverse):
C eax src high limb
C ebx
C ecx size or size-1
C edx
C esi src
C edi
C ebp divisor
movl PARAM_SIZE, %edx
cmpl $0x1000000, %ebp
L(mul_by_inverse_1c):
sbbl %ecx, %ecx
cmpl $0x10000, %ebp
sbbl $0, %ecx
cmpl $0x100, %ebp
sbbl $0, %ecx
pushl %edi FRAME_pushl()
pushl %ebx FRAME_pushl()
movl %ebp, %ebx C d
ifdef(`PIC',`
call L(here)
L(here):
popl %edi
leal 25(,%ecx,8), %ecx C 0,-1,-2,-3 -> 25,17,9,1
shrl %cl, %ebx
addl $_GLOBAL_OFFSET_TABLE_+[.-L(here)], %edi
C AGI
movl __clz_tab@GOT(%edi), %edi
addl $-34, %ecx
C AGI
movb (%ebx,%edi), %bl
',`
leal 25(,%ecx,8), %ecx C 0,-1,-2,-3 -> 25,17,9,1
shrl %cl, %ebx
addl $-34, %ecx
C AGI
movb __clz_tab(%ebx), %bl
')
movl %eax, %edi C carry -> n1
addl %ebx, %ecx C -34 + c + __clz_tab[d>>c] = -clz-1
leal -8(%esi,%edx,4), %esi C &src[size-2]
xorl $-1, %ecx C clz
movl $-1, %edx
ASSERT(e,`pushl %eax C clz calculation same as bsrl
bsrl %ebp, %eax
xorl $31, %eax
cmpl %eax, %ecx
popl %eax')
shll %cl, %ebp C d normalized
movl %ecx, VAR_NORM
subl %ebp, %edx C (b-d)-1, so edx:eax = b*(b-d)-1
movl $-1, %eax
divl %ebp C floor (b*(b-d)-1) / d
L(start_preinv):
movl %eax, VAR_INVERSE
movl %ebp, %eax C d
movl %ecx, %edx C fake high, will cancel
C For mpn_mod_1 and mpn_preinv_mod_1, the initial carry in %edi is the src
C high limb, and this may be greater than the divisor and may need one copy
C of the divisor subtracted (only one, because the divisor is normalized).
C This is accomplished by having the initial ecx:edi act as a fake previous
C n2:n10. The initial edx:eax is d, acting as a fake (q1+1)*d which is
C subtracted from ecx:edi, with the usual addback if it produces an
C underflow.
L(inverse_top):
C eax scratch (n10, n1, q1, etc)
C ebx scratch (nadj, src limit)
C ecx old n2
C edx scratch
C esi src pointer, &src[size-2] to &src[0]
C edi old n10
C ebp d
subl %eax, %edi C low n - (q1+1)*d
movl (%esi), %eax C new n10
sbbl %edx, %ecx C high n - (q1+1)*d, 0 or -1
movl %ebp, %ebx C d
sarl $31, %eax C -n1
andl %ebp, %ecx C d if underflow
addl %edi, %ecx C remainder -> n2, and possible addback
ASSERT(b,`cmpl %ebp, %ecx')
andl %eax, %ebx C -n1 & d
movl (%esi), %edi C n10
andl $1, %eax C n1
addl %ecx, %eax C n2+n1
addl %edi, %ebx C nadj = n10 + (-n1 & d), ignoring overflow
mull VAR_INVERSE C m*(n2+n1)
addl %eax, %ebx C low(m*(n2+n1) + nadj), giving carry flag
leal 1(%ecx), %eax C 1+n2
adcl %edx, %eax C 1 + high[n2<<32 + m*(n2+n1) + nadj] = q1+1
movl PARAM_SRC, %ebx
sbbl $0, %eax C use q1 if q1+1 overflows
subl $4, %esi C step src ptr
mull %ebp C (q1+1)*d
cmpl %ebx, %esi
jae L(inverse_top)
C %edi (after subtract and addback) is the remainder modulo d*2^n
C and must be reduced to 0<=r<d by calculating r*2^n mod d*2^n and
C right shifting by n.
C
C If d was already normalized on entry so that n==0 then nothing is
C needed here. This is always the case for preinv_mod_1. For mod_1
C or mod_1c the chance of n==0 is low, but about 40 cycles can be
C saved.
subl %eax, %edi C low n - (q1+1)*d
movl %ecx, %ebx C n2
sbbl %edx, %ebx C high n - (q1+1)*d, 0 or -1
xorl %esi, %esi C next n2
andl %ebp, %ebx C d if underflow
movl VAR_NORM, %ecx
addl %ebx, %edi C remainder, with possible addback
orl %ecx, %ecx
jz L(done_mul_edi)
C Here using %esi=n2 and %edi=n10, unlike the above
shldl( %cl, %edi, %esi) C n2
shll %cl, %edi C n10
movl %edi, %eax C n10
movl %edi, %ebx C n10
sarl $31, %ebx C -n1
shrl $31, %eax C n1
andl %ebp, %ebx C -n1 & d
addl %esi, %eax C n2+n1
addl %edi, %ebx C nadj = n10 + (-n1 & d), ignoring overflow
mull VAR_INVERSE C m*(n2+n1)
addl %eax, %ebx C m*(n2+n1) + nadj, low giving carry flag
leal 1(%esi), %eax C 1+n2
adcl %edx, %eax C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
sbbl $0, %eax C use q1 if q1+1 overflows
mull %ebp C (q1+1)*d
subl %eax, %edi C low n - (q1+1)*d
popl %ebx
sbbl %edx, %esi C high n - (q1+1)*d, 0 or -1
movl %edi, %eax
andl %ebp, %esi C d if underflow
popl %edi
addl %esi, %eax C addback if underflow
popl %esi
shrl %cl, %eax C denorm remainder
popl %ebp
ret
L(done_mul_edi):
movl %edi, %eax
popl %ebx
popl %edi
L(done_eax):
popl %esi
popl %ebp
ret
EPILOGUE()
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mod_1.asm CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / x86 / pentium