| << Chapter < Page | Chapter >> Page > |
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:Kn = 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-1for 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-1y(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 thedimensions 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 : columnsf = 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 matrixif 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:Ka = 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-1A = x(i*p+j+1);
for k = 0:c:c*(p-2)A = A + x(i*p+j+k+c+1);
endy(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;
endend
endy = y/p;
The permutation of
Equation 18 from Preliminaries is implemented by
pfp . It calls the function
pfp2I .
The transpose is implemented by
pfpt and it
calls
pfpt2I .
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:2a = 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-1i1 = 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:Ks = 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*sn = 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:sy(i2 + i) = x(i1 + i);
endk1 = k1 + 1;
k2 = k2 + 1;if k1>= b
k1 = k1 - b;end
if k2>= a
k2 = k2 - a;end
end
The following Matlab programs implement Rader's permutation and its transpose.They require the primitive root to be passed to them as an argument.
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) == pa = 1;
y = zeros(p,1);y(1) = x(1);
for k = 2:py(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;
breakend
endR_inv = 'Not Found';
for x = 1:Nif rem(x*R,N) == 1
R_inv = x;break
endend
Notification Switch
Would you like to follow the 'Automatic generation of prime length fft programs' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|