function [xs, Vr, Aeq, beq] = linear_constr2(f,Np,Nu)
% linear constrains on the both the Neville filters
% Aeq*[p' u'] = beq
% general solution is: x = xs + Vr*phi, where phi is an arbitrary vector
% Copyright (c) 2006 Yi Chen

format short
[d1, p_coeff] = poly2matrix_sym(f(1));
[d2, u_coeff] = poly2matrix_sym(f(2));

if Np == 0 && Nu == 0
    Aeq = [0 0 0 0]; 
    beq = 0;    
    if [size(p_coeff) size(u_coeff)] == [4 4 2 2]
        Aeq = [0 0 0 0 1 0 0 0 0 0;
            0 0 0 0 0 0 0 1 0 0];   
        beq = [0; 0];   
    end
    if [size(p_coeff) size(u_coeff)] == [4 4 4 4]
        Aeq = [0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]; 
        beq = [0; 0; 0; 0];    
    end
    if [size(p_coeff) size(u_coeff)] == [6 6 2 2]
        A_temp = zeros(6, 20);
        A_temp(1,7) = 1;
        A_temp(2,12) = 1;
        A_temp(3,13) = 1;
        A_temp(4,14) = 1;
        A_temp(5,17) = 1;
        A_temp(6,18) = 1;
        Aeq = [A_temp];
        beq = [0; 0; 0; 0; 0; 0];
    end
    if [size(p_coeff) size(u_coeff)] == [6 6 4 4]
        A_temp = zeros(8, 26);
        A_temp(1,7) = 1;
        A_temp(2,12) = 1;
        A_temp(3,13) = 1;
        A_temp(4,14) = 1;
        A_temp(5,17) = 1;
        A_temp(6,18) = 1;
        A_temp(7,23) = 1;
        A_temp(8,26) = 1;
        Aeq = A_temp;
        beq = zeros(8,1);
    end
    if [size(p_coeff) size(u_coeff)] == [6 6 6 6]
        A_temp = zeros(12, 36);
        A_temp(1,7) = 1;
        A_temp(2,12) = 1;
        A_temp(3,13) = 1;
        A_temp(4,14) = 1;
        A_temp(5,17) = 1;
        A_temp(6,18) = 1;
        A_temp(7,7+18) = 1;
        A_temp(8,12+18) = 1;
        A_temp(9,13+18) = 1;
        A_temp(10,14+18) = 1;
        A_temp(11,17+18) = 1;
        A_temp(12,18+18) = 1;
        Aeq = A_temp;
        beq = zeros(12, 1);
    end
