src/openpower/decoder/isa/remap_dct_yield.py

   1 # DCT "REMAP" scheduler
   2 #
   3 # Modifications made to create an in-place iterative DCT:
   4 # Copyright (c) 2021 Luke Kenneth Casson Leighton <lkcl@lkcl.net>
   5 #
   6 # SPDX: LGPLv3+
   7 #
   8 # Original fastdctlee.py by Nayuki:
   9 # Copyright (c) 2020 Project Nayuki. (MIT License)
  10 # https://www.nayuki.io/page/fast-discrete-cosine-transform-algorithms
  11
  12 import math
  13
  14 # bits of the integer 'val'.
  15 def reverse_bits(val, width):
  16     result = 0
  17     for _ in range(width):
  18         result = (result << 1) | (val & 1)
  19         val >>= 1
  20     return result
  21
  22
  23 # iterative version of [recursively-applied] half-rev.
  24 # relies on the list lengths being power-of-two and the fact
  25 # that bit-inversion of a list of binary numbers is the same
  26 # as reversing the order of the list
  27 # this version is dead easy to implement in hardware.
  28 # a big surprise is that the half-reversal can be done with
  29 # such a simple XOR. the inverse operation is slightly trickier
  30 def halfrev2(vec, pre_rev=True):
  31     res = []
  32     for i in range(len(vec)):
  33         if pre_rev:
  34             res.append(i ^ (i>>1))
  35         else:
  36             ri = i
  37             bl = i.bit_length()
  38             for ji in range(1, bl):
  39                 ri ^= (i >> ji)
  40             res.append(vec[ri])
  41     return res
  42
  43
  44 # python "yield" can be iterated. use this to make it clear how
  45 # the indices are generated by using natural-looking nested loops
  46 def iterate_dct_inner_butterfly_indices(SVSHAPE):
  47     # get indices to iterate over, in the required order
  48     n = SVSHAPE.lims[0]
  49     # createing lists of indices to iterate over in each dimension
  50     # has to be done dynamically, because it depends on the size
  51     # first, the size-based loop (which can be done statically)
  52     x_r = []
  53     size = 2
  54     while size <= n:
  55         x_r.append(size)
  56         size *= 2
  57     # invert order if requested
  58     if SVSHAPE.invxyz[0]:
  59         x_r.reverse()
  60
  61     if len(x_r) == 0:
  62         return
  63
  64     # reference (read/write) the in-place data in *reverse-bit-order*
  65     ri = list(range(n))
  66     if SVSHAPE.mode == 0b01:
  67         levels = n.bit_length() - 1
  68         ri = [ri[reverse_bits(i, levels)] for i in range(n)]
  69
  70     # reference list for not needing to do data-swaps, just swap what
  71     # *indices* are referenced (two levels of indirection at the moment)
  72     # pre-reverse the data-swap list so that it *ends up* in the order 0123..
  73     ji = list(range(n))
  74     inplace_mode = SVSHAPE.mode == 0b01 and SVSHAPE.skip not in [0b10, 0b11]
  75     if inplace_mode:
  76         #print ("inplace mode")
  77         ji = halfrev2(ji, True)
  78
  79     #print ("ri", ri)
  80     #print ("ji", ji)
  81
  82     # start an infinite (wrapping) loop
  83     skip = 0
  84     while True:
  85         for size in x_r:           # loop over 3rd order dimension (size)
  86             x_end = size == x_r[-1]
  87             # y_r schedule depends on size
  88             halfsize = size // 2
  89             y_r = []
  90             for i in range(0, n, size):
  91                 y_r.append(i)
  92             # invert if requested
  93             if SVSHAPE.invxyz[1]: y_r.reverse()
  94             for i in y_r:       # loop over 2nd order dimension
  95                 y_end = i == y_r[-1]
  96                 # two lists of half-range indices, e.g. j 0123, jr 7654
  97                 j = list(range(i, i + halfsize))
  98                 jr = list(range(i+halfsize, i + size))
  99                 jr.reverse()
 100                 # invert if requested
 101                 if SVSHAPE.invxyz[2]: j_r.reverse()
 102                 hz2 = halfsize // 2 # zero stops reversing 1-item lists
 103                 # if you *really* want to do the in-place swapping manually,
 104                 # this allows you to do it.  good luck...
 105                 if SVSHAPE.mode == 0b01 and not inplace_mode:
 106                     #print ("swap mode")
 107                     jr = j_r[:hz2]
 108                 #print ("xform jr", jr)
 109                 for jl, jh in zip(j, jr):   # loop over 1st order dimension
 110                     z_end = jl == j[-1]
 111                     # now depending on MODE return the index.  inner butterfly
 112                     if SVSHAPE.mode == 0b01:
 113                         if SVSHAPE.skip in [0b00, 0b10]:
 114                             result = ri[ji[jl]]        # lower half
 115                         elif SVSHAPE.skip in [0b01, 0b11]:
 116                             result = ri[ji[jh]] # upper half, reverse order
 117                     # outer butterfly
 118                     elif SVSHAPE.mode == 0b10:
 119                         if SVSHAPE.skip == 0b00:
 120                             result = ri[ji[jl]]        # lower half
 121                         elif SVSHAPE.skip == 0b01:
 122                             result = ri[ji[jl+size]]   # upper half
 123                     loopends = (z_end |
 124                                ((y_end and z_end)<<1) |
 125                                 ((y_end and x_end and z_end)<<2))
 126
 127                     yield result + SVSHAPE.offset, loopends
 128
 129                 # now in-place swap
 130                 if SVSHAPE.mode == 0b01 and inplace_mode:
 131                     for ci, (jl, jh) in enumerate(zip(j[:hz2], jr[:hz2])):
 132                         jlh = jl+halfsize
 133                         #print ("inplace swap", jh, jlh)
 134                         tmp1, tmp2 = ji[jlh], ji[jh]
 135                         ji[jlh], ji[jh] = tmp2, tmp1
 136
 137
 138 # python "yield" can be iterated. use this to make it clear how
 139 # the indices are generated by using natural-looking nested loops
 140 def iterate_dct_outer_butterfly_indices(SVSHAPE):
 141     # get indices to iterate over, in the required order
 142     n = SVSHAPE.lims[0]
 143     # createing lists of indices to iterate over in each dimension
 144     # has to be done dynamically, because it depends on the size
 145     # first, the size-based loop (which can be done statically)
 146     x_r = []
 147     size = n // 2
 148     while size >= 2:
 149         x_r.append(size)
 150         size //= 2
 151     # invert order if requested
 152     if SVSHAPE.invxyz[0]:
 153         x_r.reverse()
 154
 155     if len(x_r) == 0:
 156         return
 157
 158     #print ("outer butterfly")
 159
 160     # reference (read/write) the in-place data in *reverse-bit-order*
 161     ri = list(range(n))
 162     if SVSHAPE.mode == 0b11:
 163         levels = n.bit_length() - 1
 164         ri = [ri[reverse_bits(i, levels)] for i in range(n)]
 165
 166     # reference list for not needing to do data-swaps, just swap what
 167     # *indices* are referenced (two levels of indirection at the moment)
 168     # pre-reverse the data-swap list so that it *ends up* in the order 0123..
 169     ji = list(range(n))
 170     inplace_mode = SVSHAPE.skip in [0b10, 0b11]
 171     if inplace_mode:
 172         #print ("inplace mode", SVSHAPE.skip)
 173         ji = halfrev2(ji, True)
 174
 175     #print ("ri", ri)
 176     #print ("ji", ji)
 177
 178     # start an infinite (wrapping) loop
 179     skip = 0
 180     while True:
 181         for size in x_r:           # loop over 3rd order dimension (size)
 182             halfsize = size//2
 183             x_end = size == x_r[-1]
 184             y_r = list(range(0, halfsize))
 185             #print ("itersum", halfsize, size, y_r)
 186             # invert if requested
 187             if SVSHAPE.invxyz[1]: y_r.reverse()
 188             for i in y_r:       # loop over 2nd order dimension
 189                 y_end = i == y_r[-1]
 190                 # one list to create iterative-sum schedule
 191                 jr = list(range(i+halfsize, i+n-halfsize, size))
 192                 #print ("itersum     jr", i+halfsize, i+size, jr)
 193                 # invert if requested
 194                 if SVSHAPE.invxyz[2]: j_r.reverse()
 195                 hz2 = halfsize // 2 # zero stops reversing 1-item lists
 196                 for jh in jr:   # loop over 1st order dimension
 197                     z_end = jh == jr[-1]
 198                     #print ("     itersum", size, i, jh, jh+size)
 199                     if SVSHAPE.skip in [0b00, 0b10]:
 200                         result = ri[ji[jh]]        # lower half
 201                     elif SVSHAPE.skip in [0b01, 0b11]:
 202                         result = ri[ji[jh+size]] # upper half
 203                     loopends = (z_end |
 204                                ((y_end and z_end)<<1) |
 205                                 ((y_end and x_end and z_end)<<2))
 206
 207                     yield result + SVSHAPE.offset, loopends
 208
 209                 # now in-place swap
 210                 if SVSHAPE.mode == 0b11 and inplace_mode:
 211                     j = list(range(i, i + halfsize))
 212                     jr = list(range(i+halfsize, i + size))
 213                     jr.reverse()
 214                     for ci, (jl, jh) in enumerate(zip(j[:hz2], jr[:hz2])):
 215                         jlh = jl+halfsize
 216                         #print ("inplace swap", jh, jlh)
 217                         tmp1, tmp2 = ji[jlh], ji[jh]
 218                         ji[jlh], ji[jh] = tmp2, tmp1
 219
 220
 221 def pprint_schedule(schedule, n):
 222     size = 2
 223     idx = 0
 224     while size <= n:
 225         halfsize = size // 2
 226         tablestep = n // size
 227         print ("size %d halfsize %d tablestep %d" % \
 228                 (size, halfsize, tablestep))
 229         for i in range(0, n, size):
 230             prefix = "i %d\t" % i
 231             for j in range(i, i + halfsize):
 232                 (jl, je), (jh, he) = schedule[idx]
 233                 print ("  %-3d\t%s j=%-2d jh=%-2d "
 234                         "j[jl=%-2d] j[jh=%-2d]" % \
 235                                 (idx, prefix, j, j+halfsize,
 236                                       jl, jh,
 237                                 ),
 238                                 "end", bin(je)[2:], bin(je)[2:])
 239                 idx += 1
 240         size *= 2
 241
 242 def pprint_schedule_outer(schedule, n):
 243     size = 2
 244     idx = 0
 245     while size <= n//2:
 246         halfsize = size // 2
 247         tablestep = n // size
 248         print ("size %d halfsize %d tablestep %d" % \
 249                 (size, halfsize, tablestep))
 250         y_r = list(range(0, halfsize))
 251         for i in y_r:
 252             prefix = "i %d\t" % i
 253             jr = list(range(i+halfsize, i+n-halfsize, size))
 254             for j in jr:
 255                 (jl, je), (jh, he) = schedule[idx]
 256                 print ("  %-3d\t%s j=%-2d jh=%-2d "
 257                         "j[jl=%-2d] j[jh=%-2d]" % \
 258                                 (idx, prefix, j, j+halfsize,
 259                                       jl, jh,
 260                                 ),
 261                                 "end", bin(je)[2:], bin(je)[2:])
 262                 idx += 1
 263         size *= 2
 264
 265
 266 # totally cool *in-place* DCT algorithm using yield REMAPs
 267 def transform2(vec):
 268
 269     # Initialization
 270     n = len(vec)
 271     print ()
 272     print ("transform2", n)
 273     levels = n.bit_length() - 1
 274
 275     # reference (read/write) the in-place data in *reverse-bit-order*
 276     ri = list(range(n))
 277     ri = [ri[reverse_bits(i, levels)] for i in range(n)]
 278
 279     # and pretend we LDed data in half-swapped *and* bit-reversed order as well
 280     # TODO: merge these two
 281     vec = halfrev2(vec, False)
 282     vec = [vec[ri[i]] for i in range(n)]
 283
 284     # create a cos table: not strictly necessary but here for illustrative
 285     # purposes, to demonstrate the point that it really *is* iterative.
 286     # this table could be cached and used multiple times rather than
 287     # computed every time.
 288     ctable = []
 289     size = n
 290     while size >= 2:
 291         halfsize = size // 2
 292         for i in range(n//size):
 293             for ci in range(halfsize):
 294                 ctable.append((math.cos((ci + 0.5) * math.pi / size) * 2.0))
 295         size //= 2
 296
 297     ################
 298     # INNER butterfly
 299     ################
 300     xdim = n
 301     ydim = 0
 302     zdim = 0
 303
 304     # set up an SVSHAPE
 305     class SVSHAPE:
 306         pass
 307     # j schedule
 308     SVSHAPE0 = SVSHAPE()
 309     SVSHAPE0.lims = [xdim, ydim, zdim]
 310     SVSHAPE0.order = [0,1,2]  # experiment with different permutations, here
 311     SVSHAPE0.mode = 0b01
 312     SVSHAPE0.skip = 0b00
 313     SVSHAPE0.offset = 0       # experiment with different offset, here
 314     SVSHAPE0.invxyz = [1,0,0] # inversion if desired
 315     # j+halfstep schedule
 316     SVSHAPE1 = SVSHAPE()
 317     SVSHAPE1.lims = [xdim, ydim, zdim]
 318     SVSHAPE1.order = [0,1,2]  # experiment with different permutations, here
 319     SVSHAPE1.mode = 0b01
 320     SVSHAPE1.skip = 0b01
 321     SVSHAPE1.offset = 0       # experiment with different offset, here
 322     SVSHAPE1.invxyz = [1,0,0] # inversion if desired
 323
 324     # enumerate over the iterator function, getting new indices
 325     i0 = iterate_dct_inner_butterfly_indices(SVSHAPE0)
 326     i1 = iterate_dct_inner_butterfly_indices(SVSHAPE1)
 327     for k, ((jl, jle), (jh, jhe)) in enumerate(zip(i0, i1)):
 328         t1, t2 = vec[jl], vec[jh]
 329         coeff = ctable[k]
 330         vec[jl] = t1 + t2
 331         vec[jh] = (t1 - t2) * (1/coeff)
 332         print ("coeff", size, i, "ci", ci,
 333                 "jl", jl, "jh", jh,
 334                "i/n", (ci+0.5)/size, coeff, vec[jl],
 335                                             vec[jh],
 336                 "end", bin(jle), bin(jhe))
 337         if jle == 0b111: # all loops end
 338             break
 339
 340     print("transform2 pre-itersum", vec)
 341
 342     # now things are in the right order for the outer butterfly.
 343
 344     # j schedule
 345     SVSHAPE0 = SVSHAPE()
 346     SVSHAPE0.lims = [xdim, ydim, zdim]
 347     SVSHAPE0.order = [0,1,2]  # experiment with different permutations, here
 348     SVSHAPE0.mode = 0b10
 349     SVSHAPE0.skip = 0b00
 350     SVSHAPE0.offset = 0       # experiment with different offset, here
 351     SVSHAPE0.invxyz = [0,0,0] # inversion if desired
 352     # j+halfstep schedule
 353     SVSHAPE1 = SVSHAPE()
 354     SVSHAPE1.lims = [xdim, ydim, zdim]
 355     SVSHAPE1.order = [0,1,2]  # experiment with different permutations, here
 356     SVSHAPE1.mode = 0b10
 357     SVSHAPE1.skip = 0b01
 358     SVSHAPE1.offset = 0       # experiment with different offset, here
 359     SVSHAPE1.invxyz = [0,0,0] # inversion if desired
 360
 361     # enumerate over the iterator function, getting new indices
 362     i0 = iterate_dct_outer_butterfly_indices(SVSHAPE0)
 363     i1 = iterate_dct_outer_butterfly_indices(SVSHAPE1)
 364     for k, ((jl, jle), (jh, jhe)) in enumerate(zip(i0, i1)):
 365         print ("itersum    jr", jl, jh,
 366                 "end", bin(jle), bin(jhe))
 367         vec[jl] += vec[jh]
 368         size //= 2
 369         if jle == 0b111: # all loops end
 370             break
 371
 372     print("transform2 result", vec)
 373
 374     return vec
 375
 376
 377 def demo():
 378     # set the dimension sizes here
 379     n = 8
 380     xdim = n
 381     ydim = 0 # not needed
 382     zdim = 0 # again, not needed
 383
 384
 385     ################
 386     # INNER butterfly
 387     ################
 388
 389     # set up an SVSHAPE
 390     class SVSHAPE:
 391         pass
 392     # j schedule
 393     SVSHAPE0 = SVSHAPE()
 394     SVSHAPE0.lims = [xdim, ydim, zdim]
 395     SVSHAPE0.order = [0,1,2]  # experiment with different permutations, here
 396     SVSHAPE0.mode = 0b01
 397     SVSHAPE0.skip = 0b00
 398     SVSHAPE0.offset = 0       # experiment with different offset, here
 399     SVSHAPE0.invxyz = [0,0,0] # inversion if desired
 400     # j+halfstep schedule
 401     SVSHAPE1 = SVSHAPE()
 402     SVSHAPE1.lims = [xdim, ydim, zdim]
 403     SVSHAPE1.order = [0,1,2]  # experiment with different permutations, here
 404     SVSHAPE1.mode = 0b01
 405     SVSHAPE1.skip = 0b01
 406     SVSHAPE1.offset = 0       # experiment with different offset, here
 407     SVSHAPE1.invxyz = [0,0,0] # inversion if desired
 408
 409     # enumerate over the iterator function, getting new indices
 410     schedule = []
 411     i0 = iterate_dct_inner_butterfly_indices(SVSHAPE0)
 412     i1 = iterate_dct_inner_butterfly_indices(SVSHAPE1)
 413     for idx, (jl, jh) in enumerate(zip(i0, i1)):
 414         schedule.append((jl, jh))
 415         if jl[1] == 0b111: # end
 416             break
 417
 418     # ok now pretty-print the results, with some debug output
 419     print ("inner butterfly")
 420     pprint_schedule(schedule, n)
 421     print ("")
 422
 423     ################
 424     # outer butterfly
 425     ################
 426
 427     # j schedule
 428     SVSHAPE0 = SVSHAPE()
 429     SVSHAPE0.lims = [xdim, ydim, zdim]
 430     SVSHAPE0.order = [0,1,2]  # experiment with different permutations, here
 431     SVSHAPE0.mode = 0b10
 432     SVSHAPE0.skip = 0b10
 433     SVSHAPE0.offset = 0       # experiment with different offset, here
 434     SVSHAPE0.invxyz = [1,0,0] # inversion if desired
 435     # j+halfstep schedule
 436     SVSHAPE1 = SVSHAPE()
 437     SVSHAPE1.lims = [xdim, ydim, zdim]
 438     SVSHAPE1.order = [0,1,2]  # experiment with different permutations, here
 439     SVSHAPE1.mode = 0b10
 440     SVSHAPE1.skip = 0b11
 441     SVSHAPE1.offset = 0       # experiment with different offset, here
 442     SVSHAPE1.invxyz = [1,0,0] # inversion if desired
 443
 444     # enumerate over the iterator function, getting new indices
 445     schedule = []
 446     i0 = iterate_dct_outer_butterfly_indices(SVSHAPE0)
 447     i1 = iterate_dct_outer_butterfly_indices(SVSHAPE1)
 448     for idx, (jl, jh) in enumerate(zip(i0, i1)):
 449         schedule.append((jl, jh))
 450         if jl[1] == 0b111: # end
 451             break
 452
 453     # ok now pretty-print the results, with some debug output
 454     print ("outer butterfly")
 455     pprint_schedule_outer(schedule, n)
 456
 457 # run the demo
 458 if __name__ == '__main__':
 459     demo()