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