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

数学计算

开发平台:

Unix_Linux

  1. dnl  AMD K7 mpn_mul_basecase -- multiply two mpn numbers.
  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 K7: approx 4.42 cycles per cross product at around 20x20 limbs (16
  20. C     limbs/loop unrolling).
  21. dnl  K7 UNROLL_COUNT cycles/product (at around 20x20)
  22. dnl           8           4.67
  23. dnl          16           4.59
  24. dnl          32           4.42
  25. dnl  Maximum possible with the current code is 32.
  26. dnl
  27. dnl  At 32 the typical 13-26 limb sizes from the karatsuba code will get
  28. dnl  done with a straight run through a block of code, no inner loop.  Using
  29. dnl  32 gives 1k of code, but the k7 has a 64k L1 code cache.
  30. deflit(UNROLL_COUNT, 32)
  31. C void mpn_mul_basecase (mp_ptr wp,
  32. C                        mp_srcptr xp, mp_size_t xsize,
  33. C                        mp_srcptr yp, mp_size_t ysize);
  34. C
  35. C Calculate xp,xsize multiplied by yp,ysize, storing the result in
  36. C wp,xsize+ysize.
  37. C
  38. C This routine is essentially the same as mpn/generic/mul_basecase.c, but
  39. C it's faster because it does most of the mpn_addmul_1() startup
  40. C calculations only once.  The saving is 15-25% on typical sizes coming from
  41. C the Karatsuba multiply code.
  42. ifdef(`PIC',`
  43. deflit(UNROLL_THRESHOLD, 5)
  44. ',`
  45. deflit(UNROLL_THRESHOLD, 5)
  46. ')
  47. defframe(PARAM_YSIZE,20)
  48. defframe(PARAM_YP,   16)
  49. defframe(PARAM_XSIZE,12)
  50. defframe(PARAM_XP,   8)
  51. defframe(PARAM_WP,   4)
  52. TEXT
  53. ALIGN(32)
  54. PROLOGUE(mpn_mul_basecase)
  55. deflit(`FRAME',0)
  56. movl PARAM_XSIZE, %ecx
  57. movl PARAM_YP, %eax
  58. movl PARAM_XP, %edx
  59. movl (%eax), %eax C yp low limb
  60. cmpl $2, %ecx
  61. ja L(xsize_more_than_two)
  62. je L(two_by_something)
  63. C one limb by one limb
  64. mull (%edx)
  65. movl PARAM_WP, %ecx
  66. movl %eax, (%ecx)
  67. movl %edx, 4(%ecx)
  68. ret
  69. C -----------------------------------------------------------------------------
  70. L(two_by_something):
  71. deflit(`FRAME',0)
  72. decl PARAM_YSIZE
  73. pushl %ebx defframe_pushl(`SAVE_EBX')
  74. movl %eax, %ecx C yp low limb
  75. movl PARAM_WP, %ebx
  76. pushl %esi defframe_pushl(`SAVE_ESI')
  77. movl %edx, %esi C xp
  78. movl (%edx), %eax C xp low limb
  79. jnz L(two_by_two)
  80. C two limbs by one limb
  81. mull %ecx
  82. movl %eax, (%ebx)
  83. movl 4(%esi), %eax
  84. movl %edx, %esi C carry
  85. mull %ecx
  86. addl %eax, %esi
  87. movl %esi, 4(%ebx)
  88. movl SAVE_ESI, %esi
  89. adcl $0, %edx
  90. movl %edx, 8(%ebx)
  91. movl SAVE_EBX, %ebx
  92. addl $FRAME, %esp
  93. ret
  94. C -----------------------------------------------------------------------------
  95. C Could load yp earlier into another register.
  96. ALIGN(16)
  97. L(two_by_two):
  98. C eax xp low limb
  99. C ebx wp
  100. C ecx yp low limb
  101. C edx
  102. C esi xp
  103. C edi
  104. C ebp
  105. dnl  FRAME carries on from previous
  106. mull %ecx C xp[0] * yp[0]
  107. push %edi defframe_pushl(`SAVE_EDI')
  108. movl %edx, %edi C carry, for wp[1]
  109. movl %eax, (%ebx)
  110. movl 4(%esi), %eax
  111. mull %ecx C xp[1] * yp[0]
  112. addl %eax, %edi
  113. movl PARAM_YP, %ecx
  114. adcl $0, %edx
  115. movl 4(%ecx), %ecx C yp[1]
  116. movl %edi, 4(%ebx)
  117. movl 4(%esi), %eax C xp[1]
  118. movl %edx, %edi C carry, for wp[2]
  119. mull %ecx C xp[1] * yp[1]
  120. addl %eax, %edi
  121. adcl $0, %edx
  122. movl (%esi), %eax C xp[0]
  123. movl %edx, %esi C carry, for wp[3]
  124. mull %ecx C xp[0] * yp[1]
  125. addl %eax, 4(%ebx)
  126. adcl %edx, %edi
  127. movl %edi, 8(%ebx)
  128. adcl $0, %esi
  129. movl SAVE_EDI, %edi
  130. movl %esi, 12(%ebx)
  131. movl SAVE_ESI, %esi
  132. movl SAVE_EBX, %ebx
  133. addl $FRAME, %esp
  134. ret
  135. C -----------------------------------------------------------------------------
  136. ALIGN(16)
  137. L(xsize_more_than_two):
  138. C The first limb of yp is processed with a simple mpn_mul_1 style loop
  139. C inline.  Unrolling this doesn't seem worthwhile since it's only run once
  140. C (whereas the addmul below is run ysize-1 many times).  A call to the
  141. C actual mpn_mul_1 will be slowed down by the call and parameter pushing and
  142. C popping, and doesn't seem likely to be worthwhile on the typical 13-26
  143. C limb operations the Karatsuba code calls here with.
  144. C eax yp[0]
  145. C ebx
  146. C ecx xsize
  147. C edx xp
  148. C esi
  149. C edi
  150. C ebp
  151. dnl  FRAME doesn't carry on from previous, no pushes yet here
  152. defframe(`SAVE_EBX',-4)
  153. defframe(`SAVE_ESI',-8)
  154. defframe(`SAVE_EDI',-12)
  155. defframe(`SAVE_EBP',-16)
  156. deflit(`FRAME',0)
  157. subl $16, %esp
  158. deflit(`FRAME',16)
  159. movl %edi, SAVE_EDI
  160. movl PARAM_WP, %edi
  161. movl %ebx, SAVE_EBX
  162. movl %ebp, SAVE_EBP
  163. movl %eax, %ebp
  164. movl %esi, SAVE_ESI
  165. xorl %ebx, %ebx
  166. leal (%edx,%ecx,4), %esi C xp end
  167. leal (%edi,%ecx,4), %edi C wp end of mul1
  168. negl %ecx
  169. L(mul1):
  170. C eax scratch
  171. C ebx carry
  172. C ecx counter, negative
  173. C edx scratch
  174. C esi xp end
  175. C edi wp end of mul1
  176. C ebp multiplier
  177. movl (%esi,%ecx,4), %eax
  178. mull %ebp
  179. addl %ebx, %eax
  180. movl %eax, (%edi,%ecx,4)
  181. movl $0, %ebx
  182. adcl %edx, %ebx
  183. incl %ecx
  184. jnz L(mul1)
  185. movl PARAM_YSIZE, %edx
  186. movl PARAM_XSIZE, %ecx
  187. movl %ebx, (%edi) C final carry
  188. decl %edx
  189. jnz L(ysize_more_than_one)
  190. movl SAVE_EDI, %edi
  191. movl SAVE_EBX, %ebx
  192. movl SAVE_EBP, %ebp
  193. movl SAVE_ESI, %esi
  194. addl $FRAME, %esp
  195. ret
  196. L(ysize_more_than_one):
  197. cmpl $UNROLL_THRESHOLD, %ecx
  198. movl PARAM_YP, %eax
  199. jae L(unroll)
  200. C -----------------------------------------------------------------------------
  201. C simple addmul looping
  202. C
  203. C eax yp
  204. C ebx
  205. C ecx xsize
  206. C edx ysize-1
  207. C esi xp end
  208. C edi wp end of mul1
  209. C ebp
  210. leal 4(%eax,%edx,4), %ebp C yp end
  211. negl %ecx
  212. negl %edx
  213. movl (%esi,%ecx,4), %eax C xp low limb
  214. movl %edx, PARAM_YSIZE C -(ysize-1)
  215. incl %ecx
  216. xorl %ebx, %ebx C initial carry
  217. movl %ecx, PARAM_XSIZE C -(xsize-1)
  218. movl %ebp, PARAM_YP
  219. movl (%ebp,%edx,4), %ebp C yp second lowest limb - multiplier
  220. jmp L(simple_outer_entry)
  221. C this is offset 0x121 so close enough to aligned
  222. L(simple_outer_top):
  223. C ebp ysize counter, negative
  224. movl PARAM_YP, %edx
  225. movl PARAM_XSIZE, %ecx C -(xsize-1)
  226. xorl %ebx, %ebx C carry
  227. movl %ebp, PARAM_YSIZE
  228. addl $4, %edi C next position in wp
  229. movl (%edx,%ebp,4), %ebp C yp limb - multiplier
  230. movl -4(%esi,%ecx,4), %eax C xp low limb
  231. L(simple_outer_entry):
  232. L(simple_inner):
  233. C eax xp limb
  234. C ebx carry limb
  235. C ecx loop counter (negative)
  236. C edx scratch
  237. C esi xp end
  238. C edi wp end
  239. C ebp multiplier
  240. mull %ebp
  241. addl %eax, %ebx
  242. adcl $0, %edx
  243. addl %ebx, (%edi,%ecx,4)
  244. movl (%esi,%ecx,4), %eax
  245. adcl $0, %edx
  246. incl %ecx
  247. movl %edx, %ebx
  248. jnz L(simple_inner)
  249. mull %ebp
  250. movl PARAM_YSIZE, %ebp
  251. addl %eax, %ebx
  252. adcl $0, %edx
  253. addl %ebx, (%edi)
  254. adcl $0, %edx
  255. incl %ebp
  256. movl %edx, 4(%edi)
  257. jnz L(simple_outer_top)
  258. movl SAVE_EBX, %ebx
  259. movl SAVE_ESI, %esi
  260. movl SAVE_EDI, %edi
  261. movl SAVE_EBP, %ebp
  262. addl $FRAME, %esp
  263. ret
  264. C -----------------------------------------------------------------------------
  265. C
  266. C The unrolled loop is the same as in mpn_addmul_1(), see that code for some
  267. C comments.
  268. C
  269. C VAR_ADJUST is the negative of how many limbs the leals in the inner loop
  270. C increment xp and wp.  This is used to adjust back xp and wp, and rshifted
  271. C to given an initial VAR_COUNTER at the top of the outer loop.
  272. C
  273. C VAR_COUNTER is for the unrolled loop, running from VAR_ADJUST/UNROLL_COUNT
  274. C up to -1, inclusive.
  275. C
  276. C VAR_JMP is the computed jump into the unrolled loop.
  277. C
  278. C VAR_XP_LOW is the least significant limb of xp, which is needed at the
  279. C start of the unrolled loop.
  280. C
  281. C PARAM_YSIZE is the outer loop counter, going from -(ysize-1) up to -1,
  282. C inclusive.
  283. C
  284. C PARAM_YP is offset appropriately so that the PARAM_YSIZE counter can be
  285. C added to give the location of the next limb of yp, which is the multiplier
  286. C in the unrolled loop.
  287. C
  288. C The trick with VAR_ADJUST means it's only necessary to do one fetch in the
  289. C outer loop to take care of xp, wp and the inner loop counter.
  290. defframe(VAR_COUNTER,  -20)
  291. defframe(VAR_ADJUST,   -24)
  292. defframe(VAR_JMP,      -28)
  293. defframe(VAR_XP_LOW,   -32)
  294. deflit(VAR_EXTRA_SPACE, 16)
  295. L(unroll):
  296. C eax yp
  297. C ebx
  298. C ecx xsize
  299. C edx ysize-1
  300. C esi xp end
  301. C edi wp end of mul1
  302. C ebp
  303. movl PARAM_XP, %esi
  304. movl 4(%eax), %ebp C multiplier (yp second limb)
  305. leal 4(%eax,%edx,4), %eax C yp adjust for ysize indexing
  306. movl PARAM_WP, %edi
  307. movl %eax, PARAM_YP
  308. negl %edx
  309. movl %edx, PARAM_YSIZE
  310. leal UNROLL_COUNT-2(%ecx), %ebx C (xsize-1)+UNROLL_COUNT-1
  311. decl %ecx C xsize-1
  312. movl (%esi), %eax C xp low limb
  313. andl $-UNROLL_MASK-1, %ebx
  314. negl %ecx
  315. subl $VAR_EXTRA_SPACE, %esp
  316. deflit(`FRAME',16+VAR_EXTRA_SPACE)
  317. negl %ebx
  318. andl $UNROLL_MASK, %ecx
  319. movl %ebx, VAR_ADJUST
  320. movl %ecx, %edx
  321. shll $4, %ecx
  322. sarl $UNROLL_LOG2, %ebx
  323. C 17 code bytes per limb
  324. ifdef(`PIC',`
  325. call L(pic_calc)
  326. L(unroll_here):
  327. ',`
  328. leal L(unroll_entry) (%ecx,%edx,1), %ecx
  329. ')
  330. negl %edx
  331. movl %eax, VAR_XP_LOW
  332. movl %ecx, VAR_JMP
  333. leal 4(%edi,%edx,4), %edi C wp and xp, adjust for unrolling,
  334. leal 4(%esi,%edx,4), %esi C  and start at second limb
  335. jmp L(unroll_outer_entry)
  336. ifdef(`PIC',`
  337. L(pic_calc):
  338. C See mpn/x86/README about old gas bugs
  339. leal (%ecx,%edx,1), %ecx
  340. addl $L(unroll_entry)-L(unroll_here), %ecx
  341. addl (%esp), %ecx
  342. ret_internal
  343. ')
  344. C --------------------------------------------------------------------------
  345. ALIGN(32)
  346. L(unroll_outer_top):
  347. C ebp ysize counter, negative
  348. movl VAR_ADJUST, %ebx
  349. movl PARAM_YP, %edx
  350. movl VAR_XP_LOW, %eax
  351. movl %ebp, PARAM_YSIZE C store incremented ysize counter
  352. leal 4(%edi,%ebx,4), %edi
  353. leal (%esi,%ebx,4), %esi
  354. sarl $UNROLL_LOG2, %ebx
  355. movl (%edx,%ebp,4), %ebp C yp next multiplier
  356. movl VAR_JMP, %ecx
  357. L(unroll_outer_entry):
  358. mull %ebp
  359. testb $1, %cl C and clear carry bit
  360. movl %ebx, VAR_COUNTER
  361. movl $0, %ebx
  362. movl $0, %ecx
  363. cmovz( %eax, %ecx) C eax into low carry, zero into high carry limb
  364. cmovnz( %eax, %ebx)
  365. C Extra fetch of VAR_JMP is bad, but registers are tight
  366. jmp *VAR_JMP
  367. C -----------------------------------------------------------------------------
  368. ALIGN(32)
  369. L(unroll_top):
  370. C eax xp limb
  371. C ebx carry high
  372. C ecx carry low
  373. C edx scratch
  374. C esi xp+8
  375. C edi wp
  376. C ebp yp multiplier limb
  377. C
  378. C VAR_COUNTER  loop counter, negative
  379. C
  380. C 17 bytes each limb
  381. L(unroll_entry):
  382. deflit(CHUNK_COUNT,2)
  383. forloop(`i', 0, UNROLL_COUNT/CHUNK_COUNT-1, `
  384. deflit(`disp0', eval(i*CHUNK_COUNT*4 ifelse(UNROLL_BYTES,256,-128)))
  385. deflit(`disp1', eval(disp0 + 4))
  386. Zdisp( movl, disp0,(%esi), %eax)
  387. adcl %edx, %ebx
  388. mull %ebp
  389. Zdisp( addl, %ecx, disp0,(%edi))
  390. movl $0, %ecx
  391. adcl %eax, %ebx
  392. movl disp1(%esi), %eax
  393. adcl %edx, %ecx
  394. mull %ebp
  395. addl %ebx, disp1(%edi)
  396. movl $0, %ebx
  397. adcl %eax, %ecx
  398. ')
  399. incl VAR_COUNTER
  400. leal UNROLL_BYTES(%esi), %esi
  401. leal UNROLL_BYTES(%edi), %edi
  402. jnz L(unroll_top)
  403. C eax
  404. C ebx zero
  405. C ecx low
  406. C edx high
  407. C esi
  408. C edi wp, pointing at second last limb)
  409. C ebp
  410. C
  411. C carry flag to be added to high
  412. deflit(`disp0', ifelse(UNROLL_BYTES,256,-128))
  413. deflit(`disp1', eval(disp0-0 + 4))
  414. movl PARAM_YSIZE, %ebp
  415. adcl $0, %edx
  416. addl %ecx, disp0(%edi)
  417. adcl $0, %edx
  418. incl %ebp
  419. movl %edx, disp1(%edi)
  420. jnz L(unroll_outer_top)
  421. movl SAVE_ESI, %esi
  422. movl SAVE_EBP, %ebp
  423. movl SAVE_EDI, %edi
  424. movl SAVE_EBX, %ebx
  425. addl $FRAME, %esp
  426. ret
  427. EPILOGUE()