function [Gmin, U, mu, rowdeg, leadpos, mu_check, A_out] = minimal_encoder_from_B(B, k, A_in, opts)
% MINIMAL_ENCODER_FROM_B  Build a minimal basic encoder from the top-k rows of B,
% and compose the resulting row operations into A: A_out = U * A_in.
% J. Cioffi and ChatGPT, September 24, 2025
%
%   [Gmin,U,mu,rowdeg,leadpos,mu_check,A_out] = minimal_encoder_from_B(B,k,A_in,opts)
%
% Inputs:
%   B      : n-by-n polynomial matrix over GF(2)[D] (low->high coeffs).
%   k      : number of encoder rows (take top k rows of B).
%   A_in   : k-by-k polynomial matrix (left undo operator from Smith).
%   opts   : struct with optional fields:
%            .verify_mu (false), .max_combs (3e5), .verbose (false)
%
% Outputs:
%   Gmin   : k-by-n minimal basic encoder (row-reduced).
%   U      : k-by-k unimodular row op matrix with Gmin = U * (topkrows(B)).
%   mu     : overall constraint length = sum(rowdeg).
%   rowdeg : 1-by-k row degrees of Gmin.
%   leadpos: 1-by-k leading column positions of Gmin.
%   mu_check : [] unless verify_mu=true; then max degree among kxk minors.
%   A_out  : k-by-k composed left operator: A_out = U * A_in.
%
% Notes:
%   - If you instead want to keep the *original message basis* so that
%     c = m * (topk(B)) = (m*U^{-1}) * Gmin, ask and I’ll add U^{-1} support.

if nargin < 4, opts = struct; end
if nargin < 3 || isempty(A_in)
    error('Please provide A_in (k-by-k) from Smith decomposition.');
end

verify_mu = get_opt(opts,'verify_mu',false);
max_combs = get_opt(opts,'max_combs',3e5);
verbose   = get_opt(opts,'verbose',false);

[n, n2] = size_poly_matrix(B);
if n~=n2, error('B must be square n-by-n in polynomial sense.'); end
if k<1 || k>n, error('k must be between 1 and n.'); end

% Basic encoder from top-k rows of B
G = B(1:k,:);

% Initialize U as polynomial identity
U = eye_poly(k);

% Row-reduce to predictable-degree (Forney row-reduced) form
[Gmin, U, rowdeg, leadpos] = row_reduce_forney(G, U, verbose);

% Constraint length
mu = sum(rowdeg);

% Optional μ verification via k×k minors
mu_check = [];
if verify_mu
    mu_check = max_minor_degree(Gmin, max_combs, verbose);
end

% ---- Compose A_out = U * A_in ----
% Basic polynomial matrix multiply (GF(2)[D])
A_out = pmat_mul(U, A_in);

end

%======================== Polynomial matrix multiply =====================%
function C = pmat_mul(A, B)
% C = A * B over GF(2)[D], with entries stored as low->high coefficient row vectors
[a1, a2] = size(A); [b1, b2] = size(B);
if a2 ~= b1, error('pmat_mul: size mismatch (%dx%d)*( %dx%d )', a1,a2,b1,b2); end
C = cell(a1,b2);
for i=1:a1
    for j=1:b2
        acc = 0;
        for t=1:a2
            acc = p_add(acc, p_mul(A{i,t}, B{t,j}));
        end
        C{i,j} = p_trim(acc);
    end
end
end

%======================== Core Row-Reduction ==============================%
function [Grr, U, rowdeg, leadpos] = row_reduce_forney(G, U, verbose)
[k, n] = size_poly_matrix(G);
Grr = G;
for i=1:k, for j=1:n, Grr{i,j} = p_trim(Grr{i,j}); end, end
changed = true; iter = 0;
while changed
    iter = iter + 1;
    [rowdeg, leadpos] = row_degrees_and_leads(Grr);
    if verbose
        fprintf('[RowReduce] iter=%d; rowdeg=%s; leadpos=%s\n', iter, mat2str(rowdeg), mat2str(leadpos));
    end
    [dupCols, ~, idxG] = unique(leadpos,'stable');
    counts = histcounts(idxG, 1:(numel(dupCols)+1));
    conflicts = dupCols(counts > 1);
    if isempty(conflicts), changed = false; break; end
    for c = conflicts
        I = find(leadpos == c);
        [~, ord] = sort(rowdeg(I), 'descend'); I = I(ord);
        for t = 1:(numel(I)-1)
            i = I(t); p = I(end);
            di = poly_deg(Grr{i,c}); dp = poly_deg(Grr{p,c});
            if di < dp, continue; end
            shift = di - dp;
            for j=1:n
                Grr{i,j} = p_add(Grr{i,j}, p_shift(Grr{p,j}, shift));
            end
            for j=1:k
                U{i,j} = p_add(U{i,j}, p_shift(U{p,j}, shift));
            end
            for j=1:n, Grr{i,j} = p_trim(Grr{i,j}); end
        end
    end
