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

数学计算

开发平台:

Unix_Linux

  1. dnl  Intel Pentium MMX mpn_mul_1 -- mpn by limb multiplication.
  2. dnl  Copyright 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    cycles/limb
  20. C P5:   12.0   for 32-bit multiplier
  21. C        7.0   for 16-bit multiplier
  22. C mp_limb_t mpn_mul_1 (mp_ptr dst, mp_srcptr src, mp_size_t size,
  23. C                      mp_limb_t multiplier);
  24. C
  25. C When the multiplier is 16 bits some special case MMX code is used.  Small
  26. C multipliers might arise reasonably often from mpz_mul_ui etc.  If the size
  27. C is odd there's roughly a 5 cycle penalty, so times for say size==7 and
  28. C size==8 end up being quite close.  If src isn't aligned to an 8 byte
  29. C boundary then one limb is processed separately with roughly a 5 cycle
  30. C penalty, so in that case it's say size==8 and size==9 which are close.
  31. C
  32. C Alternatives:
  33. C
  34. C MMX is not believed to be of any use for 32-bit multipliers, since for
  35. C instance the current method would just have to be more or less duplicated
  36. C for the high and low halves of the multiplier, and would probably
  37. C therefore run at about 14 cycles, which is slower than the plain integer
  38. C at 12.
  39. C
  40. C Adding the high and low MMX products using integer code seems best.  An
  41. C attempt at using paddd and carry bit propagation with pcmpgtd didn't give
  42. C any joy.  Perhaps something could be done keeping the values signed and
  43. C thereby avoiding adjustments to make pcmpgtd into an unsigned compare, or
  44. C perhaps not.
  45. C
  46. C Future:
  47. C
  48. C An mpn_mul_1c entrypoint would need a double carry out of the low result
  49. C limb in the 16-bit code, unless it could be assumed the carry fits in 16
  50. C bits, possibly as carry<multiplier, this being true of a big calculation
  51. C done piece by piece.  But let's worry about that if/when mul_1c is
  52. C actually used.
  53. defframe(PARAM_MULTIPLIER,16)
  54. defframe(PARAM_SIZE,      12)
  55. defframe(PARAM_SRC,       8)
  56. defframe(PARAM_DST,       4)
  57. TEXT
  58. ALIGN(8)
  59. PROLOGUE(mpn_mul_1)
  60. deflit(`FRAME',0)
  61. movl PARAM_SIZE, %ecx
  62. movl PARAM_SRC, %edx
  63. cmpl $1, %ecx
  64. jne L(two_or_more)
  65. C one limb only
  66. movl PARAM_MULTIPLIER, %eax
  67. movl PARAM_DST, %ecx
  68. mull (%edx)
  69. movl %eax, (%ecx)
  70. movl %edx, %eax
  71. ret
  72. L(two_or_more):
  73. C eax size
  74. C ebx
  75. C ecx carry
  76. C edx
  77. C esi src
  78. C edi
  79. C ebp
  80. pushl %esi FRAME_pushl()
  81. pushl %edi FRAME_pushl()
  82. movl %edx, %esi C src
  83. movl PARAM_DST, %edi
  84. movl PARAM_MULTIPLIER, %eax
  85. pushl %ebx FRAME_pushl()
  86. leal (%esi,%ecx,4), %esi C src end
  87. leal (%edi,%ecx,4), %edi C dst end
  88. negl %ecx C -size
  89. pushl %ebp FRAME_pushl()
  90. cmpl $65536, %eax
  91. jb L(small)
  92. L(big):
  93. xorl %ebx, %ebx C carry limb
  94. sarl %ecx C -size/2
  95. jnc L(top) C with carry flag clear
  96. C size was odd, process one limb separately
  97. mull 4(%esi,%ecx,8) C m * src[0]
  98. movl %eax, 4(%edi,%ecx,8)
  99. incl %ecx
  100. orl %edx, %ebx C carry limb, and clear carry flag
  101. L(top):
  102. C eax
  103. C ebx carry
  104. C ecx counter, negative
  105. C edx
  106. C esi src end
  107. C edi dst end
  108. C ebp (scratch carry)
  109. adcl $0, %ebx
  110. movl (%esi,%ecx,8), %eax
  111. mull PARAM_MULTIPLIER
  112. movl %edx, %ebp
  113. addl %eax, %ebx
  114. adcl $0, %ebp
  115. movl 4(%esi,%ecx,8), %eax
  116. mull PARAM_MULTIPLIER
  117. movl %ebx, (%edi,%ecx,8)
  118. addl %ebp, %eax
  119. movl %eax, 4(%edi,%ecx,8)
  120. incl %ecx
  121. movl %edx, %ebx
  122. jnz L(top)
  123. adcl $0, %ebx
  124. popl %ebp
  125. movl %ebx, %eax
  126. popl %ebx
  127. popl %edi
  128. popl %esi
  129. ret
  130. L(small):
  131. C Special case for 16-bit multiplier.
  132. C
  133. C eax multiplier
  134. C ebx
  135. C ecx -size
  136. C edx src
  137. C esi src end
  138. C edi dst end
  139. C ebp multiplier
  140. C size<3 not supported here.  At size==3 we're already a couple of
  141. C cycles faster, so there's no threshold as such, just use the MMX
  142. C as soon as possible.
  143. cmpl $-3, %ecx
  144. ja L(big)
  145. movd %eax, %mm7 C m
  146. pxor %mm6, %mm6 C initial carry word
  147. punpcklwd %mm7, %mm7 C m replicated 2 times
  148. addl $2, %ecx C -size+2
  149. punpckldq %mm7, %mm7 C m replicated 4 times
  150. andl $4, %edx C test alignment, clear carry flag
  151. movq %mm7, %mm0 C m
  152. jz L(small_entry)
  153. C Source is unaligned, process one limb separately.
  154. C
  155. C Plain integer code is used here, since it's smaller and is about
  156. C the same 13 cycles as an mmx block would be.
  157. C
  158. C An "addl $1,%ecx" doesn't clear the carry flag when size==3, hence
  159. C the use of separate incl and orl.
  160. mull -8(%esi,%ecx,4) C m * src[0]
  161. movl %eax, -8(%edi,%ecx,4) C dst[0]
  162. incl %ecx C one limb processed
  163. movd %edx, %mm6 C initial carry
  164. orl %eax, %eax C clear carry flag
  165. jmp L(small_entry)
  166. C The scheduling here is quite tricky, since so many instructions have
  167. C pairing restrictions.  In particular the js won't pair with a movd, and
  168. C can't be paired with an adc since it wants flags from the inc, so
  169. C instructions are rotated to the top of the loop to find somewhere useful
  170. C for it.
  171. C
  172. C Trouble has been taken to avoid overlapping successive loop iterations,
  173. C since that would greatly increase the size of the startup and finishup
  174. C code.  Actually there's probably not much advantage to be had from
  175. C overlapping anyway, since the difficulties are mostly with pairing, not
  176. C with latencies as such.
  177. C
  178. C In the comments x represents the src data and m the multiplier (16
  179. C bits, but replicated 4 times).
  180. C
  181. C The m signs calculated in %mm3 are a loop invariant and could be held in
  182. C say %mm5, but that would save only one instruction and hence be no faster.
  183. L(small_top):
  184. C eax l.low, then l.high
  185. C ebx (h.low)
  186. C ecx counter, -size+2 to 0 or 1
  187. C edx (h.high)
  188. C esi &src[size]
  189. C edi &dst[size]
  190. C ebp
  191. C
  192. C %mm0 (high products)
  193. C %mm1 (low products)
  194. C %mm2 (adjust for m using x signs)
  195. C %mm3 (adjust for x using m signs)
  196. C %mm4
  197. C %mm5
  198. C %mm6 h.low, then carry
  199. C %mm7 m replicated 4 times
  200. movd %mm6, %ebx C h.low
  201. psrlq $32, %mm1 C l.high
  202. movd %mm0, %edx C h.high
  203. movq %mm0, %mm6 C new c
  204. adcl %eax, %ebx
  205. incl %ecx
  206. movd %mm1, %eax C l.high
  207. movq %mm7, %mm0
  208. adcl %eax, %edx
  209. movl %ebx, -16(%edi,%ecx,4)
  210. movl %edx, -12(%edi,%ecx,4)
  211. psrlq $32, %mm6 C c
  212. L(small_entry):
  213. pmulhw -8(%esi,%ecx,4), %mm0 C h = (x*m).high
  214. movq %mm7, %mm1
  215. pmullw -8(%esi,%ecx,4), %mm1 C l = (x*m).low
  216. movq %mm7, %mm3
  217. movq -8(%esi,%ecx,4), %mm2 C x
  218. psraw $15, %mm3 C m signs
  219. pand -8(%esi,%ecx,4), %mm3 C x selected by m signs
  220. psraw $15, %mm2 C x signs
  221. paddw %mm3, %mm0 C add x to h if m neg
  222. pand %mm7, %mm2 C m selected by x signs
  223. paddw %mm2, %mm0 C add m to h if x neg
  224. incl %ecx
  225. movd %mm1, %eax C l.low
  226. punpcklwd %mm0, %mm6 C c + h.low << 16
  227. psrlq $16, %mm0 C h.high
  228. js L(small_top)
  229. movd %mm6, %ebx C h.low
  230. psrlq $32, %mm1 C l.high
  231. adcl %eax, %ebx
  232. popl %ebp FRAME_popl()
  233. movd %mm0, %edx C h.high
  234. psrlq $32, %mm0 C l.high
  235. movd %mm1, %eax C l.high
  236. adcl %eax, %edx
  237. movl %ebx, -12(%edi,%ecx,4)
  238. movd %mm0, %eax C c
  239. adcl $0, %eax
  240. movl %edx, -8(%edi,%ecx,4)
  241. orl %ecx, %ecx
  242. jnz L(small_done) C final %ecx==1 means even, ==0 odd
  243. C Size odd, one extra limb to process.
  244. C Plain integer code is used here, since it's smaller and is about
  245. C the same speed as another mmx block would be.
  246. movl %eax, %ecx
  247. movl PARAM_MULTIPLIER, %eax
  248. mull -4(%esi)
  249. addl %ecx, %eax
  250. adcl $0, %edx
  251. movl %eax, -4(%edi)
  252. movl %edx, %eax
  253. L(small_done):
  254. popl %ebx
  255. popl %edi
  256. popl %esi
  257. emms
  258. ret
  259. EPILOGUE()