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

数学计算

开发平台:

Unix_Linux

  1. dnl  AMD K7 mpn_divrem_1, mpn_divrem_1c, mpn_preinv_divrem_1 -- mpn by limb
  2. dnl  division.
  3. dnl  Copyright 1999, 2000, 2001, 2002, 2004 Free Software Foundation, Inc.
  4. dnl
  5. dnl  This file is part of the GNU MP Library.
  6. dnl
  7. dnl  The GNU MP Library is free software; you can redistribute it and/or
  8. dnl  modify it under the terms of the GNU Lesser General Public License as
  9. dnl  published by the Free Software Foundation; either version 3 of the
  10. dnl  License, or (at your option) any later version.
  11. dnl
  12. dnl  The GNU MP Library is distributed in the hope that it will be useful,
  13. dnl  but WITHOUT ANY WARRANTY; without even the implied warranty of
  14. dnl  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15. dnl  Lesser General Public License for more details.
  16. dnl
  17. dnl  You should have received a copy of the GNU Lesser General Public License
  18. dnl  along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.
  19. include(`../config.m4')
  20. C K7: 17.0 cycles/limb integer part, 15.0 cycles/limb fraction part.
  21. C mp_limb_t mpn_divrem_1 (mp_ptr dst, mp_size_t xsize,
  22. C                         mp_srcptr src, mp_size_t size,
  23. C                         mp_limb_t divisor);
  24. C mp_limb_t mpn_divrem_1c (mp_ptr dst, mp_size_t xsize,
  25. C                          mp_srcptr src, mp_size_t size,
  26. C                          mp_limb_t divisor, mp_limb_t carry);
  27. C mp_limb_t mpn_preinv_divrem_1 (mp_ptr dst, mp_size_t xsize,
  28. C                                mp_srcptr src, mp_size_t size,
  29. C                                mp_limb_t divisor, mp_limb_t inverse,
  30. C                                unsigned shift);
  31. C
  32. C Algorithm:
  33. C
  34. C The method and nomenclature follow part 8 of "Division by Invariant
  35. C Integers using Multiplication" by Granlund and Montgomery, reference in
  36. C gmp.texi.
  37. C
  38. C The "and"s shown in the paper are done here with "cmov"s.  "m" is written
  39. C for m', and "d" for d_norm, which won't cause any confusion since it's
  40. C only the normalized divisor that's of any use in the code.  "b" is written
  41. C for 2^N, the size of a limb, N being 32 here.
  42. C
  43. C The step "sdword dr = n - 2^N*d + (2^N-1-q1) * d" is instead done as
  44. C "n-(q1+1)*d"; this rearrangement gives the same two-limb answer.  If
  45. C q1==0xFFFFFFFF, then q1+1 would overflow.  We branch to a special case
  46. C "q1_ff" if this occurs.  Since the true quotient is either q1 or q1+1 then
  47. C if q1==0xFFFFFFFF that must be the right value.
  48. C
  49. C For the last and second last steps q1==0xFFFFFFFF is instead handled by an
  50. C sbbl to go back to 0xFFFFFFFF if an overflow occurs when adding 1.  This
  51. C then goes through as normal, and finding no addback required.  sbbl costs
  52. C an extra cycle over what the main loop code does, but it keeps code size
  53. C and complexity down.
  54. C
  55. C Notes:
  56. C
  57. C mpn_divrem_1 and mpn_preinv_divrem_1 avoid one division if the src high
  58. C limb is less than the divisor.  mpn_divrem_1c doesn't check for a zero
  59. C carry, since in normal circumstances that will be a very rare event.
  60. C
  61. C The test for skipping a division is branch free (once size>=1 is tested).
  62. C The store to the destination high limb is 0 when a divide is skipped, or
  63. C if it's not skipped then a copy of the src high limb is used.  The latter
  64. C is in case src==dst.
  65. C
  66. C There's a small bias towards expecting xsize==0, by having code for
  67. C xsize==0 in a straight line and xsize!=0 under forward jumps.
  68. C
  69. C Alternatives:
  70. C
  71. C If the divisor is normalized (high bit set) then a division step can
  72. C always be skipped, since the high destination limb is always 0 or 1 in
  73. C that case.  It doesn't seem worth checking for this though, since it
  74. C probably occurs infrequently, in particular note that big_base for a
  75. C decimal mpn_get_str is not normalized in a 32-bit limb.
  76. dnl  MUL_THRESHOLD is the value of xsize+size at which the multiply by
  77. dnl  inverse method is used, rather than plain "divl"s.  Minimum value 1.
  78. dnl
  79. dnl  The inverse takes about 50 cycles to calculate, but after that the
  80. dnl  multiply is 17 c/l versus division at 42 c/l.
  81. dnl
  82. dnl  At 3 limbs the mul is a touch faster than div on the integer part, and
  83. dnl  even more so on the fractional part.
  84. deflit(MUL_THRESHOLD, 3)
  85. defframe(PARAM_PREINV_SHIFT,   28)  dnl mpn_preinv_divrem_1
  86. defframe(PARAM_PREINV_INVERSE, 24)  dnl mpn_preinv_divrem_1
  87. defframe(PARAM_CARRY,  24)          dnl mpn_divrem_1c
  88. defframe(PARAM_DIVISOR,20)
  89. defframe(PARAM_SIZE,   16)
  90. defframe(PARAM_SRC,    12)
  91. defframe(PARAM_XSIZE,  8)
  92. defframe(PARAM_DST,    4)
  93. defframe(SAVE_EBX,    -4)
  94. defframe(SAVE_ESI,    -8)
  95. defframe(SAVE_EDI,    -12)
  96. defframe(SAVE_EBP,    -16)
  97. defframe(VAR_NORM,    -20)
  98. defframe(VAR_INVERSE, -24)
  99. defframe(VAR_SRC,     -28)
  100. defframe(VAR_DST,     -32)
  101. defframe(VAR_DST_STOP,-36)
  102. deflit(STACK_SPACE, 36)
  103. TEXT
  104. ALIGN(32)
  105. PROLOGUE(mpn_preinv_divrem_1)
  106. deflit(`FRAME',0)
  107. movl PARAM_XSIZE, %ecx
  108. movl PARAM_DST, %edx
  109. subl $STACK_SPACE, %esp FRAME_subl_esp(STACK_SPACE)
  110. movl %esi, SAVE_ESI
  111. movl PARAM_SRC, %esi
  112. movl %ebx, SAVE_EBX
  113. movl PARAM_SIZE, %ebx
  114. leal 8(%edx,%ecx,4), %edx C &dst[xsize+2]
  115. movl %ebp, SAVE_EBP
  116. movl PARAM_DIVISOR, %ebp
  117. movl %edx, VAR_DST_STOP C &dst[xsize+2]
  118. movl %edi, SAVE_EDI
  119. xorl %edi, %edi C carry
  120. movl -4(%esi,%ebx,4), %eax C src high limb
  121. xor %ecx, %ecx
  122. C
  123. C
  124. cmpl %ebp, %eax C high cmp divisor
  125. cmovc( %eax, %edi) C high is carry if high<divisor
  126. cmovnc( %eax, %ecx) C 0 if skip div, src high if not
  127. C (the latter in case src==dst)
  128. movl %ecx, -12(%edx,%ebx,4) C dst high limb
  129. sbbl $0, %ebx C skip one division if high<divisor
  130. movl PARAM_PREINV_SHIFT, %ecx
  131. leal -8(%edx,%ebx,4), %edx C &dst[xsize+size]
  132. movl $32, %eax
  133. movl %edx, VAR_DST C &dst[xsize+size]
  134. shll %cl, %ebp C d normalized
  135. subl %ecx, %eax
  136. movl %ecx, VAR_NORM
  137. movd %eax, %mm7 C rshift
  138. movl PARAM_PREINV_INVERSE, %eax
  139. jmp L(start_preinv)
  140. EPILOGUE()
  141. ALIGN(16)
  142. PROLOGUE(mpn_divrem_1c)
  143. deflit(`FRAME',0)
  144. movl PARAM_CARRY, %edx
  145. movl PARAM_SIZE, %ecx
  146. subl $STACK_SPACE, %esp
  147. deflit(`FRAME',STACK_SPACE)
  148. movl %ebx, SAVE_EBX
  149. movl PARAM_XSIZE, %ebx
  150. movl %edi, SAVE_EDI
  151. movl PARAM_DST, %edi
  152. movl %ebp, SAVE_EBP
  153. movl PARAM_DIVISOR, %ebp
  154. movl %esi, SAVE_ESI
  155. movl PARAM_SRC, %esi
  156. leal -4(%edi,%ebx,4), %edi C &dst[xsize-1]
  157. jmp L(start_1c)
  158. EPILOGUE()
  159. C offset 0xa1, close enough to aligned
  160. PROLOGUE(mpn_divrem_1)
  161. deflit(`FRAME',0)
  162. movl PARAM_SIZE, %ecx
  163. movl $0, %edx C initial carry (if can't skip a div)
  164. subl $STACK_SPACE, %esp
  165. deflit(`FRAME',STACK_SPACE)
  166. movl %esi, SAVE_ESI
  167. movl PARAM_SRC, %esi
  168. movl %ebx, SAVE_EBX
  169. movl PARAM_XSIZE, %ebx
  170. movl %ebp, SAVE_EBP
  171. movl PARAM_DIVISOR, %ebp
  172. orl %ecx, %ecx C size
  173. movl %edi, SAVE_EDI
  174. movl PARAM_DST, %edi
  175. leal -4(%edi,%ebx,4), %edi C &dst[xsize-1]
  176. jz L(no_skip_div) C if size==0
  177. movl -4(%esi,%ecx,4), %eax C src high limb
  178. xorl %esi, %esi
  179. cmpl %ebp, %eax C high cmp divisor
  180. cmovc( %eax, %edx) C high is carry if high<divisor
  181. cmovnc( %eax, %esi) C 0 if skip div, src high if not
  182. movl %esi, (%edi,%ecx,4) C dst high limb
  183. sbbl $0, %ecx C size-1 if high<divisor
  184. movl PARAM_SRC, %esi C reload
  185. L(no_skip_div):
  186. L(start_1c):
  187. C eax
  188. C ebx xsize
  189. C ecx size
  190. C edx carry
  191. C esi src
  192. C edi &dst[xsize-1]
  193. C ebp divisor
  194. leal (%ebx,%ecx), %eax C size+xsize
  195. cmpl $MUL_THRESHOLD, %eax
  196. jae L(mul_by_inverse)
  197. C With MUL_THRESHOLD set to 3, the simple loops here only do 0 to 2 limbs.
  198. C It'd be possible to write them out without the looping, but no speedup
  199. C would be expected.
  200. C
  201. C Using PARAM_DIVISOR instead of %ebp measures 1 cycle/loop faster on the
  202. C integer part, but curiously not on the fractional part, where %ebp is a
  203. C (fixed) couple of cycles faster.
  204. orl %ecx, %ecx
  205. jz L(divide_no_integer)
  206. L(divide_integer):
  207. C eax scratch (quotient)
  208. C ebx xsize
  209. C ecx counter
  210. C edx scratch (remainder)
  211. C esi src
  212. C edi &dst[xsize-1]
  213. C ebp divisor
  214. movl -4(%esi,%ecx,4), %eax
  215. divl PARAM_DIVISOR
  216. movl %eax, (%edi,%ecx,4)
  217. decl %ecx
  218. jnz L(divide_integer)
  219. L(divide_no_integer):
  220. movl PARAM_DST, %edi
  221. orl %ebx, %ebx
  222. jnz L(divide_fraction)
  223. L(divide_done):
  224. movl SAVE_ESI, %esi
  225. movl SAVE_EDI, %edi
  226. movl %edx, %eax
  227. movl SAVE_EBX, %ebx
  228. movl SAVE_EBP, %ebp
  229. addl $STACK_SPACE, %esp
  230. ret
  231. L(divide_fraction):
  232. C eax scratch (quotient)
  233. C ebx counter
  234. C ecx
  235. C edx scratch (remainder)
  236. C esi
  237. C edi dst
  238. C ebp divisor
  239. movl $0, %eax
  240. divl %ebp
  241. movl %eax, -4(%edi,%ebx,4)
  242. decl %ebx
  243. jnz L(divide_fraction)
  244. jmp L(divide_done)
  245. C -----------------------------------------------------------------------------
  246. L(mul_by_inverse):
  247. C eax
  248. C ebx xsize
  249. C ecx size
  250. C edx carry
  251. C esi src
  252. C edi &dst[xsize-1]
  253. C ebp divisor
  254. bsrl %ebp, %eax C 31-l
  255. leal 12(%edi), %ebx C &dst[xsize+2], loop dst stop
  256. leal 4(%edi,%ecx,4), %edi C &dst[xsize+size]
  257. movl %edi, VAR_DST
  258. movl %ebx, VAR_DST_STOP
  259. movl %ecx, %ebx C size
  260. movl $31, %ecx
  261. movl %edx, %edi C carry
  262. movl $-1, %edx
  263. C
  264. xorl %eax, %ecx C l
  265. incl %eax C 32-l
  266. shll %cl, %ebp C d normalized
  267. movl %ecx, VAR_NORM
  268. movd %eax, %mm7
  269. movl $-1, %eax
  270. subl %ebp, %edx C (b-d)-1 giving edx:eax = b*(b-d)-1
  271. divl %ebp C floor (b*(b-d)-1) / d
  272. L(start_preinv):
  273. C eax inverse
  274. C ebx size
  275. C ecx shift
  276. C edx
  277. C esi src
  278. C edi carry
  279. C ebp divisor
  280. C
  281. C mm7 rshift
  282. orl %ebx, %ebx C size
  283. movl %eax, VAR_INVERSE
  284. leal -12(%esi,%ebx,4), %eax C &src[size-3]
  285. jz L(start_zero)
  286. movl %eax, VAR_SRC
  287. cmpl $1, %ebx
  288. movl 8(%eax), %esi C src high limb
  289. jz L(start_one)
  290. L(start_two_or_more):
  291. movl 4(%eax), %edx C src second highest limb
  292. shldl( %cl, %esi, %edi) C n2 = carry,high << l
  293. shldl( %cl, %edx, %esi) C n10 = high,second << l
  294. cmpl $2, %ebx
  295. je L(integer_two_left)
  296. jmp L(integer_top)
  297. L(start_one):
  298. shldl( %cl, %esi, %edi) C n2 = carry,high << l
  299. shll %cl, %esi C n10 = high << l
  300. movl %eax, VAR_SRC
  301. jmp L(integer_one_left)
  302. L(start_zero):
  303. C Can be here with xsize==0 if mpn_preinv_divrem_1 had size==1 and
  304. C skipped a division.
  305. shll %cl, %edi C n2 = carry << l
  306. movl %edi, %eax C return value for zero_done
  307. cmpl $0, PARAM_XSIZE
  308. je L(zero_done)
  309. jmp L(fraction_some)
  310. C -----------------------------------------------------------------------------
  311. C
  312. C The multiply by inverse loop is 17 cycles, and relies on some out-of-order
  313. C execution.  The instruction scheduling is important, with various
  314. C apparently equivalent forms running 1 to 5 cycles slower.
  315. C
  316. C A lower bound for the time would seem to be 16 cycles, based on the
  317. C following successive dependencies.
  318. C
  319. C       cycles
  320. C n2+n1 1
  321. C mul 6
  322. C q1+1 1
  323. C mul 6
  324. C sub 1
  325. C addback 1
  326. C        ---
  327. C        16
  328. C
  329. C This chain is what the loop has already, but 16 cycles isn't achieved.
  330. C K7 has enough decode, and probably enough execute (depending maybe on what
  331. C a mul actually consumes), but nothing running under 17 has been found.
  332. C
  333. C In theory n2+n1 could be done in the sub and addback stages (by
  334. C calculating both n2 and n2+n1 there), but lack of registers makes this an
  335. C unlikely proposition.
  336. C
  337. C The jz in the loop keeps the q1+1 stage to 1 cycle.  Handling an overflow
  338. C from q1+1 with an "sbbl $0, %ebx" would add a cycle to the dependent
  339. C chain, and nothing better than 18 cycles has been found when using it.
  340. C The jump is taken only when q1 is 0xFFFFFFFF, and on random data this will
  341. C be an extremely rare event.
  342. C
  343. C Branch mispredictions will hit random occurrances of q1==0xFFFFFFFF, but
  344. C if some special data is coming out with this always, the q1_ff special
  345. C case actually runs at 15 c/l.  0x2FFF...FFFD divided by 3 is a good way to
  346. C induce the q1_ff case, for speed measurements or testing.  Note that
  347. C 0xFFF...FFF divided by 1 or 2 doesn't induce it.
  348. C
  349. C The instruction groupings and empty comments show the cycles for a naive
  350. C in-order view of the code (conveniently ignoring the load latency on
  351. C VAR_INVERSE).  This shows some of where the time is going, but is nonsense
  352. C to the extent that out-of-order execution rearranges it.  In this case
  353. C there's 19 cycles shown, but it executes at 17.
  354. ALIGN(16)
  355. L(integer_top):
  356. C eax scratch
  357. C ebx scratch (nadj, q1)
  358. C ecx scratch (src, dst)
  359. C edx scratch
  360. C esi n10
  361. C edi n2
  362. C ebp divisor
  363. C
  364. C mm0 scratch (src qword)
  365. C mm7 rshift for normalization
  366. cmpl $0x80000000, %esi  C n1 as 0=c, 1=nc
  367. movl %edi, %eax         C n2
  368. movl VAR_SRC, %ecx
  369. leal (%ebp,%esi), %ebx
  370. cmovc( %esi, %ebx)    C nadj = n10 + (-n1 & d), ignoring overflow
  371. sbbl $-1, %eax          C n2+n1
  372. mull VAR_INVERSE        C m*(n2+n1)
  373. movq (%ecx), %mm0       C next limb and the one below it
  374. subl $4, %ecx
  375. movl %ecx, VAR_SRC
  376. C
  377. addl %ebx, %eax         C m*(n2+n1) + nadj, low giving carry flag
  378. leal 1(%edi), %ebx      C n2+1
  379. movl %ebp, %eax    C d
  380. C
  381. adcl %edx, %ebx         C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
  382. jz L(q1_ff)
  383. movl VAR_DST, %ecx
  384. mull %ebx    C (q1+1)*d
  385. psrlq %mm7, %mm0
  386. leal -4(%ecx), %ecx
  387. C
  388. subl %eax, %esi
  389. movl VAR_DST_STOP, %eax
  390. C
  391. sbbl %edx, %edi    C n - (q1+1)*d
  392. movl %esi, %edi    C remainder -> n2
  393. leal (%ebp,%esi), %edx
  394. movd %mm0, %esi
  395. cmovc( %edx, %edi)    C n - q1*d if underflow from using q1+1
  396. sbbl $0, %ebx    C q
  397. cmpl %eax, %ecx
  398. movl %ebx, (%ecx)
  399. movl %ecx, VAR_DST
  400. jne L(integer_top)
  401. L(integer_loop_done):
  402. C -----------------------------------------------------------------------------
  403. C
  404. C Here, and in integer_one_left below, an sbbl $0 is used rather than a jz
  405. C q1_ff special case.  This make the code a bit smaller and simpler, and
  406. C costs only 1 cycle (each).
  407. L(integer_two_left):
  408. C eax scratch
  409. C ebx scratch (nadj, q1)
  410. C ecx scratch (src, dst)
  411. C edx scratch
  412. C esi n10
  413. C edi n2
  414. C ebp divisor
  415. C
  416. C mm7 rshift
  417. cmpl $0x80000000, %esi  C n1 as 0=c, 1=nc
  418. movl %edi, %eax         C n2
  419. movl PARAM_SRC, %ecx
  420. leal (%ebp,%esi), %ebx
  421. cmovc( %esi, %ebx)    C nadj = n10 + (-n1 & d), ignoring overflow
  422. sbbl $-1, %eax          C n2+n1
  423. mull VAR_INVERSE        C m*(n2+n1)
  424. movd (%ecx), %mm0    C src low limb
  425. movl VAR_DST_STOP, %ecx
  426. C
  427. addl %ebx, %eax         C m*(n2+n1) + nadj, low giving carry flag
  428. leal 1(%edi), %ebx      C n2+1
  429. movl %ebp, %eax    C d
  430. adcl %edx, %ebx         C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
  431. sbbl $0, %ebx
  432. mull %ebx    C (q1+1)*d
  433. psllq $32, %mm0
  434. psrlq %mm7, %mm0
  435. C
  436. subl %eax, %esi
  437. C
  438. sbbl %edx, %edi    C n - (q1+1)*d
  439. movl %esi, %edi    C remainder -> n2
  440. leal (%ebp,%esi), %edx
  441. movd %mm0, %esi
  442. cmovc( %edx, %edi)    C n - q1*d if underflow from using q1+1
  443. sbbl $0, %ebx    C q
  444. movl %ebx, -4(%ecx)
  445. C -----------------------------------------------------------------------------
  446. L(integer_one_left):
  447. C eax scratch
  448. C ebx scratch (nadj, q1)
  449. C ecx dst
  450. C edx scratch
  451. C esi n10
  452. C edi n2
  453. C ebp divisor
  454. C
  455. C mm7 rshift
  456. movl VAR_DST_STOP, %ecx
  457. cmpl $0x80000000, %esi  C n1 as 0=c, 1=nc
  458. movl %edi, %eax         C n2
  459. leal (%ebp,%esi), %ebx
  460. cmovc( %esi, %ebx)    C nadj = n10 + (-n1 & d), ignoring overflow
  461. sbbl $-1, %eax          C n2+n1
  462. mull VAR_INVERSE        C m*(n2+n1)
  463. C
  464. C
  465. C
  466. addl %ebx, %eax         C m*(n2+n1) + nadj, low giving carry flag
  467. leal 1(%edi), %ebx      C n2+1
  468. movl %ebp, %eax    C d
  469. C
  470. adcl %edx, %ebx         C 1 + high(n2<<32 + m*(n2+n1) + nadj) = q1+1
  471. sbbl $0, %ebx           C q1 if q1+1 overflowed
  472. mull %ebx
  473. C
  474. C
  475. C
  476. subl %eax, %esi
  477. C
  478. sbbl %edx, %edi    C n - (q1+1)*d
  479. movl %esi, %edi    C remainder -> n2
  480. leal (%ebp,%esi), %edx
  481. cmovc( %edx, %edi)    C n - q1*d if underflow from using q1+1
  482. sbbl $0, %ebx    C q
  483. movl %ebx, -8(%ecx)
  484. subl $8, %ecx
  485. L(integer_none):
  486. cmpl $0, PARAM_XSIZE
  487. jne L(fraction_some)
  488. movl %edi, %eax
  489. L(fraction_done):
  490. movl VAR_NORM, %ecx
  491. L(zero_done):
  492. movl SAVE_EBP, %ebp
  493. movl SAVE_EDI, %edi
  494. movl SAVE_ESI, %esi
  495. movl SAVE_EBX, %ebx
  496. addl $STACK_SPACE, %esp
  497. shrl %cl, %eax
  498. emms
  499. ret
  500. C -----------------------------------------------------------------------------
  501. C
  502. C Special case for q1=0xFFFFFFFF, giving q=0xFFFFFFFF meaning the low dword
  503. C of q*d is simply -d and the remainder n-q*d = n10+d
  504. L(q1_ff):
  505. C eax (divisor)
  506. C ebx (q1+1 == 0)
  507. C ecx
  508. C edx
  509. C esi n10
  510. C edi n2
  511. C ebp divisor
  512. movl VAR_DST, %ecx
  513. movl VAR_DST_STOP, %edx
  514. subl $4, %ecx
  515. psrlq %mm7, %mm0
  516. leal (%ebp,%esi), %edi C n-q*d remainder -> next n2
  517. movl %ecx, VAR_DST
  518. movd %mm0, %esi C next n10
  519. movl $-1, (%ecx)
  520. cmpl %ecx, %edx
  521. jne L(integer_top)
  522. jmp L(integer_loop_done)
  523. C -----------------------------------------------------------------------------
  524. C
  525. C Being the fractional part, the "source" limbs are all zero, meaning
  526. C n10=0, n1=0, and hence nadj=0, leading to many instructions eliminated.
  527. C
  528. C The loop runs at 15 cycles.  The dependent chain is the same as the
  529. C general case above, but without the n2+n1 stage (due to n1==0), so 15
  530. C would seem to be the lower bound.
  531. C
  532. C A not entirely obvious simplification is that q1+1 never overflows a limb,
  533. C and so there's no need for the sbbl $0 or jz q1_ff from the general case.
  534. C q1 is the high word of m*n2+b*n2 and the following shows q1<=b-2 always.
  535. C rnd() means rounding down to a multiple of d.
  536. C
  537. C m*n2 + b*n2 <= m*(d-1) + b*(d-1)
  538. C              = m*d + b*d - m - b
  539. C              = floor((b(b-d)-1)/d)*d + b*d - m - b
  540. C              = rnd(b(b-d)-1) + b*d - m - b
  541. C              = rnd(b(b-d)-1 + b*d) - m - b
  542. C              = rnd(b*b-1) - m - b
  543. C              <= (b-2)*b
  544. C
  545. C Unchanged from the general case is that the final quotient limb q can be
  546. C either q1 or q1+1, and the q1+1 case occurs often.  This can be seen from
  547. C equation 8.4 of the paper which simplifies as follows when n1==0 and
  548. C n0==0.
  549. C
  550. C n-q1*d = (n2*k+q0*d)/b <= d + (d*d-2d)/b
  551. C
  552. C As before, the instruction groupings and empty comments show a naive
  553. C in-order view of the code, which is made a nonsense by out of order
  554. C execution.  There's 17 cycles shown, but it executes at 15.
  555. C
  556. C Rotating the store q and remainder->n2 instructions up to the top of the
  557. C loop gets the run time down from 16 to 15.
  558. ALIGN(16)
  559. L(fraction_some):
  560. C eax
  561. C ebx
  562. C ecx
  563. C edx
  564. C esi
  565. C edi carry
  566. C ebp divisor
  567. movl PARAM_DST, %esi
  568. movl VAR_DST_STOP, %ecx C &dst[xsize+2]
  569. movl %edi, %eax
  570. subl $8, %ecx C &dst[xsize]
  571. jmp L(fraction_entry)
  572. ALIGN(16)
  573. L(fraction_top):
  574. C eax n2 carry, then scratch
  575. C ebx scratch (nadj, q1)
  576. C ecx dst, decrementing
  577. C edx scratch
  578. C esi dst stop point
  579. C edi (will be n2)
  580. C ebp divisor
  581. movl %ebx, (%ecx) C previous q
  582. movl %eax, %edi C remainder->n2
  583. L(fraction_entry):
  584. mull VAR_INVERSE C m*n2
  585. movl %ebp, %eax C d
  586. subl $4, %ecx C dst
  587. leal 1(%edi), %ebx
  588. C
  589. C
  590. C
  591. C
  592. addl %edx, %ebx C 1 + high(n2<<32 + m*n2) = q1+1
  593. mull %ebx C (q1+1)*d
  594. C
  595. C
  596. C
  597. negl %eax C low of n - (q1+1)*d
  598. C
  599. sbbl %edx, %edi C high of n - (q1+1)*d, caring only about carry
  600. leal (%ebp,%eax), %edx
  601. cmovc( %edx, %eax) C n - q1*d if underflow from using q1+1
  602. sbbl $0, %ebx C q
  603. cmpl %esi, %ecx
  604. jne L(fraction_top)
  605. movl %ebx, (%ecx)
  606. jmp L(fraction_done)
  607. EPILOGUE()