Let $P\subset\mathbb{R}^n$ be a convex polyhedron described as an intersection of hyperspaces, that is, $$P:=\{\boldsymbol{x}: A\boldsymbol{x} \leq \boldsymbol{b}\}$$
Let $\boldsymbol{x} \in P$. We know from Caratheodory theorem that $\boldsymbol{x}$ can be written as a convex combination of at most $n+1$ extreme points of $P$. Is there an efficient algorithm, possibly already implemented in mathematical software such as MATLAB or SCIPY, that performs this decomposition ? That is, an algorithm that finds extreme points $\boldsymbol{p}_1,\ldots,\boldsymbol{p}_{n+1}$ of $P$ together with weights $\boldsymbol{w} \in \mathbb{R}_+^{n+1}$ satisfying $$\boldsymbol{x} = \sum_{i=1}^{n+1} w_i \boldsymbol{p}_i,\quad \sum_{i=1}^{n+1} w_i = 1.$$
Note 1
Intuitvely I think that the construction of such an algorithm is not very complicated. We can start by expressing $\boldsymbol{x}$ as a barycenter of one extreme point of $P$ (say $\boldsymbol{p}_1$), and one point situated on a face of $P$ (say $\boldsymbol{p}_2$). Then, we can continue inductively, by doing the same with $\boldsymbol{p}_2$ which lives in a face of dimension at most $(n-1)$. I wonder if such operations could be done by some kind of pivoting operations on the matrix $A$, like in the simplex algorithm ?
Note 2
For simplicity I assumed that $P$ is bounded, but the same question is of course valid for a unbounded polyhedron, with a description of $\mathbb{x}$ as "convex combination of extreme points" plus "conic combination of extreme rays".
