The reduction matrix of Equation 44 in Preliminaries is implemented
by `KRED`

which calls `RED`

.
Its transpose and inverse transpose are implemented by
`tRED`

, `tRED`

, `itKRED`

and `itRED`

.

```
function x = KRED(P,E,K,x)
% x = KRED(P,E,K,x);
% 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))) = RED(p,a,c*(p^j),x(1:a*c*(p^(j+1))));
end
end
```

```
function y = RED(p,a,c,x)
% y = RED(p,a,c,x);
y = zeros(a*c*p,1);
for i = 0:c:(a-1)*c
for j = 0:c-1
y(i+j+1) = x(i*p+j+1);
for k = 0:c:c*(p-2)
y(i+j+1) = y(i+j+1) + x(i*p+j+k+c+1);
y(i*(p-1)+j+k+a*c+1) = x(i*p+j+k+1) - x(i*p+j+c*(p-1)+1);
end
end
end
```

```
function x = tKRED(P,E,K,x)
% x = tKRED(P,E,K,x);
% (transpose)
% P : P = [P(1),...,P(K)];
% E : E = [E(K),...,E(K)];
for i = K:-1:1
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 = 0:e-1
x(1:a*c*(p^(j+1))) = tRED(p,a,c*(p^j),x(1:a*c*(p^(j+1))));
end
end
```

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