mpn/ia64/mode1o.asm

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