1dnl  AMD K7 mpn_sqr_basecase -- square an mpn number.
2
3dnl  Copyright 1999-2002 Free Software Foundation, Inc.
4
5dnl  This file is part of the GNU MP Library.
6dnl
7dnl  The GNU MP Library is free software; you can redistribute it and/or modify
8dnl  it under the terms of either:
9dnl
10dnl    * the GNU Lesser General Public License as published by the Free
11dnl      Software Foundation; either version 3 of the License, or (at your
12dnl      option) any later version.
13dnl
14dnl  or
15dnl
16dnl    * the GNU General Public License as published by the Free Software
17dnl      Foundation; either version 2 of the License, or (at your option) any
18dnl      later version.
19dnl
20dnl  or both in parallel, as here.
21dnl
22dnl  The GNU MP Library is distributed in the hope that it will be useful, but
23dnl  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
24dnl  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
25dnl  for more details.
26dnl
27dnl  You should have received copies of the GNU General Public License and the
28dnl  GNU Lesser General Public License along with the GNU MP Library.  If not,
29dnl  see https://www.gnu.org/licenses/.
30
31include(`../config.m4')
32
33
34C K7: approx 2.3 cycles/crossproduct, or 4.55 cycles/triangular product
35C     (measured on the speed difference between 25 and 50 limbs, which is
36C     roughly the Karatsuba recursing range).
37
38
39dnl  These are the same as mpn/x86/k6/sqr_basecase.asm, see that code for
40dnl  some comments.
41
42deflit(SQR_TOOM2_THRESHOLD_MAX, 66)
43
44ifdef(`SQR_TOOM2_THRESHOLD_OVERRIDE',
45`define(`SQR_TOOM2_THRESHOLD',SQR_TOOM2_THRESHOLD_OVERRIDE)')
46
47m4_config_gmp_mparam(`SQR_TOOM2_THRESHOLD')
48deflit(UNROLL_COUNT, eval(SQR_TOOM2_THRESHOLD-3))
49
50
51C void mpn_sqr_basecase (mp_ptr dst, mp_srcptr src, mp_size_t size);
52C
53C With a SQR_TOOM2_THRESHOLD around 50 this code is about 1500 bytes,
54C which is quite a bit, but is considered good value since squares big
55C enough to use most of the code will be spending quite a few cycles in it.
56
57
58defframe(PARAM_SIZE,12)
59defframe(PARAM_SRC, 8)
60defframe(PARAM_DST, 4)
61
62          TEXT
63          ALIGN(32)
64PROLOGUE(mpn_sqr_basecase)
65deflit(`FRAME',0)
66
67          movl      PARAM_SIZE, %ecx
68          movl      PARAM_SRC, %eax
69          cmpl      $2, %ecx
70
71          movl      PARAM_DST, %edx
72          je        L(two_limbs)
73          ja        L(three_or_more)
74
75
76C------------------------------------------------------------------------------
77C one limb only
78          C eax     src
79          C ecx     size
80          C edx     dst
81
82          movl      (%eax), %eax
83          movl      %edx, %ecx
84
85          mull      %eax
86
87          movl      %edx, 4(%ecx)
88          movl      %eax, (%ecx)
89          ret
90
91
92C------------------------------------------------------------------------------
93C
94C Using the read/modify/write "add"s seems to be faster than saving and
95C restoring registers.  Perhaps the loads for the first set hide under the
96C mul latency and the second gets store to load forwarding.
97
98          ALIGN(16)
99L(two_limbs):
100          C eax     src
101          C ebx
102          C ecx     size
103          C edx     dst
104deflit(`FRAME',0)
105
106          pushl     %ebx                FRAME_pushl()
107          movl      %eax, %ebx          C src
108          movl      (%eax), %eax
109
110          movl      %edx, %ecx          C dst
111
112          mull      %eax                C src[0]^2
113
114          movl      %eax, (%ecx)        C dst[0]
115          movl      4(%ebx), %eax
116
117          movl      %edx, 4(%ecx)       C dst[1]
118
119          mull      %eax                C src[1]^2
120
121          movl      %eax, 8(%ecx)       C dst[2]
122          movl      (%ebx), %eax
123
124          movl      %edx, 12(%ecx)      C dst[3]
125
126          mull      4(%ebx)             C src[0]*src[1]
127
128          popl      %ebx
129
130          addl      %eax, 4(%ecx)
131          adcl      %edx, 8(%ecx)
132          adcl      $0, 12(%ecx)
133          ASSERT(nc)
134
135          addl      %eax, 4(%ecx)
136          adcl      %edx, 8(%ecx)
137          adcl      $0, 12(%ecx)
138          ASSERT(nc)
139
140          ret
141
142
143C------------------------------------------------------------------------------
144defframe(SAVE_EBX,  -4)
145defframe(SAVE_ESI,  -8)
146defframe(SAVE_EDI, -12)
147defframe(SAVE_EBP, -16)
148deflit(STACK_SPACE, 16)
149
150L(three_or_more):
151          subl      $STACK_SPACE, %esp
152          cmpl      $4, %ecx
153          jae       L(four_or_more)
154deflit(`FRAME',STACK_SPACE)
155
156
157C------------------------------------------------------------------------------
158C Three limbs
159C
160C Writing out the loads and stores separately at the end of this code comes
161C out about 10 cycles faster than using adcls to memory.
162
163          C eax     src
164          C ecx     size
165          C edx     dst
166
167          movl      %ebx, SAVE_EBX
168          movl      %eax, %ebx          C src
169          movl      (%eax), %eax
170
171          movl      %edx, %ecx          C dst
172          movl      %esi, SAVE_ESI
173          movl      %edi, SAVE_EDI
174
175          mull      %eax                C src[0] ^ 2
176
177          movl      %eax, (%ecx)
178          movl      4(%ebx), %eax
179          movl      %edx, 4(%ecx)
180
181          mull      %eax                C src[1] ^ 2
182
183          movl      %eax, 8(%ecx)
184          movl      8(%ebx), %eax
185          movl      %edx, 12(%ecx)
186
187          mull      %eax                C src[2] ^ 2
188
189          movl      %eax, 16(%ecx)
190          movl      (%ebx), %eax
191          movl      %edx, 20(%ecx)
192
193          mull      4(%ebx)             C src[0] * src[1]
194
195          movl      %eax, %esi
196          movl      (%ebx), %eax
197          movl      %edx, %edi
198
199          mull      8(%ebx)             C src[0] * src[2]
200
201          addl      %eax, %edi
202          movl      %ebp, SAVE_EBP
203          movl      $0, %ebp
204
205          movl      4(%ebx), %eax
206          adcl      %edx, %ebp
207
208          mull      8(%ebx)             C src[1] * src[2]
209
210          xorl      %ebx, %ebx
211          addl      %eax, %ebp
212
213          adcl      $0, %edx
214
215          C eax
216          C ebx     zero, will be dst[5]
217          C ecx     dst
218          C edx     dst[4]
219          C esi     dst[1]
220          C edi     dst[2]
221          C ebp     dst[3]
222
223          adcl      $0, %edx
224          addl      %esi, %esi
225
226          adcl      %edi, %edi
227          movl      4(%ecx), %eax
228
229          adcl      %ebp, %ebp
230
231          adcl      %edx, %edx
232
233          adcl      $0, %ebx
234          addl      %eax, %esi
235          movl      8(%ecx), %eax
236
237          adcl      %eax, %edi
238          movl      12(%ecx), %eax
239          movl      %esi, 4(%ecx)
240
241          adcl      %eax, %ebp
242          movl      16(%ecx), %eax
243          movl      %edi, 8(%ecx)
244
245          movl      SAVE_ESI, %esi
246          movl      SAVE_EDI, %edi
247
248          adcl      %eax, %edx
249          movl      20(%ecx), %eax
250          movl      %ebp, 12(%ecx)
251
252          adcl      %ebx, %eax
253          ASSERT(nc)
254          movl      SAVE_EBX, %ebx
255          movl      SAVE_EBP, %ebp
256
257          movl      %edx, 16(%ecx)
258          movl      %eax, 20(%ecx)
259          addl      $FRAME, %esp
260
261          ret
262
263
264C------------------------------------------------------------------------------
265L(four_or_more):
266
267C First multiply src[0]*src[1..size-1] and store at dst[1..size].
268C Further products are added in rather than stored.
269
270          C eax     src
271          C ebx
272          C ecx     size
273          C edx     dst
274          C esi
275          C edi
276          C ebp
277
278defframe(`VAR_COUNTER',-20)
279defframe(`VAR_JMP',    -24)
280deflit(EXTRA_STACK_SPACE, 8)
281
282          movl      %ebx, SAVE_EBX
283          movl      %edi, SAVE_EDI
284          leal      (%edx,%ecx,4), %edi C &dst[size]
285
286          movl      %esi, SAVE_ESI
287          movl      %ebp, SAVE_EBP
288          leal      (%eax,%ecx,4), %esi C &src[size]
289
290          movl      (%eax), %ebp                  C multiplier
291          movl      $0, %ebx
292          decl      %ecx
293
294          negl      %ecx
295          subl      $EXTRA_STACK_SPACE, %esp
296FRAME_subl_esp(EXTRA_STACK_SPACE)
297
298L(mul_1):
299          C eax     scratch
300          C ebx     carry
301          C ecx     counter
302          C edx     scratch
303          C esi     &src[size]
304          C edi     &dst[size]
305          C ebp     multiplier
306
307          movl      (%esi,%ecx,4), %eax
308
309          mull      %ebp
310
311          addl      %ebx, %eax
312          movl      %eax, (%edi,%ecx,4)
313          movl      $0, %ebx
314
315          adcl      %edx, %ebx
316          incl      %ecx
317          jnz       L(mul_1)
318
319
320C Add products src[n]*src[n+1..size-1] at dst[2*n-1...], for each n=1..size-2.
321C
322C The last two products, which are the bottom right corner of the product
323C triangle, are left to the end.  These are src[size-3]*src[size-2,size-1]
324C and src[size-2]*src[size-1].  If size is 4 then it's only these corner
325C cases that need to be done.
326C
327C The unrolled code is the same as in mpn_addmul_1, see that routine for
328C some comments.
329C
330C VAR_COUNTER is the outer loop, running from -size+4 to -1, inclusive.
331C
332C VAR_JMP is the computed jump into the unrolled code, stepped by one code
333C chunk each outer loop.
334C
335C K7 does branch prediction on indirect jumps, which is bad since it's a
336C different target each time.  There seems no way to avoid this.
337
338dnl  This value also hard coded in some shifts and adds
339deflit(CODE_BYTES_PER_LIMB, 17)
340
341dnl  With the unmodified &src[size] and &dst[size] pointers, the
342dnl  displacements in the unrolled code fit in a byte for UNROLL_COUNT
343dnl  values up to 31, but above that an offset must be added to them.
344
345deflit(OFFSET,
346ifelse(eval(UNROLL_COUNT>31),1,
347eval((UNROLL_COUNT-31)*4),
3480))
349
350dnl  Because the last chunk of code is generated differently, a label placed
351dnl  at the end doesn't work.  Instead calculate the implied end using the
352dnl  start and how many chunks of code there are.
353
354deflit(UNROLL_INNER_END,
355`L(unroll_inner_start)+eval(UNROLL_COUNT*CODE_BYTES_PER_LIMB)')
356
357          C eax
358          C ebx     carry
359          C ecx
360          C edx
361          C esi     &src[size]
362          C edi     &dst[size]
363          C ebp
364
365          movl      PARAM_SIZE, %ecx
366          movl      %ebx, (%edi)
367
368          subl      $4, %ecx
369          jz        L(corner)
370
371          negl      %ecx
372ifelse(OFFSET,0,,`subl        $OFFSET, %edi')
373ifelse(OFFSET,0,,`subl        $OFFSET, %esi')
374
375          movl      %ecx, %edx
376          shll      $4, %ecx
377
378ifdef(`PIC',`
379          call      L(pic_calc)
380L(here):
381',`
382          leal      UNROLL_INNER_END-eval(2*CODE_BYTES_PER_LIMB)(%ecx,%edx), %ecx
383')
384
385
386          C The calculated jump mustn't come out to before the start of the
387          C code available.  This is the limit UNROLL_COUNT puts on the src
388          C operand size, but checked here directly using the jump address.
389          ASSERT(ae,
390          `movl_text_address(L(unroll_inner_start), %eax)
391          cmpl      %eax, %ecx')
392
393
394C------------------------------------------------------------------------------
395          ALIGN(16)
396L(unroll_outer_top):
397          C eax
398          C ebx     high limb to store
399          C ecx     VAR_JMP
400          C edx     VAR_COUNTER, limbs, negative
401          C esi     &src[size], constant
402          C edi     dst ptr, high of last addmul
403          C ebp
404
405          movl      -12+OFFSET(%esi,%edx,4), %ebp C next multiplier
406          movl      -8+OFFSET(%esi,%edx,4), %eax  C first of multiplicand
407
408          movl      %edx, VAR_COUNTER
409
410          mull      %ebp
411
412define(cmovX,`ifelse(eval(UNROLL_COUNT%2),0,`cmovz($@)',`cmovnz($@)')')
413
414          testb     $1, %cl
415          movl      %edx, %ebx          C high carry
416          movl      %ecx, %edx          C jump
417
418          movl      %eax, %ecx          C low carry
419          cmovX(    %ebx, %ecx)         C high carry reverse
420          cmovX(    %eax, %ebx)         C low carry reverse
421
422          leal      CODE_BYTES_PER_LIMB(%edx), %eax
423          xorl      %edx, %edx
424          leal      4(%edi), %edi
425
426          movl      %eax, VAR_JMP
427
428          jmp       *%eax
429
430
431ifdef(`PIC',`
432L(pic_calc):
433          addl      (%esp), %ecx
434          addl      $UNROLL_INNER_END-eval(2*CODE_BYTES_PER_LIMB)-L(here), %ecx
435          addl      %edx, %ecx
436          ret_internal
437')
438
439
440          C Must be an even address to preserve the significance of the low
441          C bit of the jump address indicating which way around ecx/ebx should
442          C start.
443          ALIGN(2)
444
445L(unroll_inner_start):
446          C eax     next limb
447          C ebx     carry high
448          C ecx     carry low
449          C edx     scratch
450          C esi     src
451          C edi     dst
452          C ebp     multiplier
453
454forloop(`i', UNROLL_COUNT, 1, `
455          deflit(`disp_src', eval(-i*4 + OFFSET))
456          deflit(`disp_dst', eval(disp_src - 4))
457
458          m4_assert(`disp_src>=-128 && disp_src<128')
459          m4_assert(`disp_dst>=-128 && disp_dst<128')
460
461ifelse(eval(i%2),0,`
462Zdisp(    movl,     disp_src,(%esi), %eax)
463          adcl      %edx, %ebx
464
465          mull      %ebp
466
467Zdisp(  addl,       %ecx, disp_dst,(%edi))
468          movl      $0, %ecx
469
470          adcl      %eax, %ebx
471
472',`
473          dnl  this bit comes out last
474Zdisp(  movl,       disp_src,(%esi), %eax)
475          adcl      %edx, %ecx
476
477          mull      %ebp
478
479Zdisp(    addl,     %ebx, disp_dst,(%edi))
480
481ifelse(forloop_last,0,
482`         movl      $0, %ebx')
483
484          adcl      %eax, %ecx
485')
486')
487
488          C eax     next limb
489          C ebx     carry high
490          C ecx     carry low
491          C edx     scratch
492          C esi     src
493          C edi     dst
494          C ebp     multiplier
495
496          adcl      $0, %edx
497          addl      %ecx, -4+OFFSET(%edi)
498          movl      VAR_JMP, %ecx
499
500          adcl      $0, %edx
501
502          movl      %edx, m4_empty_if_zero(OFFSET) (%edi)
503          movl      VAR_COUNTER, %edx
504
505          incl      %edx
506          jnz       L(unroll_outer_top)
507
508
509ifelse(OFFSET,0,,`
510          addl      $OFFSET, %esi
511          addl      $OFFSET, %edi
512')
513
514
515C------------------------------------------------------------------------------
516L(corner):
517          C esi     &src[size]
518          C edi     &dst[2*size-5]
519
520          movl      -12(%esi), %ebp
521          movl      -8(%esi), %eax
522          movl      %eax, %ecx
523
524          mull      %ebp
525
526          addl      %eax, -4(%edi)
527          movl      -4(%esi), %eax
528
529          adcl      $0, %edx
530          movl      %edx, %ebx
531          movl      %eax, %esi
532
533          mull      %ebp
534
535          addl      %ebx, %eax
536
537          adcl      $0, %edx
538          addl      %eax, (%edi)
539          movl      %esi, %eax
540
541          adcl      $0, %edx
542          movl      %edx, %ebx
543
544          mull      %ecx
545
546          addl      %ebx, %eax
547          movl      %eax, 4(%edi)
548
549          adcl      $0, %edx
550          movl      %edx, 8(%edi)
551
552
553
554C Left shift of dst[1..2*size-2], high bit shifted out becomes dst[2*size-1].
555
556L(lshift_start):
557          movl      PARAM_SIZE, %eax
558          movl      PARAM_DST, %edi
559          xorl      %ecx, %ecx                    C clear carry
560
561          leal      (%edi,%eax,8), %edi
562          notl      %eax                          C -size-1, preserve carry
563
564          leal      2(%eax), %eax                 C -(size-1)
565
566L(lshift):
567          C eax     counter, negative
568          C ebx
569          C ecx
570          C edx
571          C esi
572          C edi     dst, pointing just after last limb
573          C ebp
574
575          rcll      -4(%edi,%eax,8)
576          rcll      (%edi,%eax,8)
577          incl      %eax
578          jnz       L(lshift)
579
580          setc      %al
581
582          movl      PARAM_SRC, %esi
583          movl      %eax, -4(%edi)                C dst most significant limb
584
585          movl      PARAM_SIZE, %ecx
586
587
588C Now add in the squares on the diagonal, src[0]^2, src[1]^2, ...,
589C src[size-1]^2.  dst[0] hasn't yet been set at all yet, and just gets the
590C low limb of src[0]^2.
591
592          movl      (%esi), %eax                  C src[0]
593
594          mull      %eax
595
596          leal      (%esi,%ecx,4), %esi C src point just after last limb
597          negl      %ecx
598
599          movl      %eax, (%edi,%ecx,8) C dst[0]
600          incl      %ecx
601
602L(diag):
603          C eax     scratch
604          C ebx     scratch
605          C ecx     counter, negative
606          C edx     carry
607          C esi     src just after last limb
608          C edi     dst just after last limb
609          C ebp
610
611          movl      (%esi,%ecx,4), %eax
612          movl      %edx, %ebx
613
614          mull      %eax
615
616          addl      %ebx, -4(%edi,%ecx,8)
617          adcl      %eax, (%edi,%ecx,8)
618          adcl      $0, %edx
619
620          incl      %ecx
621          jnz       L(diag)
622
623
624          movl      SAVE_ESI, %esi
625          movl      SAVE_EBX, %ebx
626
627          addl      %edx, -4(%edi)                C dst most significant limb
628          movl      SAVE_EDI, %edi
629
630          movl      SAVE_EBP, %ebp
631          addl      $FRAME, %esp
632
633          ret
634
635EPILOGUE()
636