program pasfft;

{
	Copied from Apple disks 14 May, 1983 .. original written 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)));
         arg:=2.0*pi*lim2/n;
         cs:=1.0;
         si:=0.0;
         cstep:=cos(arg);
         sstep:=sin(arg);
         for m:=1 to lim1 do
            begin
               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;
               cs1:=cs*cstep-si*sstep;
               si1:=si*cstep+cs*sstep;
               cs:=cs1;
               si:=si1;
            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.