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

数学计算

开发平台:

Unix_Linux

  1. dnl  Alpha mpn_modexact_1c_odd -- mpn exact remainder
  2. dnl  Copyright 2003, 2004 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 EV4:    47
  21. C EV5:    30
  22. C EV6:    15
  23. C mp_limb_t mpn_modexact_1c_odd (mp_srcptr src, mp_size_t size, mp_limb_t d,
  24. C                                mp_limb_t c)
  25. C
  26. C This code follows the "alternate" code in mpn/generic/mode1o.c,
  27. C eliminating cbit+climb from the dependent chain.  This leaves,
  28. C
  29. C        ev4   ev5   ev6
  30. C         1     3     1    subq   y = x - h
  31. C        23    13     7    mulq   q = y * inverse
  32. C        23    14     7    umulh  h = high (q * d)
  33. C        --    --    --
  34. C        47    30    15
  35. C
  36. C In each case, the load latency, loop control, and extra carry bit handling
  37. C hide under the multiply latencies.  Those latencies are long enough that
  38. C we don't need to worry about alignment or pairing to squeeze out
  39. C performance.
  40. C
  41. C For the first limb, some of the loop code is broken out and scheduled back
  42. C since it can be done earlier.
  43. C
  44. C   - The first ldq src[0] is near the start of the routine, for maximum
  45. C     time from memory.
  46. C
  47. C   - The subq y=x-climb can be done without waiting for the inverse.
  48. C
  49. C   - The mulq y*inverse is replicated after the final subq for the inverse,
  50. C     instead of branching to the mulq in the main loop.  On ev4 a branch
  51. C     there would cost cycles, but we can hide them under the mulq latency.
  52. C
  53. C For the last limb, high<divisor is tested and if that's true a subtract
  54. C and addback is done, as per the main mpn/generic/mode1o.c code.  This is a
  55. C data-dependent branch, but we're waiting for umulh so any penalty should
  56. C hide there.  The multiplies saved would be worth the cost anyway.
  57. C
  58. C Enhancements:
  59. C
  60. C For size==1, a plain division (done bitwise say) might be faster than
  61. C calculating an inverse, the latter taking about 130 cycles on ev4 or 70 on
  62. C ev5.  A call to gcc __remqu might be a possibility.
  63. ASM_START()
  64. PROLOGUE(mpn_modexact_1c_odd,gp)
  65. C r16 src
  66. C r17 size
  67. C r18 d
  68. C r19 c
  69. LEA(r0, binvert_limb_table)
  70. srl r18, 1, r20 C d >> 1
  71. and r20, 127, r20 C idx = d>>1 & 0x7F
  72. addq r0, r20, r21 C table + idx
  73. ifelse(bwx_available_p,1,
  74. ` ldbu r20, 0(r21) C table[idx], inverse 8 bits
  75. ',`
  76. ldq_u r20, 0(r21) C table[idx] qword
  77. extbl r20, r21, r20 C table[idx], inverse 8 bits
  78. ')
  79. mull r20, r20, r7 C i*i
  80. addq r20, r20, r20 C 2*i
  81. ldq r2, 0(r16) C x = s = src[0]
  82. lda r17, -1(r17) C size--
  83. clr r0 C initial cbit=0
  84. mull r7, r18, r7 C i*i*d
  85. subq r20, r7, r20 C 2*i-i*i*d, inverse 16 bits
  86. mull r20, r20, r7 C i*i
  87. addq r20, r20, r20 C 2*i
  88. mull r7, r18, r7 C i*i*d
  89. subq r20, r7, r20 C 2*i-i*i*d, inverse 32 bits
  90. mulq r20, r20, r7 C i*i
  91. addq r20, r20, r20 C 2*i
  92. mulq r7, r18, r7 C i*i*d
  93. subq r2, r19, r3 C y = x - climb
  94. subq r20, r7, r20 C inv = 2*i-i*i*d, inverse 64 bits
  95. ASSERT(r7, C should have d*inv==1 mod 2^64
  96. ` mulq r18, r20, r7
  97. cmpeq r7, 1, r7')
  98. mulq r3, r20, r4 C first q = y * inv
  99. beq r17, L(one) C if size==1
  100. br L(entry)
  101. L(top):
  102. C r0 cbit
  103. C r16 src, incrementing
  104. C r17 size, decrementing
  105. C r18 d
  106. C r19 climb
  107. C r20 inv
  108. ldq r1, 0(r16) C s = src[i]
  109. subq r1, r0, r2 C x = s - cbit
  110. cmpult r1, r0, r0 C new cbit = s < cbit
  111. subq r2, r19, r3 C y = x - climb
  112. mulq r3, r20, r4 C q = y * inv
  113. L(entry):
  114. cmpult r2, r19, r5 C cbit2 = x < climb
  115. addq r5, r0, r0 C cbit += cbit2
  116. lda r16, 8(r16) C src++
  117. lda r17, -1(r17) C size--
  118. umulh r4, r18, r19 C climb = q * d
  119. bne r17, L(top) C while 2 or more limbs left
  120. C r0 cbit
  121. C r18 d
  122. C r19 climb
  123. C r20 inv
  124. ldq r1, 0(r16) C s = src[size-1] high limb
  125. cmpult r1, r18, r2 C test high<divisor
  126. bne r2, L(skip) C skip if so
  127. C can't skip a division, repeat loop code
  128. subq r1, r0, r2 C x = s - cbit
  129. cmpult r1, r0, r0 C new cbit = s < cbit
  130. subq r2, r19, r3 C y = x - climb
  131. mulq r3, r20, r4 C q = y * inv
  132. L(one):
  133. cmpult r2, r19, r5 C cbit2 = x < climb
  134. addq r5, r0, r0 C cbit += cbit2
  135. umulh r4, r18, r19 C climb = q * d
  136. addq r19, r0, r0 C return climb + cbit
  137. ret r31, (r26), 1
  138. ALIGN(8)
  139. L(skip):
  140. C with high<divisor, the final step can be just (cbit+climb)-s and
  141. C an addback of d if that underflows
  142. addq r19, r0, r19 C c = climb + cbit
  143. subq r19, r1, r2 C c - s
  144. cmpult r19, r1, r3 C c < s
  145. addq r2, r18, r0 C return c-s + divisor
  146. cmoveq r3, r2, r0 C return c-s if no underflow
  147. ret r31, (r26), 1
  148. EPILOGUE()
  149. ASM_END()