The following programs carry out the operation of
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
The following programs carry out the operation of
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
The following programs tell the programs for code generation relevant information about the bilinear forms for cyclotomic convolution. Specifically, they indicates the linear convolution out of which these cyclotomic convolution are composed, and the dimensions of the corresponding matrices. See the appendix Bilinear Forms for Linear Convolution.
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
Programs for Inverse Transpose Reduction Operations
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;
