diff --git a/docs/development.md b/docs/development.md
index 3ee322c5..a7221697 100644
--- a/docs/development.md
+++ b/docs/development.md
@@ -238,5 +238,5 @@ devtools::load_all('./r-phylo2vec')
 
 # A small demo
 v = sample_vector(5, FALSE)
-to_newick_from_vector(v)
+to_newick(v)
 ```
diff --git a/phylo2vec/src/matrix/ops.rs b/phylo2vec/src/matrix/ops.rs
index 5f3c0730..d4eaaca6 100644
--- a/phylo2vec/src/matrix/ops.rs
+++ b/phylo2vec/src/matrix/ops.rs
@@ -65,13 +65,13 @@ pub fn get_node_depths(m: &ArrayView2<f64>) -> Vec<f64> {
 ///     [0.0, 0.5, 0.2],
 /// ];
 /// // Root (node 4) has depth 0
-/// assert_eq!(get_node_depth(&m.view(), 4), 0.0);
+/// assert_eq!(get_node_depth(&m.view(), 4).unwrap(), 0.0);
 /// // Node 3 is 0.5 from root
-/// assert!((get_node_depth(&m.view(), 3) - 0.5).abs() < 1e-10);
+/// assert!((get_node_depth(&m.view(), 3).unwrap() - 0.5).abs() < 1e-10);
 /// // Leaf 0 is 0.5 + 0.3 = 0.8 from root
-/// assert!((get_node_depth(&m.view(), 0) - 0.8).abs() < 1e-10);
+/// assert!((get_node_depth(&m.view(), 0).unwrap() - 0.8).abs() < 1e-10);
 /// ```
