N0 := 8;
N1 := N0 + 4;
N2 := N0;
N3 := N1;

# Symmetry conditions:

# Scaling Filters
n := 'n';
for n from 0 to N0/2-1 do
        h0(N0-1-n) := h0(n);
od;

n := 'n';
for n from 0 to N1/2-1 do
        h1(N1-1-n) := h1(n);
od;


# Wavelet Filters
n := 'n';
for n from 0 to N2/2-1 do
        h2(N0-1-n) := -h2(n);
od;

n := 'n';
for n from 0 to N3/2-1 do
        h3(N1-1-n) := -h3(n);
od;


# filter H0
n := 'n';
H0 := sum(h0(n)*z^n,n=0..(N0-1));

# filter H1
n := 'n';
H1 := sum(h1(n)*z^n,n=0..(N1-1));

# filter H2
n := 'n';
H2 := sum(h2(n)*z^n,n=0..(N2-1));

# filter H3
n := 'n';
H3 := sum(h3(n)*z^n,n=0..(N3-1));


# ---------- balancing conditions ----------

P := z^3+z^2+z+1;

V1 := 1;
V2 := (3-(z^4))/2;
V3 := (15-10*(z^4)+3*(z^8))/8;
V4 := (35 - 35*(z^4) + 21*(z^8) - 5*(z^12))/16;

R1 := expand(rem(H0+V1*H1,P^1,z));
R2 := expand(rem(H0+V2*H1,P^2,z));
R3 := expand(rem(H0+V3*H1,P^3,z));
R4 := expand(rem(H0+V4*H1,P^4,z));

interface(screenwidth=300);
writeto(`eqs.B`);

lprint(`// Balancing Conditions`);
k := 'k':
for k from 0 to 20 do
	lprint(coeff(R2,z,k),`,`);
od;

# ---------- orthogonality conditions ---------- 

H0H0 := expand(collect(H0*subs(z=1/z,H0),z)):
H1H1 := expand(collect(H1*subs(z=1/z,H1),z)):
H2H2 := expand(collect(H2*subs(z=1/z,H2),z)):
H3H3 := expand(collect(H3*subs(z=1/z,H3),z)):

H0H1 := expand(collect(H0*subs(z=1/z,H1),z)):
H0H2 := expand(collect(H0*subs(z=1/z,H2),z)):
H0H3 := expand(collect(H0*subs(z=1/z,H3),z)):
H1H2 := expand(collect(H1*subs(z=1/z,H2),z)):
H1H3 := expand(collect(H1*subs(z=1/z,H3),z)):
H2H3 := expand(collect(H2*subs(z=1/z,H3),z)):

lprint(`// Self Orthogonality Conditions `);
k := 'k':
for k from 1 to 15 do
        lprint(coeff(H0H0,z,4*k),`,`);
        lprint(coeff(H1H1,z,4*k),`,`);
	lprint(coeff(H2H2,z,4*k),`,`);
        lprint(coeff(H3H3,z,4*k),`,`);
od;

lprint(`// Cross Orthogonality Conditions `);
k := 'k':
for k from -15 to 15 do
	lprint(coeff(H0H1,z,4*k),`,`);
	lprint(coeff(H0H2,z,4*k),`,`);
	lprint(coeff(H0H3,z,4*k),`,`);
	lprint(coeff(H1H2,z,4*k),`,`);
	lprint(coeff(H1H3,z,4*k),`,`);
	lprint(coeff(H2H3,z,4*k),`,`);
od;

# ---------- normalization conditions ----------

nor0 := coeff(expand(H0H0),z,0) - 1/2:
nor1 := coeff(expand(H1H1),z,0) - 1/2:
nor2 := coeff(expand(H2H2),z,0) - 1/2:
nor3 := coeff(expand(H3H3),z,0) - 1/2:

norA := subs(z=1,expand(H0)) - 1:
norB := subs(z=1,expand(H1)) - 1:

lprint(`// Normalizations`);

lprint(nor0,`,`);
lprint(nor1,`,`);
lprint(nor2,`,`);
lprint(nor3,`,`);

lprint(norA,`,`);
lprint(norB,`;`);

        # note: put a semicolon at the end of the last equation

writeto(terminal);
interface(screenwidth=80);