function rel = compute_rel(part, table)

n = sum(part);
l = length(part);

if l == 1
    rel = ones(n, n);
elseif part(1) == 1
    rel = eye(n);
else
    
    best = find_best_n(part);
    
    if best.total == n
        num = best.num;
        part2 = part(2:l);
        rel2 = Hget(part2, table);
        s2 = size(rel2, 3);
        
        rel = logical(zeros(n, n, nchoosek(n - 1, num - 1)*s2));
        vec = 1:n;
        c = 0;
        
        ind_vec = 2:num;
        while length(ind_vec)
            true_ind_vec = [1 ind_vec];
            v = vec;
            v(true_ind_vec) = 0;
            v = find(v);
            
            c = c + 1;
            rel(v, v, ((c - 1)*s2 + 1):(c*s2)) = rel2;
            rel(true_ind_vec, true_ind_vec, ((c - 1)*s2 + 1):(c*s2)) = 1;
            ind_vec = next_ind_vec(ind_vec, n);
        end
        
    else
        
        p = best.first;
        q = best.last;
        part1 = part(p:q);
        part2 = part([1:(p - 1) (q + 1):l]);
        rel1 = Hget(part1, table);
        rel2 = Hget(part2, table);
        s1 = size(rel1, 3);
        s2 = size(rel2, 3);
        s = s1*s2;
        
        rel = logical(zeros(n, n, nchoosek(n, best.total)*s));
        vec = 1:n;
        c = 0;
        
        ind_vec = 1:best.total;
        while length(ind_vec)
            v = vec;
            v(ind_vec) = 0;
            v = find(v);
            
            c = c + 1;
            for j = 1:s1
                rel(v, v, ((c - 1)*s + (j - 1)*s2 + 1):((c - 1)*s + j*s2)) = rel2;
            end
            for k = 1:s2
                rel(ind_vec, ind_vec, ((c - 1)*s + k):s2:(c*s)) = rel1;
            end
            ind_vec = next_ind_vec(ind_vec, n);
        end
    end
end