PROGRAM pasfft;

{
	Copied from Apple disks on 14 May, 1983. Original written in Sept,
    1980
}

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

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

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

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

CONST pi=3.1415926;

VAR  i,j,l,lim1,lim2,lim3,m,r:integer;
     arg,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 l:=1 TO lim2 DO
               FOR m:=1 TO lim1 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];
                     arg:=2.0*pi*((m-1)*lim2)/n;
                     a1:=c1*cos(arg)+c2*sin(arg);
                     a2:=-c1*sin(arg)+c2*cos(arg);
                     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;   (* easy  *)

BEGIN  (* main  *)
   FOR i:=1 TO max DO
      BEGIN
         x[1,i]:=10.0*sin(twopi*i/64.0);
         x[2,i]:=0.0;
      END;
   writeln('  array done');
   easy(x,max);
   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.