function [F,G,b,K] = carlbilin(F1,G1,P,Mvars,Mpos,D,Dg,K)
%CARLBILIN computes the non-redundant Carleman Bilinearization of a
%linear-analytic system
%
% [F,G,B,K] = CARLBILIN(F1,G1,P) computes the Pth order
% Carleman Bilinearization v' = Fv + Gvu + bu of the 
% linear-analytic system x'= f(x) + g(x)u from its Taylor
% expansion coeficients F1 and G1. A matrix K, which depends only on 
% m=size(F1,1) and P, is also returned and can be reused.
%
% [F,G,B,K] = CARLBILIN(F1,G1,P,MVARS,MPOS,D,DG,K) uses, for faster
% execution, pre-calculated maps that depend only on m and P.

%   CBVK Toolbox
%   JHMG, PMSB 2016


m = size(F1,1);  
  
% Number of rows of each order p 
allterms = fliplr(calcnumterms(P,m));
MP = sum(allterms);

if nargin < 4
    % Pre-computing maps of positions corresponding to each product of xis
    [Mvars, Mpos, D] = mapsvarspos(m, P);
    Dg = D(1:end-allterms(end),1:end-allterms(end));
    Dg(Dg > size(Dg,1)) = 0;
end

% Initializing output structures
F = zeros(MP,size(F1,2));
F(1:m,:) = F1;
G = zeros(MP,size(G1,2));
G(1:m,:) = G1;
if nargin < 8
    K =zeros (length(F),2);
end

% Main loop of the Carleman bilinearization
for p=2:P
    lin0 = sum(allterms(1:p-1));
    mp = allterms(p);
    Fp = zeros(mp, size(F,2));
    Gp = zeros(mp, size(G,2));
    indsF = lin0+1:size(F,2);
    indsG = lin0-allterms(p-1)+1:size(G,2);
    for linp=1:mp
        lin = lin0+linp;
        vars = Mvars{lin+1};
        [v1,v2] = splithalf(vars);
        if (nargin < 8)
            posv1 = posrecur(Mpos, v1);
            posv2 = posrecur(Mpos, v2);
            K(lin0+linp+1,:)=[posv1 posv2];
        else
            posv1=K(lin0+linp+1,1);
            posv2=K(lin0+linp+1,2);
        end
        
        dFv1 = F(posv1,:);
        dFv2 = F(posv2,:);        
        Fp(linp,indsF) = prodtaylor(dFv1,posv2,D,indsF) + ...
            prodtaylor(dFv2,posv1,D,indsF);
        dGv1 = G(posv1,:);
        dGv2 = G(posv2,:);
        Gp(linp,indsG) = prodtaylor(dGv1,posv2,Dg,indsG) + ...
            prodtaylor(dGv2,posv1,Dg,indsG);
    end;
    F(lin0+1:lin0+mp,indsF) = Fp(:,indsF);
    G(lin0+1:lin0+mp,indsG) = Gp(:,indsG);
end;    

F = F(:,2:end);
G = [G(:,2:end) zeros(MP,size(F,2)-size(G,2)+1)];
b = [G1(:,1); zeros(MP-m,1)];

    