divrem_1.asm
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:14k
源码类别:

数学计算

开发平台:

Unix_Linux

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