![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to divrem_1.asm CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / x86 / p6 / mmx|
Return to divrem_1.asm CVS log | Up to [local] / OpenXM_contrib / gmp / mpn / x86 / p6 / mmx |
1.1 maekawa 1: dnl Intel Pentium-II mpn_divrem_1 -- mpn by limb division.
2: dnl
3: dnl P6MMX: 25.0 cycles/limb integer part, 17.5 cycles/limb fraction part.
4:
5:
6: dnl Copyright (C) 1999, 2000 Free Software Foundation, Inc.
7: dnl
8: dnl This file is part of the GNU MP Library.
9: dnl
10: dnl The GNU MP Library is free software; you can redistribute it and/or
11: dnl modify it under the terms of the GNU Lesser General Public License as
12: dnl published by the Free Software Foundation; either version 2.1 of the
13: dnl License, or (at your option) any later version.
14: dnl
15: dnl The GNU MP Library is distributed in the hope that it will be useful,
16: dnl but WITHOUT ANY WARRANTY; without even the implied warranty of
17: dnl MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18: dnl Lesser General Public License for more details.
19: dnl
20: dnl You should have received a copy of the GNU Lesser General Public
21: dnl License along with the GNU MP Library; see the file COPYING.LIB. If
22: dnl not, write to the Free Software Foundation, Inc., 59 Temple Place -
23: dnl Suite 330, Boston, MA 02111-1307, USA.
24:
25:
26: include(`../config.m4')
27:
28:
29: C mp_limb_t mpn_divrem_1 (mp_ptr dst, mp_size_t xsize,
30: C mp_srcptr src, mp_size_t size,
31: C mp_limb_t divisor);
32: C mp_limb_t mpn_divrem_1c (mp_ptr dst, mp_size_t xsize,
33: C mp_srcptr src, mp_size_t size,
34: C mp_limb_t divisor, mp_limb_t carry);
35: C
36: C This code is a lightly reworked version of mpn/x86/k7/mmx/divrem_1.asm,
37: C see that file for some comments. It's likely what's here can be improved.
38:
39:
40: dnl MUL_THRESHOLD is the value of xsize+size at which the multiply by
41: dnl inverse method is used, rather than plain "divl"s. Minimum value 1.
42: dnl
43: dnl The different speeds of the integer and fraction parts means that using
44: dnl xsize+size isn't quite right. The threshold wants to be a bit higher
45: dnl for the integer part and a bit lower for the fraction part. (Or what's
46: dnl really wanted is to speed up the integer part!)
47: dnl
48: dnl The threshold is set to make the integer part right. At 4 limbs the
49: dnl div and mul are about the same there, but on the fractional part the
50: dnl mul is much faster.
51:
52: deflit(MUL_THRESHOLD, 4)
53:
54:
55: defframe(PARAM_CARRY, 24)
56: defframe(PARAM_DIVISOR,20)
57: defframe(PARAM_SIZE, 16)
58: defframe(PARAM_SRC, 12)
59: defframe(PARAM_XSIZE, 8)
60: defframe(PARAM_DST, 4)
61:
62: defframe(SAVE_EBX, -4)
63: defframe(SAVE_ESI, -8)
64: defframe(SAVE_EDI, -12)
65: defframe(SAVE_EBP, -16)
66:
67: defframe(VAR_NORM, -20)
68: defframe(VAR_INVERSE, -24)
69: defframe(VAR_SRC, -28)
70: defframe(VAR_DST, -32)
71: defframe(VAR_DST_STOP,-36)
72:
73: deflit(STACK_SPACE, 36)
74:
75: .text
76: ALIGN(16)
77:
78: PROLOGUE(mpn_divrem_1c)
79: deflit(`FRAME',0)
80: movl PARAM_CARRY, %edx
81:
82: movl PARAM_SIZE, %ecx
83: subl $STACK_SPACE, %esp
84: deflit(`FRAME',STACK_SPACE)
85:
86: movl %ebx, SAVE_EBX
87: movl PARAM_XSIZE, %ebx
88:
89: movl %edi, SAVE_EDI
90: movl PARAM_DST, %edi
91:
92: movl %ebp, SAVE_EBP
93: movl PARAM_DIVISOR, %ebp
94:
95: movl %esi, SAVE_ESI
96: movl PARAM_SRC, %esi
97:
98: leal -4(%edi,%ebx,4), %edi
99: jmp LF(mpn_divrem_1,start_1c)
100:
101: EPILOGUE()
102:
103:
104: C offset 0x31, close enough to aligned
105: PROLOGUE(mpn_divrem_1)
106: deflit(`FRAME',0)
107:
108: movl PARAM_SIZE, %ecx
109: movl $0, %edx C initial carry (if can't skip a div)
110: subl $STACK_SPACE, %esp
111: deflit(`FRAME',STACK_SPACE)
112:
113: movl %ebp, SAVE_EBP
114: movl PARAM_DIVISOR, %ebp
115:
116: movl %ebx, SAVE_EBX
117: movl PARAM_XSIZE, %ebx
118:
119: movl %esi, SAVE_ESI
120: movl PARAM_SRC, %esi
121: orl %ecx, %ecx
122:
123: movl %edi, SAVE_EDI
124: movl PARAM_DST, %edi
125:
126: leal -4(%edi,%ebx,4), %edi C &dst[xsize-1]
127: jz L(no_skip_div)
128:
129: movl -4(%esi,%ecx,4), %eax C src high limb
130: cmpl %ebp, %eax C one less div if high<divisor
131: jnb L(no_skip_div)
132:
133: movl $0, (%edi,%ecx,4) C dst high limb
134: decl %ecx C size-1
135: movl %eax, %edx C src high limb as initial carry
136: L(no_skip_div):
137:
138:
139: L(start_1c):
140: C eax
141: C ebx xsize
142: C ecx size
143: C edx carry
144: C esi src
145: C edi &dst[xsize-1]
146: C ebp divisor
147:
148: leal (%ebx,%ecx), %eax C size+xsize
149: cmpl $MUL_THRESHOLD, %eax
150: jae L(mul_by_inverse)
151:
152: orl %ecx, %ecx
153: jz L(divide_no_integer)
154:
155: L(divide_integer):
156: C eax scratch (quotient)
157: C ebx xsize
158: C ecx counter
159: C edx scratch (remainder)
160: C esi src
161: C edi &dst[xsize-1]
162: C ebp divisor
163:
164: movl -4(%esi,%ecx,4), %eax
165:
166: divl %ebp
167:
168: movl %eax, (%edi,%ecx,4)
169: decl %ecx
170: jnz L(divide_integer)
171:
172:
173: L(divide_no_integer):
174: movl PARAM_DST, %edi
175: orl %ebx, %ebx
176: jnz L(divide_fraction)
177:
178: L(divide_done):
179: movl SAVE_ESI, %esi
180:
181: movl SAVE_EDI, %edi
182:
183: movl SAVE_EBX, %ebx
184: movl %edx, %eax
185:
186: movl SAVE_EBP, %ebp
187: addl $STACK_SPACE, %esp
188:
189: ret
190:
191:
192: L(divide_fraction):
193: C eax scratch (quotient)
194: C ebx counter
195: C ecx
196: C edx scratch (remainder)
197: C esi
198: C edi dst
199: C ebp divisor
200:
201: movl $0, %eax
202:
203: divl %ebp
204:
205: movl %eax, -4(%edi,%ebx,4)
206: decl %ebx
207: jnz L(divide_fraction)
208:
209: jmp L(divide_done)
210:
211:
212:
213: C -----------------------------------------------------------------------------
214:
215: L(mul_by_inverse):
216: C eax
217: C ebx xsize
218: C ecx size
219: C edx carry
220: C esi src
221: C edi &dst[xsize-1]
222: C ebp divisor
223:
224: leal 12(%edi), %ebx
225:
226: movl %ebx, VAR_DST_STOP
227: leal 4(%edi,%ecx,4), %edi C &dst[xsize+size]
228:
229: movl %edi, VAR_DST
230: movl %ecx, %ebx C size
231:
232: bsrl %ebp, %ecx C 31-l
233: movl %edx, %edi C carry
234:
235: leal 1(%ecx), %eax C 32-l
236: xorl $31, %ecx C l
237:
238: movl %ecx, VAR_NORM
239: movl $-1, %edx
240:
241: shll %cl, %ebp C d normalized
242: movd %eax, %mm7
243:
244: movl $-1, %eax
245: subl %ebp, %edx C (b-d)-1 giving edx:eax = b*(b-d)-1
246:
247: divl %ebp C floor (b*(b-d)-1) / d
248:
249: movl %eax, VAR_INVERSE
250: orl %ebx, %ebx C size
251: leal -12(%esi,%ebx,4), %eax C &src[size-3]
252:
253: movl %eax, VAR_SRC
254: jz L(start_zero)
255:
256: movl 8(%eax), %esi C src high limb
257: cmpl $1, %ebx
258: jz L(start_one)
259:
260: L(start_two_or_more):
261: movl 4(%eax), %edx C src second highest limb
262:
263: shldl( %cl, %esi, %edi) C n2 = carry,high << l
264:
265: shldl( %cl, %edx, %esi) C n10 = high,second << l
266:
267: cmpl $2, %ebx
268: je L(integer_two_left)
269: jmp L(integer_top)
270:
271:
272: L(start_one):
273: shldl( %cl, %esi, %edi) C n2 = carry,high << l
274:
275: shll %cl, %esi C n10 = high << l
276: jmp L(integer_one_left)
277:
278:
279: L(start_zero):
280: shll %cl, %edi C n2 = carry << l
281: movl $0, %esi C n10 = 0
282:
283: C we're here because xsize+size>=MUL_THRESHOLD, so with size==0 then
284: C must have xsize!=0
285: jmp L(fraction_some)
286:
287:
288:
289: C -----------------------------------------------------------------------------
290: C
291: C This loop runs at about 25 cycles, which is probably sub-optimal, and
292: C certainly more than the dependent chain would suggest. A better loop, or
293: C a better rough analysis of what's possible, would be welcomed.
294: C
295: C In the current implementation, the following successively dependent
296: C micro-ops seem to exist.
297: C
298: C uops
299: C n2+n1 1 (addl)
300: C mul 5
301: C q1+1 3 (addl/adcl)
302: C mul 5
303: C sub 3 (subl/sbbl)
304: C addback 2 (cmov)
305: C ---
306: C 19
307: C
308: C Lack of registers hinders explicit scheduling and it might be that the
309: C normal out of order execution isn't able to hide enough under the mul
310: C latencies.
311: C
312: C Using sarl/negl to pick out n1 for the n2+n1 stage is a touch faster than
313: C cmov (and takes one uop off the dependent chain). A sarl/andl/addl
314: C combination was tried for the addback (despite the fact it would lengthen
315: C the dependent chain) but found to be no faster.
316:
317:
318: ALIGN(16)
319: L(integer_top):
320: C eax scratch
321: C ebx scratch (nadj, q1)
322: C ecx scratch (src, dst)
323: C edx scratch
324: C esi n10
325: C edi n2
326: C ebp d
327: C
328: C mm0 scratch (src qword)
329: C mm7 rshift for normalization
330:
331: movl %esi, %eax
332: movl %ebp, %ebx
333:
334: sarl $31, %eax C -n1
335: movl VAR_SRC, %ecx
336:
337: andl %eax, %ebx C -n1 & d
338: negl %eax C n1
339:
340: addl %esi, %ebx C nadj = n10 + (-n1 & d), ignoring overflow
341: addl %edi, %eax C n2+n1
342: movq (%ecx), %mm0 C next src limb and the one below it
343:
344: mull VAR_INVERSE C m*(n2+n1)
345:
346: subl $4, %ecx
347:
348: movl %ecx, VAR_SRC
349:
350: C
351:
352: C
353:
354: addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
355: movl %ebp, %eax C d
356: leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
357:
358: adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
359: jz L(q1_ff)
360:
361: mull %ebx C (q1+1)*d
362:
363: movl VAR_DST, %ecx
364: psrlq %mm7, %mm0
365:
366: C
367:
368: C
369:
370: C
371:
372: subl %eax, %esi
373: movl VAR_DST_STOP, %eax
374:
375: sbbl %edx, %edi C n - (q1+1)*d
376: movl %esi, %edi C remainder -> n2
377: leal (%ebp,%esi), %edx
378:
379: cmovc( %edx, %edi) C n - q1*d if underflow from using q1+1
380: movd %mm0, %esi
381:
382: sbbl $0, %ebx C q
383: subl $4, %ecx
384:
385: movl %ebx, (%ecx)
386: cmpl %eax, %ecx
387:
388: movl %ecx, VAR_DST
389: jne L(integer_top)
390:
391:
392: L(integer_loop_done):
393:
394:
395: C -----------------------------------------------------------------------------
396: C
397: C Here, and in integer_one_left below, an sbbl $0 is used rather than a jz
398: C q1_ff special case. This make the code a bit smaller and simpler, and
399: C costs only 2 cycles (each).
400:
401: L(integer_two_left):
402: C eax scratch
403: C ebx scratch (nadj, q1)
404: C ecx scratch (src, dst)
405: C edx scratch
406: C esi n10
407: C edi n2
408: C ebp divisor
409: C
410: C mm0 src limb, shifted
411: C mm7 rshift
412:
413:
414: movl %esi, %eax
415: movl %ebp, %ebx
416:
417: sarl $31, %eax C -n1
418: movl PARAM_SRC, %ecx
419:
420: andl %eax, %ebx C -n1 & d
421: negl %eax C n1
422:
423: addl %esi, %ebx C nadj = n10 + (-n1 & d), ignoring overflow
424: addl %edi, %eax C n2+n1
425:
426: mull VAR_INVERSE C m*(n2+n1)
427:
428: movd (%ecx), %mm0 C src low limb
429:
430: movl VAR_DST_STOP, %ecx
431:
432: C
433:
434: C
435:
436: addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
437: leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
438: movl %ebp, %eax C d
439:
440: adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
441:
442: sbbl $0, %ebx
443:
444: mull %ebx C (q1+1)*d
445:
446: psllq $32, %mm0
447:
448: psrlq %mm7, %mm0
449:
450: C
451:
452: C
453:
454: subl %eax, %esi
455:
456: sbbl %edx, %edi C n - (q1+1)*d
457: movl %esi, %edi C remainder -> n2
458: leal (%ebp,%esi), %edx
459:
460: cmovc( %edx, %edi) C n - q1*d if underflow from using q1+1
461: movd %mm0, %esi
462:
463: sbbl $0, %ebx C q
464:
465: movl %ebx, -4(%ecx)
466:
467:
468: C -----------------------------------------------------------------------------
469: L(integer_one_left):
470: C eax scratch
471: C ebx scratch (nadj, q1)
472: C ecx scratch (dst)
473: C edx scratch
474: C esi n10
475: C edi n2
476: C ebp divisor
477: C
478: C mm0 src limb, shifted
479: C mm7 rshift
480:
481:
482: movl %esi, %eax
483: movl %ebp, %ebx
484:
485: sarl $31, %eax C -n1
486: movl VAR_DST_STOP, %ecx
487:
488: andl %eax, %ebx C -n1 & d
489: negl %eax C n1
490:
491: addl %esi, %ebx C nadj = n10 + (-n1 & d), ignoring overflow
492: addl %edi, %eax C n2+n1
493:
494: mull VAR_INVERSE C m*(n2+n1)
495:
496: C
497:
498: C
499:
500: C
501:
502: addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
503: leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
504: movl %ebp, %eax C d
505:
506: C
507:
508: adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
509:
510: sbbl $0, %ebx C q1 if q1+1 overflowed
511:
512: mull %ebx
513:
514: C
515:
516: C
517:
518: C
519:
520: C
521:
522: subl %eax, %esi
523: movl PARAM_XSIZE, %eax
524:
525: sbbl %edx, %edi C n - (q1+1)*d
526: movl %esi, %edi C remainder -> n2
527: leal (%ebp,%esi), %edx
528:
529: cmovc( %edx, %edi) C n - q1*d if underflow from using q1+1
530:
531: sbbl $0, %ebx C q
532:
533: movl %ebx, -8(%ecx)
534: subl $8, %ecx
535:
536:
537:
538: orl %eax, %eax C xsize
539: jnz L(fraction_some)
540:
541: movl %edi, %eax
542: L(fraction_done):
543: movl VAR_NORM, %ecx
544: movl SAVE_EBP, %ebp
545:
546: movl SAVE_EDI, %edi
547:
548: movl SAVE_ESI, %esi
549:
550: movl SAVE_EBX, %ebx
551: addl $STACK_SPACE, %esp
552:
553: shrl %cl, %eax
554: emms
555:
556: ret
557:
558:
559: C -----------------------------------------------------------------------------
560: C
561: C Special case for q1=0xFFFFFFFF, giving q=0xFFFFFFFF meaning the low dword
562: C of q*d is simply -d and the remainder n-q*d = n10+d
563:
564: L(q1_ff):
565: C eax (divisor)
566: C ebx (q1+1 == 0)
567: C ecx
568: C edx
569: C esi n10
570: C edi n2
571: C ebp divisor
572:
573: movl VAR_DST, %ecx
574: movl VAR_DST_STOP, %edx
575: subl $4, %ecx
576:
577: movl %ecx, VAR_DST
578: psrlq %mm7, %mm0
579: leal (%ebp,%esi), %edi C n-q*d remainder -> next n2
580:
581: movl $-1, (%ecx)
582: movd %mm0, %esi C next n10
583:
584: cmpl %ecx, %edx
585: jne L(integer_top)
586:
587: jmp L(integer_loop_done)
588:
589:
590:
591: C -----------------------------------------------------------------------------
592: C
593: C In the current implementation, the following successively dependent
594: C micro-ops seem to exist.
595: C
596: C uops
597: C mul 5
598: C q1+1 1 (addl)
599: C mul 5
600: C sub 3 (negl/sbbl)
601: C addback 2 (cmov)
602: C ---
603: C 16
604: C
605: C The loop in fact runs at about 17.5 cycles. Using a sarl/andl/addl for
606: C the addback was found to be a touch slower.
607:
608:
609: ALIGN(16)
610: L(fraction_some):
611: C eax
612: C ebx
613: C ecx
614: C edx
615: C esi
616: C edi carry
617: C ebp divisor
618:
619: movl PARAM_DST, %esi
620: movl VAR_DST_STOP, %ecx
621: movl %edi, %eax
622:
623: subl $8, %ecx
624:
625:
626: ALIGN(16)
627: L(fraction_top):
628: C eax n2, then scratch
629: C ebx scratch (nadj, q1)
630: C ecx dst, decrementing
631: C edx scratch
632: C esi dst stop point
633: C edi n2
634: C ebp divisor
635:
636: mull VAR_INVERSE C m*n2
637:
638: movl %ebp, %eax C d
639: subl $4, %ecx C dst
640: leal 1(%edi), %ebx
641:
642: C
643:
644: C
645:
646: C
647:
648: addl %edx, %ebx C 1 + high(n2<<32 + m*n2) = q1+1
649:
650: mull %ebx C (q1+1)*d
651:
652: C
653:
654: C
655:
656: C
657:
658: C
659:
660: negl %eax C low of n - (q1+1)*d
661:
662: sbbl %edx, %edi C high of n - (q1+1)*d, caring only about carry
663: leal (%ebp,%eax), %edx
664:
665: cmovc( %edx, %eax) C n - q1*d if underflow from using q1+1
666:
667: sbbl $0, %ebx C q
668: movl %eax, %edi C remainder->n2
669: cmpl %esi, %ecx
670:
671: movl %ebx, (%ecx) C previous q
672: jne L(fraction_top)
673:
674:
675: jmp L(fraction_done)
676:
677: EPILOGUE()