
function c = pow(f, R, parent)
    c = zeros(size(f));
    h = zeros(numel(f), length(unique(f)));
    
    [_, s] = sort(f(:), "descend");
    
    for x = fliplr(s')
        y = x;
    
        while true
            h(y, f(x)) = h(y, f(x)) + 1;
            
            if y == parent(y)
                break;
            end
            
            y = parent(y);
        end
    end
    
    for x = s'
        for v = unique(f)
            c(x) = c(x) + h(x, v) * (f(x) - v);
        end
    end
end
