Quelle  sqr_basecase.asm   Sprache: Masm

dnl  Intel P6 mpn_sqr_basecase -- square an mpn number.

dnl  Copyright 1999, 2000, 2002 Free Software Foundation, Inc.

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 modify
dnl  it under the terms of either:
dnl
dnl    * the GNU Lesser General Public License as published by the Free
dnl      Software Foundation; either version 3 of the License, or (at your
dnl      option) any later version.
dnl
dnl  or
dnl
dnl    * the GNU General Public License as published by the Free Software
dnl      Foundation; either version 2 of the License, or (at your option) any
dnl      later version.
dnl
dnl  or both in parallel, as here.
dnl
dnl  The GNU MP Library is distributed in the hope that it will be useful, but
dnl  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
dnl  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
dnl  for more details.
dnl
dnl  You should have received copies of the GNU General Public License and the
dnl  GNU Lesser General Public License along with the GNU MP Library.  If not,
dnl  see https://www.gnu.org/licenses/.

include(`../config.m4')


C P6: approx 4.0 cycles per cross product, or 7.75 cycles per triangular
C     product (measured on the speed difference between 20 and 40 limbs,
C     which is the Karatsuba recursing range).


dnl  These are the same as in mpn/x86/k6/sqr_basecase.asm, see that file for
dnl  a description.  The only difference here is that UNROLL_COUNT can go up
dnl  to 64 (not 63) making SQR_TOOM2_THRESHOLD_MAX 67.

deflit(SQR_TOOM2_THRESHOLD_MAX, 67)

ifdef(`SQR_TOOM2_THRESHOLD_OVERRIDE',
`define(`SQR_TOOM2_THRESHOLD',SQR_TOOM2_THRESHOLD_OVERRIDE)')

