> restart;
> with(linalg):
Warning, new definition for norm
Warning, new definition for trace
> # Sande-Tukey-Algorithmus
> # Werte im Ausgabevektor in bitreversaler Reihenfolge wieder in a gespeichert!
> SandeTuckey:=proc(a)
> # a Vektor von Zweierpotenzlaenge N=2^t
> # option trace;
> local N,m,m1,w,l,r,j,s,k,t;
> N:=vectdim(a); t:=round(log[2](N));  m:=N; m1:=1;
> for k from 1 to t do
>   m:=m/2;  w:=evalf(cos(Pi/m)-I*sin(Pi/m));
>   for l from 0 to m1-1 do
>       for r from 0 to m-1 do
>          j:=r+2*l*m;
>          s:=a[j+1]+a[m+j+1];
>          a[m+j+1]:=(a[j+1]-a[m+j+1])*w^r;
>          a[j+1]:=s
>       od
>   od;
>   m1:=m1*2
> od;
> evalm(a)
> end:
> a:=vector(16,[1,2,3,4,5,6,7,8,8,7,6,5,4,3,2,1]);
>  SandeTuckey(a);
> 

        a := [1, 2, 3, 4, 5, 6, 7, 8, 8, 7, 6, 5, 4, 3, 2, 1]


[72, 0, 0, 0, 0, 0, 0, 0, -25.27414238 - 5.027339483 I,

    -.03956612 + .198912363 I, -.446462695 - .668178642 I,

    -2.239828801 + 1.496605762 I, -2.239828816 - 1.496605762 I,

    -.4464626917 + .668178630 I, -.03956613 - .198912372 I,

    -25.27414235 + 5.027339504 I]

> # Vergleich 
> a:=vector(16,[1,2,3,4,5,6,7,8,8,7,6,5,4,3,2,1]);
> dft:=proc(x)
> local n,w,j,k,r,y;
> n:=vectdim(x);
> for j from 0 to n-1 do 
>   w[j+1]:= evalf(cos( 2*Pi*j/n)- I* sin(2*Pi*j/n))
> od;
> y:=vector(n,0);
> for j from 0 to n-1 do
>   for k from 0 to n-1 do
>     r:=(k*j) mod n;
>     y[j+1]:=y[j+1]+x[k+1]*w[r+1]
>   od
> od;
> evalm(y)
> end:

        a := [1, 2, 3, 4, 5, 6, 7, 8, 8, 7, 6, 5, 4, 3, 2, 1]

> dft(a);

[72., -25.27414237 - 5.027339498 I, 0, -2.239828813 - 1.496605763 I,

    0, -.4464626915 - .6681786385 I, 0, -.0395661275 - .1989123665 I,

    0, -.0395661275 + .1989123665 I, 0, -.4464626915 + .6681786385 I,

    0, -2.239828813 + 1.496605763 I, 0, -25.27414237 + 5.027339498 I]

>