-- fft32.adb
with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types;
with Complex_Arrays; use Complex_Arrays;

procedure fft32(Z : in out Complex_Vector) is --  (0..31)
  pragma Suppress(All_Checks);
  W: Complex_Vector(0..31);  -- scratch vector, used many times
  E: Complex_Vector(0..16) :=   -- constants for FFT algorithm
   ((1.0, 0.0), (0.980785, 0.19509),
   (0.92388, 0.382683), (0.83147, 0.55557),
   (0.707107, 0.707107), (0.55557, 0.83147),
   (0.382683, 0.92388), (0.19509, 0.980785),
   (0.0, 1.0), (-0.19509, 0.980785),
   (-0.382683, 0.92388), (-0.55557, 0.83147),
   (-0.707107, 0.707107), (-0.83147, 0.55557),
   (-0.92388, 0.382683), (-0.980785, 0.19509),
   (-1.0, 0.0));
  I, J, K, L, M : Integer;
begin

   M := 16;
   L := 1;
   loop
      K := 0;
      J := L;
      I := 0;
      loop
         loop
            W(I+K) := Z(I) + Z(M+I);
            W(I+J) := E(K) * (Z(I) - Z(M+I));
            I := I+1;
            if I >= J then exit; end if;
         end loop;
         K := J;
         J := K+L;
         if J > M then exit; end if;
      end loop;
      L := L+L;
      if L > M then 
        for I in 0..2*M-1 loop
          Z(I) := W(I);
        end loop;
        return;
      end if;  -- result is in Z
                -- work back other way without copying
      K := 0;
      J := L;
      I := 0;
      loop
        loop
          Z(I+K) := W(I) + W(M+I);
          Z(I+J) := E(K) * (W(I) - W(M+I));
          I := I+1;
          if I >= J then exit; end if;
        end loop;
        K := J;
        J := K+L;
        if J > M then exit; end if;
      end loop;
      L := L+L;
      if L > M then return; end if;  -- result is in Z
   end loop;
end fft32;