-pub fn get_node_depth(m: &ArrayView2<f64>, node: usize) -> f64 {
+pub fn get_node_depth(m: &ArrayView2<f64>, node: usize) -> Result<f64, String> {
     let (v, bls) = parse_matrix(m);
     _get_node_depth(&v, Some(&bls), node)
 }
@@ -160,7 +160,7 @@ mod tests {
         #[case] node: usize,
         #[case] expected_depth: f64,
     ) {
-        let depth = get_node_depth(&m.view(), node);
+        let depth = get_node_depth(&m.view(), node).unwrap();
         assert!(
             (depth - expected_depth).abs() < 1e-10,
             "Expected depth {expected_depth}, got {depth}"
@@ -175,7 +175,7 @@ mod tests {
     fn test_get_node_depth_root_is_zero(#[case] n_leaves: usize) {
         let m = sample_matrix(n_leaves, false);
         let root = 2 * m.nrows(); // Root node index
-        let depth = get_node_depth(&m.view(), root);
+        let depth = get_node_depth(&m.view(), root).unwrap();
         assert_eq!(depth, 0.0, "Root should have depth 0");
     }
 
diff --git a/phylo2vec/src/vector/ops.rs b/phylo2vec/src/vector/ops.rs
index bb507a63..60aaaf95 100644
--- a/phylo2vec/src/vector/ops.rs
+++ b/phylo2vec/src/vector/ops.rs
@@ -160,20 +160,22 @@ fn min_common_ancestor(path1: &[usize], path2: &[usize]) -> usize {
 /// use phylo2vec::vector::ops::get_common_ancestor;
 ///
 /// let v = vec![0, 1, 2, 3, 4];
-/// let mrca = get_common_ancestor(&v, 2, 3);
+/// let mrca = get_common_ancestor(&v, 2, 3).unwrap();
 /// // Newick string of v = (0,(1,(2,(3,(4,5)6)7)8)9)10;
 /// // So 2 and 3 share the MRCA 8.
 /// assert_eq!(mrca, 8);
 /// ```
-pub fn get_common_ancestor(v: &[usize], node1: usize, node2: usize) -> usize {
+pub fn get_common_ancestor(v: &[usize], node1: usize, node2: usize) -> Result<usize, String> {
     let node_max = 2 * v.len();
 
     if node1 > node_max || node2 > node_max {
-        panic!("Node indices must be within the range of the Phylo2Vec vector. Max = {node_max}, got node1 = {node1} and node2 = {node2}");
+        return Err(format!(
+              "Node indices must be within range. Max = {node_max}, got node1 = {node1} and node2 = {node2}"
+          ));
     }
 
     if node1 == node2 {
-        return node1;
+        return Ok(node1);
     }
 
     let pairs = to_pairs(v);
@@ -181,7 +183,7 @@ pub fn get_common_ancestor(v: &[usize], node1: usize, node2: usize) -> usize {
     let path1 = get_ancestry_path_of_node(&pairs, node1);
     let path2 = get_ancestry_path_of_node(&pairs, node2);
 
-    min_common_ancestor(&path1, &path2)
+    Ok(min_common_ancestor(&path1, &path2))
 }
 
 /// Generic function to calculate the depths of all nodes in a Phylo2Vec tree.
@@ -238,20 +240,24 @@ pub fn _get_node_depths(v: &[usize], bls: Option<&Vec<[f64; 2]>>) -> Vec<f64> {
 /// The root has depth 0, and depths increase as you move toward the leaves.
 ///
 /// When `bls` is `None`, all branch lengths are assumed to be 1.0 (topological depth).
-pub fn _get_node_depth(v: &[usize], bls: Option<&Vec<[f64; 2]>>, node: usize) -> f64 {
+pub fn _get_node_depth(
+    v: &[usize],
+    bls: Option<&Vec<[f64; 2]>>,
+    node: usize,
+) -> Result<f64, String> {
     let k = v.len();
     let n_leaves = k + 1;
     let n_nodes = 2 * n_leaves - 1;
 
     if node >= n_nodes {
-        panic!(
+        return Err(format!(
             "Node index out of bounds. Max node = {}, got node = {}",
             n_nodes - 1,
             node
-        );
+        ));
     }
 
-    _get_node_depths(v, bls)[node]
+    Ok(_get_node_depths(v, bls)[node])
 }
 
 /// Get the depths of all nodes in a Phylo2Vec vector (topological).
@@ -297,15 +303,15 @@ pub fn get_node_depths(v: &[usize]) -> Vec<f64> {
 /// // Tree: (0,(1,(2,3)4)5)6
 /// let v = vec![0, 1, 2];
 /// // Root (node 6) has depth 0
-/// assert_eq!(get_node_depth(&v, 6), 0.0);
+/// assert_eq!(get_node_depth(&v, 6).unwrap(), 0.0);
 /// // Node 5 is 1 edge from root
-/// assert_eq!(get_node_depth(&v, 5), 1.0);
+/// assert_eq!(get_node_depth(&v, 5).unwrap(), 1.0);
 /// // Node 4 is 2 edges from root
-/// assert_eq!(get_node_depth(&v, 4), 2.0);
+/// assert_eq!(get_node_depth(&v, 4).unwrap(), 2.0);
 /// // Leaf 0 is 1 edge from root
-/// assert_eq!(get_node_depth(&v, 0), 1.0);
+/// assert_eq!(get_node_depth(&v, 0).unwrap(), 1.0);
 /// ```
-pub fn get_node_depth(v: &[usize], node: usize) -> f64 {
+pub fn get_node_depth(v: &[usize], node: usize) -> Result<f64, String> {
     _get_node_depth(v, None, node)
 }
 
@@ -565,7 +571,7 @@ mod tests {
         #[case] node2: usize,
         #[case] expected_mrca: usize,
     ) {
-        let mrca = get_common_ancestor(&v, node1, node2);
+        let mrca = get_common_ancestor(&v, node1, node2).unwrap();
         assert_eq!(
             mrca, expected_mrca,
             "Expected mrca of nodes {node1} and {node2} for v = {v:?} to be {expected_mrca}, but got {mrca}"
@@ -599,7 +605,7 @@ mod tests {
         #[case] node: usize,
         #[case] expected_depth: f64,
     ) {
-        let depth = get_node_depth(&v, node);
+        let depth = get_node_depth(&v, node).unwrap();
         assert!(
             (depth - expected_depth).abs() < 1e-10,
             "Expected depth {expected_depth} for node {node} in v = {v:?}, got {depth}"
@@ -616,7 +622,7 @@ mod tests {
 
         let v = sample_vector(n_leaves, false);
         let root = 2 * v.len(); // Root node index
-        let depth = get_node_depth(&v, root);
+        let depth = get_node_depth(&v, root).unwrap();
         assert_eq!(depth, 0.0, "Root should have depth 0");
     }
 
@@ -656,9 +662,9 @@ mod tests {
     #[rstest]
     #[case(vec![0, 1, 2], 7)] // Max node is 6
     #[case(vec![0], 3)] // Max node is 2
-    #[should_panic]
+    #[should_panic(expected = "Node index out of bounds")]
     fn test_get_node_depth_out_of_bounds(#[case] v: Vec<usize>, #[case] node: usize) {
-        get_node_depth(&v, node);
+        get_node_depth(&v, node).unwrap();
     }
 
     #[rstest]
diff --git a/py-phylo2vec/benchmarks/test_bench.py b/py-phylo2vec/benchmarks/test_bench.py
index 13402ce1..3027cd71 100644
--- a/py-phylo2vec/benchmarks/test_bench.py
+++ b/py-phylo2vec/benchmarks/test_bench.py
@@ -18,7 +18,7 @@ def test_vector_to_newick_ordered(benchmark, sample_size):
     sample_size : int
         Number of leaves in the vector
     """
-    benchmark(core.to_newick_from_vector, core.sample_vector(sample_size, True))
+    benchmark(core.to_newick, core.sample_vector(sample_size, True))
 
 
 @pytest.mark.parametrize("sample_size", BIG_RANGE)
@@ -32,7 +32,7 @@ def test_vector_to_newick(benchmark, sample_size):
     sample_size : int
         Number of leaves in the vector
     """
-    benchmark(core.to_newick_from_vector, core.sample_vector(sample_size, False))
+    benchmark(core.to_newick, core.sample_vector(sample_size, False))
 
 
 @pytest.mark.parametrize("sample_size", BIG_RANGE)
@@ -46,7 +46,7 @@ def test_matrix_to_newick(benchmark, sample_size):
     sample_size : int
         Number of leaves in the vector
     """
-    benchmark(core.to_newick_from_matrix, core.sample_matrix(sample_size, False))
+    benchmark(core.to_newick, core.sample_matrix(sample_size, False))
 
 
 @pytest.mark.parametrize("sample_size", BIG_RANGE)
@@ -62,7 +62,7 @@ def test_vector_from_newick(benchmark, sample_size):
     """
     benchmark(
         core.to_vector,
-        core.to_newick_from_vector(core.sample_vector(sample_size, False)),
+        core.to_newick(core.sample_vector(sample_size, False)),
     )
 
 
@@ -79,5 +79,5 @@ def test_matrix_from_newick(benchmark, sample_size):
     """
     benchmark(
         core.to_matrix,
-        core.to_newick_from_matrix(core.sample_matrix(sample_size, False)),
+        core.to_newick(core.sample_matrix(sample_size, False)),
     )
diff --git a/py-phylo2vec/phylo2vec/base/newick.py b/py-phylo2vec/phylo2vec/base/newick.py
index edaa9aef..74e9f8c2 100644
--- a/py-phylo2vec/phylo2vec/base/newick.py
+++ b/py-phylo2vec/phylo2vec/base/newick.py
@@ -39,13 +39,4 @@ def to_newick(vector_or_matrix: np.ndarray) -> str:
     newick : str
         Newick tree
     """
-    if vector_or_matrix.ndim == 2:
-        newick = core.to_newick_from_matrix(vector_or_matrix)
-    elif vector_or_matrix.ndim == 1:
-        newick = core.to_newick_from_vector(vector_or_matrix)
-    else:
-        raise ValueError(
-            "vector_or_matrix should either be a vector (ndim == 1) or matrix (ndim == 2)"
-        )
-
-    return newick
+    return core.to_newick(vector_or_matrix)
diff --git a/py-phylo2vec/phylo2vec/stats/nodewise.py b/py-phylo2vec/phylo2vec/stats/nodewise.py
index fa8f1a0b..fb851e37 100644
--- a/py-phylo2vec/phylo2vec/stats/nodewise.py
+++ b/py-phylo2vec/phylo2vec/stats/nodewise.py
@@ -7,14 +7,14 @@
 from phylo2vec.utils.vector import check_vector
 
 
-def cophenetic_distances(vector_or_matrix, unrooted=False):
+def cophenetic_distances(tree, unrooted=False):
     """
     Compute the cophenetic distance matrix of a Phylo2Vec
     vector (topological) or matrix (from branch lengths).
 
     Parameters
     ----------
-    vector_or_matrix : numpy.ndarray
+    tree : numpy.ndarray
         Phylo2Vec vector (ndim == 1)/matrix (ndim == 2)
 
     Returns
@@ -22,21 +22,13 @@ def cophenetic_distances(vector_or_matrix, unrooted=False):
     numpy.ndarray
         Cophenetic distance matrix
     """
-    if vector_or_matrix.ndim == 2:
-        coph = core.cophenetic_distances_with_bls(vector_or_matrix, unrooted=unrooted)
-    elif vector_or_matrix.ndim == 1:
-        coph = core.cophenetic_distances(vector_or_matrix, unrooted=unrooted)
-    else:
-        raise ValueError(
-            "vector_or_matrix should either be a vector (ndim == 1) or matrix (ndim == 2)"
-        )
-    return coph
+    return core.cophenetic_distances(tree, unrooted=unrooted)
 
 
 NODEWISE_DISTANCES = {"cophenetic": cophenetic_distances}
 
 
-def pairwise_distances(vector_or_matrix, metric="cophenetic"):
+def pairwise_distances(tree, metric="cophenetic"):
     """
     Compute a pairwise distance matrix
     for tree nodes from a Phylo2Vec vector.
@@ -45,7 +37,7 @@ def pairwise_distances(vector_or_matrix, metric="cophenetic"):
 
     Parameters
     ----------
-    vector_or_matrix : numpy.ndarray
+    tree : numpy.ndarray
         Phylo2Vec vector (ndim == 1)/matrix (ndim == 2)
     metric : str, optional
         Pairwise distance metric, by default "cophenetic"
@@ -55,21 +47,12 @@ def pairwise_distances(vector_or_matrix, metric="cophenetic"):
     numpy.ndarray
         Distance matrix
     """
-    if vector_or_matrix.ndim == 2:
-        check_matrix(vector_or_matrix)
-    elif vector_or_matrix.ndim == 1:
-        check_vector(vector_or_matrix)
-    else:
-        raise ValueError(
-            "vector_or_matrix should either be a vector (ndim == 1) or matrix (ndim == 2)"
-        )
-
     func = NODEWISE_DISTANCES[metric]
 
-    return func(vector_or_matrix)
+    return func(tree)
 
 
-def cov(vector_or_matrix):
+def cov(tree):
     """
     Compute the covariance matrix of a Phylo2Vec vector or matrix.
 
@@ -77,7 +60,7 @@ def cov(vector_or_matrix):
 
     Parameters
     ----------
-    vector_or_matrix : numpy.ndarray
+    tree : numpy.ndarray
         Phylo2Vec vector (ndim == 1)/matrix (ndim == 2)
 
     Returns
@@ -85,18 +68,10 @@ def cov(vector_or_matrix):
     vcv : numpy.ndarray
         Covariance matrix
     """
-    if vector_or_matrix.ndim == 2:
-        vcv = core.vcv_with_bls(vector_or_matrix)
-    elif vector_or_matrix.ndim == 1:
-        vcv = core.vcv(vector_or_matrix)
-    else:
-        raise ValueError(
-            "vector_or_matrix should either be a vector (ndim == 1) or matrix (ndim == 2)"
-        )
-    return vcv
-
-
-def precision(vector_or_matrix):
+    return core.vcv(tree)
+
+
+def precision(tree):
     """
     Compute the precision matrix of a Phylo2Vec vector or matrix.
 
@@ -104,7 +79,7 @@ def precision(vector_or_matrix):
 
     Parameters
     ----------
-    vector_or_matrix : numpy.ndarray
+    tree : numpy.ndarray
         Phylo2Vec vector (ndim == 1)/matrix (ndim == 2)
 
     Returns
@@ -112,17 +87,10 @@ def precision(vector_or_matrix):
     precision: numpy.ndarray
         Precision matrix
     """
-    if vector_or_matrix.ndim == 2:
-        precursor = core.pre_precision_with_bls(vector_or_matrix)
-    elif vector_or_matrix.ndim == 1:
-        precursor = core.pre_precision(vector_or_matrix)
-    else:
-        raise ValueError(
-            "vector_or_matrix should either be a vector (ndim == 1) or matrix (ndim == 2)"
-        )
+    precursor = core.pre_precision(tree)
 
     # Schur complement of the precursor matrix
-    n_leaves = vector_or_matrix.shape[0] + 1
+    n_leaves = tree.shape[0] + 1
     a = precursor[:n_leaves, :n_leaves]
     b = precursor[:n_leaves, n_leaves:]
     c = precursor[n_leaves:, n_leaves:]
diff --git a/py-phylo2vec/phylo2vec/stats/treewise.py b/py-phylo2vec/phylo2vec/stats/treewise.py
index 5f8064e2..134d21a6 100644
--- a/py-phylo2vec/phylo2vec/stats/treewise.py
+++ b/py-phylo2vec/phylo2vec/stats/treewise.py
@@ -55,8 +55,4 @@ def robinson_foulds(
     ete3.Tree.robinson_foulds : Reference implementation in ete3
     ape::dist.topo : Reference implementation in R's ape package
     """
-    # Extract topology (column 0) if matrix input
-    v1 = tree1[:, 0].astype(int).tolist() if tree1.ndim == 2 else tree1.tolist()
-    v2 = tree2[:, 0].astype(int).tolist() if tree2.ndim == 2 else tree2.tolist()
-
-    return core.robinson_foulds(v1, v2, normalize)
+    return core.robinson_foulds(tree1, tree2, normalize)
diff --git a/py-phylo2vec/phylo2vec/utils/vector.py b/py-phylo2vec/phylo2vec/utils/vector.py
index c31d309b..217258c1 100644
--- a/py-phylo2vec/phylo2vec/utils/vector.py
+++ b/py-phylo2vec/phylo2vec/utils/vector.py
@@ -251,22 +251,10 @@ def get_common_ancestor(vector_or_matrix, node1, node2):
 
     Returns
     -------
-    mrca : int
+    int
         Most recent common ancestor node between node1 and node2
     """
-    if vector_or_matrix.ndim == 2:
-        v = vector_or_matrix[:, 0].astype(int)
-    elif vector_or_matrix.ndim == 1:
-        v = vector_or_matrix
-    else:
-        raise ValueError(
-            "vector_or_matrix should either be a vector (ndim == 1) or matrix (ndim == 2)"
-        )
-
-    if not (0 <= node1 <= 2 * len(v) and 0 <= node2 <= 2 * len(v)):
-        raise ValueError("Nodes must be in the range [0, 2 * n_leaves]")
-
-    return core.get_common_ancestor(v.tolist(), node1, node2)
+    return core.get_common_ancestor(vector_or_matrix, node1, node2)
 
 
 def get_node_depth(vector_or_matrix, node: int) -> float:
@@ -303,19 +291,7 @@ def get_node_depth(vector_or_matrix, node: int) -> float:
     >>> get_node_depth(v, 0)  # Leaf 0 is 1 edge from root
     1.0
     """
-    n_leaves = vector_or_matrix.shape[0] + 1
-    max_node = 2 * n_leaves - 1
-    if not 0 <= node < max_node:
-        raise ValueError(f"node must be in the range [0, {max_node})")
-
-    if vector_or_matrix.ndim == 2:
-        return core.get_node_depth_with_bls(vector_or_matrix, node)
-    if vector_or_matrix.ndim == 1:
-        return core.get_node_depth(vector_or_matrix.tolist(), node)
-
-    raise ValueError(
-        "vector_or_matrix should either be a vector (ndim == 1) or matrix (ndim == 2)"
-    )
+    return core.get_node_depth(vector_or_matrix, node)
 
 
 def get_node_depths(vector_or_matrix) -> np.ndarray:
@@ -352,14 +328,7 @@ def get_node_depths(vector_or_matrix) -> np.ndarray:
     >>> depths[0]  # Leaf 0 is 1 edge from root
     1.0
     """
-    if vector_or_matrix.ndim == 2:
-        return np.asarray(core.get_node_depths_with_bls(vector_or_matrix))
-    if vector_or_matrix.ndim == 1:
-        return np.asarray(core.get_node_depths(vector_or_matrix.tolist()))
-
-    raise ValueError(
-        "vector_or_matrix should either be a vector (ndim == 1) or matrix (ndim == 2)"
-    )
+    return np.asarray(core.get_node_depths(vector_or_matrix))
 
 
 def reroot(vector_or_matrix, node) -> np.ndarray:
diff --git a/py-phylo2vec/src/lib.rs b/py-phylo2vec/src/lib.rs
index 13168ae6..f603ec5b 100644
--- a/py-phylo2vec/src/lib.rs
+++ b/py-phylo2vec/src/lib.rs
@@ -17,6 +17,12 @@ use phylo2vec::vector::distance as vdist;
 use phylo2vec::vector::graph as vgraph;
 use phylo2vec::vector::ops as vops;
 
+#[derive(FromPyObject)]
+enum PyTree<'py> {
+    Matrix(PyReadonlyArray2<'py, f64>),
+    Vector(Vec<usize>),
+}
+
 #[pyfunction]
 fn sample_vector(n_leaves: isize, ordered: bool) -> PyResult<Vec<usize>> {
     if n_leaves < 2 {
@@ -82,16 +88,11 @@ fn check_m(input_matrix: PyReadonlyArray2<f64>) -> PyResult<()> {
 }
 
 #[pyfunction]
-fn to_newick_from_vector(input_vector: Vec<usize>) -> PyResult<String> {
-    let newick = vconvert::to_newick(&input_vector);
-    Ok(newick)
-}
-
-#[pyfunction]
-fn to_newick_from_matrix(input_matrix: PyReadonlyArray2<f64>) -> PyResult<String> {
-    let arr = input_matrix.as_array();
-    let newick = mconvert::to_newick(&arr);
-    Ok(newick)
+fn to_newick(tree: PyTree) -> String {
+    match tree {
+        PyTree::Matrix(m) => mconvert::to_newick(&m.as_array()),
+        PyTree::Vector(v) => vconvert::to_newick(&v),
+    }
 }
 
 #[pyfunction]
@@ -141,50 +142,33 @@ fn build_newick(input_pairs: Vec<(usize, usize)>) -> String {
 }
 
 #[pyfunction]
-fn cophenetic_distances(
-    py: Python<'_>,
-    input_vector: Vec<usize>,
-    unrooted: bool,
-) -> Bound<'_, PyArray2<f64>> {
-    vgraph::cophenetic_distances(&input_vector, unrooted).into_pyarray(py)
-}
-
-#[pyfunction]
-fn cophenetic_distances_with_bls<'py>(
+fn cophenetic_distances<'py>(
     py: Python<'py>,
-    input_matrix: PyReadonlyArray2<f64>,
+    tree: PyTree<'py>,
     unrooted: bool,
-) -> Bound<'py, PyArray2<f64>> {
-    let m = input_matrix.as_array();
-    mgraph::cophenetic_distances(&m, unrooted).into_pyarray(py)
-}
-
-#[pyfunction]
-fn pre_precision(py: Python<'_>, input_vector: Vec<usize>) -> PyResult<Bound<'_, PyArray2<f64>>> {
-    Ok(vgraph::pre_precision(&input_vector).into_pyarray(py))
-}
-
-#[pyfunction]
-fn pre_precision_with_bls<'py>(
-    py: Python<'py>,
-    input_matrix: PyReadonlyArray2<f64>,
 ) -> PyResult<Bound<'py, PyArray2<f64>>> {
-    let m = input_matrix.as_array();
-    Ok(mgraph::pre_precision(&m).into_pyarray(py))
+    match tree {
+        PyTree::Matrix(m) => {
+            Ok(mgraph::cophenetic_distances(&m.as_array(), unrooted).into_pyarray(py))
+        }
+        PyTree::Vector(v) => Ok(vgraph::cophenetic_distances(&v, unrooted).into_pyarray(py)),
+    }
 }
 
 #[pyfunction]
-fn vcv(py: Python<'_>, input_vector: Vec<usize>) -> PyResult<Bound<'_, PyArray2<f64>>> {
-    Ok(vgraph::vcv(&input_vector).into_pyarray(py))
+fn pre_precision<'py>(py: Python<'py>, tree: PyTree<'py>) -> PyResult<Bound<'py, PyArray2<f64>>> {
+    match tree {
+        PyTree::Matrix(m) => Ok(mgraph::pre_precision(&m.as_array()).into_pyarray(py)),
+        PyTree::Vector(v) => Ok(vgraph::pre_precision(&v).into_pyarray(py)),
+    }
 }
 
 #[pyfunction]
-fn vcv_with_bls<'py>(
-    py: Python<'py>,
-    input_matrix: PyReadonlyArray2<f64>,
-) -> PyResult<Bound<'py, PyArray2<f64>>> {
-    let m = input_matrix.as_array();
-    Ok(mgraph::vcv(&m).into_pyarray(py))
+fn vcv<'py>(py: Python<'py>, tree: PyTree<'py>) -> PyResult<Bound<'py, PyArray2<f64>>> {
+    match tree {
+        PyTree::Matrix(m) => Ok(mgraph::vcv(&m.as_array()).into_pyarray(py)),
+        PyTree::Vector(v) => Ok(vgraph::vcv(&v).into_pyarray(py)),
+    }
 }
 
 #[pyfunction]
@@ -263,50 +247,43 @@ fn create_label_mapping(newick: String) -> PyResult<(String, HashMap<usize, Stri
 }
 
 #[pyfunction]
-fn get_common_ancestor(v: Vec<usize>, node1: usize, node2: usize) -> usize {
-    vops::get_common_ancestor(&v, node1, node2)
-}
-
-#[pyfunction]
-fn get_node_depth(v: Vec<usize>, node: isize) -> PyResult<f64> {
-    if node < 0 {
-        return Err(PyValueError::new_err("node must be non-negative"));
-    }
-    let node = node as usize;
-    let n_nodes = 2 * v.len() + 1;
-    if node >= n_nodes {
-        return Err(PyValueError::new_err(format!(
-            "node must be less than {n_nodes}"
-        )));
+fn get_common_ancestor(tree: PyTree, node1: isize, node2: isize) -> PyResult<usize> {
+    if node1 < 0 || node2 < 0 {
+        return Err(PyValueError::new_err("nodes must be non-negative"));
     }
-    Ok(vops::get_node_depth(&v, node))
+    let node1 = node1 as usize;
+    let node2 = node2 as usize;
+    let v = match tree {
+        PyTree::Matrix(m) => {
+            let (v, _) = mbase::parse_matrix(&m.as_array());
+            v
+        }
+        PyTree::Vector(v) => v,
+    };
+
+    vops::get_common_ancestor(&v, node1, node2).map_err(PyValueError::new_err)
 }
 
 #[pyfunction]
-fn get_node_depth_with_bls(input_matrix: PyReadonlyArray2<f64>, node: isize) -> PyResult<f64> {
+fn get_node_depth(tree: PyTree, node: isize) -> PyResult<f64> {
     if node < 0 {
         return Err(PyValueError::new_err("node must be non-negative"));
     }
     let node = node as usize;
-    let m = input_matrix.as_array();
-    let n_nodes = 2 * m.nrows() + 1;
-    if node >= n_nodes {
-        return Err(PyValueError::new_err(format!(
-            "node must be less than {n_nodes}"
-        )));
-    }
-    Ok(mops::get_node_depth(&m, node))
-}
+    let depth = match tree {
+        PyTree::Matrix(m) => mops::get_node_depth(&m.as_array(), node),
+        PyTree::Vector(v) => vops::get_node_depth(&v, node),
+    };
 
-#[pyfunction]
-fn get_node_depths(v: Vec<usize>) -> Vec<f64> {
-    vops::get_node_depths(&v)
+    depth.map_err(PyValueError::new_err)
 }
 
 #[pyfunction]
-fn get_node_depths_with_bls(input_matrix: PyReadonlyArray2<f64>) -> Vec<f64> {
-    let m = input_matrix.as_array();
-    mops::get_node_depths(&m)
+fn get_node_depths(tree: PyTree) -> Vec<f64> {
+    match tree {
+        PyTree::Matrix(m) => mops::get_node_depths(&m.as_array()),
+        PyTree::Vector(v) => vops::get_node_depths(&v),
+    }
 }
 
 #[pyfunction]
@@ -315,8 +292,19 @@ fn queue_shuffle(v: Vec<usize>, shuffle_cherries: bool) -> (Vec<usize>, Vec<usiz
 }
 
 #[pyfunction]
-#[pyo3(signature = (v1, v2, normalize=false))]
-fn robinson_foulds(v1: Vec<usize>, v2: Vec<usize>, normalize: bool) -> f64 {
+fn robinson_foulds(tree1: PyTree, tree2: PyTree, normalize: bool) -> f64 {
+    let (v1, v2) = match (tree1, tree2) {
+        (PyTree::Matrix(m1), PyTree::Matrix(m2)) => (
+            mbase::parse_matrix(&m1.as_array()).0,
+            mbase::parse_matrix(&m2.as_array()).0,
+        ),
+        (PyTree::Vector(v1), PyTree::Vector(v2)) => (v1, v2),
+        (PyTree::Matrix(m), PyTree::Vector(v)) | (PyTree::Vector(v), PyTree::Matrix(m)) => {
+            let v_from_m = mbase::parse_matrix(&m.as_array()).0;
+            (v_from_m, v)
+        }
+    };
+
     vdist::robinson_foulds(&v1, &v2, normalize)
 }
 
@@ -348,7 +336,6 @@ fn _phylo2vec_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(check_m, m)?)?;
     m.add_function(wrap_pyfunction!(check_v, m)?)?;
     m.add_function(wrap_pyfunction!(cophenetic_distances, m)?)?;
-    m.add_function(wrap_pyfunction!(cophenetic_distances_with_bls, m)?)?;
     m.add_function(wrap_pyfunction!(create_label_mapping, m)?)?;
     m.add_function(wrap_pyfunction!(find_num_leaves, m)?)?;
     m.add_function(wrap_pyfunction!(from_ancestry, m)?)?;
@@ -358,9 +345,7 @@ fn _phylo2vec_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(get_common_ancestor, m)?)?;
     m.add_function(wrap_pyfunction!(get_edges, m)?)?;
     m.add_function(wrap_pyfunction!(get_node_depth, m)?)?;
-    m.add_function(wrap_pyfunction!(get_node_depth_with_bls, m)?)?;
     m.add_function(wrap_pyfunction!(get_node_depths, m)?)?;
-    m.add_function(wrap_pyfunction!(get_node_depths_with_bls, m)?)?;
     m.add_function(wrap_pyfunction!(get_pairs, m)?)?;
     m.add_function(wrap_pyfunction!(has_branch_lengths, m)?)?;
     m.add_function(wrap_pyfunction!(has_parents, m)?)?;
@@ -370,7 +355,6 @@ fn _phylo2vec_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(incidence_dense, m)?)?;
     m.add_function(wrap_pyfunction!(leaf_depth_variance, m)?)?;
     m.add_function(wrap_pyfunction!(pre_precision, m)?)?;
-    m.add_function(wrap_pyfunction!(pre_precision_with_bls, m)?)?;
     m.add_function(wrap_pyfunction!(queue_shuffle, m)?)?;
     m.add_function(wrap_pyfunction!(remove_branch_lengths, m)?)?;
     m.add_function(wrap_pyfunction!(remove_leaf, m)?)?;
@@ -379,12 +363,10 @@ fn _phylo2vec_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(sackin, m)?)?;
     m.add_function(wrap_pyfunction!(sample_matrix, m)?)?;
     m.add_function(wrap_pyfunction!(sample_vector, m)?)?;
-    m.add_function(wrap_pyfunction!(to_newick_from_vector, m)?)?;
-    m.add_function(wrap_pyfunction!(to_newick_from_matrix, m)?)?;
     m.add_function(wrap_pyfunction!(to_matrix, m)?)?;
+    m.add_function(wrap_pyfunction!(to_newick, m)?)?;
     m.add_function(wrap_pyfunction!(to_vector, m)?)?;
     m.add_function(wrap_pyfunction!(vcv, m)?)?;
-    m.add_function(wrap_pyfunction!(vcv_with_bls, m)?)?;
     // Metadata about the package bindings
     m.add("__version__", env!("CARGO_PKG_VERSION"))?;
     Ok(())
diff --git a/py-phylo2vec/tests/test_base.py b/py-phylo2vec/tests/test_base.py
index c7e96e2d..925337e8 100644
--- a/py-phylo2vec/tests/test_base.py
+++ b/py-phylo2vec/tests/test_base.py
@@ -35,17 +35,17 @@ def test_v2newick2v(n_leaves):
 
 class TestToNewickEdgeCases:
     def test_to_newick_empty(self):
-        # dim 0
+        # Empty array
         v0 = np.array(0)
-        # Check that we raise a ValueError
-        with pytest.raises(ValueError):
+        # Check that we raise a TypeError
+        with pytest.raises(TypeError):
             to_newick(v0)
 
     def test_to_newick_ndim3(self):
         # array with 3 dimensions
         t = np.zeros((MIN_N_LEAVES, 3, 1))
-        # Check that we raise a ValueError
-        with pytest.raises(ValueError):
+        # Check that we raise a TypeError
+        with pytest.raises(TypeError):
             to_newick(t)
 
 
diff --git a/r-phylo2vec/DESCRIPTION b/r-phylo2vec/DESCRIPTION
index b9e5d5b1..9a1cc54d 100644
--- a/r-phylo2vec/DESCRIPTION
+++ b/r-phylo2vec/DESCRIPTION
@@ -7,7 +7,7 @@ Description: A tool for working with binary (phylogenetic) trees
 License: GPL-3 + file LICENSE
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Config/rextendr/version: 0.4.2
 Suggests:
     testthat (>= 3.0.0),
diff --git a/r-phylo2vec/NAMESPACE b/r-phylo2vec/NAMESPACE
index 2b9d6005..67144aa5 100644
--- a/r-phylo2vec/NAMESPACE
+++ b/r-phylo2vec/NAMESPACE
@@ -18,6 +18,7 @@ export(has_branch_lengths)
 export(incidence)
 export(load_newick)
 export(load_p2v)
+export(pre_precision)
 export(precision)
 export(queue_shuffle)
 export(remove_branch_lengths)
diff --git a/r-phylo2vec/R/base_newick.R b/r-phylo2vec/R/base_newick.R
index d201d60b..210d4113 100644
--- a/r-phylo2vec/R/base_newick.R
+++ b/r-phylo2vec/R/base_newick.R
@@ -1,21 +1,3 @@
-#' Convert a phylo2vec vector or matrix to Newick format
-#'
-#' @param vector_or_matrix phylo2vec matrix (tree with branch
-#' lengths) or a vector (tree topology only)
-#' @return Newick string representation of the tree
-#' @export
-to_newick <- function(vector_or_matrix) {
-  if (is.vector(vector_or_matrix, "integer")) {
-    # If the input is a vector, call the C function for vector
-    .Call(wrap__to_newick_from_vector, vector_or_matrix)
-  } else if (is.matrix(vector_or_matrix)) {
-    # If the input is a matrix, call the C function for matrix
-    .Call(wrap__to_newick_from_matrix, vector_or_matrix)
-  } else {
-    stop("Input must be either an integer vector or a 2D matrix.")
-  }
-}
-
 #' Convert a Newick string to a phylo2vec vector or matrix
 #'
 #' @param newick Newick string representation of the tree
diff --git a/r-phylo2vec/R/extendr-wrappers.R b/r-phylo2vec/R/extendr-wrappers.R
index 84c43f93..54108357 100644
--- a/r-phylo2vec/R/extendr-wrappers.R
+++ b/r-phylo2vec/R/extendr-wrappers.R
@@ -52,9 +52,16 @@ check_m <- function(vector) invisible(.Call(wrap__check_m, vector))
 #' @export
 check_v <- function(vector) invisible(.Call(wrap__check_v, vector))
 
-cophenetic_from_matrix <- function(matrix, unrooted) .Call(wrap__cophenetic_from_matrix, matrix, unrooted)
-
-cophenetic_from_vector <- function(vector, unrooted) .Call(wrap__cophenetic_from_vector, vector, unrooted)
+#' Get the (topological) cophenetic distance matrix of a phylo2vec tree
+#'
+#' The cophenetic distance between two leaves is the distance from each leaf to their most recent common ancestor.
+#' For phylo2vec vectors, this is the topological distance. For phylo2vec matrices, this is the distance with branch lengths.
+#'
+#' @param tree A phylo2vec tree
+#' @param unrooted If true, the distance is calculated as the distance between each leaf and their most recent common ancestor, multiplied by 2. If false, the distance is calculated as the distance from each leaf to their most recent common ancestor.
+#' @return Cophenetic distance matrix (shape: (n_leaves, n_leaves))
+#' @export
+cophenetic_distances <- function(tree, unrooted) .Call(wrap__cophenetic_distances, tree, unrooted)
 
 #' Create an integer-taxon label mapping (label_mapping)
 #' from a string-based newick (where leaves are strings)
@@ -100,15 +107,35 @@ from_edges <- function(edges) .Call(wrap__from_edges, edges)
 #' @export
 from_pairs <- function(pairs) .Call(wrap__from_pairs, pairs)
 
-get_common_ancestor_from_vector <- function(vector, node1, node2) .Call(wrap__get_common_ancestor, vector, node1, node2)
-
-get_node_depth_from_matrix <- function(matrix, node) .Call(wrap__get_node_depth_from_matrix, matrix, node)
-
-get_node_depth_from_vector <- function(vector, node) .Call(wrap__get_node_depth_from_vector, vector, node)
+#' Get the first recent common ancestor between two nodes in a phylo2vec tree
+#' node1 and node2 can be leaf nodes (0 to n_leaves) or internal nodes (n_leaves to 2*(n_leaves-1)).
+#' Similar to ape's `getMRCA` function in R (for leaf nodes)
+#' and ETE's `get_common_ancestor` in Python (for all nodes), but for phylo2vec vectors.
+#'
+#' @param vector phylo2vec vector representation of a tree topology
+#' @param node1 The first node (0-indexed)
+#' @param node2 The second node (0-indexed)
+#' @return The common ancestor node (0-indexed)
+#' @export
+get_common_ancestor <- function(tree, node1, node2) .Call(wrap__get_common_ancestor, tree, node1, node2)
 
-get_node_depths_from_matrix <- function(matrix) .Call(wrap__get_node_depths_from_matrix, matrix)
+#' Get the depth of a node in a phylo2vec tree
+#' Depth is the distance from root to that node. Root has depth 0.
+#' For vectors, this is the topological depth. For matrices, this is the depth with branch lengths.
+#'
+#' @param tree A phylo2vec tree
+#' @param node The node to get the depth of (0-indexed)
+#' @return The depth of the node
+#' @export
+get_node_depth <- function(tree, node) .Call(wrap__get_node_depth, tree, node)
 
-get_node_depths_from_vector <- function(vector) .Call(wrap__get_node_depths_from_vector, vector)
+#' Get the depths of all nodes in a phylo2vec tree
+#' Depth is the distance from root to each node. Root has depth 0.
+#'
+#' @param tree A phylo2vec tree
+#' @return A vector of depths for all nodes (length: 2*n_leaves - 1)
+#' @export
+get_node_depths <- function(tree) .Call(wrap__get_node_depths, tree)
 
 #' Check if a newick string has branch lengths
 #'
@@ -125,10 +152,6 @@ incidence_csr <- function(input_vector) .Call(wrap__incidence_csr, input_vector)
 
 incidence_dense <- function(input_vector) .Call(wrap__incidence_dense, input_vector)
 
-pre_precision_from_matrix <- function(matrix) .Call(wrap__pre_precision_from_matrix, matrix)
-
-pre_precision_from_vector <- function(vector) .Call(wrap__pre_precision_from_vector, vector)
-
 #' Produce an ordered version (i.e., birth-death process version)
 #' of a phylo2vec vector using the Queue Shuffle algorithm.
 #'
@@ -159,6 +182,16 @@ remove_branch_lengths <- function(newick) .Call(wrap__remove_branch_lengths, new
 #' @export
 remove_parent_labels <- function(newick) .Call(wrap__remove_parent_labels, newick)
 
+#' Get a precursor of the precision matrix of a phylo2vec tree
+#' The precision matrix is the inverse of the variance-covariance matrix.
+#' The precision matrix can be obtained using Schur's complement on this precursor.
+#' Output shape: (2 * (n_leaves - 1), 2 * (n_leaves- 1)]
+#'
+#' @param tree A phylo2vec tree
+#' @return A precursor of the precision matrix (shape: (2 * (n_leaves - 1), 2 * (n_leaves- 1)))
+#' @export
+pre_precision <- function(tree) .Call(wrap__pre_precision, tree)
+
 #' Remove a leaf from a phylo2vec vector
 #'
 #' @param vector phylo2vec vector representation of a tree topology
@@ -176,7 +209,7 @@ remove_leaf <- function(vector, leaf) .Call(wrap__remove_leaf, vector, leaf)
 #'
 #' @param v1 First tree as phylo2vec vector
 #' @param v2 Second tree as phylo2vec vector
-#' @param normalize If TRUE, return normalized distance in range [0.0, 1.0]
+#' @param normalize If TRUE, return normalized distance in range `[0.0, 1.0]`
 #' @return RF distance (numeric)
 #' @export
 robinson_foulds <- function(v1, v2, normalize) .Call(wrap__robinson_foulds, v1, v2, normalize)
@@ -214,9 +247,12 @@ to_ancestry <- function(vector) .Call(wrap__to_ancestry, vector)
 #' @export
 to_edges <- function(vector) .Call(wrap__to_edges, vector)
 
-to_newick_from_matrix <- function(matrix) .Call(wrap__to_newick_from_matrix, matrix)
-
-to_newick_from_vector <- function(vector) .Call(wrap__to_newick_from_vector, vector)
+#' Recover a rooted tree (in Newick format) from a phylo2vec tree
+#'
+#' @param tree A phylo2vec tree
+#' @return Newick string representation of the tree
+#' @export
+to_newick <- function(tree) .Call(wrap__to_newick, tree)
 
 #' Get pairs from a phylo2vec vector
 #'
@@ -229,9 +265,12 @@ to_matrix <- function(newick) .Call(wrap__to_matrix, newick)
 
 to_vector <- function(newick) .Call(wrap__to_vector, newick)
 
-vcov_from_matrix <- function(matrix) .Call(wrap__vcov_from_matrix, matrix)
-
-vcov_from_vector <- function(vector) .Call(wrap__vcov_from_vector, vector)
+#' Get the variance-covariance matrix of a phylo2vec tree
+#'
+#' @param tree A phylo2vec tree
+#' @return Variance-covariance matrix (shape: (n_leaves, n_leaves))
+#' @export
+vcovp <- function(tree) .Call(wrap__vcovp, tree)
 
 
 # nolint end
diff --git a/r-phylo2vec/R/stats.R b/r-phylo2vec/R/stats.R
index cb987c7f..134f13f4 100644
--- a/r-phylo2vec/R/stats.R
+++ b/r-phylo2vec/R/stats.R
@@ -1,41 +1,3 @@
-#' Compute the cophenetic distance matrix of a phylo2vec
-#' vector (tree topology) or matrix (topology + branch lengths).
-#'
-#' Output shape: (n_leaves, n_leaves)
-#'
-#' @param vector_or_matrix phylo2vec vector (ndim == 1)/matrix (ndim == 2)
-#' @param unrooted Logical indicating whether applying the cophenetic distance
-#'  to an unrooted tree or not (see ete3). Default is FALSE.
-#' @return Cophenetic distance matrix
-#' @export
-cophenetic_distances <- function(vector_or_matrix, unrooted = FALSE) {
-  if (is.vector(vector_or_matrix, "integer")) {
-    .Call(wrap__cophenetic_from_vector, vector_or_matrix, unrooted)
-  } else if (is.matrix(vector_or_matrix)) {
-    .Call(wrap__cophenetic_from_matrix, vector_or_matrix, unrooted)
-  } else {
-    stop("Input must be either an integer vector or a 2D matrix.")
-  }
-}
-
-#' Compute the variance-covariance matrix of a phylo2vec
-#' vector (tree topology) or matrix (topology + branch lengths).
-#'
-#' Output shape: (n_leaves, n_leaves)
-#'
-#' @param vector_or_matrix phylo2vec vector (ndim == 1)/matrix (ndim == 2)
-#' @return Variance-covariance matrix
-#' @export
-vcovp <- function(vector_or_matrix) {
-  if (is.vector(vector_or_matrix, "integer")) {
-    .Call(wrap__vcov_from_vector, vector_or_matrix)
-  } else if (is.matrix(vector_or_matrix)) {
-    .Call(wrap__vcov_from_matrix, vector_or_matrix)
-  } else {
-    stop("Input must be either an integer vector or a 2D matrix.")
-  }
-}
-
 #' Compute the precision matrix of a phylo2vec
 #' vector (tree topology) or matrix (topology + branch lengths).
 #'
@@ -47,14 +9,7 @@ vcovp <- function(vector_or_matrix) {
 #' @return Precision matrix
 #' @export
 precision <- function(vector_or_matrix) {
-  # Rust panics if len(v) == 0
-  if (is.vector(vector_or_matrix, "integer")) {
-    precursor <- .Call(wrap__pre_precision_from_vector, vector_or_matrix)
-  } else if (is.matrix(vector_or_matrix)) {
-    precursor <- .Call(wrap__pre_precision_from_matrix, vector_or_matrix)
-  } else {
-    stop("Input must be either an integer vector or a 2D matrix.")
-  }
+  precursor <- pre_precision(vector_or_matrix)
 
   # nrow(precursor) = ncol(precursor) = 2*k where k = n_leaves - 1
   n <- nrow(precursor)
@@ -72,6 +27,21 @@ precision <- function(vector_or_matrix) {
   a - b %*% solve(c, d)
 }
 
+#' Compute the cophenetic distance matrix of a phylo2vec tree.
+#'
+#' The cophenetic distance between two leaves is the distance from each leaf
+#' to their most recent common ancestor.
+#' For vectors, this is the topological distance.
+#' For matrices, this uses branch lengths.
+#'
+#' @param tree phylo2vec vector (1D) or matrix (2D)
+#' @param unrooted If TRUE, compute unrooted distances. Default is FALSE.
+#' @return Cophenetic distance matrix (shape: (n_leaves, n_leaves))
+#' @export
+cophenetic_distances <- function(tree, unrooted = FALSE) {
+  .Call(wrap__cophenetic_distances, tree, unrooted)
+}
+
 #' Compute the Robinson-Foulds distance between two trees.
 #'
 #' RF distance counts the number of bipartitions (splits) that differ
@@ -81,10 +51,10 @@ precision <- function(vector_or_matrix) {
 #'   Only topology is used; branch lengths are ignored.
 #' @param tree2 Second tree as phylo2vec vector (1D) or matrix (2D).
 #'   Only topology is used; branch lengths are ignored.
-#' @param normalize If TRUE, return normalized distance in range [0.0, 1.0].
+#' @param normalize If TRUE, return normalized distance in range `[0.0, 1.0]`.
 #'   Default is FALSE.
 #' @return RF distance (numeric). Integer value if normalize=FALSE,
-#'   float in [0,1] otherwise.
+#'   float in `[0, 1]` otherwise.
 #' @examples
 #' v1 <- sample_tree(10, topology_only = TRUE)
 #' v2 <- sample_tree(10, topology_only = TRUE)
diff --git a/r-phylo2vec/R/utils.R b/r-phylo2vec/R/utils.R
index e2207133..6c2810f3 100644
--- a/r-phylo2vec/R/utils.R
+++ b/r-phylo2vec/R/utils.R
@@ -15,78 +15,3 @@ sample_tree <- function(n_leaves, topology_only = FALSE, ordered = FALSE) {
     .Call(wrap__sample_matrix, n_leaves, ordered)
   }
 }
-
-#' Get the most recent common ancestor between two nodes in a phylo2vec tree.
-#'
-#' `node1` and `node2` can be leaf nodes (0 to n_leaves - 1)
-#' or internal nodes (n_leaves to 2 * n_leaves - 2).
-#'
-#' Similar to ape's `getMRCA` function (for leaf nodes)
-#' and ETE's `get_common_ancestor` (for all nodes), but for phylo2vec.
-#'
-#' Note: The MRCA is purely topological, so for matrices, only the vector
-#' component (column 1) is used.
-#'
-#' @param vector_or_matrix phylo2vec vector (integer) or matrix (numeric)
-#' @param node1 A node in the tree (0-indexed)
-#' @param node2 A node in the tree (0-indexed)
-#' @return The most recent common ancestor node (0-indexed)
-#' @export
-get_common_ancestor <- function(vector_or_matrix, node1, node2) {
-  if (is.vector(vector_or_matrix, "integer")) {
-    get_common_ancestor_from_vector(vector_or_matrix, node1, node2)
-  } else if (is.matrix(vector_or_matrix)) {
-    v <- as.integer(vector_or_matrix[, 1])
-    get_common_ancestor_from_vector(v, node1, node2)
-  } else {
-    stop("Input must be either an integer vector or a 2D matrix.")
-  }
-}
-
-#' Get the depth of a node in a phylo2vec tree.
-#'
-#' The depth of a node is the length of the path from the root to that node
-#' (i.e., distance from root). This follows the BEAST/ETE convention.
-#'
-#' The root has depth 0, and depths increase as you move toward the leaves.
-#'
-#' For vectors, topological depth is returned (all branch lengths = 1).
-#' For matrices, actual branch lengths are used.
-#'
-#' @param vector_or_matrix phylo2vec vector (ndim == 1)/matrix (ndim == 2)
-#' @param node A node in the tree (0-indexed, from 0 to 2*n_leaves - 2)
-#' @return Depth of the node (distance from root)
-#' @export
-get_node_depth <- function(vector_or_matrix, node) {
-  if (is.vector(vector_or_matrix, "integer")) {
-    get_node_depth_from_vector(vector_or_matrix, node)
-  } else if (is.matrix(vector_or_matrix)) {
-    get_node_depth_from_matrix(vector_or_matrix, node)
-  } else {
-    stop("Input must be either an integer vector or a 2D matrix.")
-  }
-}
-
-#' Get the depths of all nodes in a phylo2vec tree.
-#'
-#' The depth of a node is the length of the path from the root to that node
-#' (i.e., distance from root). This follows the BEAST/ETE convention.
-#'
-#' The root has depth 0, and depths increase as you move toward the leaves.
-#'
-#' For vectors, topological depth is returned (all branch lengths = 1).
-#' For matrices, actual branch lengths are used.
-#'
-#' @param vector_or_matrix phylo2vec vector (ndim == 1)/matrix (ndim == 2)
-#' @return Vector of depths for all nodes (length = 2 * n_leaves - 1).
-#'   Index i+1 contains the depth of node i (R is 1-indexed).
-#' @export
-get_node_depths <- function(vector_or_matrix) {
-  if (is.vector(vector_or_matrix, "integer")) {
-    get_node_depths_from_vector(vector_or_matrix)
-  } else if (is.matrix(vector_or_matrix)) {
-    get_node_depths_from_matrix(vector_or_matrix)
-  } else {
-    stop("Input must be either an integer vector or a 2D matrix.")
-  }
-}
diff --git a/r-phylo2vec/man/cophenetic_distances.Rd b/r-phylo2vec/man/cophenetic_distances.Rd
index f0e6adcf..2e940973 100644
--- a/r-phylo2vec/man/cophenetic_distances.Rd
+++ b/r-phylo2vec/man/cophenetic_distances.Rd
@@ -1,21 +1,29 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/stats.R
+% Please edit documentation in R/extendr-wrappers.R, R/stats.R
 \name{cophenetic_distances}
 \alias{cophenetic_distances}
-\title{Compute the cophenetic distance matrix of a phylo2vec
-vector (tree topology) or matrix (topology + branch lengths).}
+\title{Get the (topological) cophenetic distance matrix of a phylo2vec tree}
 \usage{
-cophenetic_distances(vector_or_matrix, unrooted = FALSE)
+cophenetic_distances(tree, unrooted = FALSE)
+
+cophenetic_distances(tree, unrooted = FALSE)
 }
 \arguments{
-\item{vector_or_matrix}{phylo2vec vector (ndim == 1)/matrix (ndim == 2)}
+\item{tree}{phylo2vec vector (1D) or matrix (2D)}
 
-\item{unrooted}{Logical indicating whether applying the cophenetic distance
-to an unrooted tree or not (see ete3). Default is FALSE.}
+\item{unrooted}{If TRUE, compute unrooted distances. Default is FALSE.}
 }
 \value{
-Cophenetic distance matrix
+Cophenetic distance matrix (shape: (n_leaves, n_leaves))
+
+Cophenetic distance matrix (shape: (n_leaves, n_leaves))
 }
 \description{
-Output shape: (n_leaves, n_leaves)
+The cophenetic distance between two leaves is the distance from each leaf to their most recent common ancestor.
+For phylo2vec vectors, this is the topological distance. For phylo2vec matrices, this is the distance with branch lengths.
+
+The cophenetic distance between two leaves is the distance from each leaf
+to their most recent common ancestor.
+For vectors, this is the topological distance.
+For matrices, this uses branch lengths.
 }
diff --git a/r-phylo2vec/man/get_common_ancestor.Rd b/r-phylo2vec/man/get_common_ancestor.Rd
index 19963c40..2cfe87a7 100644
--- a/r-phylo2vec/man/get_common_ancestor.Rd
+++ b/r-phylo2vec/man/get_common_ancestor.Rd
@@ -1,29 +1,27 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils.R
+% Please edit documentation in R/extendr-wrappers.R
 \name{get_common_ancestor}
 \alias{get_common_ancestor}
-\title{Get the most recent common ancestor between two nodes in a phylo2vec tree.}
+\title{Get the first recent common ancestor between two nodes in a phylo2vec tree
+node1 and node2 can be leaf nodes (0 to n_leaves) or internal nodes (n_leaves to 2*(n_leaves-1)).
+Similar to ape's \code{getMRCA} function in R (for leaf nodes)
+and ETE's \code{get_common_ancestor} in Python (for all nodes), but for phylo2vec vectors.}
 \usage{
-get_common_ancestor(vector_or_matrix, node1, node2)
+get_common_ancestor(tree, node1, node2)
 }
 \arguments{
-\item{vector_or_matrix}{phylo2vec vector (integer) or matrix (numeric)}
+\item{node1}{The first node (0-indexed)}
 
-\item{node1}{A node in the tree (0-indexed)}
+\item{node2}{The second node (0-indexed)}
 
-\item{node2}{A node in the tree (0-indexed)}
+\item{vector}{phylo2vec vector representation of a tree topology}
 }
 \value{
-The most recent common ancestor node (0-indexed)
+The common ancestor node (0-indexed)
 }
 \description{
-\code{node1} and \code{node2} can be leaf nodes (0 to n_leaves - 1)
-or internal nodes (n_leaves to 2 * n_leaves - 2).
-}
-\details{
-Similar to ape's \code{getMRCA} function (for leaf nodes)
-and ETE's \code{get_common_ancestor} (for all nodes), but for phylo2vec.
-
-Note: The MRCA is purely topological, so for matrices, only the vector
-component (column 1) is used.
+Get the first recent common ancestor between two nodes in a phylo2vec tree
+node1 and node2 can be leaf nodes (0 to n_leaves) or internal nodes (n_leaves to 2*(n_leaves-1)).
+Similar to ape's \code{getMRCA} function in R (for leaf nodes)
+and ETE's \code{get_common_ancestor} in Python (for all nodes), but for phylo2vec vectors.
 }
diff --git a/r-phylo2vec/man/get_node_depth.Rd b/r-phylo2vec/man/get_node_depth.Rd
index 4b732113..794e96a7 100644
--- a/r-phylo2vec/man/get_node_depth.Rd
+++ b/r-phylo2vec/man/get_node_depth.Rd
@@ -1,26 +1,23 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils.R
+% Please edit documentation in R/extendr-wrappers.R
 \name{get_node_depth}
 \alias{get_node_depth}
-\title{Get the depth of a node in a phylo2vec tree.}
+\title{Get the depth of a node in a phylo2vec tree
+Depth is the distance from root to that node. Root has depth 0.
+For vectors, this is the topological depth. For matrices, this is the depth with branch lengths.}
 \usage{
-get_node_depth(vector_or_matrix, node)
+get_node_depth(tree, node)
 }
 \arguments{
-\item{vector_or_matrix}{phylo2vec vector (ndim == 1)/matrix (ndim == 2)}
+\item{tree}{A phylo2vec tree}
 
-\item{node}{A node in the tree (0-indexed, from 0 to 2*n_leaves - 2)}
+\item{node}{The node to get the depth of (0-indexed)}
 }
 \value{
-Depth of the node (distance from root)
+The depth of the node
 }
 \description{
-The depth of a node is the length of the path from the root to that node
-(i.e., distance from root). This follows the BEAST/ETE convention.
-}
-\details{
-The root has depth 0, and depths increase as you move toward the leaves.
-
-For vectors, topological depth is returned (all branch lengths = 1).
-For matrices, actual branch lengths are used.
+Get the depth of a node in a phylo2vec tree
+Depth is the distance from root to that node. Root has depth 0.
+For vectors, this is the topological depth. For matrices, this is the depth with branch lengths.
 }
diff --git a/r-phylo2vec/man/get_node_depths.Rd b/r-phylo2vec/man/get_node_depths.Rd
index 91d19ded..c51034c6 100644
--- a/r-phylo2vec/man/get_node_depths.Rd
+++ b/r-phylo2vec/man/get_node_depths.Rd
@@ -1,25 +1,19 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils.R
+% Please edit documentation in R/extendr-wrappers.R
 \name{get_node_depths}
 \alias{get_node_depths}
-\title{Get the depths of all nodes in a phylo2vec tree.}
+\title{Get the depths of all nodes in a phylo2vec tree
+Depth is the distance from root to each node. Root has depth 0.}
 \usage{
-get_node_depths(vector_or_matrix)
+get_node_depths(tree)
 }
 \arguments{
-\item{vector_or_matrix}{phylo2vec vector (ndim == 1)/matrix (ndim == 2)}
+\item{tree}{A phylo2vec tree}
 }
 \value{
-Vector of depths for all nodes (length = 2 * n_leaves - 1).
-Index i+1 contains the depth of node i (R is 1-indexed).
+A vector of depths for all nodes (length: 2*n_leaves - 1)
 }
 \description{
-The depth of a node is the length of the path from the root to that node
-(i.e., distance from root). This follows the BEAST/ETE convention.
-}
-\details{
-The root has depth 0, and depths increase as you move toward the leaves.
-
-For vectors, topological depth is returned (all branch lengths = 1).
-For matrices, actual branch lengths are used.
+Get the depths of all nodes in a phylo2vec tree
+Depth is the distance from root to each node. Root has depth 0.
 }
diff --git a/r-phylo2vec/man/pre_precision.Rd b/r-phylo2vec/man/pre_precision.Rd
new file mode 100644
index 00000000..11aa400b
--- /dev/null
+++ b/r-phylo2vec/man/pre_precision.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/extendr-wrappers.R
+\name{pre_precision}
+\alias{pre_precision}
+\title{Get a precursor of the precision matrix of a phylo2vec tree
+The precision matrix is the inverse of the variance-covariance matrix.
+The precision matrix can be obtained using Schur's complement on this precursor.
+Output shape: (2 * (n_leaves - 1), 2 * (n_leaves- 1)]}
+\usage{
+pre_precision(tree)
+}
+\arguments{
+\item{tree}{A phylo2vec tree}
+}
+\value{
+A precursor of the precision matrix (shape: (2 * (n_leaves - 1), 2 * (n_leaves- 1)))
+}
+\description{
+Get a precursor of the precision matrix of a phylo2vec tree
+The precision matrix is the inverse of the variance-covariance matrix.
+The precision matrix can be obtained using Schur's complement on this precursor.
+Output shape: (2 * (n_leaves - 1), 2 * (n_leaves- 1)]
+}
diff --git a/r-phylo2vec/man/robinson_foulds.Rd b/r-phylo2vec/man/robinson_foulds.Rd
index f89487a9..411874ca 100644
--- a/r-phylo2vec/man/robinson_foulds.Rd
+++ b/r-phylo2vec/man/robinson_foulds.Rd
@@ -15,7 +15,7 @@ Only topology is used; branch lengths are ignored.}
 \item{tree2}{Second tree as phylo2vec vector (1D) or matrix (2D).
 Only topology is used; branch lengths are ignored.}
 
-\item{normalize}{If TRUE, return normalized distance in range \link{0.0, 1.0}.
+\item{normalize}{If TRUE, return normalized distance in range \verb{[0.0, 1.0]}.
 Default is FALSE.}
 
 \item{v1}{First tree as phylo2vec vector}
@@ -26,7 +26,7 @@ Default is FALSE.}
 RF distance (numeric)
 
 RF distance (numeric). Integer value if normalize=FALSE,
-float in \link{0,1} otherwise.
+float in \verb{[0, 1]} otherwise.
 }
 \description{
 RF distance counts the number of bipartitions (splits) that differ
diff --git a/r-phylo2vec/man/to_newick.Rd b/r-phylo2vec/man/to_newick.Rd
index fdc57d2c..d8808758 100644
--- a/r-phylo2vec/man/to_newick.Rd
+++ b/r-phylo2vec/man/to_newick.Rd
@@ -1,18 +1,17 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/base_newick.R
+% Please edit documentation in R/extendr-wrappers.R
 \name{to_newick}
 \alias{to_newick}
-\title{Convert a phylo2vec vector or matrix to Newick format}
+\title{Recover a rooted tree (in Newick format) from a phylo2vec tree}
 \usage{
-to_newick(vector_or_matrix)
+to_newick(tree)
 }
 \arguments{
-\item{vector_or_matrix}{phylo2vec matrix (tree with branch
-lengths) or a vector (tree topology only)}
+\item{tree}{A phylo2vec tree}
 }
 \value{
 Newick string representation of the tree
 }
 \description{
-Convert a phylo2vec vector or matrix to Newick format
+Recover a rooted tree (in Newick format) from a phylo2vec tree
 }
diff --git a/r-phylo2vec/man/vcovp.Rd b/r-phylo2vec/man/vcovp.Rd
index 1bbc4c62..a403d914 100644
--- a/r-phylo2vec/man/vcovp.Rd
+++ b/r-phylo2vec/man/vcovp.Rd
@@ -1,18 +1,17 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/stats.R
+% Please edit documentation in R/extendr-wrappers.R
 \name{vcovp}
 \alias{vcovp}
-\title{Compute the variance-covariance matrix of a phylo2vec
-vector (tree topology) or matrix (topology + branch lengths).}
+\title{Get the variance-covariance matrix of a phylo2vec tree}
 \usage{
-vcovp(vector_or_matrix)
+vcovp(tree)
 }
 \arguments{
-\item{vector_or_matrix}{phylo2vec vector (ndim == 1)/matrix (ndim == 2)}
+\item{tree}{A phylo2vec tree}
 }
 \value{
-Variance-covariance matrix
+Variance-covariance matrix (shape: (n_leaves, n_leaves))
 }
 \description{
-Output shape: (n_leaves, n_leaves)
+Get the variance-covariance matrix of a phylo2vec tree
 }
diff --git a/r-phylo2vec/src/rust/src/lib.rs b/r-phylo2vec/src/rust/src/lib.rs
index 6b8ec66d..f31e92e3 100644
--- a/r-phylo2vec/src/rust/src/lib.rs
+++ b/r-phylo2vec/src/rust/src/lib.rs
@@ -2,7 +2,7 @@ use std::collections::HashMap;
 use std::result::Result;
 
 use extendr_api::prelude::*;
-use ndarray::{Array2, ArrayView2};
+use ndarray::Array2;
 
 use phylo2vec::matrix::base as mbase;
 use phylo2vec::matrix::convert as mconvert;
@@ -15,6 +15,23 @@ use phylo2vec::vector::distance as vdist;
 use phylo2vec::vector::graph as vgraph;
 use phylo2vec::vector::ops as vops;
 
+enum RTree {
+    Matrix(RMatrix<f64>),
+    Vector(Vec<i32>),
+}
+
+impl TryFrom<Robj> for RTree {
+    type Error = Error;
+
+    fn try_from(robj: Robj) -> Result<Self, Self::Error> {
+        if robj.is_matrix() {
+            Ok(RTree::Matrix(robj.try_into()?))
+        } else {
+            Ok(RTree::Vector(robj.try_into()?))
+        }
+    }
+}
+
 // Convert Vec<i32> to Vec<usize>
 fn as_usize(v: Vec<i32>) -> Vec<usize> {
     v.iter().map(|&x| x as usize).collect()
@@ -83,18 +100,23 @@ fn sample_matrix(n_leaves: isize, ordered: bool) -> Result<RMatrix<f64>, Error>
     }
 }
 
-// Recover a rooted tree (in Newick format) from a phylo2vec vector
-#[extendr]
-fn to_newick_from_vector(vector: Vec<i32>) -> String {
-    let v_usize = as_usize(vector);
-    vconvert::to_newick(&v_usize)
-}
-
-// Recover a rooted tree (in Newick format) from a phylo2vec matrix
+/// Recover a rooted tree (in Newick format) from a phylo2vec tree
+///
+/// @param tree A phylo2vec tree
+/// @return Newick string representation of the tree
+/// @export
 #[extendr]
-fn to_newick_from_matrix(matrix: RMatrix<f64>) -> String {
-    let matrix = convert_from_rmatrix(&matrix).unwrap();
-    mconvert::to_newick(&matrix.view())
+fn to_newick(tree: RTree) -> String {
+    match tree {
+        RTree::Matrix(m) => {
+            let matrix = convert_from_rmatrix(&m).unwrap();
+            mconvert::to_newick(&matrix.view())
+        }
+        RTree::Vector(v) => {
+            let v_usize = as_usize(v);
+            vconvert::to_newick(&v_usize)
+        }
+    }
 }
 
 // Convert a newick string to a phylo2vec vector
@@ -133,7 +155,6 @@ fn check_v(vector: Vec<i32>) {
 #[extendr]
 fn check_m(vector: RMatrix<f64>) {
     let matrix = convert_from_rmatrix(&vector).unwrap();
-
     mbase::check_m(&matrix.view());
 }
 
@@ -280,48 +301,72 @@ fn remove_leaf(vector: Vec<i32>, leaf: i32) -> Robj {
 /// Similar to ape's `getMRCA` function in R (for leaf nodes)
 /// and ETE's `get_common_ancestor` in Python (for all nodes), but for phylo2vec vectors.
 ///
-/// @param vector phylo2vec vector representation of a tree topology
+/// @param tree A phylo2vec tree
 /// @param node1 The first node (0-indexed)
 /// @param node2 The second node (0-indexed)
 /// @return The common ancestor node (0-indexed)
 /// @export
 #[extendr]
-fn get_common_ancestor(vector: Vec<i32>, node1: i32, node2: i32) -> i32 {
-    let v_usize: Vec<usize> = as_usize(vector);
-    let common_ancestor = vops::get_common_ancestor(&v_usize, node1 as usize, node2 as usize);
-    common_ancestor as i32
-}
-
-// Get the topological depth of a node in a phylo2vec vector
-// Depth is the distance from root to that node. Root has depth 0.
-#[extendr]
-fn get_node_depth_from_vector(vector: Vec<i32>, node: i32) -> f64 {
-    let v_usize: Vec<usize> = as_usize(vector);
-    vops::get_node_depth(&v_usize, node as usize)
-}
-
-// Get the depth of a node in a phylo2vec matrix
-// Depth is the distance from root to that node. Root has depth 0.
-#[extendr]
-fn get_node_depth_from_matrix(matrix: RMatrix<f64>, node: i32) -> f64 {
-    let matrix_rs = convert_from_rmatrix(&matrix).unwrap();
-    mops::get_node_depth(&matrix_rs.view(), node as usize)
+fn get_common_ancestor(tree: RTree, node1: i32, node2: i32) -> Result<i32, Error> {
+    let node1 = node1 as usize;
+    let node2 = node2 as usize;
+    let mrca = match tree {
+        RTree::Matrix(m) => {
+            let matrix = convert_from_rmatrix(&m).unwrap();
+            let (v, _) = mbase::parse_matrix(&matrix.view());
+            vops::get_common_ancestor(&v, node1, node2)
+        }
+        RTree::Vector(v) => {
+            let v_usize: Vec<usize> = as_usize(v);
+            vops::get_common_ancestor(&v_usize, node1, node2)
+        }
+    };
+    mrca.map(|node| node as i32)
+        .map_err(|e| Error::OutOfRange(e.into()))
 }
 
-// Get the depths of all nodes in a phylo2vec vector (topological)
-// Depth is the distance from root to each node. Root has depth 0.
+/// Get the depth of a node in a phylo2vec tree
+/// Depth is the distance from root to that node. Root has depth 0.
+/// For vectors, this is the topological depth. For matrices, this is the depth with branch lengths.
+///
+/// @param tree A phylo2vec tree
+/// @param node The node to get the depth of (0-indexed)
+/// @return The depth of the node
+/// @export
 #[extendr]
-fn get_node_depths_from_vector(vector: Vec<i32>) -> Vec<f64> {
-    let v_usize: Vec<usize> = as_usize(vector);
-    vops::get_node_depths(&v_usize)
+fn get_node_depth(tree: RTree, node: i32) -> Result<f64, Error> {
+    let node_usize = node as usize;
+    let depth = match tree {
+        RTree::Matrix(m) => {
+            let matrix = convert_from_rmatrix(&m).unwrap();
+            mops::get_node_depth(&matrix.view(), node_usize)
+        }
+        RTree::Vector(v) => {
+            let v_usize: Vec<usize> = as_usize(v);
+            vops::get_node_depth(&v_usize, node_usize)
+        }
+    };
+    depth.map_err(|e| Error::OutOfRange(e.into()))
 }
 
-// Get the depths of all nodes in a phylo2vec matrix
-// Depth is the distance from root to each node. Root has depth 0.
+/// Get the depths of all nodes in a phylo2vec tree
+/// Depth is the distance from root to each node. Root has depth 0.
+///
+/// @param tree A phylo2vec tree
+/// @return A vector of depths for all nodes (length: 2*n_leaves - 1)
+/// @export
 #[extendr]
-fn get_node_depths_from_matrix(matrix: RMatrix<f64>) -> Vec<f64> {
-    let matrix_rs = convert_from_rmatrix(&matrix).unwrap();
-    mops::get_node_depths(&matrix_rs.view())
+fn get_node_depths(tree: RTree) -> Vec<f64> {
+    match tree {
+        RTree::Matrix(m) => {
+            let matrix = convert_from_rmatrix(&m).unwrap();
+            mops::get_node_depths(&matrix.view())
+        }
+        RTree::Vector(v) => {
+            let v_usize: Vec<usize> = as_usize(v);
+            vops::get_node_depths(&v_usize)
+        }
+    }
 }
 
 /// Produce an ordered version (i.e., birth-death process version)
@@ -439,23 +484,27 @@ fn remove_parent_labels(newick: &str) -> String {
     newick::remove_parent_labels(newick)
 }
 
-// Get the topological cophenetic distance matrix of a phylo2vec vector
-#[extendr]
-fn cophenetic_from_vector(vector: Vec<i32>, unrooted: bool) -> RMatrix<i32> {
-    let v_usize: Vec<usize> = as_usize(vector);
-    let k = v_usize.len();
-    let coph_rs = vgraph::cophenetic_distances(&v_usize, unrooted);
-    let mut coph_r = RMatrix::new_matrix(k + 1, k + 1, |r, c| coph_rs[[r, c]] as i32);
-    let dimnames = (0..=k).map(|x| x as i32).collect::<Vec<i32>>();
-    coph_r.set_dimnames(List::from_values(vec![dimnames.clone(), dimnames]));
-
-    coph_r
-}
-
-// Get the cophenetic distance matrix of a phylo2vec matrix
+/// Get the (topological) cophenetic distance matrix of a phylo2vec tree
+///
+/// The cophenetic distance between two leaves is the distance from each leaf to their most recent common ancestor.
+/// For phylo2vec vectors, this is the topological distance. For phylo2vec matrices, this is the distance with branch lengths.
+///
+/// @param tree A phylo2vec tree
+/// @param unrooted If true, the distance is calculated as the distance between each leaf and their most recent common ancestor, multiplied by 2. If false, the distance is calculated as the distance from each leaf to their most recent common ancestor.
+/// @return Cophenetic distance matrix (shape: (n_leaves, n_leaves))
+/// @export
 #[extendr]
-fn cophenetic_from_matrix(matrix: ArrayView2<f64>, unrooted: bool) -> RMatrix<f64> {
-    let coph_rs = mgraph::cophenetic_distances(&matrix.view(), unrooted);
+fn cophenetic_distances(tree: RTree, #[default = "FALSE"] unrooted: bool) -> RMatrix<f64> {
+    let coph_rs = match tree {
+        RTree::Matrix(m) => {
+            let matrix = convert_from_rmatrix(&m).unwrap();
+            mgraph::cophenetic_distances(&matrix.view(), unrooted)
+        }
+        RTree::Vector(v) => {
+            let v_usize: Vec<usize> = as_usize(v);
+            vgraph::cophenetic_distances(&v_usize, unrooted)
+        }
+    };
     let n_leaves = coph_rs.shape()[0];
     let mut coph_r = RMatrix::new_matrix(n_leaves, n_leaves, |r, c| coph_rs[[r, c]]);
     let dimnames = (0..n_leaves).map(|x| x as i32).collect::<Vec<i32>>();
@@ -464,14 +513,26 @@ fn cophenetic_from_matrix(matrix: ArrayView2<f64>, unrooted: bool) -> RMatrix<f6
     coph_r
 }
 
-// Get a precursor of the precision matrix of a phylo2vec vector
-// The precision matrix is the inverse of the variance-covariance matrix.
-// The precision matrix can be obtained using Schur's complement on this precursor.
-// Output shape: (2 * (n_leaves - 1), 2 * (n_leaves- 1)]
+/// Get a precursor of the precision matrix of a phylo2vec tree
+/// The precision matrix is the inverse of the variance-covariance matrix.
+/// The precision matrix can be obtained using Schur's complement on this precursor.
+/// Output shape: (2 * (n_leaves - 1), 2 * (n_leaves- 1)]
+///
+/// @param tree A phylo2vec tree
+/// @return A precursor of the precision matrix (shape: (2 * (n_leaves - 1), 2 * (n_leaves- 1)))
+/// @export
 #[extendr]
-fn pre_precision_from_vector(vector: Vec<i32>) -> RMatrix<f64> {
-    let v_usize: Vec<usize> = as_usize(vector);
-    let pre_precision_rs = vgraph::pre_precision(&v_usize);
+fn pre_precision(tree: RTree) -> RMatrix<f64> {
+    let pre_precision_rs = match tree {
+        RTree::Matrix(m) => {
+            let matrix = convert_from_rmatrix(&m).unwrap();
+            mgraph::pre_precision(&matrix.view())
+        }
+        RTree::Vector(v) => {
+            let v_usize: Vec<usize> = as_usize(v);
+            vgraph::pre_precision(&v_usize)
+        }
+    };
     let n = pre_precision_rs.shape()[0];
     let mut preprecision_r = RMatrix::new_matrix(n, n, |r, c| pre_precision_rs[[r, c]]);
     let dimnames = (0..n).map(|x| x as i32).collect::<Vec<i32>>();
@@ -480,27 +541,23 @@ fn pre_precision_from_vector(vector: Vec<i32>) -> RMatrix<f64> {
     preprecision_r
 }
 
-// Get the precision matrix of a phylo2vec matrix
-// The precision matrix is the inverse of the variance-covariance matrix.
-// The precision matrix can be obtained using Schur's complement on this precursor.
-// Output shape: (2 * (n_leaves - 1), 2 * (n_leaves- 1)]
-#[extendr]
-fn pre_precision_from_matrix(matrix: RMatrix<f64>) -> RMatrix<f64> {
-    let matrix_rs = convert_from_rmatrix(&matrix).unwrap();
-    let pre_precision_rs = mgraph::pre_precision(&matrix_rs.view());
-    let n = pre_precision_rs.shape()[0];
-    let mut pre_precision_r = RMatrix::new_matrix(n, n, |r, c| pre_precision_rs[[r, c]]);
-    let dimnames = (0..n).map(|x| x as i32).collect::<Vec<i32>>();
-    pre_precision_r.set_dimnames(List::from_values(vec![dimnames.clone(), dimnames]));
-
-    pre_precision_r
-}
-
-// Get the variance-covariance matrix of a phylo2vec vector
+/// Get the variance-covariance matrix of a phylo2vec tree
+///
+/// @param tree A phylo2vec tree
+/// @return Variance-covariance matrix (shape: (n_leaves, n_leaves))
+/// @export
 #[extendr]
-fn vcov_from_vector(vector: Vec<i32>) -> RMatrix<f64> {
-    let v_usize: Vec<usize> = as_usize(vector);
-    let vcv_rs = vgraph::vcv(&v_usize);
+fn vcovp(tree: RTree) -> RMatrix<f64> {
+    let vcv_rs = match tree {
+        RTree::Matrix(m) => {
+            let matrix = convert_from_rmatrix(&m).unwrap();
+            mgraph::vcv(&matrix.view())
+        }
+        RTree::Vector(v) => {
+            let v_usize: Vec<usize> = as_usize(v);
+            vgraph::vcv(&v_usize)
+        }
+    };
     let n_leaves = vcv_rs.shape()[0];
     let mut vcv_r = RMatrix::new_matrix(n_leaves, n_leaves, |r, c| vcv_rs[[r, c]]);
     let dimnames = (0..n_leaves).map(|x| x as i32).collect::<Vec<i32>>();
@@ -509,19 +566,6 @@ fn vcov_from_vector(vector: Vec<i32>) -> RMatrix<f64> {
     vcv_r
 }
 
-// Get the variance-covariance matrix of a phylo2vec matrix
-#[extendr]
-fn vcov_from_matrix(matrix: RMatrix<f64>) -> RMatrix<f64> {
-    let matrix_rs = convert_from_rmatrix(&matrix).unwrap();
-    let vcv_rs = mgraph::vcv(&matrix_rs.view());
-    let n_leaves = vcv_rs.shape()[0];
-    let mut vcv_matrix = RMatrix::new_matrix(n_leaves, n_leaves, |r, c| vcv_rs[[r, c]]);
-    let dimnames = (0..n_leaves).map(|x| x as i32).collect::<Vec<i32>>();
-    vcv_matrix.set_dimnames(List::from_values(vec![dimnames.clone(), dimnames]));
-
-    vcv_matrix
-}
-
 // Get the oriented incidence matrix of a phylo2vec vector in dense format
 #[extendr]
 fn incidence_dense(input_vector: Vec<i32>) -> RMatrix<i32> {
@@ -580,7 +624,7 @@ fn incidence_csr(input_vector: Vec<i32>) -> List {
 ///
 /// @param v1 First tree as phylo2vec vector
 /// @param v2 Second tree as phylo2vec vector
-/// @param normalize If TRUE, return normalized distance in range [0.0, 1.0]
+/// @param normalize If TRUE, return normalized distance in range `[0.0, 1.0]`
 /// @return RF distance (numeric)
 /// @export
 #[extendr]
@@ -599,39 +643,33 @@ extendr_module! {
     fn apply_label_mapping;
     fn check_m;
     fn check_v;
-    fn cophenetic_from_matrix;
-    fn cophenetic_from_vector;
+    fn cophenetic_distances;
     fn create_label_mapping;
     fn find_num_leaves;
     fn from_ancestry;
     fn from_edges;
     fn from_pairs;
     fn get_common_ancestor;
-    fn get_node_depth_from_matrix;
-    fn get_node_depth_from_vector;
-    fn get_node_depths_from_matrix;
-    fn get_node_depths_from_vector;
+    fn get_node_depth;
+    fn get_node_depths;
     fn has_branch_lengths;
     fn incidence_coo;
     fn incidence_csc;
     fn incidence_csr;
     fn incidence_dense;
-    fn pre_precision_from_matrix;
-    fn pre_precision_from_vector;
     fn queue_shuffle;
     fn remove_branch_lengths;
     fn remove_parent_labels;
+    fn pre_precision;
     fn remove_leaf;
     fn robinson_foulds;
     fn sample_matrix;
     fn sample_vector;
     fn to_ancestry;
     fn to_edges;
-    fn to_newick_from_matrix;
-    fn to_newick_from_vector;
+    fn to_newick;
     fn to_pairs;
     fn to_matrix;
     fn to_vector;
-    fn vcov_from_matrix;
-    fn vcov_from_vector;
+    fn vcovp;
 }
diff --git a/r-phylo2vec/tests/testthat/test_stats.R b/r-phylo2vec/tests/testthat/test_stats.R
index d04dd12b..b315e866 100644
--- a/r-phylo2vec/tests/testthat/test_stats.R
+++ b/r-phylo2vec/tests/testthat/test_stats.R
@@ -107,7 +107,7 @@ test_that(desc = "vcv_and_precision_vector", {
 
       # Test that the variance-covariance and precision matrices are inverses
       # vcov_from_vector = same as vcovp for vectors
-      vcv_p2v <- vcov_from_vector(v)
+      vcv_p2v <- vcovp(v)
       prec_p2v <- precision(v)
       identity <- diag(nrow(vcv_p2v))
 
@@ -124,7 +124,7 @@ test_that(desc = "vcv_and_precision_matrix", {
 
       # Test that the variance-covariance and precision matrices are inverses
       # vcov_from_vector = same as vcovp for matrices
-      vcv_p2v <- vcov_from_matrix(m)
+      vcv_p2v <- vcovp(m)
       prec_p2v <- precision(m)
       identity <- diag(nrow(vcv_p2v))