--- /dev/null
+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 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 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()
projects / msp430-gcc.git / blobdiff