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

数学计算

开发平台:

Unix_Linux

  1. dnl  Intel P6 mpn_sqr_basecase -- square an mpn number.
  2. dnl  Copyright 1999, 2000, 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 P6: approx 4.0 cycles per cross product, or 7.75 cycles per triangular
  20. C     product (measured on the speed difference between 20 and 40 limbs,
  21. C     which is the Karatsuba recursing range).
  22. dnl  These are the same as in mpn/x86/k6/sqr_basecase.asm, see that file for
  23. dnl  a description.  The only difference here is that UNROLL_COUNT can go up
  24. dnl  to 64 (not 63) making SQR_TOOM2_THRESHOLD_MAX 67.
  25. deflit(SQR_TOOM2_THRESHOLD_MAX, 67)
  26. ifdef(`SQR_TOOM2_THRESHOLD_OVERRIDE',
  27. `define(`SQR_TOOM2_THRESHOLD',SQR_TOOM2_THRESHOLD_OVERRIDE)')
  28. m4_config_gmp_mparam(`SQR_TOOM2_THRESHOLD')
  29. deflit(UNROLL_COUNT, eval(SQR_TOOM2_THRESHOLD-3))
  30. C void mpn_sqr_basecase (mp_ptr dst, mp_srcptr src, mp_size_t size);
  31. C
  32. C The algorithm is basically the same as mpn/generic/sqr_basecase.c, but a
  33. C lot of function call overheads are avoided, especially when the given size
  34. C is small.
  35. C
  36. C The code size might look a bit excessive, but not all of it is executed so
  37. C it won't all get into the code cache.  The 1x1, 2x2 and 3x3 special cases
  38. C clearly apply only to those sizes; mid sizes like 10x10 only need part of
  39. C the unrolled addmul; and big sizes like 40x40 that do use the full
  40. C unrolling will least be making good use of it, because 40x40 will take
  41. C something like 7000 cycles.
  42. defframe(PARAM_SIZE,12)
  43. defframe(PARAM_SRC, 8)
  44. defframe(PARAM_DST, 4)
  45. TEXT
  46. ALIGN(32)
  47. PROLOGUE(mpn_sqr_basecase)
  48. deflit(`FRAME',0)
  49. movl PARAM_SIZE, %edx
  50. movl PARAM_SRC, %eax
  51. cmpl $2, %edx
  52. movl PARAM_DST, %ecx
  53. je L(two_limbs)
  54. movl (%eax), %eax
  55. ja L(three_or_more)
  56. C -----------------------------------------------------------------------------
  57. C one limb only
  58. C eax src limb
  59. C ebx
  60. C ecx dst
  61. C edx
  62. mull %eax
  63. movl %eax, (%ecx)
  64. movl %edx, 4(%ecx)
  65. ret
  66. C -----------------------------------------------------------------------------
  67. L(two_limbs):
  68. C eax src
  69. C ebx
  70. C ecx dst
  71. C edx
  72. defframe(SAVE_ESI, -4)
  73. defframe(SAVE_EBX, -8)
  74. defframe(SAVE_EDI, -12)
  75. defframe(SAVE_EBP, -16)
  76. deflit(`STACK_SPACE',16)
  77. subl $STACK_SPACE, %esp
  78. deflit(`FRAME',STACK_SPACE)
  79. movl %esi, SAVE_ESI
  80. movl %eax, %esi
  81. movl (%eax), %eax
  82. mull %eax C src[0]^2
  83. movl %eax, (%ecx) C dst[0]
  84. movl 4(%esi), %eax
  85. movl %ebx, SAVE_EBX
  86. movl %edx, %ebx C dst[1]
  87. mull %eax C src[1]^2
  88. movl %edi, SAVE_EDI
  89. movl %eax, %edi C dst[2]
  90. movl (%esi), %eax
  91. movl %ebp, SAVE_EBP
  92. movl %edx, %ebp C dst[3]
  93. mull 4(%esi) C src[0]*src[1]
  94. addl %eax, %ebx
  95. movl SAVE_ESI, %esi
  96. adcl %edx, %edi
  97. adcl $0, %ebp
  98. addl %ebx, %eax
  99. movl SAVE_EBX, %ebx
  100. adcl %edi, %edx
  101. movl SAVE_EDI, %edi
  102. adcl $0, %ebp
  103. movl %eax, 4(%ecx)
  104. movl %ebp, 12(%ecx)
  105. movl SAVE_EBP, %ebp
  106. movl %edx, 8(%ecx)
  107. addl $FRAME, %esp
  108. ret
  109. C -----------------------------------------------------------------------------
  110. L(three_or_more):
  111. C eax src low limb
  112. C ebx
  113. C ecx dst
  114. C edx size
  115. deflit(`FRAME',0)
  116. pushl %esi defframe_pushl(`SAVE_ESI')
  117. cmpl $4, %edx
  118. movl PARAM_SRC, %esi
  119. jae L(four_or_more)
  120. C -----------------------------------------------------------------------------
  121. C three limbs
  122. C eax src low limb
  123. C ebx
  124. C ecx dst
  125. C edx
  126. C esi src
  127. C edi
  128. C ebp
  129. pushl %ebp defframe_pushl(`SAVE_EBP')
  130. pushl %edi defframe_pushl(`SAVE_EDI')
  131. mull %eax C src[0] ^ 2
  132. movl %eax, (%ecx)
  133. movl %edx, 4(%ecx)
  134. movl 4(%esi), %eax
  135. xorl %ebp, %ebp
  136. mull %eax C src[1] ^ 2
  137. movl %eax, 8(%ecx)
  138. movl %edx, 12(%ecx)
  139. movl 8(%esi), %eax
  140. pushl %ebx defframe_pushl(`SAVE_EBX')
  141. mull %eax C src[2] ^ 2
  142. movl %eax, 16(%ecx)
  143. movl %edx, 20(%ecx)
  144. movl (%esi), %eax
  145. mull 4(%esi) C src[0] * src[1]
  146. movl %eax, %ebx
  147. movl %edx, %edi
  148. movl (%esi), %eax
  149. mull 8(%esi) C src[0] * src[2]
  150. addl %eax, %edi
  151. movl %edx, %ebp
  152. adcl $0, %ebp
  153. movl 4(%esi), %eax
  154. mull 8(%esi) C src[1] * src[2]
  155. xorl %esi, %esi
  156. addl %eax, %ebp
  157. C eax
  158. C ebx dst[1]
  159. C ecx dst
  160. C edx dst[4]
  161. C esi zero, will be dst[5]
  162. C edi dst[2]
  163. C ebp dst[3]
  164. adcl $0, %edx
  165. addl %ebx, %ebx
  166. adcl %edi, %edi
  167. adcl %ebp, %ebp
  168. adcl %edx, %edx
  169. movl 4(%ecx), %eax
  170. adcl $0, %esi
  171. addl %ebx, %eax
  172. movl %eax, 4(%ecx)
  173. movl 8(%ecx), %eax
  174. adcl %edi, %eax
  175. movl 12(%ecx), %ebx
  176. adcl %ebp, %ebx
  177. movl 16(%ecx), %edi
  178. movl %eax, 8(%ecx)
  179. movl SAVE_EBP, %ebp
  180. movl %ebx, 12(%ecx)
  181. movl SAVE_EBX, %ebx
  182. adcl %edx, %edi
  183. movl 20(%ecx), %eax
  184. movl %edi, 16(%ecx)
  185. movl SAVE_EDI, %edi
  186. adcl %esi, %eax C no carry out of this
  187. movl SAVE_ESI, %esi
  188. movl %eax, 20(%ecx)
  189. addl $FRAME, %esp
  190. ret
  191. C -----------------------------------------------------------------------------
  192. defframe(VAR_COUNTER,-20)
  193. defframe(VAR_JMP,    -24)
  194. deflit(`STACK_SPACE',24)
  195. L(four_or_more):
  196. C eax src low limb
  197. C ebx
  198. C ecx
  199. C edx size
  200. C esi src
  201. C edi
  202. C ebp
  203. deflit(`FRAME',4)  dnl  %esi already pushed
  204. C First multiply src[0]*src[1..size-1] and store at dst[1..size].
  205. subl $STACK_SPACE-FRAME, %esp
  206. deflit(`FRAME',STACK_SPACE)
  207. movl $1, %ecx
  208. movl %edi, SAVE_EDI
  209. movl PARAM_DST, %edi
  210. movl %ebx, SAVE_EBX
  211. subl %edx, %ecx C -(size-1)
  212. movl %ebp, SAVE_EBP
  213. movl $0, %ebx C initial carry
  214. leal (%esi,%edx,4), %esi C &src[size]
  215. movl %eax, %ebp C multiplier
  216. leal -4(%edi,%edx,4), %edi C &dst[size-1]
  217. C This loop runs at just over 6 c/l.
  218. L(mul_1):
  219. C eax scratch
  220. C ebx carry
  221. C ecx counter, limbs, negative, -(size-1) to -1
  222. C edx scratch
  223. C esi &src[size]
  224. C edi &dst[size-1]
  225. C ebp multiplier
  226. movl %ebp, %eax
  227. mull (%esi,%ecx,4)
  228. addl %ebx, %eax
  229. movl $0, %ebx
  230. adcl %edx, %ebx
  231. movl %eax, 4(%edi,%ecx,4)
  232. incl %ecx
  233. jnz L(mul_1)
  234. movl %ebx, 4(%edi)
  235. C Addmul src[n]*src[n+1..size-1] at dst[2*n-1...], for each n=1..size-2.
  236. C
  237. C The last two addmuls, which are the bottom right corner of the product
  238. C triangle, are left to the end.  These are src[size-3]*src[size-2,size-1]
  239. C and src[size-2]*src[size-1].  If size is 4 then it's only these corner
  240. C cases that need to be done.
  241. C
  242. C The unrolled code is the same as mpn_addmul_1(), see that routine for some
  243. C comments.
  244. C
  245. C VAR_COUNTER is the outer loop, running from -(size-4) to -1, inclusive.
  246. C
  247. C VAR_JMP is the computed jump into the unrolled code, stepped by one code
  248. C chunk each outer loop.
  249. dnl  This is also hard-coded in the address calculation below.
  250. deflit(CODE_BYTES_PER_LIMB, 15)
  251. dnl  With &src[size] and &dst[size-1] pointers, the displacements in the
  252. dnl  unrolled code fit in a byte for UNROLL_COUNT values up to 32, but above
  253. dnl  that an offset must be added to them.
  254. deflit(OFFSET,
  255. ifelse(eval(UNROLL_COUNT>32),1,
  256. eval((UNROLL_COUNT-32)*4),
  257. 0))
  258. C eax
  259. C ebx carry
  260. C ecx
  261. C edx
  262. C esi &src[size]
  263. C edi &dst[size-1]
  264. C ebp
  265. movl PARAM_SIZE, %ecx
  266. subl $4, %ecx
  267. jz L(corner)
  268. movl %ecx, %edx
  269. negl %ecx
  270. shll $4, %ecx
  271. ifelse(OFFSET,0,,`subl $OFFSET, %esi')
  272. ifdef(`PIC',`
  273. call L(pic_calc)
  274. L(here):
  275. ',`
  276. leal L(unroll_inner_end)-eval(2*CODE_BYTES_PER_LIMB)(%ecx,%edx), %ecx
  277. ')
  278. negl %edx
  279. ifelse(OFFSET,0,,`subl $OFFSET, %edi')
  280. C The calculated jump mustn't be before the start of the available
  281. C code.  This is the limit that UNROLL_COUNT puts on the src operand
  282. C size, but checked here using the jump address directly.
  283. ASSERT(ae,
  284. `movl_text_address( L(unroll_inner_start), %eax)
  285. cmpl %eax, %ecx')
  286. C -----------------------------------------------------------------------------
  287. ALIGN(16)
  288. L(unroll_outer_top):
  289. C eax
  290. C ebx high limb to store
  291. C ecx VAR_JMP
  292. C edx VAR_COUNTER, limbs, negative
  293. C esi &src[size], constant
  294. C edi dst ptr, second highest limb of last addmul
  295. C ebp
  296. movl -12+OFFSET(%esi,%edx,4), %ebp C multiplier
  297. movl %edx, VAR_COUNTER
  298. movl -8+OFFSET(%esi,%edx,4), %eax C first limb of multiplicand
  299. mull %ebp
  300. define(cmovX,`ifelse(eval(UNROLL_COUNT%2),1,`cmovz($@)',`cmovnz($@)')')
  301. testb $1, %cl
  302. movl %edx, %ebx C high carry
  303. leal 4(%edi), %edi
  304. movl %ecx, %edx C jump
  305. movl %eax, %ecx C low carry
  306. leal CODE_BYTES_PER_LIMB(%edx), %edx
  307. cmovX( %ebx, %ecx) C high carry reverse
  308. cmovX( %eax, %ebx) C low carry reverse
  309. movl %edx, VAR_JMP
  310. jmp *%edx
  311. C Must be on an even address here so the low bit of the jump address
  312. C will indicate which way around ecx/ebx should start.
  313. ALIGN(2)
  314. L(unroll_inner_start):
  315. C eax scratch
  316. C ebx carry high
  317. C ecx carry low
  318. C edx scratch
  319. C esi src pointer
  320. C edi dst pointer
  321. C ebp multiplier
  322. C
  323. C 15 code bytes each limb
  324. C ecx/ebx reversed on each chunk
  325. forloop(`i', UNROLL_COUNT, 1, `
  326. deflit(`disp_src', eval(-i*4 + OFFSET))
  327. deflit(`disp_dst', eval(disp_src))
  328. m4_assert(`disp_src>=-128 && disp_src<128')
  329. m4_assert(`disp_dst>=-128 && disp_dst<128')
  330. ifelse(eval(i%2),0,`
  331. Zdisp( movl, disp_src,(%esi), %eax)
  332. mull %ebp
  333. Zdisp( addl, %ebx, disp_dst,(%edi))
  334. adcl %eax, %ecx
  335. movl %edx, %ebx
  336. adcl $0, %ebx
  337. ',`
  338. dnl  this one comes out last
  339. Zdisp( movl, disp_src,(%esi), %eax)
  340. mull %ebp
  341. Zdisp( addl, %ecx, disp_dst,(%edi))
  342. adcl %eax, %ebx
  343. movl %edx, %ecx
  344. adcl $0, %ecx
  345. ')
  346. ')
  347. L(unroll_inner_end):
  348. addl %ebx, m4_empty_if_zero(OFFSET)(%edi)
  349. movl VAR_COUNTER, %edx
  350. adcl $0, %ecx
  351. movl %ecx, m4_empty_if_zero(OFFSET+4)(%edi)
  352. movl VAR_JMP, %ecx
  353. incl %edx
  354. jnz L(unroll_outer_top)
  355. ifelse(OFFSET,0,,`
  356. addl $OFFSET, %esi
  357. addl $OFFSET, %edi
  358. ')
  359. C -----------------------------------------------------------------------------
  360. ALIGN(16)
  361. L(corner):
  362. C eax
  363. C ebx
  364. C ecx
  365. C edx
  366. C esi &src[size]
  367. C edi &dst[2*size-5]
  368. C ebp
  369. movl -12(%esi), %eax
  370. mull -8(%esi)
  371. addl %eax, (%edi)
  372. movl -12(%esi), %eax
  373. movl $0, %ebx
  374. adcl %edx, %ebx
  375. mull -4(%esi)
  376. addl %eax, %ebx
  377. movl -8(%esi), %eax
  378. adcl $0, %edx
  379. addl %ebx, 4(%edi)
  380. movl $0, %ebx
  381. adcl %edx, %ebx
  382. mull -4(%esi)
  383. movl PARAM_SIZE, %ecx
  384. addl %ebx, %eax
  385. adcl $0, %edx
  386. movl %eax, 8(%edi)
  387. movl %edx, 12(%edi)
  388. movl PARAM_DST, %edi
  389. C Left shift of dst[1..2*size-2], the bit shifted out becomes dst[2*size-1].
  390. subl $1, %ecx C size-1
  391. xorl %eax, %eax C ready for final adcl, and clear carry
  392. movl %ecx, %edx
  393. movl PARAM_SRC, %esi
  394. L(lshift):
  395. C eax
  396. C ebx
  397. C ecx counter, size-1 to 1
  398. C edx size-1 (for later use)
  399. C esi src (for later use)
  400. C edi dst, incrementing
  401. C ebp
  402. rcll 4(%edi)
  403. rcll 8(%edi)
  404. leal 8(%edi), %edi
  405. decl %ecx
  406. jnz L(lshift)
  407. adcl %eax, %eax
  408. movl %eax, 4(%edi) C dst most significant limb
  409. movl (%esi), %eax C src[0]
  410. leal 4(%esi,%edx,4), %esi C &src[size]
  411. subl %edx, %ecx C -(size-1)
  412. C Now add in the squares on the diagonal, src[0]^2, src[1]^2, ...,
  413. C src[size-1]^2.  dst[0] hasn't yet been set at all yet, and just gets the
  414. C low limb of src[0]^2.
  415. mull %eax
  416. movl %eax, (%edi,%ecx,8) C dst[0]
  417. L(diag):
  418. C eax scratch
  419. C ebx scratch
  420. C ecx counter, negative
  421. C edx carry
  422. C esi &src[size]
  423. C edi dst[2*size-2]
  424. C ebp
  425. movl (%esi,%ecx,4), %eax
  426. movl %edx, %ebx
  427. mull %eax
  428. addl %ebx, 4(%edi,%ecx,8)
  429. adcl %eax, 8(%edi,%ecx,8)
  430. adcl $0, %edx
  431. incl %ecx
  432. jnz L(diag)
  433. movl SAVE_ESI, %esi
  434. movl SAVE_EBX, %ebx
  435. addl %edx, 4(%edi) C dst most significant limb
  436. movl SAVE_EDI, %edi
  437. movl SAVE_EBP, %ebp
  438. addl $FRAME, %esp
  439. ret
  440. C -----------------------------------------------------------------------------
  441. ifdef(`PIC',`
  442. L(pic_calc):
  443. addl (%esp), %ecx
  444. addl $L(unroll_inner_end)-L(here)-eval(2*CODE_BYTES_PER_LIMB), %ecx
  445. addl %edx, %ecx
  446. ret_internal
  447. ')
  448. EPILOGUE()