function G = recombine_AGB(A, Gamma, B)
% RECOMBINE_AGB  Reconstruct G = A*Gamma*B over F2[D] (cell-array polynomial matrices)
% J. Cioffi and ChatGPT - September 25, 2025
% function G = recombine_AGB(A, Gamma, B)
%
% Representation:
%   - Matrices are cell arrays.
%   - Each entry is a row vector of 0/1 coefficients [c0 c1 ...], low->high in D.
%   - [] denotes the zero polynomial.
%
% Usage:
%   G = recombine_AGB(A, Gamma, B);
%
% If A/Gamma/B are numeric (e.g., gf arrays), this falls back to ordinary *.
% ===============================================================================

    if ~(iscell(A) || iscell(Gamma) || iscell(B))
        % numeric/gf case
        G = A * (Gamma * B);
        return
    end

    % Normalize inputs to cell-array polynomial matrices
    A     = ensure_cell_pmat(A);
    Gamma = ensure_cell_pmat(Gamma);
    B     = ensure_cell_pmat(B);

    % Dimension checks
    [kA, rA] = size(A);
    [kG, nG] = size(Gamma);
    [kB, nB] = size(B);
    if rA ~= kG
        error('Dimension mismatch: size(A)=%dx%d, size(Gamma)=%dx%d.', kA,rA,kG,nG);
    end
    if nG ~= kB
        error('Dimension mismatch: size(Gamma)=%dx%d, size(B)=%dx%d.', kG,nG,kB,nB);
    end

    % Multiply over F2[D]: (A*Gamma)*B
    T = pmat_mul(A, Gamma);
    G = pmat_mul(T, B);
end

% ---------- helpers ----------

function Xc = ensure_cell_pmat(X)
% Ensure X is a cell-array polynomial matrix with entries as 0/1 row vectors.
    if iscell(X)
        Xc = X;
        for i = 1:numel(Xc)
            Xi = Xc{i};
            if isempty(Xi)
                Xc{i} = [];
            elseif isnumeric(Xi) || islogical(Xi)
                v = double(Xi(:).');              % row
                if numel(v) == 1
                    v = mod(v,2);
                    if v == 0, v = []; else, v = 1; end
                else
                    v = mod(v,2);
                end
                Xc{i} = trim_poly(v);
            else
                error('Unsupported cell entry class: %s', class(Xi));
            end
        end
    elseif isnumeric(X) || islogical(X)
        % Wrap numeric matrix as constant polynomials
        [m,n] = size(X);
        Xc = cell(m,n);
        for i = 1:m
            for j = 1:n
                v = mod(double(X(i,j)),2);
                if v == 0
                    Xc{i,j} = [];
                else
                    Xc{i,j} = 1; % constant "1"
                end
            end
        end
    else
        error('Unsupported matrix class: %s (expected cell or numeric).', class(X));
    end
end

function C = pmat_mul(A, B)
% Polynomial matrix multiply over F2[D].
    [m,r]  = size(A);
    [r2,n] = size(B);
    if r ~= r2
        error('Inner dimensions must agree for polynomial matrix multiply.');
    end
    C = cell(m,n);
    for i = 1:m
        for j = 1:n
            s = []; % zero polynomial
            for t = 1:r
                s = p_add(s, p_mul(A{i,t}, B{t,j}));
            end
            C{i,j} = s;
        end
    end
end

function r = p_add(p, q)
% Add over F2 (XOR), low->high vectors.
    if isempty(p), r = trim_poly(q); return; end
    if isempty(q), r = trim_poly(p); return; end
    lp = numel(p); lq = numel(q);
    L  = max(lp,lq);
    pp = zeros(1,L); qq = zeros(1,L);
    pp(1:lp) = p; qq(1:lq) = q;
    r = mod(pp + qq, 2);
    r = trim_poly(r);
end

function r = p_mul(p, q)
% Convolution over F2.
    if isempty(p) || isempty(q), r = []; return; end
    r = conv(double(p), double(q));
    r = mod(r, 2);
    r = trim_poly(r);
end

function r = trim_poly(p)
% Canonicalize: row vector, drop trailing zeros; [] for zero.
    p = double(p(:).');
    idx = find(p ~= 0, 1, 'last');
    if isempty(idx), r = []; else, r = p(1:idx); end
end
