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

数学计算

开发平台:

Unix_Linux

  1. dnl  Itanium-2 mpn_gcd_1 -- mpn by 1 gcd.
  2. dnl  Copyright 2002, 2003, 2004, 2005 Free Software Foundation, Inc.
  3. dnl  This file is part of the GNU MP Library.
  4. dnl  The GNU MP Library is free software; you can redistribute it and/or modify
  5. dnl  it under the terms of the GNU Lesser General Public License as published
  6. dnl  by the Free Software Foundation; either version 3 of the License, or (at
  7. dnl  your option) any later version.
  8. dnl  The GNU MP Library is distributed in the hope that it will be useful, but
  9. dnl  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  10. dnl  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  11. dnl  License for more details.
  12. dnl  You should have received a copy of the GNU Lesser General Public License
  13. dnl  along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.
  14. include(`../config.m4')
  15. C           cycles/bitpair (1x1 gcd)
  16. C Itanium:      14 (approx)
  17. C Itanium 2:     6.3
  18. C mpn_gcd_1 (mp_srcptr xp, mp_size_t xsize, mp_limb_t y);
  19. C
  20. C The entry sequence is designed to expect xsize>1 and hence a modexact
  21. C call.  This ought to be more common than a 1x1 operation.  Our critical
  22. C path is thus stripping factors of 2 from y, calling modexact, then
  23. C stripping factors of 2 from the x remainder returned.
  24. C
  25. C The common factors of 2 between x and y must be determined using the
  26. C original x, not the remainder from the modexact.  This is done with
  27. C x_orig which is xp[0].  There's plenty of time to do this while the rest
  28. C of the modexact etc is happening.
  29. C
  30. C It's possible xp[0] is zero.  In this case the trailing zeros calculation
  31. C popc((x-1)&~x) gives 63, and that's clearly no less than what y will
  32. C have, making min(x_twos,y_twos) == y_twos.
  33. C
  34. C The main loop consists of transforming x,y to abs(x-y),min(x,y), and then
  35. C stripping factors of 2 from abs(x-y).  Those factors of two are
  36. C determined from just y-x, without the abs(), since there's the same
  37. C number of trailing zeros on n or -n in twos complement.  That makes the
  38. C dependent chain
  39. C
  40. C cycles
  41. C   1    sub     x-y and x-y-1
  42. C   3    andcm   (x-y-1)&~(x-y)
  43. C   2    popcnt  trailing zeros
  44. C   3    shr.u   strip abs(x-y)
  45. C  ---
  46. C   9
  47. C
  48. C The selection of x-y versus y-x for abs(x-y), and the selection of the
  49. C minimum of x and y, is done in parallel with the above.
  50. C
  51. C The algorithm takes about 0.68 iterations per bit (two N bit operands) on
  52. C average, hence the final 6.3 cycles/bitpair.
  53. C
  54. C The loop is not as fast as one might hope, since there's extra latency
  55. C from andcm going across to the `multimedia' popcnt, and vice versa from
  56. C multimedia shr.u back to the integer sub.
  57. C
  58. C The loop branch is .sptk.clr since we usually expect a good number of
  59. C iterations, and the iterations are data dependent so it's unlikely past
  60. C results will predict anything much about the future.
  61. C
  62. C Not done:
  63. C
  64. C An alternate algorithm which didn't strip all twos, but instead applied
  65. C tbit and predicated extr on x, and then y, was attempted.  The loop was 6
  66. C cycles, but the algorithm is an average 1.25 iterations per bitpair for a
  67. C total 7.25 c/bp, which is slower than the current approach.
  68. C
  69. C Alternatives:
  70. C
  71. C Perhaps we could do something tricky by extracting a few high bits and a
  72. C few low bits from the operands, and looking up a table which would give a
  73. C set of predicates to control some shifts or subtracts or whatever.  That
  74. C could knock off multiple bits per iteration.
  75. C
  76. C The right shifts are a bit of a bottleneck (shr at 2 or 3 cycles, or extr
  77. C only going down I0), perhaps it'd be possible to shift left instead,
  78. C using add.  That would mean keeping track of the lowest not-yet-zeroed
  79. C bit, using some sort of mask.
  80. C
  81. C Itanium-1:
  82. C
  83. C This code is not designed for itanium-1 and in fact doesn't run well on
  84. C that chip.  The loop seems to be about 21 cycles, probably because we end
  85. C up with a 10 cycle replay for not forcibly scheduling the shr.u latency.
  86. C Lack of branch hints might introduce a couple of bubbles too.
  87. C
  88. ASM_START()
  89. .explicit C What does this mean?
  90. C HP's assembler requires these declarations for importing mpn_modexact_1c_odd
  91. .global mpn_modexact_1c_odd
  92. .type mpn_modexact_1c_odd,@function
  93. PROLOGUE(mpn_gcd_1)
  94. C r32 xp
  95. C r33 xsize
  96. C r34 y
  97. define(x,           r8)
  98. define(xp_orig,     r32)
  99. define(xsize,       r33)
  100. define(y,           r34)  define(inputs, 3)
  101. define(save_rp,     r35)
  102. define(save_pfs,    r36)
  103. define(x_orig,      r37)
  104. define(x_orig_one,  r38)
  105. define(y_twos,      r39)  define(locals, 5)
  106. define(out_xp,      r40)
  107. define(out_xsize,   r41)
  108. define(out_divisor, r42)
  109. define(out_carry,   r43)  define(outputs, 4)
  110. .prologue
  111. { .mmi;
  112. ifdef(`HAVE_ABI_32',
  113. ` addp4 r9 = 0, xp_orig   define(xp,r9)', C M0
  114. `   define(xp,xp_orig)')
  115. .save ar.pfs, save_pfs
  116. alloc save_pfs = ar.pfs, inputs, locals, outputs, 0 C M2
  117. .save rp, save_rp
  118. mov save_rp = b0 C I0
  119. }{ .body
  120. add r10 = -1, y C M3  y-1
  121. } ;;
  122. { .mmi; ld8 x = [xp] C M0  x = xp[0] if no modexact
  123. ld8 x_orig = [xp] C M1  orig x for common twos
  124. cmp.ne p6,p0 = 1, xsize C I0
  125. }{ .mmi; andcm y_twos = r10, y C M2  (y-1)&~y
  126. mov out_xp = xp_orig C M3
  127. mov out_xsize = xsize C I1
  128. } ;;
  129. mov out_carry = 0
  130. C
  131. popcnt y_twos = y_twos C I0  y twos
  132. ;;
  133. C
  134. { .mmi; add x_orig_one = -1, x_orig C M0  orig x-1
  135. shr.u out_divisor = y, y_twos C I0  y without twos
  136. }{ shr.u y = y, y_twos C I1  y without twos
  137. (p6) br.call.sptk.many b0 = mpn_modexact_1c_odd  C if xsize>1
  138. } ;;
  139. C modexact can leave x==0
  140. { .mmi; cmp.eq p6,p0 = 0, x C M0  if {xp,xsize} % y == 0
  141. andcm x_orig = x_orig_one, x_orig C M1  orig (x-1)&~x
  142. add r9 = -1, x C I0  x-1
  143. } ;;
  144. { .mmi; andcm r9 = r9, x C M0  (x-1)&~x
  145. mov b0 = save_rp C I0
  146. } ;;
  147. C
  148. popcnt x_orig = x_orig C I0  orig x twos
  149. popcnt r9 = r9 C I0  x twos
  150. ;;
  151. C
  152. { cmp.lt p7,p0 = x_orig, y_twos C M0  orig x_twos < y_twos
  153. shr.u x = x, r9 C I0  x odd
  154. } ;;
  155. { (p7) mov y_twos = x_orig C M0  common twos
  156. add r10 = -1, y C I0  y-1
  157. (p6) br.dpnt.few .Ldone_y C B0  x%y==0 then result y
  158. } ;;
  159. C
  160. C No noticable difference in speed for the loop aligned to
  161. C 32 or just 16.
  162. .Ltop:
  163. C r8 x
  164. C r10  y-1
  165. C r34 y
  166. C r38 common twos, for use at end
  167. {  .mmi; cmp.gtu p8,p9 = x, y C M0  x>y
  168. cmp.ne p10,p0 = x, y C M1  x==y
  169. sub r9 = y, x C I0  d = y - x
  170. }{ .mmi; sub r10 = r10, x C M2  d-1 = y - x - 1
  171. } ;;
  172. { .mmi; .pred.rel "mutex", p8, p9
  173. (p8) sub x = x, y C M0  x>y  use x=x-y, y unchanged
  174. (p9) mov y = x C M1  y>=x use y=x
  175. (p9) mov x = r9 C I0  y>=x use x=y-x
  176. }{ .mmi; andcm r9 = r10, r9 C M2  (d-1)&~d
  177. ;;
  178. add r10 = -1, y C M0  new y-1
  179. popcnt r9 = r9 C I0  twos on x-y
  180. } ;;
  181. { shr.u x = x, r9 C I0   new x without twos
  182. (p10) br.sptk.few.clr .Ltop
  183. } ;;
  184. C result is y
  185. .Ldone_y:
  186. shl r8 = y, y_twos C I   common factors of 2
  187. ;;
  188. mov ar.pfs = save_pfs C I0
  189. br.ret.sptk.many b0
  190. EPILOGUE()