NTT数论变换（四）

上一篇：https://blog.csdn.net/wxkhturfun/article/details/111937428

Stockham算法（采用GS蝶形运算单元）

关于该算法，可以详细见链接：http://wwwa.pikara.ne.jp/okojisan/otfft-en/stockham1.html
Stockham和FFT的区别在于是否需要额外进行码位变换，在基于库里-图基的FFT(或NTT)中，需要首先进行Rader变换，而基于桑德-图基的FFT（或NTT）中，最后需要进行Rader变换，而存储数组的下标是始终在变化的。以库里-图基为例，当进行8点FFT(或NTT)时，其变换图如https://blog.csdn.net/wxkhturfun/article/details/111186205中的图如示，变换为：

0 <==> 4        :wn[ 1 ]^ 0
*****小循环*****
2 <==> 6        :wn[ 1 ]^ 0
*****小循环*****
1 <==> 5        :wn[ 1 ]^ 0
*****小循环*****
3 <==> 7        :wn[ 1 ]^ 0
*****小循环*****
============大循环============
0 <==> 2        :wn[ 2 ]^ 0
4 <==> 6        :wn[ 2 ]^ 1
*****小循环*****
1 <==> 3        :wn[ 2 ]^ 0
5 <==> 7        :wn[ 2 ]^ 1
*****小循环*****
============大循环============
0 <==> 1        :wn[ 3 ]^ 0
4 <==> 5        :wn[ 3 ]^ 1
2 <==> 3        :wn[ 3 ]^ 2
6 <==> 7        :wn[ 3 ]^ 3
*****小循环*****
============大循环============

也就是说在大循环中，数组下标为非0，非7时，在三层运算（log8 = 3）时，其对应的数组下标都是在变化的，零和七是个例外。Stockham就是在每层大循环时，将数组有目的地写入到对应的位置，保证其下一次大循环时需要进行蝶形运算操作的双方的数组下标不变。这种有目的操作看起来就像是一个trick，但是在硬件实现时，这更易对其进行流水实现！

def Stockham(x,order):
    '''
        Input: n-length signal x = (x0, x1, · · · , xn−1), where n = 2k. ω, n-th
        primitive root of unity. An external n-length continuous memory
        y.
        Output: FFT result of vector x
    '''
    k=order
    n=2**order
    J = 1
    #w=xxxxx
    w=1
    y=x
    for  i in range(k):
        J=2**i
        for s in range(int(n/(2*J))):
            rc =w**(J*s)
            for j in range(J):
                var=s*J 
                y[2*var+j] = x[var+j] + x[var +j+ int(n/2)]
                y[(2*s + 1)*J + j] = rc*(x[var+j]- x[var +j+ int(n/2)])
                print('y[',str(2*var+j),']=x[',str(var+j),']+x[',str(var +j+ int(n/2)),']')
                print('y[',str((2*s + 1)*J + j),']=w^(',str(J*s),')*(x[',str(var+j),']-x[',str(var +j+ int(n/2)),'])')
            x,y = y,x
        print('==============',i,'====================')
    return x
if __name__=='__main__':
    x=[]
    order=3
    for i in range(2**order):
        x.append(i)
    Stockham(x,order)

打印的结果如下 ：

y[ 0 ]=x[ 0 ]+x[ 4 ]
y[ 1 ]=w^( 0 )*(x[ 0 ]-x[ 4 ])       
y[ 2 ]=x[ 1 ]+x[ 5 ]
y[ 3 ]=w^( 1 )*(x[ 1 ]-x[ 5 ])       
y[ 4 ]=x[ 2 ]+x[ 6 ]
y[ 5 ]=w^( 2 )*(x[ 2 ]-x[ 6 ])       
y[ 6 ]=x[ 3 ]+x[ 7 ]
y[ 7 ]=w^( 3 )*(x[ 3 ]-x[ 7 ])       
============== 0 ====================
y[ 0 ]=x[ 0 ]+x[ 4 ]
y[ 2 ]=w^( 0 )*(x[ 0 ]-x[ 4 ])       
y[ 1 ]=x[ 1 ]+x[ 5 ]
y[ 3 ]=w^( 0 )*(x[ 1 ]-x[ 5 ])       
y[ 4 ]=x[ 2 ]+x[ 6 ]
y[ 6 ]=w^( 2 )*(x[ 2 ]-x[ 6 ])       
y[ 5 ]=x[ 3 ]+x[ 7 ]
y[ 7 ]=w^( 2 )*(x[ 3 ]-x[ 7 ])       
============== 1 ====================
y[ 0 ]=x[ 0 ]+x[ 4 ]
y[ 4 ]=w^( 0 )*(x[ 0 ]-x[ 4 ])       
y[ 1 ]=x[ 1 ]+x[ 5 ]
y[ 5 ]=w^( 0 )*(x[ 1 ]-x[ 5 ])       
y[ 2 ]=x[ 2 ]+x[ 6 ]
y[ 6 ]=w^( 0 )*(x[ 2 ]-x[ 6 ])       
y[ 3 ]=x[ 3 ]+x[ 7 ]
y[ 7 ]=w^( 0 )*(x[ 3 ]-x[ 7 ])       
============== 2 ====================

可以见到x[0]始终与x[4]进行蝶形运算，依次类推！