else    
    %     linear constraints
    %     predict filter    
    sizep = size(p_coeff);
    [cp, bp] = solv_quin(sizep(2)/2, sizep(1)/2, Np);
    
    %     update filter    
    sizeu = size(u_coeff);
    [cu, bu] = solv_quin(sizeu(2)/2, sizeu(1)/2, Nu);
    
    %     combine them together
    cp1 = [cp zeros(length(cp(:,1)),length(cu(1,:)))];
    cu1 = [zeros(length(cu(:,1)),length(cp(1,:))) cu];
    %     linear equality
    beq = [-bp; 0.5*bu];
    Aeq = [cp1; cu1];
    
    %     for those have zeros in the lifting filters
    [d, p] = poly2matrix_sym(f(1));
    [d, u] = poly2matrix_sym(f(2));
    pr = reshape(p, [1, prod(size(p))]);
    pr = pr(length(pr)/2+1:length(pr));
    ur = reshape(u, [1, prod(size(u))]);
    ur = ur(length(ur)/2+1:length(ur));
    Aeq = Aeq(:, [find(pr ~= '0'), find(ur ~= '0') + prod(size(p))/2]);
    
    if Np == 2 && Nu == 2
        Aeq = Aeq([1, 4],:);
        beq = beq([1, 4]);
    elseif Np == 4 && Nu == 2
        Aeq = Aeq([4,5,6,10,11],:);
        beq = beq([4,5,6,10,11]);
    elseif Np == 4 && Nu == 4
        Aeq = Aeq([4,5,6,10,14 15 16 20],:);
        beq = beq([4,5,6,10,14 15 16 20]);
    elseif Np == 6 && Nu == 2
        Aeq = Aeq([1 4 6 8 11 14 16 20 21 22],:);
        beq = beq([1 4 6 8 11 14 16 20 21 22]);    
    elseif Np == 6 && Nu == 4
        Aeq = Aeq([1 4 6 8 11 14 16 20 21 25 26 27 31],:);
        beq = beq([1 4 6 8 11 14 16 20 21 25 26 27 31]);
    elseif Np == 6 && Nu == 6
        Aeq = Aeq([1 4 6 8 11 14 16 20 21 22 25 27 29 32 35 37 41 42],:);
        beq = beq([1 4 6 8 11 14 16 20 21 22 25 27 29 32 35 37 41 42]);        
    end
    
    if [size(p_coeff) size(u_coeff)] == [4 4 2 2]
        Aeq = [Aeq;
            0 0 0 0 1 0 0 0 0 0;
            0 0 0 0 0 0 0 1 0 0]; 
        beq = [beq; 0; 0];
    end
    if [size(p_coeff) size(u_coeff)] == [4 4 4 4]
        Aeq = [Aeq;
            0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]; 
        beq = [beq; 0; 0; 0; 0];
    end
    if [size(p_coeff) size(u_coeff)] == [6 6 2 2]
        A_temp = zeros(6, 20);
        A_temp(1,7) = 1;
        A_temp(2,12) = 1;
        A_temp(3,13) = 1;
        A_temp(4,14) = 1;
        A_temp(5,17) = 1;
        A_temp(6,18) = 1;
        Aeq = [Aeq; A_temp];
        beq = [beq; 0; 0; 0; 0; 0; 0];
    end
    if [size(p_coeff) size(u_coeff)] == [6 6 4 4]
        A_temp = zeros(8, 26);
        A_temp(1,7) = 1;
        A_temp(2,12) = 1;
        A_temp(3,13) = 1;
        A_temp(4,14) = 1;
        A_temp(5,17) = 1;
        A_temp(6,18) = 1;
        A_temp(7,23) = 1;
        A_temp(8,26) = 1;
        Aeq = [Aeq; A_temp];
        beq = [beq; zeros(8,1)];
    end
    if [size(p_coeff) size(u_coeff)] == [6 6 6 6]
        A_temp = zeros(12, 36);
        A_temp(1,7) = 1;
        A_temp(2,12) = 1;
        A_temp(3,13) = 1;
        A_temp(4,14) = 1;
        A_temp(5,17) = 1;
        A_temp(6,18) = 1;
        A_temp(7,7+18) = 1;
        A_temp(8,12+18) = 1;
        A_temp(9,13+18) = 1;
        A_temp(10,14+18) = 1;
        A_temp(11,17+18) = 1;
        A_temp(12,18+18) = 1;
        Aeq = [Aeq; A_temp];
        beq = [beq; zeros(12, 1)];
    end
    if size(p_coeff) == [8 8]
        A_tempp = zeros(10, 32);
        A_tempp(1,9) = 1;
        A_tempp(2,16) = 1;
        A_tempp(3,17) = 1;
        A_tempp(4,24) = 1;
        A_tempp(5,25) = 1;
        A_tempp(6,26) = 1;
        A_tempp(7,27) = 1;
        A_tempp(8,30) = 1;
        A_tempp(9,31) = 1;
        A_tempp(10,32) = 1;
        if size(u_coeff) == [2 2]
            A_temp = zeros(10,34);
            A_temp(1:10, 1:32) = A_tempp;
        end
        if size(u_coeff) == [4 4]
            A_temp = zeros(12,40);
            A_temp(11,37) = 1;
            A_temp(12,40) = 1;
            A_temp(1:10, 1:32) = A_tempp;
        end
        if size(u_coeff) == [6 6]
            A_temp = zeros(16,50);
            A_temp(11,7+32) = 1;
            A_temp(12,12+32) = 1;
            A_temp(13,13+32) = 1;
            A_temp(14,14+32) = 1;
            A_temp(15,17+32) = 1;
            A_temp(16,18+32) = 1;
            A_temp(1:10, 1:32) = A_tempp;
        end
        if size(u_coeff) == [8 8]
            A_temp = zeros(20,64);
            A_temp(1:10, 1:32) = A_tempp;
            A_temp(11:20, 33:64) = A_tempp;
        end
        Aeq = [Aeq; A_temp];
        beq = [beq; zeros(size(A_temp,1),1)];
    end
end

r = rank(Aeq);
[U,S,V] = svd(Aeq);
% last n-r columns of V
Vr = V(:, r+1:length(V(1,:)));

xs = pinv(Aeq)*beq;