function [rl, U] = takagi(N, tol, rounding)
    % Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.
    %
    % Args:
    %   N (complex matrix): square, symmetric matrix N
    %   tol (double): tolerance for symmetry check (default: 1e-13)
    %   rounding (integer): decimal places for rounding singular values (default: 13)
    %
    % Returns:
    %   rl (vector): rounded singular values
    %   U (matrix): Takagi unitary, such that N = U * diag(rl) * U.'

    if nargin < 2
        tol = 1e-13;
    end
    if nargin < 3
        rounding = 13;
    end

    [n, m] = size(N);

    if n ~= m
        error('The input matrix must be square');
    end

    if norm(N - N.') >= tol
        error('The input matrix is not symmetric');
    end

    if all(abs(imag(N(:))) < eps)
        N = real(N);
    end

    if all(N(:) == 0)
        rl = zeros(n, 1);
        U = eye(n);
        return;
    end

    if isreal(N)
        % If the matrix N is real, use its eigendecomposition
        [U, L] = eig(N);
        l = diag(L);
        vals = abs(l);
        phases = sqrt(complex(sign(l)));
        Uc = U * diag(phases);

        [sorted_l, idx] = sort(vals, 'descend');
        Uc = Uc(:, idx);

        rl = sorted_l;
        U = Uc;
        return;
    end

    [v, l, ws] = svd(N);
    w = ws';
    rl = round(diag(l), rounding);

    % Group degenerate singular values
    [~, ~, ic] = unique(rl);
    result = accumarray(ic, (1:numel(rl))', [], @(x) {x});

    % Generate the matrices qs of the degenerate subspaces
    qs = cell(1, length(result));
    for i = 1:length(result)
        idx = result{i};
        vas = v(:, idx);
        was = w(:, idx);
        qs{i} = sqrtm(vas' * was);
    end

    % Construct the Takagi unitary
    qb = blkdiag(qs{:});
    U = v * conj(qb);
end