end
[rowdeg, leadpos] = row_degrees_and_leads(Grr);
end

%======================== μ via k×k minors (optional) =====================%
function mu = max_minor_degree(G, max_combs, verbose)
[k, n] = size_poly_matrix(G);
if k > n, mu = -inf; return; end
numC = nchoosek(n, k);
if numC > max_combs
    if verbose, fprintf('[mu-check] C(n,k)=%g too large; random-sampling %g combos.\n', numC, max_combs); end
    rng(0);
    picks = zeros(max_combs, k);
    for s=1:max_combs, picks(s,:) = sort(randsample(n, k)); end
else
    picks = nchoosek(1:n, k);
end
mu_best = -inf;
for r = 1:size(picks, 1)
    Sub = G(:, picks(r,:));
    detp = poly_det(Sub);
    mu_best = max(mu_best, poly_deg(detp));
end
mu = max(0, mu_best);
end

%======================== Degrees/leads, det, and poly ops ===============%
function [rowdeg, leadpos] = row_degrees_and_leads(G)
[k, n] = size_poly_matrix(G);
rowdeg = -inf(1,k); leadpos = zeros(1,k);
for i=1:k
    bestDeg = -inf; bestJ = 1;
    for j=1:n
        d = poly_deg(G{i,j});
        if d > bestDeg, bestDeg = d; bestJ = j; end
    end
    rowdeg(i) = max(bestDeg, 0);
    leadpos(i) = bestJ;
end
end

function detp = poly_det(A)
[m, n] = size_poly_matrix(A);
if m ~= n, error('poly_det requires square.'); end
M = A; detp = 1;
for col = 1:n
    piv = 0; pivdeg = -inf;
    for r = col:n
        d = poly_deg(M{r,col});
        if d > pivdeg, pivdeg = d; piv = r; end
    end
    if pivdeg < 0, detp = 0; return; end
    if piv ~= col, tmp = M(col,:); M(col,:) = M(piv,:); M(piv,:) = tmp; end
    for r = col+1:n
        if poly_deg(M{r,col}) >= 0
            [q, ~] = p_divmod(M{r,col}, M{col,col});
            for c = col:n
                M{r,c} = p_add(M{r,c}, p_mul(q, M{col,c}));
            end
        end
    end
    detp = p_mul(detp, M{col,col});
end
detp = p_trim(detp);
end

function [m, n] = size_poly_matrix(M), m = size(M,1); n = size(M,2); end
function I = eye_poly(k)
I = cell(k,k); for i=1:k, for j=1:k, I{i,j} = double(i==j); end, end
end
function c = p_add(a,b)
la = length(a); lb = length(b); L = max(la, lb);
a = [a, zeros(1,L-la)]; b = [b, zeros(1,L-lb)];
c = mod(a + b, 2); c = p_trim(c);
end
function c = p_shift(a, t), if isempty(a)||all(a==0), c = 0; else, c = [zeros(1,t), a]; end, end
function c = p_mul(a,b)
if (isempty(a)||all(a==0)) || (isempty(b)||all(b==0)), c = 0; return; end
la = length(a); lb = length(b); c = zeros(1, la+lb-1);
for i=1:la, if a(i)==0, continue; end, c(i:(i+lb-1)) = mod(c(i:(i+lb-1)) + b, 2); end
c = p_trim(c);
end
function [q,r] = p_divmod(a,b)
a = p_trim(a); b = p_trim(b);
if all(b==0), error('Division by zero polynomial.'); end
q = 0; r = a; db = poly_deg(b);
while poly_deg(r) >= db && ~all(r==0)
    dr = poly_deg(r); t = dr - db; qt = zeros(1, t+1); qt(end) = 1;
    q = p_add(q, qt); r = p_add(r, p_shift(b, t));
end
q = p_trim(q); r = p_trim(r);
end
function d = poly_deg(a)
a = p_trim(a);
if isempty(a) || all(a==0), d = -inf; else, d = find(a~=0, 1, 'last') - 1; end
end
function a = p_trim(a)
if isempty(a), a = 0; return; end
if ~any(a), a = 0; return; end
i = find(a~=0,1,'last'); a = a(1:i);
end
function v = get_opt(s, name, default)
if isfield(s,name) && ~isempty(s.(name)), v = s.(name); else, v = default; end
end
