数学计算

开发平台：
Unix_Linux

dive_1.asm：源码内容
							dnl  AMD K6 mpn_divexact_1 -- mpn by limb exact division.
dnl  Copyright 2000, 2001, 2002, 2007 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 3 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 License
dnl  along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.
include(`../config.m4')
C         divisor
C       odd   even
C K6:   10.0  12.0  cycles/limb
C K6-2: 10.0  11.5
C void mpn_divexact_1 (mp_ptr dst, mp_srcptr src, mp_size_t size,
C                      mp_limb_t divisor);
C
C A simple divl is used for size==1.  This is about 10 cycles faster for an
C odd divisor or 20 cycles for an even divisor.
C
C The loops are quite sensitive to code alignment, speeds should be
C rechecked (odd and even divisor, pic and non-pic) if contemplating
C changing anything.
defframe(PARAM_DIVISOR,16)
defframe(PARAM_SIZE,   12)
defframe(PARAM_SRC,    8)
defframe(PARAM_DST,    4)
dnl  re-use parameter space
define(VAR_INVERSE,`PARAM_DST')
	TEXT
	ALIGN(32)
PROLOGUE(mpn_divexact_1)
deflit(`FRAME',0)
	movl	PARAM_SIZE, %ecx
	movl	PARAM_SRC, %eax
	xorl	%edx, %edx
	cmpl	$1, %ecx
	jnz	L(two_or_more)
	movl	(%eax), %eax
	divl	PARAM_DIVISOR
	movl	PARAM_DST, %ecx
	movl	%eax, (%ecx)
	ret
L(two_or_more):
	movl	PARAM_DIVISOR, %eax
	pushl	%ebx		FRAME_pushl()
	movl	PARAM_SRC, %ebx
	pushl	%ebp		FRAME_pushl()
L(strip_twos):
	shrl	%eax
	incl	%edx			C will get shift+1
	jnc	L(strip_twos)
	pushl	%esi		FRAME_pushl()
	leal	1(%eax,%eax), %esi	C d without twos
	andl	$127, %eax		C d/2, 7 bits
ifdef(`PIC',`
	LEA(	binvert_limb_table, %ebp)
Zdisp(	movzbl,	0,(%eax,%ebp), %eax)
',`
	movzbl	binvert_limb_table(%eax), %eax	C inv 8 bits
')
	pushl	%edi		FRAME_pushl()
	leal	(%eax,%eax), %ebp	C 2*inv
	imull	%eax, %eax		C inv*inv
	movl	PARAM_DST, %edi
	imull	%esi, %eax		C inv*inv*d
	subl	%eax, %ebp		C inv = 2*inv - inv*inv*d
	leal	(%ebp,%ebp), %eax	C 2*inv
	imull	%ebp, %ebp		C inv*inv
	movl	%esi, PARAM_DIVISOR	C d without twos
	leal	(%ebx,%ecx,4), %ebx	C src end
	imull	%esi, %ebp		C inv*inv*d
	leal	(%edi,%ecx,4), %edi	C dst end
	negl	%ecx			C -size
	subl	%ebp, %eax		C inv = 2*inv - inv*inv*d
	subl	$1, %edx		C shift amount, and clear carry
	ASSERT(e,`	C expect d*inv == 1 mod 2^GMP_LIMB_BITS
	pushl	%eax	FRAME_pushl()
	imull	PARAM_DIVISOR, %eax
	cmpl	$1, %eax
	popl	%eax	FRAME_popl()')
	movl	%eax, VAR_INVERSE
	jnz	L(even)
	movl	(%ebx,%ecx,4), %esi	C src low limb
	jmp	L(odd_entry)
	ALIGN(16)
	nop	C code alignment
L(odd_top):
	C eax	scratch
	C ebx	src end
	C ecx	counter, limbs, negative
	C edx	inverse
	C esi	next limb, adjusted for carry
	C edi	dst end
	C ebp	carry bit, 0 or -1
	imull	%edx, %esi
	movl	PARAM_DIVISOR, %eax
	movl	%esi, -4(%edi,%ecx,4)
	mull	%esi			C carry limb in edx
	subl	%ebp, %edx		C apply carry bit
	movl	(%ebx,%ecx,4), %esi
L(odd_entry):
	subl	%edx, %esi		C apply carry limb
	movl	VAR_INVERSE, %edx
	sbbl	%ebp, %ebp		C 0 or -1
	incl	%ecx
	jnz	L(odd_top)
	imull	%edx, %esi
	movl	%esi, -4(%edi,%ecx,4)
	popl	%edi
	popl	%esi
	popl	%ebp
	popl	%ebx
	ret
L(even):
	C eax
	C ebx	src end
	C ecx	-size
	C edx	twos
	C esi
	C edi	dst end
	C ebp
	xorl	%ebp, %ebp
Zdisp(	movq,	0,(%ebx,%ecx,4), %mm0)	C src[0,1]
	movd	%edx, %mm7
	movl	VAR_INVERSE, %edx
	addl	$2, %ecx
	psrlq	%mm7, %mm0
	movd	%mm0, %esi
	jz	L(even_two)		C if only two limbs
C Out-of-order execution is good enough to hide the load/rshift/movd
C latency.  Having imul at the top of the loop gives 11.5 c/l instead of 12,
C on K6-2.  In fact there's only 11 of decode, but nothing running at 11 has
C been found.  Maybe the fact every second movq is unaligned costs the extra
C 0.5.
L(even_top):
	C eax	scratch
	C ebx	src end
	C ecx	counter, limbs, negative
	C edx	inverse
	C esi	next limb, adjusted for carry
	C edi	dst end
	C ebp	carry bit, 0 or -1
	C
	C mm0	scratch, source limbs
	C mm7	twos
	imull	%edx, %esi
	movl	%esi, -8(%edi,%ecx,4)
	movl	PARAM_DIVISOR, %eax
	mull	%esi			C carry limb in edx
	movq	-4(%ebx,%ecx,4), %mm0
	psrlq	%mm7, %mm0
	movd	%mm0, %esi
	subl	%ebp, %edx		C apply carry bit
	subl	%edx, %esi		C apply carry limb
	movl	VAR_INVERSE, %edx
	sbbl	%ebp, %ebp		C 0 or -1
	incl	%ecx
	jnz	L(even_top)
L(even_two):
	movd	-4(%ebx), %mm0		C src high limb
	psrlq	%mm7, %mm0
	imull	%edx, %esi
	movl	%esi, -8(%edi)
	movl	PARAM_DIVISOR, %eax
	mull	%esi			C carry limb in edx
	movd	%mm0, %esi
	subl	%ebp, %edx		C apply carry bit
	movl	VAR_INVERSE, %eax
	subl	%edx, %esi		C apply carry limb
	imull	%eax, %esi
	movl	%esi, -4(%edi)
	popl	%edi
	popl	%esi
	popl	%ebp
	popl	%ebx
	emms_or_femms
	ret
EPILOGUE()

						