return res
+# python "yield" can be iterated. use this to make it clear how
+# the indices are generated by using natural-looking nested loops
+def iterate_dct_inner_costable_indices(SVSHAPE):
+ # get indices to iterate over, in the required order
+ n = SVSHAPE.lims[0]
+ mode = SVSHAPE.lims[1]
+ print ("inner costable", mode)
+ # creating lists of indices to iterate over in each dimension
+ # has to be done dynamically, because it depends on the size
+ # first, the size-based loop (which can be done statically)
+ x_r = []
+ size = 2
+ while size <= n:
+ x_r.append(size)
+ size *= 2
+ # invert order if requested
+ if SVSHAPE.invxyz[0]:
+ x_r.reverse()
+
+ if len(x_r) == 0:
+ return
+
+ #print ("ri", ri)
+ #print ("ji", ji)
+
+ # start an infinite (wrapping) loop
+ skip = 0
+ z_end = 1 # doesn't exist in this, only 2 loops
+ k = 0
+ while True:
+ for size in x_r: # loop over 3rd order dimension (size)
+ x_end = size == x_r[-1]
+ # y_r schedule depends on size
+ halfsize = size // 2
+ y_r = []
+ for i in range(0, n, size):
+ y_r.append(i)
+ # invert if requested
+ if SVSHAPE.invxyz[1]: y_r.reverse()
+ # two lists of half-range indices, e.g. j 0123, jr 7654
+ j = list(range(0, halfsize))
+ # invert if requested
+ if SVSHAPE.invxyz[2]: j_r.reverse()
+ #print ("xform jr", jr)
+ # loop over 1st order dimension
+ for ci, jl in enumerate(j):
+ y_end = jl == j[-1]
+ # now depending on MODE return the index. inner butterfly
+ if SVSHAPE.skip == 0b00: # in [0b00, 0b10]:
+ result = k # offset into COS table
+ elif SVSHAPE.skip == 0b10: #
+ result = ci # coefficient helper
+ elif SVSHAPE.skip == 0b11: #
+ result = size # coefficient helper
+ loopends = (z_end |
+ ((y_end and z_end)<<1) |
+ ((y_end and x_end and z_end)<<2))
+
+ yield result + SVSHAPE.offset, loopends
+ k += 1
+
# python "yield" can be iterated. use this to make it clear how
# the indices are generated by using natural-looking nested loops
def iterate_dct_inner_butterfly_indices(SVSHAPE):
# start an infinite (wrapping) loop
skip = 0
+ k = 0
+ k_start = 0
while True:
for size in x_r: # loop over 3rd order dimension (size)
halfsize = size//2
# invert if requested
if SVSHAPE.invxyz[2]: j_r.reverse()
hz2 = halfsize // 2 # zero stops reversing 1-item lists
+ k = k_start
for ci, jh in enumerate(jr): # loop over 1st order dimension
z_end = jh == jr[-1]
#print (" itersum", size, i, jh, jh+size)
- if SVSHAPE.skip == 0b00: # in [0b00, 0b10]:
- result = ri[ji[jh]] # lower half
- elif SVSHAPE.skip == 0b01: # in [0b01, 0b11]:
- result = ri[ji[jh+size]] # upper half
- elif SVSHAPE.skip == 0b10: #
- result = ci # coefficient helper
- elif SVSHAPE.skip == 0b11: #
- result = size # coefficient helper
+ if mode == 4:
+ # COS table pre-generated mode
+ if SVSHAPE.skip == 0b00: # in [0b00, 0b10]:
+ result = ri[ji[jh]] # lower half
+ elif SVSHAPE.skip == 0b01: # in [0b01, 0b11]:
+ result = ri[ji[jh+size]] # upper half
+ elif SVSHAPE.skip == 0b10: #
+ result = k # cos table offset
+ else:
+ # COS table generated on-demand ("Vertical-First") mode
+ if SVSHAPE.skip == 0b00: # in [0b00, 0b10]:
+ result = ri[ji[jh]] # lower half
+ elif SVSHAPE.skip == 0b01: # in [0b01, 0b11]:
+ result = ri[ji[jh+size]] # upper half
+ elif SVSHAPE.skip == 0b10: #
+ result = ci # coefficient helper
+ elif SVSHAPE.skip == 0b11: #
+ result = size # coefficient helper
loopends = (z_end |
((y_end and z_end)<<1) |
((y_end and x_end and z_end)<<2))
yield result + SVSHAPE.offset, loopends
+ k += 1
# now in-place swap
if SVSHAPE.submode2 == 0b11 and inplace_mode:
tmp1, tmp2 = ji[jlh], ji[jh]
ji[jlh], ji[jh] = tmp2, tmp1
+ # new k_start point for cos tables( runs inside x_r loop NOT i loop)
+ k_start += halfsize
+
def pprint_schedule(schedule, n):
size = 2
print ("transform2", n)
levels = n.bit_length() - 1
+ # set up dims
+ xdim = n
+
# reference (read/write) the in-place data in *reverse-bit-order*
ri = list(range(n))
ri = [ri[reverse_bits(i, levels)] for i in range(n)]
size = n
while size >= 2:
halfsize = size // 2
- for i in range(n//size):
- for ci in range(halfsize):
- ctable.append((math.cos((ci + 0.5) * math.pi / size) * 2.0))
+ for ci in range(halfsize):
+ coeff = (math.cos((ci + 0.5) * math.pi / size) * 2.0)
+ ctable.append(coeff)
+ print ("coeff", size, "ci", ci, "k", len(ctable)-1,
+ "i/n", (ci+0.5)/size, coeff)
size //= 2
+ # set up an SVSHAPE
+ class SVSHAPE:
+ pass
+ # ci schedule
+ SVSHAPE0 = SVSHAPE()
+ SVSHAPE0.lims = [xdim, 4, 0]
+ SVSHAPE0.mode = 0b01
+ SVSHAPE0.submode2 = 0b01
+ SVSHAPE0.skip = 0b10
+ SVSHAPE0.offset = 0 # experiment with different offset, here
+ SVSHAPE0.invxyz = [1,0,0] # inversion if desired
+ # size schedule
+ SVSHAPE1 = SVSHAPE()
+ SVSHAPE1.lims = [xdim, 4, 0]
+ SVSHAPE1.mode = 0b01
+ SVSHAPE1.submode2 = 0b01
+ SVSHAPE1.skip = 0b11
+ SVSHAPE1.offset = 0 # experiment with different offset, here
+ SVSHAPE1.invxyz = [1,0,0] # inversion if desired
+ # k schedule
+ SVSHAPE2 = SVSHAPE()
+ SVSHAPE2.lims = [xdim, 4, 0]
+ SVSHAPE2.mode = 0b01
+ SVSHAPE2.submode2 = 0b01
+ SVSHAPE2.skip = 0b00
+ SVSHAPE2.offset = 0 # experiment with different offset, here
+ SVSHAPE2.invxyz = [1,0,0] # inversion if desired
+
+ # enumerate over the iterator function, getting new indices
+ i0 = iterate_dct_inner_costable_indices(SVSHAPE0)
+ i1 = iterate_dct_inner_costable_indices(SVSHAPE1)
+ i2 = iterate_dct_inner_costable_indices(SVSHAPE2)
+ for ((ci, cie), (size, sze), (k, ke)) in \
+ zip(i0, i1, i2):
+ print ("xform2 cos", ci, size, k)
+ coeff = (math.cos((ci + 0.5) * math.pi / size) * 2.0)
+ assert coeff == ctable[k]
+ print ("coeff", size, "ci", ci, "k", k,
+ "i/n", (ci+0.5)/size, coeff,
+ "end", bin(cie), bin(sze), bin(ke))
+ if cie == 0b111: # all loops end
+ break
+
################
# INNER butterfly
################
- xdim = n
- ydim = 0
- zdim = 0
- # set up an SVSHAPE
- class SVSHAPE:
- pass
# j schedule
SVSHAPE0 = SVSHAPE()
- SVSHAPE0.lims = [xdim, 0b000001, zdim]
+ SVSHAPE0.lims = [xdim, 0b000001, 0]
SVSHAPE0.mode = 0b01
SVSHAPE0.submode2 = 0b01
SVSHAPE0.skip = 0b00
SVSHAPE0.invxyz = [1,0,0] # inversion if desired
# j+halfstep schedule
SVSHAPE1 = SVSHAPE()
- SVSHAPE1.lims = [xdim, 0b000001, zdim]
+ SVSHAPE1.lims = [xdim, 0b000001, 0]
SVSHAPE1.mode = 0b01
SVSHAPE1.submode2 = 0b01
SVSHAPE1.skip = 0b01
SVSHAPE1.invxyz = [1,0,0] # inversion if desired
# ci schedule
SVSHAPE2 = SVSHAPE()
- SVSHAPE2.lims = [xdim, 0b000001, zdim]
+ SVSHAPE2.lims = [xdim, 0b000001, 0]
SVSHAPE2.mode = 0b01
SVSHAPE2.submode2 = 0b01
SVSHAPE2.skip = 0b10
SVSHAPE2.invxyz = [1,0,0] # inversion if desired
# size schedule
SVSHAPE3 = SVSHAPE()
- SVSHAPE3.lims = [xdim, 0b000001, zdim]
+ SVSHAPE3.lims = [xdim, 0b000001, 0]
SVSHAPE3.mode = 0b01
SVSHAPE3.submode2 = 0b01
SVSHAPE3.skip = 0b11
t1, t2 = vec[jl], vec[jh]
print ("xform2", jl, jh, ci, size)
coeff = (math.cos((ci + 0.5) * math.pi / size) * 2.0)
- assert coeff == ctable[k]
+ #assert coeff == ctable[k]
vec[jl] = t1 + t2
vec[jh] = (t1 - t2) * (1/coeff)
- print ("coeff", size, i, "ci", ci,
+ print ("coeff", size, "ci", ci,
"jl", jl, "jh", jh,
"i/n", (ci+0.5)/size, coeff, vec[jl],
vec[jh],
# j schedule
SVSHAPE0 = SVSHAPE()
- SVSHAPE0.lims = [xdim, 0b0000010, zdim]
+ SVSHAPE0.lims = [xdim, 0b0000010, 0]
SVSHAPE0.submode2 = 0b100
SVSHAPE0.mode = 0b01
SVSHAPE0.skip = 0b00
SVSHAPE0.invxyz = [0,0,0] # inversion if desired
# j+halfstep schedule
SVSHAPE1 = SVSHAPE()
- SVSHAPE1.lims = [xdim, 0b0000010, zdim]
+ SVSHAPE1.lims = [xdim, 0b0000010, 0]
SVSHAPE1.mode = 0b01
SVSHAPE1.submode2 = 0b100
SVSHAPE1.skip = 0b01