1 dnl AMD K7 mpn_sqr_basecase -- square an mpn number.
3 dnl K7: approx 2.3 cycles/crossproduct, or 4.55 cycles/triangular product
4 dnl (measured on the speed difference between 25 and 50 limbs, which is
5 dnl roughly the Karatsuba recursing range).
8 dnl Copyright (C) 1999, 2000 Free Software Foundation, Inc.
10 dnl This file is part of the GNU MP Library.
12 dnl The GNU MP Library is free software; you can redistribute it and/or
13 dnl modify it under the terms of the GNU Lesser General Public License as
14 dnl published by the Free Software Foundation; either version 2.1 of the
15 dnl License, or (at your option) any later version.
17 dnl The GNU MP Library is distributed in the hope that it will be useful,
18 dnl but WITHOUT ANY WARRANTY; without even the implied warranty of
19 dnl MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 dnl Lesser General Public License for more details.
22 dnl You should have received a copy of the GNU Lesser General Public
23 dnl License along with the GNU MP Library; see the file COPYING.LIB. If
24 dnl not, write to the Free Software Foundation, Inc., 59 Temple Place -
25 dnl Suite 330, Boston, MA 02111-1307, USA.
28 include(`../config.m4')
31 dnl These are the same as mpn/x86/k6/sqr_basecase.asm, see that code for
34 deflit(KARATSUBA_SQR_THRESHOLD_MAX, 66)
36 ifdef(`KARATSUBA_SQR_THRESHOLD_OVERRIDE',
37 `define(`KARATSUBA_SQR_THRESHOLD',KARATSUBA_SQR_THRESHOLD_OVERRIDE)')
39 m4_config_gmp_mparam(`KARATSUBA_SQR_THRESHOLD')
40 deflit(UNROLL_COUNT, eval(KARATSUBA_SQR_THRESHOLD-3))
43 C void mpn_sqr_basecase (mp_ptr dst, mp_srcptr src, mp_size_t size);
45 C With a KARATSUBA_SQR_THRESHOLD around 50 this code is about 1500 bytes,
46 C which is quite a bit, but is considered good value since squares big
47 C enough to use most of the code will be spending quite a few cycles in it.
50 defframe(PARAM_SIZE,12)
51 defframe(PARAM_SRC, 8)
52 defframe(PARAM_DST, 4)
56 PROLOGUE(mpn_sqr_basecase)
68 C------------------------------------------------------------------------------
84 C------------------------------------------------------------------------------
86 C Using the read/modify/write "add"s seems to be faster than saving and
87 C restoring registers. Perhaps the loads for the first set hide under the
88 C mul latency and the second gets store to load forwarding.
98 pushl %ebx FRAME_pushl()
102 movl %edx, %ecx C dst
106 movl %eax, (%ecx) C dst[0]
109 movl %edx, 4(%ecx) C dst[1]
113 movl %eax, 8(%ecx) C dst[2]
116 movl %edx, 12(%ecx) C dst[3]
118 mull 4(%ebx) C src[0]*src[1]
135 C------------------------------------------------------------------------------
136 defframe(SAVE_EBX, -4)
137 defframe(SAVE_ESI, -8)
138 defframe(SAVE_EDI, -12)
139 defframe(SAVE_EBP, -16)
140 deflit(STACK_SPACE, 16)
143 subl $STACK_SPACE, %esp
146 deflit(`FRAME',STACK_SPACE)
149 C------------------------------------------------------------------------------
152 C Writing out the loads and stores separately at the end of this code comes
153 C out about 10 cycles faster than using adcls to memory.
160 movl %eax, %ebx C src
163 movl %edx, %ecx C dst
167 mull %eax C src[0] ^ 2
173 mull %eax C src[1] ^ 2
179 mull %eax C src[2] ^ 2
185 mull 4(%ebx) C src[0] * src[1]
191 mull 8(%ebx) C src[0] * src[2]
200 mull 8(%ebx) C src[1] * src[2]
208 C ebx zero, will be dst[5]
256 C------------------------------------------------------------------------------
259 C First multiply src[0]*src[1..size-1] and store at dst[1..size].
260 C Further products are added in rather than stored.
270 defframe(`VAR_COUNTER',-20)
271 defframe(`VAR_JMP', -24)
272 deflit(EXTRA_STACK_SPACE, 8)
276 leal (%edx,%ecx,4), %edi C &dst[size]
280 leal (%eax,%ecx,4), %esi C &src[size]
282 movl (%eax), %ebp C multiplier
287 subl $EXTRA_STACK_SPACE, %esp
288 FRAME_subl_esp(EXTRA_STACK_SPACE)
299 movl (%esi,%ecx,4), %eax
304 movl %eax, (%edi,%ecx,4)
312 C Add products src[n]*src[n+1..size-1] at dst[2*n-1...], for each n=1..size-2.
314 C The last two products, which are the bottom right corner of the product
315 C triangle, are left to the end. These are src[size-3]*src[size-2,size-1]
316 C and src[size-2]*src[size-1]. If size is 4 then it's only these corner
317 C cases that need to be done.
319 C The unrolled code is the same as in mpn_addmul_1, see that routine for
322 C VAR_COUNTER is the outer loop, running from -size+4 to -1, inclusive.
324 C VAR_JMP is the computed jump into the unrolled code, stepped by one code
325 C chunk each outer loop.
327 C K7 does branch prediction on indirect jumps, which is bad since it's a
328 C different target each time. There seems no way to avoid this.
330 dnl This value also hard coded in some shifts and adds
331 deflit(CODE_BYTES_PER_LIMB, 17)
333 dnl With the unmodified &src[size] and &dst[size] pointers, the
334 dnl displacements in the unrolled code fit in a byte for UNROLL_COUNT
335 dnl values up to 31, but above that an offset must be added to them.
338 ifelse(eval(UNROLL_COUNT>31),1,
339 eval((UNROLL_COUNT-31)*4),
342 dnl Because the last chunk of code is generated differently, a label placed
343 dnl at the end doesn't work. Instead calculate the implied end using the
344 dnl start and how many chunks of code there are.
346 deflit(UNROLL_INNER_END,
347 `L(unroll_inner_start)+eval(UNROLL_COUNT*CODE_BYTES_PER_LIMB)')
357 movl PARAM_SIZE, %ecx
364 ifelse(OFFSET,0,,`subl $OFFSET, %edi')
365 ifelse(OFFSET,0,,`subl $OFFSET, %esi')
374 leal UNROLL_INNER_END-eval(2*CODE_BYTES_PER_LIMB)(%ecx,%edx), %ecx
378 C The calculated jump mustn't come out to before the start of the
379 C code available. This is the limit UNROLL_COUNT puts on the src
380 C operand size, but checked here directly using the jump address.
382 `movl_text_address(L(unroll_inner_start), %eax)
386 C------------------------------------------------------------------------------
390 C ebx high limb to store
392 C edx VAR_COUNTER, limbs, negative
393 C esi &src[size], constant
394 C edi dst ptr, high of last addmul
397 movl -12+OFFSET(%esi,%edx,4), %ebp C next multiplier
398 movl -8+OFFSET(%esi,%edx,4), %eax C first of multiplicand
400 movl %edx, VAR_COUNTER
404 define(cmovX,`ifelse(eval(UNROLL_COUNT%2),0,`cmovz($@)',`cmovnz($@)')')
407 movl %edx, %ebx C high carry
408 movl %ecx, %edx C jump
410 movl %eax, %ecx C low carry
411 cmovX( %ebx, %ecx) C high carry reverse
412 cmovX( %eax, %ebx) C low carry reverse
414 leal CODE_BYTES_PER_LIMB(%edx), %eax
426 addl $UNROLL_INNER_END-eval(2*CODE_BYTES_PER_LIMB)-L(here), %ecx
432 C Must be an even address to preserve the significance of the low
433 C bit of the jump address indicating which way around ecx/ebx should
437 L(unroll_inner_start):
446 forloop(`i', UNROLL_COUNT, 1, `
447 deflit(`disp_src', eval(-i*4 + OFFSET))
448 deflit(`disp_dst', eval(disp_src - 4))
450 m4_assert(`disp_src>=-128 && disp_src<128')
451 m4_assert(`disp_dst>=-128 && disp_dst<128')
454 Zdisp( movl, disp_src,(%esi), %eax)
459 Zdisp( addl, %ecx, disp_dst,(%edi))
465 dnl this bit comes out last
466 Zdisp( movl, disp_src,(%esi), %eax)
471 dnl Zdisp( addl %ebx, disp_src,(%edi))
472 addl %ebx, disp_dst(%edi)
473 ifelse(forloop_last,0,
489 addl %ecx, -4+OFFSET(%edi)
494 movl %edx, m4_empty_if_zero(OFFSET) (%edi)
495 movl VAR_COUNTER, %edx
498 jnz L(unroll_outer_top)
507 C------------------------------------------------------------------------------
546 C Left shift of dst[1..2*size-2], high bit shifted out becomes dst[2*size-1].
549 movl PARAM_SIZE, %eax
551 xorl %ecx, %ecx C clear carry
553 leal (%edi,%eax,8), %edi
554 notl %eax C -size-1, preserve carry
556 leal 2(%eax), %eax C -(size-1)
559 C eax counter, negative
564 C edi dst, pointing just after last limb
575 movl %eax, -4(%edi) C dst most significant limb
577 movl PARAM_SIZE, %ecx
580 C Now add in the squares on the diagonal, src[0]^2, src[1]^2, ...,
581 C src[size-1]^2. dst[0] hasn't yet been set at all yet, and just gets the
582 C low limb of src[0]^2.
584 movl (%esi), %eax C src[0]
588 leal (%esi,%ecx,4), %esi C src point just after last limb
591 movl %eax, (%edi,%ecx,8) C dst[0]
597 C ecx counter, negative
599 C esi src just after last limb
600 C edi dst just after last limb
603 movl (%esi,%ecx,4), %eax
608 addl %ebx, -4(%edi,%ecx,8)
609 adcl %eax, (%edi,%ecx,8)
619 addl %edx, -4(%edi) C dst most significant limb