function [u,ip,op,ADDS,MULTS] = ff(p,e);
% [u,ip,op,ADDS,MULTS] = ff(p,e);
% u  : multiplicative constants
% ip : input permutation
% op : output permutation

K = length(p);
N = prod(p.^e);
P = N + 1;
[pr, ipr] = primitive_root(P);
Red_Adds = 2 * N * (K - sum(1./(p.^e)) );
ADDS = 2 * Red_Adds;

FS = sprintf('fft%d.m',P);
fid = fopen(FS,'w');

fprintf(fid,'function y = fft%d(x,u,ip,op)\n',P);
fprintf(fid,'%% y = fft%d(x,u,ip,op)\n',P);
fprintf(fid,'%% y  : the %d point DFT of x \n',P);
fprintf(fid,'%% u  : a vector of precomputed multiplicative constants\n');
fprintf(fid,'%% ip : input permutation\n');
fprintf(fid,'%% op : ouput permutation\n');

Pstr = sprintf('[%d',p(1));
for k = 2:K, Pstr = [Pstr, sprintf(',%d',p(k))]; end
Pstr = [Pstr,']'];
Estr = sprintf('[%d',e(1));
for k = 2:K, Estr = [Estr, sprintf(',%d',e(k))]; end
Estr = [Estr,']'];
PEstr = sprintf('[%d',p(1)^e(1));
for k = 2:K, PEstr = [PEstr, sprintf(',%d',p(k)^e(k))]; end
PEstr = [PEstr,']'];

fprintf(fid,'\n');
S = sprintf('y = zeros(%d,1);\n',P);
fprintf(fid,S);
S1 = sprintf('x = x(ip);');
S2 = sprintf('%% input permutation\n');
fprintf(fid,'%-50s%s',S1,S2);
S1 = sprintf(['x(2:%d) = KRED(',Pstr,',',Estr,',%d,x(2:%d));'],P,K,P);
S2 = sprintf('%% reduction operations\n');
fprintf(fid,'%-50s%s',S1,S2);

