File: [local] / OpenXM_contrib / gmp / mpn / x86 / p6 / mmx / Attic / mod_1.asm (download)
Revision 1.1.1.1 (vendor branch), Sat Sep 9 14:12:44 2000 UTC (23 years, 10 months ago) by maekawa
Branch: GMP
CVS Tags: maekawa-ipv6, VERSION_3_1_1, VERSION_3_1, RELEASE_1_2_2, RELEASE_1_2_1, RELEASE_1_1_3
Changes since 1.1: +0 -0
lines
Import gmp 3.1
|
dnl Intel Pentium-II mpn_mod_1 -- mpn by limb remainder.
dnl
dnl P6MMX: 24.0 cycles/limb.
dnl Copyright (C) 1999, 2000 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 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
C The code here very similar to mpn_divrem_1, but with the quotient
C discarded. What's here probably isn't optimal.
C
C See mpn/x86/p6/mmx/divrem_1.c and mpn/x86/k7/mmx/mod_1.asm for some
C comments.
dnl MUL_THRESHOLD is the size at which the multiply by inverse method is
dnl used, rather than plain "divl"s. Minimum value 2.
deflit(MUL_THRESHOLD, 4)
defframe(PARAM_CARRY, 16)
defframe(PARAM_DIVISOR,12)
defframe(PARAM_SIZE, 8)
defframe(PARAM_SRC, 4)
defframe(SAVE_EBX, -4)
defframe(SAVE_ESI, -8)
defframe(SAVE_EDI, -12)
defframe(SAVE_EBP, -16)
defframe(VAR_NORM, -20)
defframe(VAR_INVERSE, -24)
defframe(VAR_SRC_STOP,-28)
deflit(STACK_SPACE, 28)
.text
ALIGN(16)
PROLOGUE(mpn_mod_1c)
deflit(`FRAME',0)
movl PARAM_CARRY, %edx
movl PARAM_SIZE, %ecx
subl $STACK_SPACE, %esp
deflit(`FRAME',STACK_SPACE)
movl %ebp, SAVE_EBP
movl PARAM_DIVISOR, %ebp
movl %esi, SAVE_ESI
movl PARAM_SRC, %esi
jmp LF(mpn_mod_1,start_1c)
EPILOGUE()
ALIGN(16)
PROLOGUE(mpn_mod_1)
deflit(`FRAME',0)
movl $0, %edx C initial carry (if can't skip a div)
movl PARAM_SIZE, %ecx
subl $STACK_SPACE, %esp
deflit(`FRAME',STACK_SPACE)
movl %esi, SAVE_ESI
movl PARAM_SRC, %esi
movl %ebp, SAVE_EBP
movl PARAM_DIVISOR, %ebp
orl %ecx, %ecx
jz L(divide_done)
movl -4(%esi,%ecx,4), %eax C src high limb
cmpl %ebp, %eax C carry flag if high<divisor
cmovc( %eax, %edx) C src high limb as initial carry
sbbl $0, %ecx C size-1 to skip one div
jz L(divide_done)
ALIGN(16)
L(start_1c):
C eax
C ebx
C ecx size
C edx carry
C esi src
C edi
C ebp divisor
cmpl $MUL_THRESHOLD, %ecx
jae L(mul_by_inverse)
orl %ecx, %ecx
jz L(divide_done)
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
movl -4(%esi,%ecx,4), %eax
divl %ebp
decl %ecx
jnz L(divide_top)
L(divide_done):
movl SAVE_ESI, %esi
movl %edx, %eax
movl SAVE_EBP, %ebp
addl $STACK_SPACE, %esp
ret
C -----------------------------------------------------------------------------
L(mul_by_inverse):
C eax
C ebx
C ecx size
C edx carry
C esi src
C edi
C ebp divisor
movl %ebx, SAVE_EBX
leal -4(%esi), %ebx
movl %ebx, VAR_SRC_STOP
movl %ecx, %ebx C size
movl %edi, SAVE_EDI
movl %edx, %edi C carry
bsrl %ebp, %ecx C 31-l
movl $-1, %edx
leal 1(%ecx), %eax C 32-l
xorl $31, %ecx C l
movl %ecx, VAR_NORM
shll %cl, %ebp C d normalized
movd %eax, %mm7
movl $-1, %eax
subl %ebp, %edx C (b-d)-1 so edx:eax = b*(b-d)-1
divl %ebp C floor (b*(b-d)-1) / d
C
movl %eax, VAR_INVERSE
leal -12(%esi,%ebx,4), %eax C &src[size-3]
movl 8(%eax), %esi C src high limb
movl 4(%eax), %edx C src second highest limb
shldl( %cl, %esi, %edi) C n2 = carry,high << l
shldl( %cl, %edx, %esi) C n10 = high,second << l
movl %eax, %ecx C &src[size-3]
ifelse(MUL_THRESHOLD,2,`
cmpl $2, %ebx
je L(inverse_two_left)
')
C The dependent chain here is the same as in mpn_divrem_1, but a few
C instructions are saved by not needing to store the quotient limbs. This
C gets it down to 24 c/l, which is still a bit away from a theoretical 19
C c/l.
ALIGN(16)
L(inverse_top):
C eax scratch
C ebx scratch (nadj, q1)
C ecx src pointer, decrementing
C edx scratch
C esi n10
C edi n2
C ebp divisor
C
C mm0 scratch (src qword)
C mm7 rshift for normalization
movl %esi, %eax
movl %ebp, %ebx
sarl $31, %eax C -n1
andl %eax, %ebx C -n1 & d
negl %eax C n1
addl %esi, %ebx C nadj = n10 + (-n1 & d), ignoring overflow
addl %edi, %eax C n2+n1
mull VAR_INVERSE C m*(n2+n1)
movq (%ecx), %mm0 C next src limb and the one below it
subl $4, %ecx
C
C
C
addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
movl %ebp, %eax C d
adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
jz L(q1_ff)
mull %ebx C (q1+1)*d
psrlq %mm7, %mm0
movl VAR_SRC_STOP, %ebx
C
C
C
subl %eax, %esi
sbbl %edx, %edi C n - (q1+1)*d
movl %esi, %edi C remainder -> n2
leal (%ebp,%esi), %edx
cmovc( %edx, %edi) C n - q1*d if underflow from using q1+1
movd %mm0, %esi
cmpl %ebx, %ecx
jne L(inverse_top)
L(inverse_loop_done):
C -----------------------------------------------------------------------------
L(inverse_two_left):
C eax scratch
C ebx scratch (nadj, q1)
C ecx &src[-1]
C edx scratch
C esi n10
C edi n2
C ebp divisor
C
C mm0 scratch (src dword)
C mm7 rshift
movl %esi, %eax
movl %ebp, %ebx
sarl $31, %eax C -n1
andl %eax, %ebx C -n1 & d
negl %eax C n1
addl %esi, %ebx C nadj = n10 + (-n1 & d), ignoring overflow
addl %edi, %eax C n2+n1
mull VAR_INVERSE C m*(n2+n1)
movd 4(%ecx), %mm0 C src low limb
C
C
C
addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
sbbl $0, %ebx
movl %ebp, %eax C d
mull %ebx C (q1+1)*d
psllq $32, %mm0
psrlq %mm7, %mm0
C
C
subl %eax, %esi
sbbl %edx, %edi C n - (q1+1)*d
movl %esi, %edi C remainder -> n2
leal (%ebp,%esi), %edx
cmovc( %edx, %edi) C n - q1*d if underflow from using q1+1
movd %mm0, %esi
C One limb left
C eax scratch
C ebx scratch (nadj, q1)
C ecx
C edx scratch
C esi n10
C edi n2
C ebp divisor
C
C mm0 src limb, shifted
C mm7 rshift
movl %esi, %eax
movl %ebp, %ebx
sarl $31, %eax C -n1
andl %eax, %ebx C -n1 & d
negl %eax C n1
addl %esi, %ebx C nadj = n10 + (-n1 & d), ignoring overflow
addl %edi, %eax C n2+n1
mull VAR_INVERSE C m*(n2+n1)
movl VAR_NORM, %ecx C for final denorm
C
C
C
addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
sbbl $0, %ebx
movl %ebp, %eax C d
mull %ebx C (q1+1)*d
movl SAVE_EBX, %ebx
C
C
C
subl %eax, %esi
sbbl %edx, %edi C n - (q1+1)*d
leal (%ebp,%esi), %edx
movl SAVE_EBP, %ebp
movl %esi, %eax C remainder
movl SAVE_ESI, %esi
cmovc( %edx, %eax) C n - q1*d if underflow from using q1+1
movl SAVE_EDI, %edi
shrl %cl, %eax C denorm remainder
addl $STACK_SPACE, %esp
emms
ret
C -----------------------------------------------------------------------------
C
C Special case for q1=0xFFFFFFFF, giving q=0xFFFFFFFF meaning the low dword
C of q*d is simply -d and the remainder n-q*d = n10+d
L(q1_ff):
C eax (divisor)
C ebx (q1+1 == 0)
C ecx src pointer
C edx
C esi n10
C edi (n2)
C ebp divisor
leal (%ebp,%esi), %edi C n-q*d remainder -> next n2
movl VAR_SRC_STOP, %edx
psrlq %mm7, %mm0
movd %mm0, %esi C next n10
cmpl %ecx, %edx
jne L(inverse_top)
jmp L(inverse_loop_done)
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 / p6 / mmx