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

数学计算

开发平台:

Unix_Linux

  1. dnl  Itanium-2 mpn_modexact_1c_odd -- mpn by 1 exact remainder.
  2. dnl  Copyright 2003, 2004, 2005 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 Itanium:      15
  21. C Itanium 2:     8
  22. dnl  Usage: ABI32(`code')
  23. dnl
  24. dnl  Emit the given code only under HAVE_ABI_32.
  25. dnl
  26. define(ABI32,
  27. m4_assert_onearg()
  28. `ifdef(`HAVE_ABI_32',`$1')')
  29. C mp_limb_t mpn_modexact_1c_odd (mp_srcptr src, mp_size_t size,
  30. C                                mp_limb_t divisor, mp_limb_t carry);
  31. C
  32. C The modexact algorithm is usually conceived as a dependent chain
  33. C
  34. C l = src[i] - c
  35. C q = low(l * inverse)
  36. C c = high(q*divisor) + (src[i]<c)
  37. C
  38. C but we can work the src[i]-c into an xma by calculating si=src[i]*inverse
  39. C separately (off the dependent chain) and using
  40. C
  41. C q = low(c * inverse + si)
  42. C c = high(q*divisor + c)
  43. C
  44. C This means the dependent chain is simply xma.l followed by xma.hu, for a
  45. C total 8 cycles/limb on itanium-2.
  46. C
  47. C The reason xma.hu works for the new c is that the low of q*divisor is
  48. C src[i]-c (being the whole purpose of the q generated, and it can be
  49. C verified algebraically).  If there was an underflow from src[i]-c, then
  50. C there will be an overflow from (src-c)+c, thereby adding 1 to the new c
  51. C the same as the borrow bit (src[i]<c) gives in the first style shown.
  52. C
  53. C Incidentally, fcmp is not an option for treating src[i]-c, since it
  54. C apparently traps to the kernel for unnormalized operands like those used
  55. C and generated by ldf8 and xma.  On one GNU/Linux system it took about 1200
  56. C cycles.
  57. C
  58. C
  59. C First Limb:
  60. C
  61. C The first limb uses q = (src[0]-c) * inverse shown in the first style.
  62. C This lets us get the first q as soon as the inverse is ready, without
  63. C going through si=s*inverse.  Basically at the start we have c and can use
  64. C it while waiting for the inverse, whereas for the second and subsequent
  65. C limbs it's the other way around, ie. we have the inverse and are waiting
  66. C for c.
  67. C
  68. C At .Lentry the first two instructions in the loop have been done already.
  69. C The load of f11=src[1] at the start (predicated on size>=2), and the
  70. C calculation of q by the initial different scheme.
  71. C
  72. C
  73. C Entry Sequence:
  74. C
  75. C In the entry sequence, the critical path is the calculation of the
  76. C inverse, so this is begun first and optimized.  Apart from that, ar.lc is
  77. C established nice and early so the br.cloop's should predict perfectly.
  78. C And the load for the low limbs src[0] and src[1] can be initiated long
  79. C ahead of where they're needed.
  80. C
  81. C
  82. C Inverse Calculation:
  83. C
  84. C The initial 8-bit inverse is calculated using a table lookup.  If it hits
  85. C L1 (which is likely if we're called several times) then it should take a
  86. C total 4 cycles, otherwise hopefully L2 for 9 cycles.  This is considered
  87. C the best approach, on balance.  It could be done bitwise, but that would
  88. C probably be about 14 cycles (2 per bit beyond the first couple).  Or it
  89. C could be taken from 4 bits to 8 with xmpy doubling as used beyond 8 bits,
  90. C but that would be about 11 cycles.
  91. C
  92. C The table is not the same as binvert_limb_table, instead it's 256 bytes,
  93. C designed to be indexed by the low byte of the divisor.  The divisor is
  94. C always odd, so the relevant data is every second byte in the table.  The
  95. C padding lets us use zxt1 instead of extr.u, the latter would cost an extra
  96. C cycle because it must go down I0, and we're using the first I0 slot to get
  97. C ip.  The extra 128 bytes of padding should be insignificant compared to
  98. C typical ia64 code bloat.
  99. C
  100. C Having the table in .text allows us to use IP-relative addressing,
  101. C avoiding a fetch from ltoff.  .rodata is apparently not suitable for use
  102. C IP-relative, it gets a linker relocation overflow on GNU/Linux.
  103. C
  104. C
  105. C Load Scheduling:
  106. C
  107. C In the main loop, the data loads are scheduled for an L2 hit, which means
  108. C 6 cycles for the data ready to use.  In fact we end up 7 cycles ahead.  In
  109. C any case that scheduling is achieved simply by doing the load (and xmpy.l
  110. C for "si") in the immediately preceding iteration.
  111. C
  112. C The main loop requires size >= 2, and we handle size==1 by an initial
  113. C br.cloop to enter the loop only if size>1.  Since ar.lc is established
  114. C early, this should predict perfectly.
  115. C
  116. C
  117. C Not done:
  118. C
  119. C Consideration was given to using a plain "(src[0]-c) % divisor" for
  120. C size==1, but cycle counting suggests about 50 for the sort of approach
  121. C taken by gcc __umodsi3, versus about 47 for the modexact.  (Both assuming
  122. C L1 hits for their respective fetching.)
  123. C
  124. C Consideration was given to a test for high<divisor and replacing the last
  125. C loop iteration with instead c-=src[size-1] followed by c+=d if underflow.
  126. C Branching on high<divisor wouldn't be good since a mispredict would cost
  127. C more than the loop iteration saved, and the condition is of course data
  128. C dependent.  So the theory would be to shorten the loop count if
  129. C high<divisor, and predicate extra operations at the end.  That would mean
  130. C a gain of 6 when high<divisor, or a cost of 2 if not.
  131. C
  132. C Whether such a tradeoff is a win on average depends on assumptions about
  133. C how many bits in the high and the divisor.  If both are uniformly
  134. C distributed then high<divisor about 50% of the time.  But smallish
  135. C divisors (less chance of high<divisor) might be more likely from
  136. C applications (mpz_divisible_ui, mpz_gcd_ui, etc).  Though biggish divisors
  137. C would be normal internally from say mpn/generic/perfsqr.c.  On balance,
  138. C for the moment, it's felt the gain is not really enough to be worth the
  139. C trouble.
  140. C
  141. C
  142. C Enhancement:
  143. C
  144. C Process two source limbs per iteration using a two-limb inverse and a
  145. C sequence like
  146. C
  147. C ql  = low (c * il + sil) quotient low limb
  148. C qlc = high(c * il + sil)
  149. C qh1 = low (c * ih + sih) quotient high, partial
  150. C
  151. C cl = high (ql * d + c) carry out of low
  152. C qh = low (qlc * 1 + qh1) quotient high limb
  153. C
  154. C new c = high (qh * d + cl) carry out of high
  155. C
  156. C This would be 13 cycles/iteration, giving 6.5 cycles/limb.  The two limb
  157. C s*inverse as sih:sil = sh:sl * ih:il would be calculated off the dependent
  158. C chain with 4 multiplies.  The bigger inverse would take extra time to
  159. C calculate, but a one limb iteration to handle an odd size could be done as
  160. C soon as 64-bits of inverse were ready.
  161. C
  162. C Perhaps this could even extend to a 3 limb inverse, which might promise 17
  163. C or 18 cycles for 3 limbs, giving 5.66 or 6.0 cycles/limb.
  164. C
  165. ASM_START()
  166. .explicit
  167. .text
  168. .align 32
  169. .Ltable:
  170. data1 0,0x01, 0,0xAB, 0,0xCD, 0,0xB7, 0,0x39, 0,0xA3, 0,0xC5, 0,0xEF
  171. data1 0,0xF1, 0,0x1B, 0,0x3D, 0,0xA7, 0,0x29, 0,0x13, 0,0x35, 0,0xDF
  172. data1 0,0xE1, 0,0x8B, 0,0xAD, 0,0x97, 0,0x19, 0,0x83, 0,0xA5, 0,0xCF
  173. data1 0,0xD1, 0,0xFB, 0,0x1D, 0,0x87, 0,0x09, 0,0xF3, 0,0x15, 0,0xBF
  174. data1 0,0xC1, 0,0x6B, 0,0x8D, 0,0x77, 0,0xF9, 0,0x63, 0,0x85, 0,0xAF
  175. data1 0,0xB1, 0,0xDB, 0,0xFD, 0,0x67, 0,0xE9, 0,0xD3, 0,0xF5, 0,0x9F
  176. data1 0,0xA1, 0,0x4B, 0,0x6D, 0,0x57, 0,0xD9, 0,0x43, 0,0x65, 0,0x8F
  177. data1 0,0x91, 0,0xBB, 0,0xDD, 0,0x47, 0,0xC9, 0,0xB3, 0,0xD5, 0,0x7F
  178. data1 0,0x81, 0,0x2B, 0,0x4D, 0,0x37, 0,0xB9, 0,0x23, 0,0x45, 0,0x6F
  179. data1 0,0x71, 0,0x9B, 0,0xBD, 0,0x27, 0,0xA9, 0,0x93, 0,0xB5, 0,0x5F
  180. data1 0,0x61, 0,0x0B, 0,0x2D, 0,0x17, 0,0x99, 0,0x03, 0,0x25, 0,0x4F
  181. data1 0,0x51, 0,0x7B, 0,0x9D, 0,0x07, 0,0x89, 0,0x73, 0,0x95, 0,0x3F
  182. data1 0,0x41, 0,0xEB, 0,0x0D, 0,0xF7, 0,0x79, 0,0xE3, 0,0x05, 0,0x2F
  183. data1 0,0x31, 0,0x5B, 0,0x7D, 0,0xE7, 0,0x69, 0,0x53, 0,0x75, 0,0x1F
  184. data1 0,0x21, 0,0xCB, 0,0xED, 0,0xD7, 0,0x59, 0,0xC3, 0,0xE5, 0,0x0F
  185. data1 0,0x11, 0,0x3B, 0,0x5D, 0,0xC7, 0,0x49, 0,0x33, 0,0x55, 0,0xFF
  186. PROLOGUE(mpn_modexact_1c_odd)
  187. C r32 src
  188. C r33 size
  189. C r34 divisor
  190. C r35 carry
  191. .prologue
  192. .Lhere:
  193. { .mmi; add r33 = -1, r33 C M0  size-1
  194. mov r14 = 2 C M1  2
  195. mov r15 = ip C I0  .Lhere
  196. }{.mmi; setf.sig f6 = r34 C M2  divisor
  197. setf.sig f9 = r35 C M3  carry
  198. zxt1 r3 = r34 C I1  divisor low byte
  199. } ;;
  200. { .mmi; add r3 = .Ltable-.Lhere, r3 C M0  table offset ip and index
  201. sub r16 = 0, r34 C M1  -divisor
  202. .save ar.lc, r2
  203. mov r2 = ar.lc C I0
  204. }{.mmi; .body
  205. setf.sig f13 = r14 C M2  2 in significand
  206. mov r17 = -1 C M3  -1
  207. ABI32(` zxt4 r33 = r33') C I1  size extend
  208. } ;;
  209. { .mmi; add r3 = r3, r15 C M0  table entry address
  210. ABI32(` addp4 r32 = 0, r32') C M1  src extend
  211. mov ar.lc = r33 C I0  size-1 loop count
  212. }{.mmi; setf.sig f12 = r16 C M2  -divisor
  213. setf.sig f8 = r17 C M3  -1
  214. } ;;
  215. { .mmi; ld1 r3 = [r3] C M0  inverse, 8 bits
  216. ldf8 f10 = [r32], 8 C M1  src[0]
  217. cmp.ne p6,p0 = 0, r33 C I0  test size!=1
  218. } ;;
  219. C Wait for table load.
  220. C Hope for an L1 hit of 1 cycles to ALU, but could be more.
  221. setf.sig f7 = r3 C M2  inverse, 8 bits
  222. (p6) ldf8 f11 = [r32], 8 C M1  src[1], if size!=1
  223. ;;
  224. C 5 cycles
  225. C f6 divisor
  226. C f7 inverse, being calculated
  227. C f8 -1, will be -inverse
  228. C f9 carry
  229. C f10 src[0]
  230. C f11 src[1]
  231. C f12 -divisor
  232. C f13 2
  233. C f14 scratch
  234. xmpy.l f14 = f13, f7 C 2*i
  235. xmpy.l f7 = f7, f7 C i*i
  236. ;;
  237. xma.l f7 = f7, f12, f14 C i*i*-d + 2*i, inverse 16 bits
  238. ;;
  239. xmpy.l f14 = f13, f7 C 2*i
  240. xmpy.l f7 = f7, f7 C i*i
  241. ;;
  242. xma.l f7 = f7, f12, f14 C i*i*-d + 2*i, inverse 32 bits
  243. ;;
  244. xmpy.l f14 = f13, f7 C 2*i
  245. xmpy.l f7 = f7, f7 C i*i
  246. ;;
  247. xma.l f7 = f7, f12, f14 C i*i*-d + 2*i, inverse 64 bits
  248. xma.l f10 = f9, f8, f10 C sc = c * -1 + src[0]
  249. ;;
  250. ASSERT(p6, `
  251. xmpy.l f15 = f6, f7 ;; C divisor*inverse
  252. getf.sig r31 = f15 ;;
  253. cmp.eq p6,p0 = 1, r31 C should == 1
  254. ')
  255. xmpy.l f10 = f10, f7 C q = sc * inverse
  256. xmpy.l f8 = f7, f8 C -inverse = inverse * -1
  257. br.cloop.sptk.few.clr .Lentry C main loop, if size > 1
  258. ;;
  259. C size==1, finish up now
  260. xma.hu f9 = f10, f6, f9 C c = high(q * divisor + c)
  261. mov ar.lc = r2 C I0
  262. ;;
  263. getf.sig r8 = f9 C M2  return c
  264. br.ret.sptk.many b0
  265. .Ltop:
  266. C r2 saved ar.lc
  267. C f6 divisor
  268. C f7 inverse
  269. C f8 -inverse
  270. C f9 carry
  271. C f10 src[i] * inverse
  272. C f11 scratch src[i+1]
  273. add r16 = 160, r32
  274. ldf8 f11 = [r32], 8 C src[i+1]
  275. ;;
  276. C 2 cycles
  277. lfetch [r16]
  278. xma.l f10 = f9, f8, f10 C q = c * -inverse + si
  279. ;;
  280. C 3 cycles
  281. .Lentry:
  282. xma.hu f9 = f10, f6, f9 C c = high(q * divisor + c)
  283. xmpy.l f10 = f11, f7 C si = src[i] * inverse
  284. br.cloop.sptk.few.clr .Ltop
  285. ;;
  286. xma.l f10 = f9, f8, f10 C q = c * -inverse + si
  287. mov ar.lc = r2 C I0
  288. ;;
  289. xma.hu f9 = f10, f6, f9 C c = high(q * divisor + c)
  290. ;;
  291. getf.sig r8 = f9 C M2  return c
  292. br.ret.sptk.many b0
  293. EPILOGUE()