e_table = [0:e(1)]';
a = e(1)+1;
for i = 2:K
   e_table = [kron(ones(e(i)+1,1),e_table), kron([0:e(i)]',ones(a,1))];
   a = a * (e(i)+1);
end
R = prod(e+1);


% ------------------------ MULTIPLICATIVE CONSTANTS ------------------------
k = rp(P,ipr,0:N);
I = sqrt(-1);
W = exp(-I*2*pi*k/P);
h = W(2:P);
h = h(N:-1:1);
h = pfp(p.^e,K,h);
h = itKRED(p,e,K,h);
u = h(1);

S1 = sprintf('y(1) = x(1)+x(2);');
S2 = sprintf('%% DC term calculation\n');
fprintf(fid,'%-50s%s',S1,S2);

DC_ADDS = 2;
ADDS = ADDS + DC_ADDS;

SLINE = '--------------------------------------------------------------------------------';
SB = ' block : 1 ';
SC = SLINE;
BL = 21;
SC(BL:BL-1+length(SB)) = SB;
fprintf(fid,'%% %s\n',SC);
S1 = sprintf('y(2) = x(2)*u(1);');
fprintf(fid,'%-40s\n',S1);
a = 1;
MULTS = 1;
for i = 2:R
   v = e_table(i,:);
   f = find(v>0);
   q = p(f);
   t = v(f);
   L = prod(q-1)*prod(q.^(t-1));

   B = prod(q.^t);
   bs = sprintf('%d',q(1)^t(1));
   for k = 2:length(q), bs = [bs, sprintf(' * %d',q(k)^t(k))]; end
   if length(q) > 1
	SB = sprintf(' block : %d = %s ',B,bs);
	SC = SLINE;
	SC(BL:BL-1+length(SB)) = SB;
	fprintf(fid,'%% %s\n',SC);
   else
	SB = sprintf(' block : %d ',B);
        SC = SLINE;
        SC(BL:BL-1+length(SB)) = SB;
        fprintf(fid,'%% %s\n',SC);
   end
   if prod(q.^t) == 2
      S1 = sprintf('y(%d) = x(%d)*u(%d);',a+2,a+2,MULTS+1);
      fprintf(fid,'%-40s\n',S1);
      Mk = 1;
   else
      d = []; r = []; c = []; Q = []; Qt = [];
      for j = 1:length(q)
         [dk,rk,ck,Qk,Qtk] = A_data(q(j)^t(j));
         if dk > 1
            d = [d dk]; r = [r rk]; c = [c ck]; Q = [Q Qk]; Qt = [Qt Qtk];
         end
      end
      [g,C1] = cgc(Q,r,c,length(Q));
      ADDS = ADDS + C1;
      Mk = prod(r);
      BEG = int2str(a+2); FIN = int2str(a+1+L);
      XX = ['x(',BEG,':',FIN,')']; YY = 'v';
      kpi(d,g,r,c,length(Q),YY,XX,fid);
      S1 = ['v = v.*u(',int2str(MULTS+1),':',int2str(MULTS+Mk),');'];
      fprintf(fid,'%-40s\n',S1);
      [g,C2] = cgc(Qt,c,r,length(Q));
      ADDS = ADDS + C2;
      XX = 'v'; YY = ['y(',BEG,':',FIN,')'];
      kpit(d,g,c,r,length(Q),YY,XX,fid);
   end

   c = [];
   r = [];
   lq = length(q);
   for j = 1:lq
      [fk,rk,ck] = C_data(q(j),t(j));
      r = [r rk]; c = [c ck];
   end
   f = (q-1).*(q.^(t-1));
   temp = Kcrot(q,t,lq,h(a+1:a+L));
   temp = KFt(f,r,c,temp);
   u = [u; temp(:)];
   a = a + L;
   MULTS = MULTS + Mk;
end

u(1) = u(1)-1;
fprintf(fid,'%% %s\n',SLINE);
S1 = sprintf('y(2) = y(1)+y(2);');
S2 = sprintf('%% DC term calculation\n');
fprintf(fid,'%-50s%s',S1,S2);
S1 = sprintf(['y(2:%d) = tKRED(',Pstr,',',Estr,',%d,y(2:%d));'],P,K,P);
S2 = sprintf('%% transpose reduction operations\n');
fprintf(fid,'%-50s%s',S1,S2);
S1 = sprintf('y = y(op);');
S2 = sprintf('%% output permutation\n');
fprintf(fid,'%-50s%s',S1,S2);
fprintf(fid,'\n');

MULTS = 2 * MULTS;
ADDS = 2* ADDS;
fprintf(fid,'%% For complex data - \n');
fprintf(fid,'%% Total Number of Real Multiplications : %d\n',MULTS);
fprintf(fid,'%% Total Number of Real Additions: %d\n\n',ADDS);
fclose(fid);

%%%%%%%%%%%%%%%%%%%% COMPUTE INPUT AND OUTPUT PERMUTATIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%

id = 1:P;   % identity permutation
ip = rp(P,pr,id);
ip(2:P) = pfp(p.^e,K,ip(2:P));

op = id;
op(2:P) = pfpt(p.^e,K,op(2:P));
op(2:P) = op(P:-1:2);
op = rpt(P,ipr,op);

%%%%%%%%%%%%%%%%% PUT MULTIPLICATIVE CONSTANTS AND PERMUTATIONS IN A FILE %%%%%%%%%%%%%%

CFS = sprintf('cap%d.m',P);
fid = fopen(CFS,'w');
fprintf(fid,'\n%% The multiplicative constants for the %d point FFT\n\n',P);
fprintf(fid,'I = sqrt(-1);\n');
fprintf(fid,'u = [\n');
for k = 1:MULTS/2
   if abs(real(u(k))) < 0.000001
      fprintf(fid,'%25.15f*I\n',imag(u(k)));
   elseif abs(imag(u(k))) < 0.00001
      fprintf(fid,'%25.15f\n',real(u(k)));
   else
      fprintf(fid,'%25.15f + %25.15f*I\n',real(u(k)),imag(u(k)));
   end
end
fprintf(fid,'];\n\n');
fprintf(fid,'\n%% The input permutation for the %d point FFT\n\n',P);
fprintf(fid,'ip = [\n');
for k = 1:P
	fprintf(fid,'   %d\n',ip(k));
end
fprintf(fid,'];\n\n');
fprintf(fid,'\n%% The output permutation for the %d point FFT\n\n',P);
fprintf(fid,'op = [\n');
for k = 1:P
        fprintf(fid,'   %d\n',op(k));
end
fprintf(fid,'];\n\n');
fclose(fid);

The following programs print the program statements that carry out the operation I⊗D_k⊗I and I⊗D_k^t⊗I. They are modeled after kpi in the text.

function kpi(d,g,r,c,n,Y,X,fid)
% kpi(d,g,r,c,n,Y,X,fid);
% Kronecker Product : A(d(1)) kron ... kron A(d(n))
% g : permutation of 1,...,n 
% r : [r(1),...,r(n)]
% c : [c(1),..,c(n)]
% r(i) : rows of A(d(i))
% c(i) : columns of A(d(i))
% n : number of terms
for i = 1:n
   a = 1;
   for k = 1:(g(i)-1)
      if i > find(g==k) 
         a = a * r(k);
      else
         a = a * c(k);
      end
   end
   b = 1;
   for k = (g(i)+1):n
      if i > find(g==k)
         b = b * r(k);
      else
         b = b * c(k);
      end
   end
   % Y = (I(a) kron A(d(g(i))) kron I(b)) * X;
   if i == 1
      S1 = sprintf([Y,' = ID%dI(%d,%d,',X,');    '],d(g(i)),a,b);
      S2 = sprintf(['%% ',Y,' = (I(%d) kron D%d kron I(%d)) * ',X],a,d(g(i)),b);
      fprintf(fid,'%-35s%s\n',S1,S2);
   elseif d(g(i)) ~= 1
      S1 = sprintf([Y,' = ID%dI(%d,%d,',Y,');    '],d(g(i)),a,b);
      S2 = sprintf(['%% ',Y,' = (I(%d) kron D%d kron I(%d)) * ',Y],a,d(g(i)),b);
      fprintf(fid,'%-35s%s\n',S1,S2);
   end
end

function kpit(d,g,r,c,n,Y,X,fid)
% kpit(g,r,c,n,Y,X,fid);
% (transpose)
% Kronecker Product : A(d(1))' kron ... kron A(d(n))'
% g : permutation of 1,...,n 
% r : [r(1),...,r(n)]
% c : [c(1),..,c(n)]
% r(i) : rows of A(d(i))'
% c(i) : columns of A(d(i))'
% n : number of terms

for i = 1:n
   a = 1;
   for k = 1:(g(i)-1)
      if i > find(g==k)
         a = a * r(k);
      else
         a = a * c(k);
      end
   end
   b = 1;
   for k = (g(i)+1):n
      if i > find(g==k)
         b = b * r(k);
      else
         b = b * c(k);
      end
   end
   % x = (I(a) kron A(d(g(i)))'' kron I(b)) * x;
   if i == n
      S1 = sprintf([Y,' = ID%dtI(%d,%d,',X,');    '],d(g(i)),a,b);
      S2 = sprintf(['%% ',Y,' = (I(%d) kron D%d'' kron I(%d)) * ',X],a,d(g(i)),b);
      fprintf(fid,'%-35s%s\n',S1,S2);
   elseif d(g(i)) ~= 1
      S1 = sprintf([X,' = ID%dtI(%d,%d,',X,');    '],d(g(i)),a,b);
      S2 = sprintf(['%% ',X,' = (I(%d) kron D%d'' kron I(%d)) * ',X],a,d(g(i)),b);
      fprintf(fid,'%-35s%s\n',S1,S2);
   end
end

function u = KFt(f,r,c,u)
% u = (F^t kron ... kron F^t)*u
% (transpose)
% f = [f(1),...,f(K)]
% r : r(i) = rows of F(i)
% c : c(i) = columns of F(i)
% u : length(u) = prod(c);
K = length(f);
for i = 1:K
   m = prod(c(1:i-1));
   n = prod(r(i+1:K));
   u = IFtI(f(i),r(i),c(i),m,n,u);
end

function y = IFtI(s,r,c,m,n,x);
% y = (I(m) kron F(s)^t kron I(n))*x
% (transpose)
% r : rows of F(s)
% c : columns of F(s)
v = 0:n:n*(c-1);
u = 0:n:n*(r-1);
for i = 0:m-1
   for j = 0:n-1
      y(v+i*c*n+j+1) = Ftop(s,x(u+i*r*n+j+1));
   end
end

function y = Ftop(k,x)
if k == 1, y = x;
elseif k == 2, y = F2t(x);
elseif k == 3, y = F3t(x);
elseif k == 4, y = F4t(x);
elseif k == 6, y = F6t(x);
elseif k == 8, y = F8t(x);
elseif k == 18, y = F18t(x);
end

function x = Kcrot(p,e,K,x)
% Kronecker product of Cyclotomic Reduction Operations.
% x = (G(p(1)^e(1)) kron ... kron G(p(K)^(K)))^t*x
% (transpose)
% p : p = [p(1),...,p(K)];
% e : e = [e(1),...,e(K)];
a = (p-1).*((p).^(e-1));
r = a;		% r(i) = number of rows of G(i)
c = 2*a-1;	% c(i) = number of columns of G(i)
m = 1;
n = prod(r);
for i = 1:K
   n = n / r(i);
   x = IcrotI(p(i),e(i),m,n,x);
   m = m * c(i);
end

function y = IcrotI(p,e,m,n,x)
%  y = (eye(m) kron G(p^e)^t kron eye(n))*x
%  (transpose)
a = (p-1)*(p^(e-1));
c = a;
r = 2*a-1;
y = zeros(r*m*n,1);
v = 0:n:(r-1)*n;
u = 0:n:(c-1)*n;
for i = 0:m-1
   for j = 0:n-1
      y(v+i*r*n+j+1) = crot(p,e,x(u+i*c*n+j+1));
   end
end

function y = crot(p,e,x)
% y = crot(p,x)
% cyclotomic reduction matrix (transpose)
% length(x) == 2*n-1
% length(y) == n
% where n = (p-1)*(p^(e-1))
n = (p-1)*(p^(e-1));
y = zeros(2*n-1,1); 
if p == 2
  n = p^(e-1);
  y(1:n) = x;
  y(n+1:2*n-1) = -x(1:n-1);
else
  y(1:n) = x;
  L = p^(e-1);
  y(n+1:n+L) = -x(1:L);
  a = L;
  for k = 2:p-1
     y(n+1:n+L) = y(n+1:n+L) - x(a+1:a+L);
     a = a + L;
  end
  b = 2*n-1 - p*(p^(e-1));
  y(p*L+1:p*L+b) = x(1:b);
end

function [d,r,c,Q,Qt] = A_data(n)
% A : A matrix in bilinear form for cyclotomic convolution
% d : linear convolution modules used
% r : rows
% c : columns
% Q : Q(i) = cost associated with D(d(i))
% Qt : Qt(i) = cost associated with D(d(i))'
if n == 2, d = [1];
elseif n == 4, d = [2];
elseif n == 8, d = [2 2];
elseif n == 16, d = [2 2 2];
elseif n == 3, d = [2];
elseif n == 9, d = [2 3];
elseif n == 27, d = [2 3 3];
elseif n == 5, d = [2 2];
elseif n == 7, d = [2 3];
end
r = []; c = []; Q = []; Qt = [];
for k = 1:length(d)
   [rk, ck, Qk, Qtk] = D_data(d(k));
   r = [r rk]; c = [c ck]; Q = [Q Qk]; Qt = [Qt Qtk];
end

function [r,c,Q,Qt] = D_data(d);
% D : D matrix in bilinear form for linear convolution
% r : rows
% c : columns
% Q : cost associated with D(d)
% Qt : cost associated with D(d)'
if d == 1, r = 1; c = 1; Q = 0; Qt = 0;
elseif d == 2, r = 3; c = 2; Q = 1; Qt = 2;
elseif d == 3, r = 5; c = 3; Q = 7; Qt = 9;
end

function [f,r,c] = C_data(p,e)
% f : length of linear convolution
% r : rows
% c : columns
f = prod((p-1).*(p.^(e-1)));
% (Euler Totient Function)
r = 2*f-1;
c = F_data(f);

function c = F_data(n)
% c : columns of F matrix
if n == 1, c = 1;
elseif n == 2, c = 3;
elseif n == 4, c = 9;
elseif n == 8, c = 27;
elseif n == 3, c = 5;
elseif n == 6, c = 15;
elseif n == 18, c = 75;
end

function x = itKRED(P,E,K,x)
% x = itKRED(P,E,K,x);
% (inverse transpose)
% P : P = [P(1),...,P(K)];
% E : E = [E(K),...,E(K)];
for i = 1:K
   a = prod(P(1:i-1).^E(1:i-1));
   c = prod(P(i+1:K).^E(i+1:K));
   p = P(i);
   e = E(i);
   for j = e-1:-1:0
      x(1:a*c*(p^(j+1))) = itRED(p,a,c*(p^j),x(1:a*c*(p^(j+1))));
   end
end

function y = itRED(p,a,c,x)
% y = itRED(p,a,c,x);
% (inverse transpose)
y = zeros(a*c*p,1);
for i = 0:c:(a-1)*c
   for j = 0:c-1
      A = x(i*p+j+1);
      for k = 0:c:c*(p-2)
         A = A + x(i*p+j+k+c+1);
      end
      y(i+j+1) = A;
      for k = 0:c:c*(p-2)
         y(i*(p-1)+j+k+a*c+1) = p*x(i*p+j+k+1) - A;
      end
   end
end
y = y/p;

function x = pfp(n,K,x)
% x = P(n(1),...,n(K)) * x
% n = [n(1),...,n(K)];
% length(x) = prod(n(1),...,n(K))
a = prod(n);
s = 1;
for i = K:-1:2
  a = a / n(i);
  x = pfp2I(a,n(i),s,x);
  s = s * n(i);
end

function y = pfp2I(a,b,s,x)
% y = kron(P(a,b),I(s)) * x;
% length(x) = a*b*s
n = a * b;
y = zeros(n*s,1);
k1 = 0;
k2 = 0;
for k = 0:n-1
  i1 = s * (k1 + b * k2);
  i2 = s * k;
  for i = 1:s
    y(i1 + i) = x(i2 + i);
  end
  k1 = k1 + 1;
  k2 = k2 + 1;
  if k1 >= b
    k1 = k1 - b;
  end
  if k2 >= a
    k2 = k2 - a;
  end
end

function x = pfpt(n,K,x)
% x = P(n(1),...,n(K))' * x
% (tanspose)
% n = [n(1),...,n(K)];
% length(x) = prod(n(1),...,n(K))
% a = prod(n);
a = n(1);
s = prod(n(2:K));
for i = 2:K
  s = s / n(i);
  x = pfpt2I(a,n(i),s,x);
  a = a * n(i);
end

function y = pfpt2I(a,b,s,x)
% y = P(a,b)' kron I(s) * x;
% (transpose)
% length(x) = a*b*s
n = a * b;
y = zeros(n*s,1);
k1 = 0;
k2 = 0;
for k = 0:n-1
  i1 = s * (k1 + b * k2);
  i2 = s * k;
  for i = 1:s
    y(i2 + i) = x(i1 + i);
  end
  k1 = k1 + 1;
  k2 = k2 + 1;
  if k1 >= b
    k1 = k1 - b;
  end
  if k2 >= a
    k2 = k2 - a;
  end
end

function y = rp(p,r,x)
% Rader's Permutation
% p : prime
% r : a primitive root of p
% x : length(x) == p
a = 1;
y = zeros(p,1);
y(1) = x(1);
for k = 2:p
   y(k) = x(a+1);
   a = rem(a*r,p);
end

function y = rpt(p,r,x)
% Rader's Permutation
% (transpose)
% p : prime
% r : a primitive root of p
% x : length(x) == p
a = 1;
y = zeros(p,1);
y(1) = x(1);
for k = 2:p
   y(a+1) = x(k);
   a = rem(a*r,p);
end

function [R, R_inv] = primitive_root(N)

%  function [R, R_inv] = primitive_root(N)
%  Ivan Selesnick
%  N is assumed to be prime.  This function returns R,
%  the smallest	primitive root of N, and R_inv, the 
%  inverse of R modulo N.

R = 'Not Found';
m = 0:(N-2);
for x = 1:(N-1)
   if ( 1:(N-1) == sort(rem2(x,m,N)) )
      R = x;
      break
   end
end
R_inv = 'Not Found';
for x = 1:N
   if rem(x*R,N) == 1
      R_inv = x;
      break
   end
end