PROGRAM pasfft;

{
	Copied from Apple disk 14 May, 1983 .. orig. written in Sept, 1980
}

CONST  twopi=6.283185;
       max=256;
       maxless=255;
       
TYPE data=ARRAY[1..2,1..max] OF real;
     table=ARRAY[0..maxless] OF real;

VAR  x:data;
     c:real;
     i:integer;
     costab:table;
     
(*----------------------------------------------------------------------*)

PROCEDURE easy(VAR a:data;n:integer;costab:table);

(*   fft procedure with a two dimensional data array and n is the size of
the array  *)

CONST pi=3.1415926;

VAR  arg,i,j,l,lim1,lim2,lim3,m,r:integer;
     a1,a2,b1,b2,c1,c2,cs,cs1,cstep,si,si1,sstep:real;
     
BEGIN
   r:=round(ln(n)/ln(2.0));
   lim1:=n-1;
   lim2:=n DIV 2;
   j:=1;                          (* begin shuffle *)
   FOR i:=1 TO lim1 DO
      BEGIN
         IF i<j THEN
            BEGIN
               a1:=a[1,j];
               a[1,j]:=a[1,i];
               a[2,j]:=0.0;
               a[1,i]:=a1;
               a[2,i]:=0.0;
            END;
         l:=lim2;
         WHILE l<j DO
            BEGIN
               j:=j-l;
               l:=l DIV 2;
            END;
         j:=j+l;
      END;               (* shuffle ends, decimation begins  *)
   FOR i:=1 TO r DO
      BEGIN
         lim1:=round(exp((i-1)*ln(2.0)));
         lim2:=round(exp((r-i)*ln(2.0)));
         FOR m:=1 TO lim1 DO
            BEGIN
               arg:=((m-1)*lim2) MOD max;
               j:=(max DIV 4)-arg;
               IF j<0 THEN j:=j+max;
               si:=costab[j];
               cs:=costab[arg];
               FOR l:=1 TO lim2 DO
                  BEGIN
                     lim3:=1+(m-1)+(l-1)*2*lim1;
                     b1:=a[1,lim3];
                     b2:=a[2,lim3];
                     c1:=a[1,lim3+lim1];
                     c2:=a[2,lim3+lim1];
                     a1:=c1*cs+c2*si;
                     a2:=-c1*si+c2*cs;
                     a[1,lim3]:=b1+a1;
                     a[2,lim3]:=b2+a2;
                     a[1,lim3+lim1]:=b1-a1;
                     a[2,lim3+lim1]:=b2-a2;
                  END;
            END;
      END;
END;   (* easy  *)

BEGIN  (* main  *)
   FOR i:=1 TO max DO
      BEGIN
         costab[i-1]:=cos(twopi*(i-1)/max);
         x[1,i]:=10.0*sin(twopi*4.02/i);
         x[2,i]:=0.0;
      END;
   writeln('  array done');
   easy(x,max,costab);
   FOR i:=1 TO max DO
      BEGIN
         c:=sqrt(sqr(x[1,i])+sqr(x[2,i]));
         writeln(i:10,c:10:2);
      END;
END.