m4_config_gmp_mparam(`SQR_TOOM2_THRESHOLD')
deflit(UNROLL_COUNT, eval(SQR_TOOM2_THRESHOLD-3))


C void mpn_sqr_basecase (mp_ptr dst, mp_srcptr src, mp_size_t size);
C
C The algorithm is basically the same as mpn/generic/sqr_basecase.c, but a
C lot of function call overheads are avoided, especially when the given size
C is small.
C
C The code size might look a bit excessive, but not all of it is executed so
C it won't all get into the code cache. The 1x1, 2x2 and 3x3 special cases
C clearly apply only to those sizes; mid sizes like 10x10 only need part of
C the unrolled addmul; and big sizes like 40x40 that do use the full
C unrolling will least be making good use of it, because 40x40 will take
C something like 7000 cycles.

defframe(PARAM_SIZE,12)
defframe(PARAM_SRC, 8)
defframe(PARAM_DST, 4)

 TEXT
 ALIGN(32)
PROLOGUE(mpn_sqr_basecase)
deflit(`FRAME',0)

 movl PARAM_SIZE, %edx

 movl PARAM_SRC, %eax

 cmpl $2, %edx
 movl PARAM_DST, %ecx
 je L(two_limbs)

 movl (%eax), %eax
 ja L(three_or_more)


C -----------------------------------------------------------------------------
C one limb only
 C eax src limb
 C ebx
 C ecx dst
 C edx

 mull %eax

 movl %eax, (%ecx)
 movl %edx, 4(%ecx)

 ret


C -----------------------------------------------------------------------------
L(two_limbs):
 C eax src
 C ebx
 C ecx dst
 C edx

defframe(SAVE_ESI, -4)
defframe(SAVE_EBX, -8)
defframe(SAVE_EDI, -12)
defframe(SAVE_EBP, -16)
deflit(`STACK_SPACE',16)

 subl $STACK_SPACE, %esp
deflit(`FRAME',STACK_SPACE)

 movl %esi, SAVE_ESI
 movl %eax, %esi
 movl (%eax), %eax

 mull %eax  C src[0]^2

 movl %eax, (%ecx) C dst[0]
 movl 4(%esi), %eax

 movl %ebx, SAVE_EBX
 movl %edx, %ebx C dst[1]

 mull %eax  C src[1]^2

 movl %edi, SAVE_EDI
 movl %eax, %edi C dst[2]
 movl (%esi), %eax

 movl %ebp, SAVE_EBP
 movl %edx, %ebp C dst[3]

 mull 4(%esi)  C src[0]*src[1]

 addl %eax, %ebx
 movl SAVE_ESI, %esi

 adcl %edx, %edi

 adcl $0, %ebp
 addl %ebx, %eax
 movl SAVE_EBX, %ebx

 adcl %edi, %edx
 movl SAVE_EDI, %edi

 adcl $0, %ebp

 movl %eax, 4(%ecx)

 movl %ebp, 12(%ecx)
 movl SAVE_EBP, %ebp

 movl %edx, 8(%ecx)
 addl $FRAME, %esp

 ret


C -----------------------------------------------------------------------------
L(three_or_more):
 C eax src low limb
 C ebx
 C ecx dst
 C edx size
deflit(`FRAME',0)

 pushl %esi defframe_pushl(`SAVE_ESI')
 cmpl $4, %edx

 movl PARAM_SRC, %esi
 jae L(four_or_more)


C -----------------------------------------------------------------------------
C three limbs

 C eax src low limb
 C ebx
 C ecx dst
 C edx
 C esi src
 C edi
 C ebp

 pushl %ebp defframe_pushl(`SAVE_EBP')
 pushl %edi defframe_pushl(`SAVE_EDI')

 mull %eax  C src[0] ^ 2

 movl %eax, (%ecx)
 movl %edx, 4(%ecx)

 movl 4(%esi), %eax
 xorl %ebp, %ebp

 mull %eax  C src[1] ^ 2

 movl %eax, 8(%ecx)
 movl %edx, 12(%ecx)
 movl 8(%esi), %eax

 pushl %ebx defframe_pushl(`SAVE_EBX')

 mull %eax  C src[2] ^ 2

 movl %eax, 16(%ecx)
 movl %edx, 20(%ecx)

 movl (%esi), %eax

 mull 4(%esi)  C src[0] * src[1]

 movl %eax, %ebx
 movl %edx, %edi

 movl (%esi), %eax

 mull 8(%esi)  C src[0] * src[2]

 addl %eax, %edi
 movl %edx, %ebp

 adcl $0, %ebp
 movl 4(%esi), %eax

 mull 8(%esi)  C src[1] * src[2]

 xorl %esi, %esi
 addl %eax, %ebp

 C eax
 C ebx dst[1]
 C ecx dst
 C edx dst[4]
 C esi zero, will be dst[5]
 C edi dst[2]
 C ebp dst[3]

 adcl $0, %edx
 addl %ebx, %ebx

 adcl %edi, %edi

 adcl %ebp, %ebp

 adcl %edx, %edx
 movl 4(%ecx), %eax

 adcl $0, %esi
 addl %ebx, %eax

 movl %eax, 4(%ecx)
 movl 8(%ecx), %eax

 adcl %edi, %eax
 movl 12(%ecx), %ebx

 adcl %ebp, %ebx
 movl 16(%ecx), %edi

 movl %eax, 8(%ecx)
 movl SAVE_EBP, %ebp

 movl %ebx, 12(%ecx)
 movl SAVE_EBX, %ebx

 adcl %edx, %edi
 movl 20(%ecx), %eax

 movl %edi, 16(%ecx)
 movl SAVE_EDI, %edi

 adcl %esi, %eax C no carry out of this
 movl SAVE_ESI, %esi

 movl %eax, 20(%ecx)
 addl $FRAME, %esp

 ret



C -----------------------------------------------------------------------------
defframe(VAR_COUNTER,-20)
defframe(VAR_JMP,    -24)
deflit(`STACK_SPACE',24)

L(four_or_more):
 C eax src low limb
 C ebx
 C ecx
 C edx size
 C esi src
 C edi
 C ebp
deflit(`FRAME',4) dnl %esi already pushed

C First multiply src[0]*src[1..size-1] and store at dst[1..size].

 subl $STACK_SPACE-FRAME, %esp
deflit(`FRAME',STACK_SPACE)
 movl $1, %ecx

 movl %edi, SAVE_EDI
 movl PARAM_DST, %edi

 movl %ebx, SAVE_EBX
 subl %edx, %ecx  C -(size-1)

 movl %ebp, SAVE_EBP
 movl $0, %ebx  C initial carry

 leal (%esi,%edx,4), %esi C &src[size]
 movl %eax, %ebp  C multiplier

 leal -4(%edi,%edx,4), %edi C &dst[size-1]


C This loop runs at just over 6 c/l.

L(mul_1):
 C eax scratch
 C ebx carry
 C ecx counter, limbs, negative, -(size-1) to -1
 C edx scratch
 C esi &src[size]
 C edi &dst[size-1]
 C ebp multiplier

 movl %ebp, %eax

 mull (%esi,%ecx,4)

 addl %ebx, %eax
 movl $0, %ebx

 adcl %edx, %ebx
 movl %eax, 4(%edi,%ecx,4)

 incl %ecx
 jnz L(mul_1)


 movl %ebx, 4(%edi)


C Addmul src[n]*src[n+1..size-1] at dst[2*n-1...], for each n=1..size-2.
C
C The last two addmuls, which are the bottom right corner of the product
C triangle, are left to the end.  These are src[size-3]*src[size-2,size-1]
C and src[size-2]*src[size-1].  If size is 4 then it's only these corner
C cases that need to be done.
C
C The unrolled code is the same as mpn_addmul_1(), see that routine for some
C comments.
C
C VAR_COUNTER is the outer loop, running from -(size-4) to -1, inclusive.
C
C VAR_JMP is the computed jump into the unrolled code, stepped by one code
C chunk each outer loop.

dnl  This is also hard-coded in the address calculation below.
deflit(CODE_BYTES_PER_LIMB, 15)

dnl  With &src[size] and &dst[size-1] pointers, the displacements in the
dnl  unrolled code fit in a byte for UNROLL_COUNT values up to 32, but above
dnl  that an offset must be added to them.
deflit(OFFSET,
ifelse(eval(UNROLL_COUNT>32),1,
eval((UNROLL_COUNT-32)*4),
0))

 C eax
 C ebx carry
 C ecx
 C edx
 C esi &src[size]
 C edi &dst[size-1]
 C ebp

 movl PARAM_SIZE, %ecx

 subl $4, %ecx
 jz L(corner)

 movl %ecx, %edx
 negl %ecx

 shll $4, %ecx
ifelse(OFFSET,0,,`subl $OFFSET, %esi')

ifdef(`PIC',`
 call L(pic_calc)
L(here):
',`
 leal L(unroll_inner_end)-eval(2*CODE_BYTES_PER_LIMB)(%ecx,%edx), %ecx
')
 negl %edx

ifelse(OFFSET,0,,`subl $OFFSET, %edi')

 C The calculated jump mustn't be before the start of the available
 C code.  This is the limit that UNROLL_COUNT puts on the src operand
 C size, but checked here using the jump address directly.

 ASSERT(ae,
 `movl_text_address( L(unroll_inner_start), %eax)
 cmpl %eax, %ecx')


C -----------------------------------------------------------------------------
 ALIGN(16)
L(unroll_outer_top):
 C eax
 C ebx high limb to store
 C ecx VAR_JMP
 C edx VAR_COUNTER, limbs, negative
 C esi &src[size], constant
 C edi dst ptr, second highest limb of last addmul
 C ebp

 movl -12+OFFSET(%esi,%edx,4), %ebp C multiplier
 movl %edx, VAR_COUNTER

 movl -8+OFFSET(%esi,%edx,4), %eax C first limb of multiplicand

 mull %ebp

define(cmovX,`ifelse(eval(UNROLL_COUNT%2),1,`cmovz($@)',`cmovnz($@)')')

 testb $1, %cl

 movl %edx, %ebx C high carry
 leal 4(%edi), %edi

 movl %ecx, %edx C jump

 movl %eax, %ecx C low carry
 leal CODE_BYTES_PER_LIMB(%edx), %edx

 cmovX( %ebx, %ecx) C high carry reverse
 cmovX( %eax, %ebx) C low carry reverse
 movl %edx, VAR_JMP
 jmp *%edx


 C Must be on an even address here so the low bit of the jump address
 C will indicate which way around ecx/ebx should start.

 ALIGN(2)

L(unroll_inner_start):
 C eax scratch
 C ebx carry high
 C ecx carry low
 C edx scratch
 C esi src pointer
 C edi dst pointer
 C ebp multiplier
 C
 C 15 code bytes each limb
 C ecx/ebx reversed on each chunk

forloop(`i', UNROLL_COUNT, 1, `
 deflit(`disp_src', eval(-i*4 + OFFSET))
 deflit(`disp_dst', eval(disp_src))

 m4_assert(`disp_src>=-128 && disp_src<128')
 m4_assert(`disp_dst>=-128 && disp_dst<128')

ifelse(eval(i%2),0,`
Zdisp( movl, disp_src,(%esi), %eax)
 mull %ebp
Zdisp( addl, %ebx, disp_dst,(%edi))
 adcl %eax, %ecx
 movl %edx, %ebx
 adcl $0, %ebx
',`
 dnl  this one comes out last
Zdisp( movl, disp_src,(%esi), %eax)
 mull %ebp
Zdisp( addl, %ecx, disp_dst,(%edi))
 adcl %eax, %ebx
 movl %edx, %ecx
 adcl $0, %ecx
')
')
L(unroll_inner_end):

 addl %ebx, m4_empty_if_zero(OFFSET)(%edi)

 movl VAR_COUNTER, %edx
 adcl $0, %ecx

 movl %ecx, m4_empty_if_zero(OFFSET+4)(%edi)
 movl VAR_JMP, %ecx

 incl %edx
 jnz L(unroll_outer_top)


ifelse(OFFSET,0,,`
 addl $OFFSET, %esi
 addl $OFFSET, %edi
')


C -----------------------------------------------------------------------------
 ALIGN(16)
L(corner):
 C eax
 C ebx
 C ecx
 C edx
 C esi &src[size]
 C edi &dst[2*size-5]
 C ebp

 movl -12(%esi), %eax

 mull -8(%esi)

 addl %eax, (%edi)
 movl -12(%esi), %eax
 movl $0, %ebx

 adcl %edx, %ebx

 mull -4(%esi)

 addl %eax, %ebx
 movl -8(%esi), %eax

 adcl $0, %edx

 addl %ebx, 4(%edi)
 movl $0, %ebx

 adcl %edx, %ebx

 mull -4(%esi)

 movl PARAM_SIZE, %ecx
 addl %ebx, %eax

 adcl $0, %edx

 movl %eax, 8(%edi)

 movl %edx, 12(%edi)
 movl PARAM_DST, %edi


C Left shift of dst[1..2*size-2], the bit shifted out becomes dst[2*size-1].

 subl $1, %ecx  C size-1
 xorl %eax, %eax  C ready for final adcl, and clear carry

 movl %ecx, %edx
 movl PARAM_SRC, %esi


L(lshift):
 C eax
 C ebx
 C ecx counter, size-1 to 1
 C edx size-1 (for later use)
 C esi src (for later use)
 C edi dst, incrementing
 C ebp

 rcll 4(%edi)
 rcll 8(%edi)

 leal 8(%edi), %edi
 decl %ecx
 jnz L(lshift)


 adcl %eax, %eax

 movl %eax, 4(%edi)  C dst most significant limb
 movl (%esi), %eax  C src[0]

 leal 4(%esi,%edx,4), %esi C &src[size]
 subl %edx, %ecx  C -(size-1)


C Now add in the squares on the diagonal, src[0]^2, src[1]^2, ...,
C src[size-1]^2.  dst[0] hasn't yet been set at all yet, and just gets the
C low limb of src[0]^2.


 mull %eax

 movl %eax, (%edi,%ecx,8) C dst[0]


L(diag):
 C eax scratch
 C ebx scratch
 C ecx counter, negative
 C edx carry
 C esi &src[size]
 C edi dst[2*size-2]
 C ebp

 movl (%esi,%ecx,4), %eax
 movl %edx, %ebx

 mull %eax

 addl %ebx, 4(%edi,%ecx,8)
 adcl %eax, 8(%edi,%ecx,8)
 adcl $0, %edx

 incl %ecx
 jnz L(diag)


 movl SAVE_ESI, %esi
 movl SAVE_EBX, %ebx

 addl %edx, 4(%edi)  C dst most significant limb

 movl SAVE_EDI, %edi
 movl SAVE_EBP, %ebp
 addl $FRAME, %esp
 ret



C -----------------------------------------------------------------------------
ifdef(`PIC',`
L(pic_calc):
 addl (%esp), %ecx
 addl $L(unroll_inner_end)-L(here)-eval(2*CODE_BYTES_PER_LIMB), %ecx
 addl %edx, %ecx
 ret_internal
')


EPILOGUE()