![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mod_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 mod_1.asm CVS log | Up to [local] / OpenXM_contrib / gmp / mpn / x86 / p6 / mmx |
1.1 maekawa 1: dnl Intel Pentium-II mpn_mod_1 -- mpn by limb remainder.
2: dnl
3: dnl P6MMX: 24.0 cycles/limb.
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_mod_1 (mp_srcptr src, mp_size_t size, mp_limb_t divisor);
30: C mp_limb_t mpn_mod_1c (mp_srcptr src, mp_size_t size, mp_limb_t divisor,
31: C mp_limb_t carry);
32: C
33: C The code here very similar to mpn_divrem_1, but with the quotient
34: C discarded. What's here probably isn't optimal.
35: C
36: C See mpn/x86/p6/mmx/divrem_1.c and mpn/x86/k7/mmx/mod_1.asm for some
37: C comments.
38:
39:
40: dnl MUL_THRESHOLD is the size at which the multiply by inverse method is
41: dnl used, rather than plain "divl"s. Minimum value 2.
42:
43: deflit(MUL_THRESHOLD, 4)
44:
45:
46: defframe(PARAM_CARRY, 16)
47: defframe(PARAM_DIVISOR,12)
48: defframe(PARAM_SIZE, 8)
49: defframe(PARAM_SRC, 4)
50:
51: defframe(SAVE_EBX, -4)
52: defframe(SAVE_ESI, -8)
53: defframe(SAVE_EDI, -12)
54: defframe(SAVE_EBP, -16)
55:
56: defframe(VAR_NORM, -20)
57: defframe(VAR_INVERSE, -24)
58: defframe(VAR_SRC_STOP,-28)
59:
60: deflit(STACK_SPACE, 28)
61:
62: .text
63: ALIGN(16)
64:
65: PROLOGUE(mpn_mod_1c)
66: deflit(`FRAME',0)
67: movl PARAM_CARRY, %edx
68: movl PARAM_SIZE, %ecx
69: subl $STACK_SPACE, %esp
70: deflit(`FRAME',STACK_SPACE)
71:
72: movl %ebp, SAVE_EBP
73: movl PARAM_DIVISOR, %ebp
74:
75: movl %esi, SAVE_ESI
76: movl PARAM_SRC, %esi
77: jmp LF(mpn_mod_1,start_1c)
78:
79: EPILOGUE()
80:
81:
82: ALIGN(16)
83: PROLOGUE(mpn_mod_1)
84: deflit(`FRAME',0)
85:
86: movl $0, %edx C initial carry (if can't skip a div)
87: movl PARAM_SIZE, %ecx
88: subl $STACK_SPACE, %esp
89: deflit(`FRAME',STACK_SPACE)
90:
91: movl %esi, SAVE_ESI
92: movl PARAM_SRC, %esi
93:
94: movl %ebp, SAVE_EBP
95: movl PARAM_DIVISOR, %ebp
96:
97: orl %ecx, %ecx
98: jz L(divide_done)
99:
100: movl -4(%esi,%ecx,4), %eax C src high limb
101:
102: cmpl %ebp, %eax C carry flag if high<divisor
103:
104: cmovc( %eax, %edx) C src high limb as initial carry
105: sbbl $0, %ecx C size-1 to skip one div
106: jz L(divide_done)
107:
108:
109: ALIGN(16)
110: L(start_1c):
111: C eax
112: C ebx
113: C ecx size
114: C edx carry
115: C esi src
116: C edi
117: C ebp divisor
118:
119: cmpl $MUL_THRESHOLD, %ecx
120: jae L(mul_by_inverse)
121:
122:
123: orl %ecx, %ecx
124: jz L(divide_done)
125:
126:
127: L(divide_top):
128: C eax scratch (quotient)
129: C ebx
130: C ecx counter, limbs, decrementing
131: C edx scratch (remainder)
132: C esi src
133: C edi
134: C ebp
135:
136: movl -4(%esi,%ecx,4), %eax
137:
138: divl %ebp
139:
140: decl %ecx
141: jnz L(divide_top)
142:
143:
144: L(divide_done):
145: movl SAVE_ESI, %esi
146: movl %edx, %eax
147:
148: movl SAVE_EBP, %ebp
149: addl $STACK_SPACE, %esp
150:
151: ret
152:
153:
154:
155: C -----------------------------------------------------------------------------
156:
157: L(mul_by_inverse):
158: C eax
159: C ebx
160: C ecx size
161: C edx carry
162: C esi src
163: C edi
164: C ebp divisor
165:
166: movl %ebx, SAVE_EBX
167: leal -4(%esi), %ebx
168:
169: movl %ebx, VAR_SRC_STOP
170: movl %ecx, %ebx C size
171:
172: movl %edi, SAVE_EDI
173: movl %edx, %edi C carry
174:
175: bsrl %ebp, %ecx C 31-l
176: movl $-1, %edx
177:
178: leal 1(%ecx), %eax C 32-l
179: xorl $31, %ecx C l
180:
181: movl %ecx, VAR_NORM
182: shll %cl, %ebp C d normalized
183:
184: movd %eax, %mm7
185: movl $-1, %eax
186: subl %ebp, %edx C (b-d)-1 so edx:eax = b*(b-d)-1
187:
188: divl %ebp C floor (b*(b-d)-1) / d
189:
190: C
191:
192: movl %eax, VAR_INVERSE
193: leal -12(%esi,%ebx,4), %eax C &src[size-3]
194:
195: movl 8(%eax), %esi C src high limb
196: movl 4(%eax), %edx C src second highest limb
197:
198: shldl( %cl, %esi, %edi) C n2 = carry,high << l
199:
200: shldl( %cl, %edx, %esi) C n10 = high,second << l
201:
202: movl %eax, %ecx C &src[size-3]
203:
204:
205: ifelse(MUL_THRESHOLD,2,`
206: cmpl $2, %ebx
207: je L(inverse_two_left)
208: ')
209:
210:
211: C The dependent chain here is the same as in mpn_divrem_1, but a few
212: C instructions are saved by not needing to store the quotient limbs. This
213: C gets it down to 24 c/l, which is still a bit away from a theoretical 19
214: C c/l.
215:
216: ALIGN(16)
217: L(inverse_top):
218: C eax scratch
219: C ebx scratch (nadj, q1)
220: C ecx src pointer, decrementing
221: C edx scratch
222: C esi n10
223: C edi n2
224: C ebp divisor
225: C
226: C mm0 scratch (src qword)
227: C mm7 rshift for normalization
228:
229:
230: movl %esi, %eax
231: movl %ebp, %ebx
232:
233: sarl $31, %eax C -n1
234:
235: andl %eax, %ebx C -n1 & d
236: negl %eax C n1
237:
238: addl %esi, %ebx C nadj = n10 + (-n1 & d), ignoring overflow
239: addl %edi, %eax C n2+n1
240:
241: mull VAR_INVERSE C m*(n2+n1)
242:
243: movq (%ecx), %mm0 C next src limb and the one below it
244: subl $4, %ecx
245:
246: C
247:
248: C
249:
250: C
251:
252: addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
253: leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
254: movl %ebp, %eax C d
255:
256: adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
257: jz L(q1_ff)
258:
259: mull %ebx C (q1+1)*d
260:
261: psrlq %mm7, %mm0
262: movl VAR_SRC_STOP, %ebx
263:
264: C
265:
266: C
267:
268: C
269:
270: subl %eax, %esi
271:
272: sbbl %edx, %edi C n - (q1+1)*d
273: movl %esi, %edi C remainder -> n2
274: leal (%ebp,%esi), %edx
275:
276: cmovc( %edx, %edi) C n - q1*d if underflow from using q1+1
277: movd %mm0, %esi
278: cmpl %ebx, %ecx
279:
280: jne L(inverse_top)
281:
282:
283: L(inverse_loop_done):
284:
285:
286: C -----------------------------------------------------------------------------
287:
288: L(inverse_two_left):
289: C eax scratch
290: C ebx scratch (nadj, q1)
291: C ecx &src[-1]
292: C edx scratch
293: C esi n10
294: C edi n2
295: C ebp divisor
296: C
297: C mm0 scratch (src dword)
298: C mm7 rshift
299:
300: movl %esi, %eax
301: movl %ebp, %ebx
302:
303: sarl $31, %eax C -n1
304:
305: andl %eax, %ebx C -n1 & d
306: negl %eax C n1
307:
308: addl %esi, %ebx C nadj = n10 + (-n1 & d), ignoring overflow
309: addl %edi, %eax C n2+n1
310:
311: mull VAR_INVERSE C m*(n2+n1)
312:
313: movd 4(%ecx), %mm0 C src low limb
314:
315: C
316:
317: C
318:
319: C
320:
321: addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
322: leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
323:
324: adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
325:
326: sbbl $0, %ebx
327: movl %ebp, %eax C d
328:
329: mull %ebx C (q1+1)*d
330:
331: psllq $32, %mm0
332:
333: psrlq %mm7, %mm0
334:
335: C
336:
337: C
338:
339: subl %eax, %esi
340:
341: sbbl %edx, %edi C n - (q1+1)*d
342: movl %esi, %edi C remainder -> n2
343: leal (%ebp,%esi), %edx
344:
345: cmovc( %edx, %edi) C n - q1*d if underflow from using q1+1
346: movd %mm0, %esi
347:
348:
349: C One limb left
350:
351: C eax scratch
352: C ebx scratch (nadj, q1)
353: C ecx
354: C edx scratch
355: C esi n10
356: C edi n2
357: C ebp divisor
358: C
359: C mm0 src limb, shifted
360: C mm7 rshift
361:
362: movl %esi, %eax
363: movl %ebp, %ebx
364:
365: sarl $31, %eax C -n1
366:
367: andl %eax, %ebx C -n1 & d
368: negl %eax C n1
369:
370: addl %esi, %ebx C nadj = n10 + (-n1 & d), ignoring overflow
371: addl %edi, %eax C n2+n1
372:
373: mull VAR_INVERSE C m*(n2+n1)
374:
375: movl VAR_NORM, %ecx C for final denorm
376:
377: C
378:
379: C
380:
381: C
382:
383: addl %ebx, %eax C m*(n2+n1) + nadj, low giving carry flag
384: leal 1(%edi), %ebx C n2<<32 + m*(n2+n1))
385:
386: adcl %edx, %ebx C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
387:
388: sbbl $0, %ebx
389: movl %ebp, %eax C d
390:
391: mull %ebx C (q1+1)*d
392:
393: movl SAVE_EBX, %ebx
394:
395: C
396:
397: C
398:
399: C
400:
401: subl %eax, %esi
402:
403: sbbl %edx, %edi C n - (q1+1)*d
404: leal (%ebp,%esi), %edx
405: movl SAVE_EBP, %ebp
406:
407: movl %esi, %eax C remainder
408: movl SAVE_ESI, %esi
409:
410: cmovc( %edx, %eax) C n - q1*d if underflow from using q1+1
411: movl SAVE_EDI, %edi
412:
413: shrl %cl, %eax C denorm remainder
414: addl $STACK_SPACE, %esp
415: emms
416:
417: ret
418:
419:
420: C -----------------------------------------------------------------------------
421: C
422: C Special case for q1=0xFFFFFFFF, giving q=0xFFFFFFFF meaning the low dword
423: C of q*d is simply -d and the remainder n-q*d = n10+d
424:
425: L(q1_ff):
426: C eax (divisor)
427: C ebx (q1+1 == 0)
428: C ecx src pointer
429: C edx
430: C esi n10
431: C edi (n2)
432: C ebp divisor
433:
434: leal (%ebp,%esi), %edi C n-q*d remainder -> next n2
435: movl VAR_SRC_STOP, %edx
436: psrlq %mm7, %mm0
437:
438: movd %mm0, %esi C next n10
439: cmpl %ecx, %edx
440: jne L(inverse_top)
441:
442: jmp L(inverse_loop_done)
443:
444: EPILOGUE()