diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 3bb325a..886a077 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -19,6 +19,7 @@
 ### Improvements
 
 * The `README` has been improved to describe the expanding capabilities of `haarpy` [(#39)](https://github.com/polyquantique/haarpy/pull/39).
+* Added references and examples to the docstrings. Slight modification to the docstring format  [(#47)](https://github.com/polyquantique/haarpy/pull/47).
 
 ### Bug fixes
 
diff --git a/haarpy/circular_ensembles.py b/haarpy/circular_ensembles.py
index 3c9eee2..a668d5a 100644
--- a/haarpy/circular_ensembles.py
+++ b/haarpy/circular_ensembles.py
@@ -13,6 +13,11 @@
 # limitations under the License.
 """
 Circular ensembles Python interface
+
+References
+----------
+    [1] Matsumoto, S. (2013). Weingarten calculus for matrix ensembles associated with compact
+        symmetric spaces. arXiv preprint arXiv:1301.5401.
 """
 
 from math import prod
@@ -36,14 +41,35 @@ def weingarten_circular_orthogonal(
     permutation: Union[Permutation, tuple[int]],
     coe_dimension: Symbol,
 ) -> Expr:
-    """Returns the circular orthogonal ensembles Weingarten functions
-
-    Args:
-        permutation (Permutation): A permutation of S_2k or its coset-type
-        coe_dimension (int): The dimension of the COE
-
-    Returns:
-        Expr: The Weingarten function
+    """Returns the circular orthogonal ensemble's Weingarten functions
+
+    Parameters
+    ----------
+        permutation (Permutation) : A permutation of S_2k or its coset-type
+        coe_dimension (Symbol) : The dimension of the COE
+
+    Returns
+    -------
+        Expr : The Weingarten function
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from sympy.combinatorics import Permutation
+        >>> from haarpy import weingarten_circular_orthogonal
+        >>> d = Symbol("d")
+        >>> weingarten_circular_orthogonal(Permutation(3)(0,1), 4)
+        Fraction(3, 70)
+        >>> weingarten_circular_orthogonal(Permutation(3)(0,1), d)
+        (d + 2)/(d*(d + 1)*(d + 3))
+        >>> weingarten_circular_orthogonal((1,1), d)
+        (d + 2)/(d*(d + 1)*(d + 3))
+
+        Where (1,1) is the coset-type of Permutation(3)(0,1)
+
+    See Also
+    --------
+        coset_type, weingarten_orthogonal
     """
     return weingarten_orthogonal(permutation, coe_dimension + 1)
 
@@ -52,12 +78,29 @@ def weingarten_circular_orthogonal(
 def weingarten_circular_symplectic(permutation: Permutation, cse_dimension: Symbol) -> Expr:
     """Returns the circular symplectic ensembles Weingarten functions
 
-    Args:
-        permutation (Permutation): A permutation of the symmetric group S_2k
-        cse_dimension (int): The dimension of the CSE
-
-    Returns:
-        Expr: The Weingarten function
+    Parameters
+    ----------
+        permutation (Permutation) : A permutation of the symmetric group S_2k
+        cse_dimension (int) : The dimension of the CSE
+
+    Returns
+    -------
+        Expr : The Weingarten function
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from sympy.combinatorics import Permutation
+        >>> from haarpy import weingarten_circular_symplectic
+        >>> d = Symbol("d")
+        >>> weingarten_circular_symplectic(Permutation(3)(0,1), 4)
+        Fraction(-3, 140)
+        >>> weingarten_circular_symplectic(Permutation(3)(0,1), d)
+        (1 - d)/(d*(2*d - 3)*(2*d - 1))
+
+    See Also
+    --------
+        weingarten_symplectic
     """
     symplectic_dimension = (
         (2 * cse_dimension - 1) / 2
@@ -74,16 +117,34 @@ def haar_integral_circular_orthogonal(
     """Returns integral over circular orthogonal ensemble polynomial
     sampled at random from the Haar measure
 
-    Args:
+    Parameters
+    ----------
         sequences (tuple[tuple[int]]) : Indices of matrix elements
         group_dimension (Symbol) : Dimension of the orthogonal group
 
-    Returns:
+    Returns
+    -------
         Expr : Integral under the Haar measure
 
-    Raise:
-        ValueError : If sequences doesn't contain 2 tuples
-        ValueError : If tuples i and j are of odd size
+    Raise
+    -----
+        ValueError : if sequences do not contain 2 tuples
+        ValueError : if tuples i and j are of odd size
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from haarpy import haar_integral_circular_orthogonal
+        >>> d = Symbol("d")
+        >>> seq_i, seq_j = (0, 0, 1, 2), (1, 0, 0, 2)
+        >>> haar_integral_circular_orthogonal((seq_i, seq_j), 7)
+        Fraction(-1, 280)
+        >>> haar_integral_circular_orthogonal((seq_i, seq_j), d)
+        -2/(d*(d + 1)*(d + 3))
+
+    See Also
+    --------
+        coset_type, stabilizer_coset, weingarten_circular_orthogonal
     """
     if len(sequences) != 2:
         raise ValueError("Wrong tuple format")
@@ -114,21 +175,40 @@ def haar_integral_circular_symplectic(sequences: tuple[tuple[Expr]], half_dimens
     """Returns integral over circular symplectic ensemble polynomial
     sampled at random from the Haar measure
 
-    Args:
+    Parameters
+    ----------
         sequences (tuple[tuple[int]]) : Indices of matrix elements
         half_dimension (Symbol) : Half the dimension of the unitary group
 
-    Returns:
+    Returns
+    -------
         Expr : Integral under the Haar measure
 
-    Raise:
-        ValueError : If sequences doesn't contain 2 tuples
-        ValueError : If tuples i and j are of odd size
-        TypeError: If dimension is int and sequence is not
-        TypeError: If the half_dimension is not int nor Symbol
-        ValueError: If all sequence indices are not between 0 and 2*dimension - 1
-        TypeError: If sequence containt something else than Expr
-        TypeError: If symbolic sequences have the wrong format
+    Raise
+    -----
+        ValueError : if sequences do not contain 2 tuples
+        ValueError : if tuples i and j are of odd size
+        TypeError : if dimension is int and sequence is not
+        TypeError : if the half_dimension is neither int nor Symbol
+        ValueError : if all sequence indices are not between 0 and 2*dimension - 1
+        TypeError : if sequence contains something other than Expr
+        TypeError : if symbolic sequences have the wrong format
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from haarpy import haar_integral_circular_symplectic
+        >>> d = Symbol("d")
+        >>> seq_i_num, seq_j_num = (0, 3, 2, 1), (0, 1, 2, 3)
+        >>> haar_integral_circular_symplectic((seq_i_num, seq_j_num), 2)
+        Fraction(1, 6)
+        >>> seq_i_symb, seq_j_symb = (0, d+1, d, 1), (0, 1, d, d + 1)
+        >>> haar_integral_circular_symplectic((seq_i_symb, seq_j_symb), d)
+        -1/(2*d*(2*d - 3)*(2*d - 1))
+
+    See Also
+    --------
+        weingarten_circular_symplectic
     """
     if len(sequences) != 2:
         raise ValueError("Wrong sequence format")
diff --git a/haarpy/orthogonal.py b/haarpy/orthogonal.py
index 2ac540c..d9c86ab 100644
--- a/haarpy/orthogonal.py
+++ b/haarpy/orthogonal.py
@@ -13,6 +13,14 @@
 # limitations under the License.
 """
 Orthogonal group Python interface
+
+References
+----------
+    [1] Collins, B., & Śniady, P. (2006). Integration with respect to the Haar measure on unitary,
+    orthogonal and symplectic group. Communications in Mathematical Physics, 264(3), 773-795.
+    [2] Matsumoto, S. (2013). Weingarten calculus for matrix ensembles associated with compact
+    symmetric spaces. arXiv preprint arXiv:1301.5401.
+    [3] Macdonald, I. G. (1998). Symmetric functions and Hall polynomials. Oxford university press.
 """
 
 from math import prod
@@ -36,20 +44,33 @@
 
 
 @lru_cache
-def zonal_spherical_function(permutation: Permutation, partition: tuple[int]) -> float:
+def zonal_spherical_function(permutation: Permutation, partition: tuple[int]) -> Fraction:
     """Returns the zonal spherical function of the Gelfand pair (S_2k, H_k)
-    as seen in Macdonald's "Symmetric Functions and Hall Polynomials" chapter VII
-
-    Args:
-        permutation (Permutation): A permutation of the symmetric group S_2k
-        partition (tuple[int]): A partition of k
 
-    Returns:
-        (float): The zonal spherical function of the given permutation
-
-    Raise:
-        TypeError: If partition argument is not a tuple
-        TypeError: If permutation argument is not a permutation
+    Parameters
+    ----------
+        permutation (Permutation) : a permutation of the symmetric group S_2k
+        partition (tuple[int]) : a partition of k
+
+    Returns
+    -------
+        Fraction : the zonal spherical function of the given permutation
+
+    Raise
+    -----
+        TypeError : if partition argument is not a tuple
+        TypeError : if permutation argument is not a permutation
+
+    Examples
+    --------
+        >>> from sympy.combinatorics import Permutation
+        >>> from haarpy import zonal_spherical_function
+        >>> zonal_spherical_function(Permutation(5)(0,1,2), (2,1))
+        Fraction(1, 6)
+
+    See Also
+    --------
+        HyperoctahedralGroup, murn_naka_rule
     """
     if not isinstance(partition, tuple):
         raise TypeError
@@ -79,17 +100,39 @@ def weingarten_orthogonal(
 ) -> Expr:
     """Returns the orthogonal Weingarten function
 
-    Args:
-        permutation (Permutation, tuple[int]): A permutation of S_2k or its coset-type
-        orthogonal_dimension (int): Dimension of the orthogonal group
-
-    Returns:
-        Symbol: The Weingarten function
-
-    Raise:
-        TypeError: if unitary_dimension has the wrong type
-        TypeError: if permutation has the wrong type
-        ValueError: if the degree 2k of the symmetric group S_2k is not a factor of 2
+    Parameters
+    ----------
+        permutation (Permutation, tuple[int]) : a permutation of S_2k or its coset-type
+        orthogonal_dimension (int): dimension of the orthogonal group
+
+    Returns
+    -------
+        Symbol : the Weingarten function
+
+    Raise
+    -----
+        TypeError : if unitary_dimension has the wrong type
+        TypeError : if permutation has the wrong type
+        ValueError : if the degree 2k of the symmetric group S_2k is not a factor of 2
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from sympy.combinatorics import Permutation
+        >>> from haarpy import weingarten_orthogonal
+        >>> d = Symbol("d")
+        >>> weingarten_orthogonal(Permutation(5)(0,1,2), 6)
+        Fraction(-1, 1200)
+        >>> weingarten_orthogonal(Permutation(5)(0,1,2), d)
+        -1/(d*(d - 2)*(d - 1)*(d + 4))
+        >>> weingarten_orthogonal((2,1), d)
+        -1/(d*(d - 2)*(d - 1)*(d + 4))
+
+        Where (2,1) is the coset-type of Permutation(5)(0,1,2)
+
+    See Also
+    --------
+        zonal_spherical_function, coset_type_representative
     """
     if not isinstance(orthogonal_dimension, (Expr, int)):
         raise TypeError("orthogonal_dimension must be an instance of int or sympy.Symbol")
@@ -161,18 +204,36 @@ def weingarten_orthogonal(
 
 @lru_cache
 def haar_integral_orthogonal(sequences: tuple[tuple[int]], orthogonal_dimension: Symbol) -> Expr:
-    """Returns integral over orthogonal group polynomial sampled at random from the Haar measure
-
-    Args:
-        sequences (tuple[tuple[int]]): Indices of matrix elements
-        orthogonal_dimension (int): Dimension of the orthogonal group
-
-    Returns:
-        Expr: Integral under the Haar measure
-
-    Raise:
-        ValueError: If sequences doesn't contain 2 tuples
-        ValueError: If tuples i and j are of different length
+    """Returns the integral over orthogonal group polynomial sampled at random from the Haar measure
+
+    Parameters
+    ----------
+        sequences (tuple[tuple[int]]) : indices of matrix elements
+        orthogonal_dimension (int) : dimension of the orthogonal group
+
+    Returns
+    -------
+        Expr : integral under the Haar measure
+
+    Raise
+    -----
+        ValueError : if sequences do not contain 2 tuples
+        ValueError : if tuples i and j are of different length
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from haarpy import haar_integral_orthogonal
+        >>> d = Symbol("d")
+        >>> row_indices, column_indices = (0, 0, 1, 1, 2, 2), (0, 2, 2, 1, 1, 0)
+        >>> haar_integral_orthogonal((row_indices, column_indices), 4)
+        Fraction(1, 576)
+        >>> haar_integral_orthogonal((row_indices, column_indices), d)
+        2/(d*(d - 2)*(d - 1)*(d + 2)*(d + 4))
+
+    See Also
+    --------
+        hyperoctahedral_transversal, coset_type, weingarten_orthogonal
     """
     if len(sequences) != 2:
         raise ValueError("Wrong tuple format")
diff --git a/haarpy/partition.py b/haarpy/partition.py
index 6e80829..42bd208 100644
--- a/haarpy/partition.py
+++ b/haarpy/partition.py
@@ -13,6 +13,13 @@
 # limitations under the License.
 """
 Partition Python interface
+
+References
+----------
+    [1] Collins, B., & Nagatsu, M. (2025). Weingarten calculus for centered random permutation
+    matrices. arXiv preprint arXiv:2503.18453.
+    [2] Matsumoto, S. (2013). Weingarten calculus for matrix ensembles associated with compact
+    symmetric spaces. arXiv preprint arXiv:1301.5401.
 """
 
 from functools import lru_cache
@@ -24,14 +31,28 @@ def set_partitions(collection: tuple) -> Generator[tuple[tuple], None, None]:
     """Returns the partitionning of a given collection (set) of objects
     into non-empty subsets.
 
-    Args:
-        collection (tuple): An indexable iterable to be partitionned
-
-    Returns:
-        generator(tuple[tuple]): all partitions of the input collection
-
-    Raise:
-        ValueError: if the collection is not a tuple
+    Parameters
+    ----------
+        collection (tuple) : an indexable iterable to be partitionned
+
+    Returns
+    -------
+        generator(tuple[tuple]) : all partitions of the input collection
+
+    Raise
+    -----
+        ValueError : if the collection is not a tuple
+
+    Examples
+    --------
+        >>> from haarpy import set_partitions
+        >>> for partition in set_partitions((0,1,2)):
+        >>>     print(partition)
+        ((0, 1, 2),)
+        ((0,), (1, 2))
+        ((0, 1), (2,))
+        ((0, 2), (1,))
+        ((0,), (1,), (2,))
     """
     if not isinstance(collection, tuple):
         raise TypeError("collection must be a tuple")
@@ -52,16 +73,27 @@ def perfect_matchings(
 ) -> Generator[tuple[tuple[int]], None, None]:
     """Returns the partitions of a tuple in terms of perfect matchings.
 
-    Args:
-        seed (tuple[int]): a tuple representing the (multi-)set that will be partitioned.
-            Note that it must hold that ``len(s) >= 2``.
-
-    Returns:
-        generator: a generators that goes through all the single-double
-        partitions of the tuple
-
-    Raise:
-        TypeError: if the seed is not a tuple
+    Parameters
+    ----------
+        seed (tuple[int]) : a tuple representing the (multi-)set that will be partitioned.
+            Note that it must hold that ``len(s) >= 2``
+
+    Returns
+    -------
+        generator : all the single-double partitions of the tuple
+
+    Raise
+    -----
+        TypeError : if the seed is not a tuple
+
+    Examples
+    --------
+        >>> from haarpy import perfect_matchings
+        >>> for matching in perfect_matchings((0,1,2,3)):
+        >>>     print(matching)
+        ((0, 1), (2, 3))
+        ((0, 2), (1, 3))
+        ((0, 3), (1, 2))
     """
     if not isinstance(seed, tuple):
         raise TypeError("seed must be a tuple")
@@ -79,21 +111,31 @@ def perfect_matchings(
 
 @lru_cache
 def partial_order(partition_1: tuple[tuple[int]], partition_2: tuple[tuple[int]]) -> bool:
-    """Returns True if parition_1 <= partition_2 in terms of partial order
+    """Checks if parition_1 <= partition_2 in terms of partial order
 
     For parition_1 and partition_2, two partitions of the same set, we call
-    parition_1 <= partition_2 if and only if each block of parition_1 is
+    partition_1 <= partition_2 if and only if each block of partition_1 is
     contained in some block of partition_2
 
-    Ex.
-    ((0,1), (2,3), (4,)) < ((0,1), (2,3,4))
-
-    Args:
-        partition_1 tuple[tuple[int]]: The partition of lower order
-        partition_2 tuple[tuple[int]]: The partition of higher order
-
-    Returns:
-        bool: True if parition_1 <= partition_2
+    Parameters
+    ----------
+        partition_1 (tuple[tuple[int]]) : the partition of lower order
+        partition_2 (tuple[tuple[int]]) : the partition of higher order
+
+    Returns
+    -------
+        bool : True if partition_1 <= partition_2
+
+    Examples
+    --------
+        >>> from haarpy import partial_order
+        >>> partition_1 = ((0, 1), (2, 3), (4,))
+        >>> partition_2 = ((0, 1), (2, 3, 4))
+        >>> partition_3 = ((0, 4), (1, 2, 3))
+        >>> partial_order(partition_1, partition_2)
+        True
+        >>> partial_order(partition_1, partition_3)
+        False
     """
     for part in partition_1:
         if all(not set(part).issubset(bigger_part) for bigger_part in partition_2):
@@ -109,15 +151,25 @@ def meet_operation(partition_1: tuple[tuple[int]], partition_2: tuple[tuple[int]
     For parition_1 and partition_2, two partitions of the same set,
     the meet operation yields the greatest lower bound of both partitions
 
-    Ex.
-    ((0,1), (2,3), (4,)) ∧ ((0,1,2), (3,4)) = ((0,1), (2,), (3,), (4,))
-
-    Args:
-        partition_1 (tuple[tuple[int]]): partition of a set
-        partition_2 (tuple[tuple[int]]): partition of a set
-
-    Return:
-        tuple[tuple]: Greatest lower bound
+    The meet operation symbol is ∧, for instance
+    ((0, 1), (2, 3), (4,)) ∧ ((0, 1, 2), (3, 4)) = ((0, 1), (2,), (3,), (4,))
+
+    Parameters
+    ----------
+        partition_1 (tuple[tuple[int]]) : partition of a set
+        partition_2 (tuple[tuple[int]]) : partition of a set
+
+    Returns
+    -------
+        tuple[tuple] : greatest lower bound
+
+    Examples
+    --------
+        >>> from haarpy import meet_operation
+        >>> partition_1 = ((0, 1), (2, 3), (4,))
+        >>> partition_2 = ((0, 1, 2), (3, 4))
+        >>> meet_operation(partition_1, partition_2)
+        ((0, 1), (2,), (3,), (4,))
     """
     partition_1 = tuple(set(part) for part in partition_1)
     partition_2 = tuple(set(part) for part in partition_2)
@@ -138,15 +190,25 @@ def join_operation(partition_1: tuple[tuple[int]], partition_2: tuple[tuple[int]
     For parition_1 and partition_2, two partitions of the same set,
     the join operation yields the least upper bound of both partitions
 
-    Ex.
-    ((0,1), (2,), (3,4)) ∨ ((0,2), (1,), (3,), (4,)) = ((0,1,2), (3,4))
-
-    Args:
-        partition_1 (tuple[tuple[int]]): partition of a set
-        partition_2 (tuple[tuple[int]]): partition of a set
-
-    Return:
-        tuple[tuple[int]]: Least upper bound
+    The join operation symbol is ∨, for instance
+    ((0, 1), (2,), (3, 4)) ∨ ((0, 2), (1,), (3,), (4,)) = ((0, 1, 2), (3, 4))
+
+    Parameters
+    ----------
+        partition_1 (tuple[tuple[int]]) : partition of a set
+        partition_2 (tuple[tuple[int]]) : partition of a set
+
+    Returns
+    ------
+        tuple[tuple[int]] : least upper bound
+
+    Examples
+    --------
+        >>> from haarpy import join_operation
+        >>> partition_1 = ((0, 1), (2,), (3, 4))
+        >>> partition_2 = ((0, 2), (1,), (3,), (4,))
+        >>> join_operation(partition_1, partition_2)
+        ((0, 1, 2), (3, 4))
     """
     parent = [
         {index for value in block1 for index, block2 in enumerate(partition_2) if value in block2}
@@ -171,18 +233,23 @@ def join_operation(partition_1: tuple[tuple[int]], partition_2: tuple[tuple[int]
 
 @lru_cache
 def is_crossing_partition(partition: tuple[tuple[int]]) -> bool:
-    """
-    Checks if the partition is a crossing partition
-
-    Ex.
-    Crossing partition : ((0,2,4), (1,3))
-    Non crossing partition : ((0,3,4), (1,2))
-
-    Args:
-        partition (tuple[tuple[int]])): partition of a set
-
-    Returns:
-        bool: True if the partition is crossing, False otherwise
+    """Checks if a given partition is crossing
+
+    Parameters
+    ----------
+        partition (tuple[tuple[int]])) : partition of a set
+
+    Returns
+    -------
+        bool : True if the partition is crossing, False otherwise
+
+    Examples
+    --------
+        >>> from haarpy import is_crossing_partition
+        >>> is_crossing_partition(((0,2,4), (1,3)))
+        True
+        >>> is_crossing_partition(((0,3,4), (1,2)))
+        False
     """
     filtered_partition = tuple(
         block for block in partition if len(block) != 1 and block[0] + 1 != block[-1]
diff --git a/haarpy/permutation.py b/haarpy/permutation.py
index 05b711a..1bed0b0 100644
--- a/haarpy/permutation.py
+++ b/haarpy/permutation.py
@@ -13,28 +13,42 @@
 # limitations under the License.
 """
 Permutation matrices Python interface
+
+References
+----------
+    [1] Collins, B., & Nagatsu, M. (2025). Weingarten calculus for centered random permutation
+    matrices. arXiv preprint arXiv:2503.18453.
 """
 
 from math import factorial, prod
 from functools import lru_cache
 from itertools import product
 from collections.abc import Sequence
+from fractions import Fraction
 from sympy import Symbol, simplify, binomial, factor, fraction
 from haarpy import set_partitions, meet_operation, join_operation, partial_order
 
 
 @lru_cache
 def mobius_function(partition_1: tuple[tuple[int]], partition_2: tuple[tuple[int]]) -> int:
-    """Return the Möbius function
-    as seen in `Collins and Nagatsu. Weingarten Calculus for Centered Random
-    Permutation Matrices <https://arxiv.org/abs/2503.18453>`_
-
-    Args:
-        partition_1 (tuple[tuple[int]): The intersected partition
-        partition_2 (tuple[tuple[int]]): The partition summed over
-
-    Returns:
-        int: The value of the Möbius function
+    """Returns the Möbius function of two given partitions
+
+    Parameters
+    ----------
+        partition_1 (tuple[tuple[int]) : the intersected partition
+        partition_2 (tuple[tuple[int]]) : the partition summed over
+
+    Returns
+    -------
+        int : the value of the Möbius function
+
+    Examples
+    --------
+        >>> from haarpy import mobius_function
+        >>> partition_1 = ((0, 3), (1, 2), (4,), (5,))
+        >>> partition_2 = ((0, 1, 2), (3, 4, 5))
+        >>> mobius_function(partition_1, partition_2)
+        -2
     """
     partition_set_1 = tuple(set(block) for block in partition_1)
     partition_set_2 = tuple(set(block) for block in partition_2)
@@ -56,16 +70,32 @@ def weingarten_permutation(
     dimension: Symbol,
 ) -> Symbol:
     """Returns the Weingarten function for random permutation matrices
-    as seen in `Collins and Nagatsu. Weingarten Calculus for Centered Random
-    Permutation Matrices <https://arxiv.org/abs/2503.18453>`_
-
-    Args:
-        first_partition (tuple(tuple(int))): a set partition of integer k
-        second_partition (tuple(tuple(int))): a set partition of integer k
-        dimension (Symbol): Dimension of the random permutation matrices
 
-    Returns:
-        Symbol : The Weingarten function
+    Parameters
+    ----------
+        first_partition (tuple(tuple(int))) : a set partition of integer k
+        second_partition (tuple(tuple(int))) : a set partition of integer k
+        dimension (Symbol) : the dimension of the random permutation matrices
+
+    Returns
+    -------
+        Symbol : the Weingarten function
+
+    Example
+    -------
+        >>> from sympy import Symbol
+        >>> from haarpy import weingarten_permutation
+        >>> d = Symbol('d')
+        >>> partition_1 = ((0, 1, 2), (3,))
+        >>> partition_2 = ((0,), (1, 2, 3))
+        >>> weingarten_permutation(partition_1, partition_2, 4)
+        Fraction(5, 24)
+        >>> weingarten_permutation(partition_1, partition_2, d)
+        (d + 1)/(d*(d - 3)*(d - 2)*(d - 1))
+
+    See Also
+    --------
+        meet_operation, mobius_function
     """
     disjoint_partition_tuple = tuple(
         (partition for partition in set_partitions(block))
@@ -77,14 +107,26 @@ def weingarten_permutation(
         for partition_tuple in product(*disjoint_partition_tuple)
     )
 
-    weingarten = sum(
-        mobius_function(partition, first_partition)
-        * mobius_function(partition, second_partition)
-        / prod(dimension - i for i, _ in enumerate(partition))
-        for partition in inferieur_partition_tuple
-    )
+    if isinstance(dimension, int):
+        weingarten = sum(
+            Fraction(
+                mobius_function(partition, first_partition)
+                * mobius_function(partition, second_partition),
+                prod(dimension - i for i, _ in enumerate(partition)),
+            )
+            for partition in inferieur_partition_tuple
+        )
+    else:
+        weingarten = sum(
+            mobius_function(partition, first_partition)
+            * mobius_function(partition, second_partition)
+            / prod(dimension - i for i, _ in enumerate(partition))
+            for partition in inferieur_partition_tuple
+        )
+        numerator, denominator = fraction(simplify(weingarten))
+        weingarten = factor(numerator) / factor(denominator)
 
-    return weingarten if isinstance(dimension, int) else simplify(weingarten)
+    return weingarten
 
 
 @lru_cache
@@ -94,16 +136,32 @@ def weingarten_centered_permutation(
     dimension: Symbol,
 ) -> Symbol:
     """Returns the Weingarten function for centered random permutation matrices
-    as seen in `Collins and Nagatsu. Weingarten Calculus for Centered Random
-    Permutation Matrices <https://arxiv.org/abs/2503.18453>`_
-
-    Args:
-        first_partition (tuple(tuple(int))): a set partition of integer k
-        second_partition (tuple(tuple(int))): a set partition of integer k
-        dimension (Symbol): Dimension of the centered random permutation matrices
 
-    Returns:
-        Symbol : The Weingarten function
+    Parameters
+    ----------
+        first_partition (tuple(tuple(int))) : a set partition of integer k
+        second_partition (tuple(tuple(int))) : a set partition of integer k
+        dimension (Symbol) : the dimension of the centered random permutation matrices
+
+    Returns
+    -------
+        Symbol : the Weingarten function
+
+    Example
+    -------
+        >>> from sympy import Symbol
+        >>> from haarpy import weingarten_centered_permutation
+        >>> d = Symbol('d')
+        >>> partition_1 = ((0,), (1,), (2,))
+        >>> partition_2 = ((0, 2), (1,))
+        >>> weingarten_centered_permutation(partition_1, partition_2, 3)
+        Fraction(-1, 9)
+        >>> weingarten_centered_permutation(partition_1, partition_2, d)
+        -2/(d**2*(d - 2)*(d - 1))
+
+    See Also
+    --------
+        join_operation, meet_operation, mobius_function
     """
     singleton_set_size = sum(
         1 for block in join_operation(first_partition, second_partition) if len(block) == 1
@@ -119,20 +177,34 @@ def weingarten_centered_permutation(
         for partition_tuple in product(*disjoint_partition_tuple)
     )
 
-    weingarten = sum(
-        (-1) ** i
-        * binomial(singleton_set_size, i)
-        * sum(
-            mobius_function(partition, first_partition)
-            * mobius_function(partition, second_partition)
-            / dimension**i
-            / prod(dimension - j for j in range(len(partition) - i))
-            for partition in inferieur_partition_tuple
+    if isinstance(dimension, int):
+        weingarten = sum(
+            (-1) ** i
+            * Fraction(binomial(singleton_set_size, i))
+            * sum(
+                Fraction(
+                    mobius_function(partition, first_partition)
+                    * mobius_function(partition, second_partition),
+                    dimension**i * prod(dimension - j for j in range(len(partition) - i)),
+                )
+                for partition in inferieur_partition_tuple
+            )
+            for i in range(singleton_set_size + 1)
+        )
+    else:
+        weingarten = sum(
+            (-1) ** i
+            * binomial(singleton_set_size, i)
+            * sum(
+                mobius_function(partition, first_partition)
+                * mobius_function(partition, second_partition)
+                / dimension**i
+                / prod(dimension - j for j in range(len(partition) - i))
+                for partition in inferieur_partition_tuple
+            )
+            for i in range(singleton_set_size + 1)
         )
-        for i in range(singleton_set_size + 1)
-    )
 
-    if isinstance(dimension, Symbol):
         num, denum = fraction(simplify(weingarten))
         weingarten = factor(num) / factor(denum)
 
@@ -146,19 +218,32 @@ def haar_integral_permutation(
     dimension: Symbol,
 ) -> Symbol:
     """Returns the integral over Haar random permutation matrices
-    as seen in `Collins and Nagatsu. Weingarten Calculus for Centered Random
-    Permutation Matrices <https://arxiv.org/abs/2503.18453>`_
 
-    Args:
+    Parameters
+    ----------
         row_indices (tuple(int)) : sequence of row indices
         column_indices (tuple(int)) : sequence of column indices
 
-    Returns:
-        Symbol : Integral under the Haar measure
-
-    Raise:
-        TypeError : If row_indices and column_indices are not Sequence
-        ValueError : If row_indices and column_indices are of different length
+    Returns
+    -------
+        Symbol : integral under the Haar measure
+
+    Raise
+    -----
+        TypeError : if row_indices and column_indices are not Sequence
+        ValueError : if row_indices and column_indices are of different length
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from haarpy import haar_integral_permutation
+        >>> d = Symbol('d')
+        >>> row_indices = (0, 1, 2, 2)
+        >>> column_indices = (1, 0, 2, 2)
+        >>> haar_integral_permutation(row_indices, column_indices, 3)
+        Fraction(1, 6)
+        >>> haar_integral_permutation(row_indices, column_indices, d)
+        1/(d*(d - 2)*(d - 1))
     """
     if not (isinstance(row_indices, Sequence) and isinstance(column_indices, Sequence)):
         raise TypeError
@@ -176,9 +261,13 @@ def sequence_to_partition(sequence: tuple) -> tuple[tuple[int]]:
     column_partition = sequence_to_partition(column_indices)
 
     return (
-        1 / prod(dimension - i for i, _ in enumerate(row_partition))
-        if row_partition == column_partition
-        else 0
+        Fraction(1, prod(dimension - i for i, _ in enumerate(row_partition)))
+        if row_partition == column_partition and isinstance(dimension, int)
+        else (
+            1 / prod(dimension - i for i, _ in enumerate(row_partition))
+            if row_partition == column_partition
+            else 0
+        )
     )
 
 
@@ -189,19 +278,33 @@ def haar_integral_centered_permutation(
     dimension: Symbol,
 ) -> Symbol:
     """Returns the integral over Haar random centered permutation matrices
-    as seen in `Collins and Nagatsu. Weingarten Calculus for Centered Random
-    Permutation Matrices <https://arxiv.org/abs/2503.18453>`_
 
     Args:
         row_indices (tuple(int)) : sequence of row indices
         column_indices (tuple(int)) : sequence of column indices
 
     Returns:
-        Symbol : Integral under the Haar measure
+        Symbol : integral under the Haar measure
 
     Raise:
-        TypeError : If row_indices and column_indices are not Sequence
-        ValueError : If row_indices and column_indices are of different length
+        TypeError : if row_indices and column_indices are not Sequence
+        ValueError : if row_indices and column_indices are of different length
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from haarpy import haar_integral_centered_permutation
+        >>> d = Symbol('d')
+        >>> row_indices = (0, 1, 2)
+        >>> column_indices = (0, 1, 1)
+        >>> haar_integral_centered_permutation(row_indices, column_indices, 3)
+        Fraction(-1, 27)
+        >>> haar_integral_centered_permutation(row_indices, column_indices, d)
+        -2/(d**3*(d - 1))
+
+    See Also
+    --------
+        partial_order, weingarten_centered_permutation
     """
     if not (isinstance(row_indices, Sequence) and isinstance(column_indices, Sequence)):
         raise TypeError
diff --git a/haarpy/symmetric.py b/haarpy/symmetric.py
index 28cfea5..a1c579f 100644
--- a/haarpy/symmetric.py
+++ b/haarpy/symmetric.py
@@ -13,6 +13,15 @@
 # limitations under the License.
 """
 Symmetric group Python interface
+
+References
+----------
+    [1] Collins, B. (2003). Moments and cumulants of polynomial random variables on unitarygroups,
+    the Itzykson-Zuber integral, and free probability. International Mathematics Research Notices,
+    2003(17), 953-982.
+    [2] Matsumoto, S. (2013). Weingarten calculus for matrix ensembles associated with compact
+    symmetric spaces. arXiv preprint arXiv:1301.5401.
+    [3] Macdonald, I. G. (1998). Symmetric functions and Hall polynomials. Oxford university press.
 """
 
 from math import factorial, prod
@@ -32,18 +41,28 @@
 def get_conjugacy_class(perm: Permutation, degree: int) -> tuple:
     """Returns the conjugacy class of an element of the symmetric group Sp
 
-    Args:
-        perm (Permutation): Permutation cycle from the symmetric group
-        degree (integer): Order of the symmetric group
+    Parameters
+    ----------
+        perm (Permutation) : permutation cycle from the symmetric group
+        degree (integer) : order of the symmetric group
 
-    Returns:
-        tuple[int]: the conjugacy class in partition form
+    Returns
+    -------
+        tuple[int] : the conjugacy class in partition form
 
-    Raise:
+    Raise
+    -----
         TypeError : order must be of type int
-        ValueError : If order is not an integer greaten than 1
+        ValueError : if order is not an integer greaten than 1
         TypeError : cycle must be of type sympy.combinatorics.permutations.Permutation
-        ValueError : Incompatible degree and permutation cycle
+        ValueError : incompatible degree and permutation cycle
+
+    Examples
+    --------
+        >>> from sympy.combinatorics import Permutation
+        >>> from haarpy import get_conjugacy_class
+        >>> get_conjugacy_class(Permutation(5)(0,1,3), 6)
+        (3, 1, 1, 1)
     """
     if not isinstance(degree, int):
         raise TypeError("degree must be of type int")
@@ -74,13 +93,15 @@ def derivative_tableaux(
     """Takes a single tableau and adds the selected number to its contents
     in a way that keeps it semi-standard. All possible tableaux are yielded
 
-    Args:
-        tableau (tuple[tuple[int]]): An incomplet Young tableau
-        increment (int): Selected number to be added
+    Parameters
+    ----------
+        tableau (tuple[tuple[int]]) : an incomplete Young tableau
+        increment (int) : selected number to be added
         partition (tuple[int]) : partition characterizing an irrep of Sp
 
-    Yields:
-        tuple[tuple[int]]: Modified tableaux
+    Yields
+    ------
+        tuple[tuple[int]] : modified tableaux
     """
     # empty tableau
     if not tableau[0]:
@@ -117,12 +138,22 @@ def semi_standard_young_tableaux(
 ) -> set[tuple[tuple[int]]]:
     """all eligible semi-standard young tableaux based of the partition
 
-    Args:
+    Parameters
+    ----------
         partition (tuple[int]) : partition characterizing an irrep of Sp
-        conjugacy_class (tuple[int]) : A conjugacy class, in partition form, of Sp
-
-    Returns:
-        set[tuple[tuple[int]]]: all eligible semi-standard young tableaux
+        conjugacy_class (tuple[int]) : a conjugacy class, in partition form, of Sp
+
+    Returns
+    -------
+        set[tuple[tuple[int]]] : all eligible semi-standard young tableaux
+
+    Examples
+    --------
+        >>> from haarpy import semi_standard_young_tableaux
+        >>> semi_standard_young_tableaux((3, 1), (2, 1, 1))
+        {((0, 1, 2), (0,)), ((0, 0, 2), (1,)), ((0, 0, 1), (2,))}
+        >>> semi_standard_young_tableaux((3,2),(4,1))
+        {((0, 0, 1), (0, 0)), ((0, 0, 0), (0, 1))}
     """
     tableaux = (tuple(() for _ in partition),)
     cell_values = (i for i, m in enumerate(conjugacy_class) for _ in range(m))
@@ -140,12 +171,22 @@ def semi_standard_young_tableaux(
 def proper_border_strip(tableau: tuple[tuple[int]], conjugacy_class: tuple[int]) -> bool:
     """Returns True if input Young tableau is a valid border-strip tableau
 
-    Args:
-        tableau (tuple[tuple[int]]) : A semi-standard Young tableau
-        conjugacy_class (tuple[int]) : A conjugacy class, in partition form, of Sp
+    Parameters
+    ----------
+        tableau (tuple[tuple[int]]) : a semi-standard Young tableau
+        conjugacy_class (tuple[int]) : a conjugacy class, in partition form, of Sp
 
-    Returns:
+    Returns
+    -------
         bool : True if the tableau is a border-strip
+
+    Examples
+    --------
+        >>> from haarpy import proper_border_strip
+        >>> ap.proper_border_strip(((0, 0, 0), (0, 1)), (4,1))
+        True
+        >>> ap.proper_border_strip(((0, 0, 1), (0, 0)), (4,1))
+        False
     """
     if len(tableau) == 1:
         return True
@@ -179,12 +220,26 @@ def murn_naka_rule(partition: tuple[int], conjugacy_class: tuple[int]) -> int:
     """Implementation of the Murnaghan-Nakayama rule for the characters irreducible
     representations of the symmetric group Sp
 
-    Args:
+    Parameters
+    ----------
         partition (tuple[int]) : partition characterizing an irrep of Sp
-        conjugacy_class (tuple[int]) : A conjugacy class, in partition form, of Sp
-
-    Returns:
-        int : Character of the elements in class mu of the irrep of the symmetric group
+        conjugacy_class (tuple[int]) : a conjugacy class, in partition form, of Sp
+
+    Returns
+    -------
+        int : character of the elements in class mu of the irrep of the symmetric group
+
+    Examples
+    --------
+        >>> from haarpy import murn_naka_rule
+        >>> murn_naka_rule((4, 1, 1), (3, 2, 1))
+        -1
+        >>> murn_naka_rule((3, 1), (1, 1, 1, 1))
+        3
+
+    See Also
+    --------
+        semi_standard_young_tableaux, proper_border_strip
     """
     if sum(partition) != sum(conjugacy_class):
         return 0
@@ -206,11 +261,19 @@ def murn_naka_rule(partition: tuple[int], conjugacy_class: tuple[int]) -> int:
 def irrep_dimension(partition: tuple[int]) -> int:
     """Returns the dimension of the irrep of the symmetric group Sp labelled by the input partition
 
-    Args:
-        partition (tuple[int]) : A partition labelling an irrep of Sp
+    Parameters
+    ----------
+        partition (tuple[int]) : a partition labelling an irrep of Sp
+
+    Returns
+    -------
+        int : the dimension of the irrep
 
-    Returns:
-        int : The dimension of the irrep
+    Examples
+    --------
+        >>> from haarpy import irrep_dimension
+        >>> irrep_dimension((2, 1, 1))
+        3
     """
     numerator = prod(
         part_i - part_j + j + 1
@@ -227,15 +290,34 @@ def irrep_dimension(partition: tuple[int]) -> int:
 def sorting_permutation(*sequence: tuple[int]) -> Permutation:
     """Returns the sorting permutation of a given sequence
 
-    Args:
-        sequence (tuple[int]): a sequence of unorderd elements
-
-    Returns:
-        Permutation: the sorting permutation
-
-    Raise:
-        ValueError: for incompatible sequence inputs
-        TypeError: if more than two sequences are passed as arguments
+    If two sequences, sequence_1 and sequence_2, are given as parameters,
+    the function returns the permutation such that
+    permutation(sequence_1) = sequence_2
+
+    Note that the function return the first permutation found if the sequences
+    contain multiplicity
+
+    Parameters
+    ----------
+        *sequence (tuple[int]) : one or two sequences of unorderd elements
+
+    Returns
+    -------
+        Permutation : the sorting permutation
+
+    Raise
+    -----
+        ValueError : for incompatible sequence inputs
+        TypeError : if more than two sequences are passed as arguments
+
+    Examples
+    --------
+        >>> from haarpy import sorting_permutation
+        >>> sequence_1, sequence_2 = (2, 1, 2, 1, 3), (3, 2, 2, 1, 1)
+        >>> sorting_permutation(sequence_1)
+        Permutation(4)(0, 1, 3, 2)
+        >>> sorting_permutation(sequence_1, sequence_2)
+        Permutation(0, 4, 3, 1)
     """
     if len(sequence) == 1:
         return Permutation(sorted(range(len(sequence[0])), key=lambda k: sequence[0][k]))
@@ -247,19 +329,33 @@ def sorting_permutation(*sequence: tuple[int]) -> Permutation:
     raise TypeError
 
 
+# pylint: disable=invalid-name
 def YoungSubgroup(partition: tuple[int]) -> PermutationGroup:
     """Returns the Young subgroup of a given input partition
-    See `<https://en.wikipedia.org/wiki/YoungSubgroup>`_
-
-    Args:
-        partition (tuple[int]): A partition
 
-    Returns:
-        PermutationGroup: the associated Young subgroup
-
-    Raise:
-        TypeError: if partition is not a tuple or a list
-        TypeError: if partition is not made of positive integers
+    Parameters
+    ----------
+        partition (tuple[int]) : a partition
+
+    Returns
+    -------
+        PermutationGroup : the associated Young subgroup
+
+    Raise
+    -----
+        TypeError : if partition is not a tuple or a list
+        TypeError : if partition is not made of positive integers
+
+    Examples
+    --------
+        >>> from haarpy import YoungSubgroup
+        >>> young = YoungSubgroup((2, 2))
+        >>> list(young.generate())
+        [Permutation(3), Permutation(3)(0, 1), Permutation(2, 3), Permutation(0, 1)(2, 3)]
+
+    See Also
+    --------
+        sympy.combinatorics.DirectProduct, sympy.combinatorics.SymmetricGroup
     """
     if not isinstance(partition, (tuple, list)):
         raise TypeError
@@ -271,14 +367,34 @@ def YoungSubgroup(partition: tuple[int]) -> PermutationGroup:
 def stabilizer_coset(*sequence: tuple) -> Generator[Permutation, None, None]:
     """Returns all permutations that, when acting on sequence[0], return sequence[1]
 
-    Args:
-        *sequence (tuple): the sequences acted upon
-
-    Returns:
-        Generator[Permutation]: permutations that, when acting on sequence[0], return sequence[1]
-
-    Raise:
-        TypeError: if the sequence argument contains more than two sequences
+    For a single input, the function returns the stabilizer group with respect to the sequence.
+    For two inputs, it returns the stabilizer of the first sequence with respect to the second,
+    that is, all permutations such that permutation(sequence[0]) == sequence[1].
+
+    Parameters
+    ----------
+        *sequence (tuple) : the sequences acted upon
+
+    Returns
+    -------
+        Generator[Permutation] : permutations that, when acting on sequence[0], return sequence[1]
+
+    Raise
+    -----
+        TypeError : if the sequence argument contains more than two sequences
+
+    Examples
+    --------
+        >>> from haarpy import stabilizer_coset
+        >>> sequence_1, sequence_2 = (2, 2, 1, 3), (3, 2, 2, 1)
+        >>> list(stabilizer_coset(sequence_1))
+        [Permutation(3), Permutation(3)(0, 1)]
+        >>> list(stabilizer_coset(sequence_1, sequence_2))
+        [Permutation(0, 3, 2, 1), Permutation(0, 3, 2)]
+
+    See Also
+    --------
+        sorting_permutation, YoungSubgroup
     """
     if len(sequence) == 1:
         sequence = tuple(sequence[0] for _ in range(2))
@@ -296,18 +412,33 @@ def stabilizer_coset(*sequence: tuple) -> Generator[Permutation, None, None]:
     )
 
 
+# pylint: disable=invalid-name
 @lru_cache
 def HyperoctahedralGroup(degree: int) -> PermutationGroup:
     """Return the hyperoctahedral group
 
-    Args:
-        degree (int): The degree k of the hyperoctahedral group H_k
+    Parameters
+    ----------
+        degree (int) : the degree k of the hyperoctahedral group H_k
+
+    Returns
+    -------
+        PermutationGroup : the hyperoctahedral group
+
+    Raise
+    -----
+        TypeError : if degree is not of type int
 
-    Returns:
-        (PermutationGroup): The hyperoctahedral group
+    Examples
+    --------
+        >>> from haarpy import HyperoctahedralGroup
+        >>> hyperoctahedral = ap.HyperoctahedralGroup(3)
+        >>> hyperoctahedral.order()
+        48
 
-    Raise:
-        TypeError: If degree is not of type int
+    See Also
+    --------
+        sympy.combinatorics.PermutationGroup
     """
     if not isinstance(degree, int):
         raise TypeError
@@ -322,14 +453,26 @@ def HyperoctahedralGroup(degree: int) -> PermutationGroup:
 
 def hyperoctahedral_transversal(degree: int) -> Generator[Permutation, None, None]:
     """Returns a generator with the permutations of M_2k, the complete set of coset
-    representatives of S_2k/H_k as seen in Macdonald's "Symmetric Functions and Hall
-    Polynomials" chapter VII
-
-    Args:
-        degree (int): Degree 2k of the set M_2k
-
-    Returns:
-        (Generator[Permutation]): The permutations of M_2k
+    representatives of S_2k/H_k
+
+    Parameters
+    ----------
+        degree (int) : degree 2k of the set M_2k
+
+    Returns
+    -------
+        Generator[Permutation] : the permutations of M_2k
+
+    Examples
+    --------
+        >>> from haarpy import hyperoctahedral_transversal
+        >>> transversal = ap.hyperoctahedral_transversal(4)
+        >>> list(transversal)
+        [Permutation(3), Permutation(3)(1, 2), Permutation(1, 3, 2)]
+
+    See Also
+    --------
+        perfect_matchings
     """
     if degree % 2:
         raise ValueError("degree should be a factor of 2")
@@ -343,19 +486,31 @@ def hyperoctahedral_transversal(degree: int) -> Generator[Permutation, None, Non
 
 @lru_cache
 def coset_type(permutation: Permutation) -> tuple[int]:
-    """Returns the coset-type of a given permutation of S_2k as seen
-    `Matsumoto. Weingarten calculus for matrix ensembles associated with
-    compact symmetric spaces <https://arxiv.org/pdf/1301.5401>`_
-
-    Args:
-        permutation (Permutation): A permutation of the symmetric group S_2k
-
-    Returns:
-        tuple[int]: The associated coset-type as a partition of k
-
-    Raise:
-        TypeError: If partition is not a Permutation
-        ValueError: If the symmetric group is of odd degree
+    """Returns the coset-type of a given permutation of S_2k
+
+    Parameters
+    ----------
+        permutation (Permutation) : a permutation of the symmetric group S_2k
+
+    Returns
+    -------
+        tuple[int] : the associated coset-type as a partition of k
+
+    Raise
+    -----
+        TypeError : if partition is not a Permutation
+        ValueError : if the symmetric group is of odd degree
+
+    Examples
+    --------
+        >>> from sympy.combinatorics import Permutation
+        >>> from haarpy import coset_type
+        >>> coset_type(Permutation(0, 5, 4, 2, 1, 3))
+        (2, 1)
+
+    See Also
+    --------
+        join_operation
     """
     if not isinstance(permutation, Permutation):
         raise TypeError
@@ -380,18 +535,25 @@ def coset_type(permutation: Permutation) -> tuple[int]:
 @lru_cache
 def coset_type_representative(partition: tuple[int]) -> Permutation:
     """Returns a representative permutation of S_2k for a given
-    input coset-type (partition of k) as seen
-    `Matsumoto. Weingarten calculus for matrix ensembles associated with
-    compact symmetric spaces <https://arxiv.org/pdf/1301.5401>`_
+    input coset-type (partition of k)
+
+    Parameters
+    ----------
+        partition (tuple[int]) : the coset-type (partition of k)
 
-    Args:
-        partition (tuple[int]): The coset-type (partition of k)
+    Returns
+    -------
+        Permutation : the associated permutation of S_2k
 
-    Returns:
-        (Permutation): The associated permutation of S_2k
+    Raise
+    -----
+        TypeError : if partition is not a tuple
 
-    Raise:
-        TypeError: If partition is not a tuple
+    Examples
+    --------
+        >>> from haarpy import coset_type_representative
+        >>> coset_type_representative((2, 1))
+        Permutation(5)(1, 3, 2)
     """
     if not isinstance(partition, tuple):
         raise TypeError
diff --git a/haarpy/symplectic.py b/haarpy/symplectic.py
index b7f1e64..a25ea34 100644
--- a/haarpy/symplectic.py
+++ b/haarpy/symplectic.py
@@ -13,6 +13,14 @@
 # limitations under the License.
 """
 Symplectic group Python interface
+
+References
+----------
+    [1] Collins, B., & Śniady, P. (2006). Integration with respect to the Haar measure on unitary,
+    orthogonal and symplectic group. Communications in Mathematical Physics, 264(3), 773-795.
+    [2] Matsumoto, S. (2013). Weingarten calculus for matrix ensembles associated with compact
+    symmetric spaces. arXiv preprint arXiv:1301.5401.
+    [3] Macdonald, I. G. (1998). Symmetric functions and Hall polynomials. Oxford university press.
 """
 
 from math import prod
@@ -33,22 +41,35 @@
 
 
 @lru_cache
-def twisted_spherical_function(permutation: Permutation, partition: tuple[int]) -> float:
+def twisted_spherical_function(permutation: Permutation, partition: tuple[int]) -> Fraction:
     """Returns the twisted spherical function of the Gelfand pair (S_2k, H_k)
-    as seen in Macdonald's "Symmetric Functions and Hall Polynomials" chapter VII
-
-    Args:
-        permutation (Permutation): A permutation of the symmetric group S_2k
-        partition (tuple[int]): A partition of k
-
-    Returns:
-        (float): The twisted spherical function of the given permutation
 
-    Raise:
-        TypeError: If partition is not a tuple
-        TypeError: If permutation argument is not a permutation.
-        ValueError: If the degree of the permutation is not a factor of 2
-        ValueError: If the degree of the partition and the permutation are incompatible
+    Parameters
+    ----------
+        permutation (Permutation) : a permutation of the symmetric group S_2k
+        partition (tuple[int]) : a partition of k
+
+    Returns
+    -------
+        Fraction : the twisted spherical function of the given permutation
+
+    Raise
+    -----
+        TypeError : if partition is not a tuple
+        TypeError : if permutation argument is not a permutation.
+        ValueError : if the degree of the permutation is not a factor of 2
+        ValueError : if the degree of the partition and the permutation are incompatible
+
+    Examples
+    --------
+        >>> from sympy.combinatorics import Permutation
+        >>> from haarpy import twisted_spherical_function
+        >>> twisted_spherical_function(Permutation(5, 0, 1, 2), (2, 1))
+        Fraction(1, 4)
+
+    See Also
+    --------
+        murn_naka_rule
     """
     if not isinstance(partition, tuple):
         raise TypeError
@@ -77,15 +98,33 @@ def twisted_spherical_function(permutation: Permutation, partition: tuple[int])
 def weingarten_symplectic(permutation: Permutation, half_dimension: Symbol) -> Expr:
     """Returns the symplectic Weingarten function
 
-    Args:
-        permutation (Permutation): A permutation of the symmetric group S_2k
-        half_dimension (Symbol): Half the dimension of the symplectic group
-
-    Returns:
-        Expr : The Weingarten function
-
-    Raise:
-        ValueError : If the degree 2k of the symmetric group S_2k is not a factor of 2
+    Parameters
+    ----------
+        permutation (Permutation) : a permutation of the symmetric group S_2k
+        half_dimension (Symbol) : half the dimension of the symplectic group
+
+    Returns
+    -------
+        Expr : the Weingarten function
+
+    Raise
+    -----
+        ValueError : if the degree 2k of the symmetric group S_2k is not a factor of 2
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from sympy.combinatorics import Permutation
+        >>> from haarpy import weingarten_symplectic
+        >>> d = Symbol("d")
+        >>> weingarten_symplectic(Permutation(0, 1, 2, 3, 4, 5), 7)
+        Fraction(-1, 201600)
+        >>> weingarten_symplectic(Permutation(0, 1, 2, 3, 4, 5), d)
+        -1/(8*d*(d - 2)*(d - 1)*(d + 1)*(2*d + 1))
+
+    See Also
+    --------
+        twisted_spherical_function, sympy.utilities.iterables.partitions
     """
     degree = permutation.size
     if degree % 2:
@@ -111,7 +150,7 @@ def weingarten_symplectic(permutation: Permutation, half_dimension: Symbol) -> E
         for partition in partition_tuple
     )
 
-    if isinstance(half_dimension, int):
+    if isinstance(half_dimension, (int, Fraction)):
         weingarten = sum(
             Fraction(
                 irrep_dim * zonal_spherical,
@@ -151,21 +190,43 @@ def haar_integral_symplectic(
 ) -> Expr:
     """Returns integral over symplectic group polynomial sampled at random from the Haar measure
 
-    Args:
-        sequences (tuple[tuple[Expr]]): Indices of matrix elements
-        half_dimension (Symbol): Half the dimension of the symplectic group
-
-    Returns:
-        Expr: Integral under the Haar measure
-
-    Raise:
-        ValueError: If sequences don't contain 2 tuples
-        ValueError: If tuples i and j are of different length
-        TypeError: If the half_dimension is not int nor Symbol
-        TypeError: If dimension is int and sequence is not
-        ValueError: If all sequence indices are not between 0 and 2*dimension - 1
-        TypeError: If sequence containt something else than Expr
-        TypeError: If symbolic sequences have the wrong format
+    Parameters
+    ----------
+        sequences (tuple[tuple[Expr]]) : indices of matrix elements
+        half_dimension (Symbol) : half the dimension of the symplectic group
+
+    Returns
+    -------
+        Expr : integral under the Haar measure
+
+    Raise
+    -----
+        ValueError : if sequences do not contain 2 tuples
+        ValueError : if tuples i and j are of different length
+        TypeError : if the half_dimension is not int nor Symbol
+        TypeError : if dimension is int and sequence is not
+        ValueError : if all sequence indices are not between 0 and 2*dimension - 1
+        TypeError : if sequence containt something else than Expr
+        TypeError : if symbolic sequences have the wrong format
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from haarpy import haar_integral_symplectic
+
+        >>> d = Symbol("d")
+        >>> sequence_1, sequence_2 = (0, 1, 2, d, d+1, d+2), (0, 0, 0, d, d, d)
+        >>> haar_integral_symplectic((sequence_1, sequence_2), d)
+        1/(4*d*(d + 1)*(2*d + 1))
+
+        >>> d = 4
+        >>> sequence_1, sequence_2 = (0, 1, 2, d, d+1, d+2), (0, 0, 0, d, d, d)
+        >>> haar_integral_symplectic((sequence_1, sequence_2), d)
+        Fraction(1, 720)
+
+    See Also
+    --------
+        hyperoctahedral_transversal, weingarten_symplectic
     """
     if len(sequences) != 2:
         raise ValueError("Wrong sequence format")
diff --git a/haarpy/tests/test_circular_ensembles.py b/haarpy/tests/test_circular_ensembles.py
index bca6111..75cda0f 100644
--- a/haarpy/tests/test_circular_ensembles.py
+++ b/haarpy/tests/test_circular_ensembles.py
@@ -16,13 +16,14 @@
 """
 
 from math import prod
-from random import randint
+from random import seed, randint
 from fractions import Fraction
 from sympy import Symbol, simplify, factorial, factorial2
 from sympy.combinatorics import SymmetricGroup
 import pytest
 import haarpy as ap
 
+seed(137)
 d = Symbol("d")
 
 # The values in the dictionary below were verified against Monte Carlo simulations
diff --git a/haarpy/tests/test_permutation.py b/haarpy/tests/test_permutation.py
index 2efcb05..a435efc 100644
--- a/haarpy/tests/test_permutation.py
+++ b/haarpy/tests/test_permutation.py
@@ -17,10 +17,13 @@
 
 import pytest
 from math import factorial
+from random import seed, randint
 from itertools import product
+from fractions import Fraction
 from sympy import Symbol, simplify, factor, fraction
 import haarpy as ap
 
+seed(137)
 d = Symbol("d")
 
 hand_calculated_weingarten = {
@@ -130,9 +133,13 @@ def test_mobius_inversion_formula(size):
 )
 def test_weingarten_centered_permutation_hand_calculated(partition1, partition2, result_key):
     "Test Weingarten centered permutation function against hand calculated cases"
-    assert (
-        ap.weingarten_centered_permutation(partition1, partition2, d)
-        == hand_calculated_weingarten[result_key]
+    dimension = randint(10, 99)
+    assert ap.weingarten_centered_permutation(
+        partition1, partition2, d
+    ) == hand_calculated_weingarten[result_key] and ap.weingarten_centered_permutation(
+        partition1, partition2, dimension
+    ) == Fraction(
+        hand_calculated_weingarten[result_key].subs(d, dimension)
     )
 
 
@@ -148,11 +155,14 @@ def test_weingarten_centered_permutation_hand_calculated(partition1, partition2,
         ((3, 3, 2, 2, 3, 2), (3, 3, 2, 2, 1, 2)),
     ],
 )
-def test_haar_integral_permutation_weingarten(row_indices, column_indices):
+@pytest.mark.parametrize("symbolic", [True, False])
+def test_haar_integral_permutation_weingarten(row_indices, column_indices, symbolic):
     """Test haar integral for permutation matrices against the Weingarten
     sum as seen in Eq.(2.2) and (2.4) of `Collins and Nagatsu. Weingarten
     Calculus for Centered Random Permutation Matrices <https://arxiv.org/abs/2503.18453>`_
     """
+    dimension = d if symbolic else randint(10, 99)
+
     partition_row = tuple(
         tuple(index for index, value in enumerate(row_indices) if value == unique)
         for unique in set(row_indices)
@@ -166,7 +176,7 @@ def test_haar_integral_permutation_weingarten(row_indices, column_indices):
         ap.weingarten_permutation(
             partition_sigma,
             partition_tau,
-            d,
+            dimension,
         )
         for partition_sigma, partition_tau in product(
             ap.set_partitions(tuple(i for i, _ in enumerate(row_indices))),
@@ -176,10 +186,11 @@ def test_haar_integral_permutation_weingarten(row_indices, column_indices):
         and ap.partial_order(partition_tau, partition_column)
     )
 
-    num, denum = fraction(simplify(weingarten_integral))
-    weingarten_integral = factor(num) / factor(denum)
+    if symbolic:
+        num, denum = fraction(simplify(weingarten_integral))
+        weingarten_integral = factor(num) / factor(denum)
 
-    assert ap.haar_integral_permutation(row_indices, column_indices, d) == weingarten_integral
+    assert ap.haar_integral_permutation(row_indices, column_indices, dimension) == weingarten_integral
 
 
 @pytest.mark.parametrize(
diff --git a/haarpy/unitary.py b/haarpy/unitary.py
index 04fd050..a7ac36b 100644
--- a/haarpy/unitary.py
+++ b/haarpy/unitary.py
@@ -1,4 +1,4 @@
-# Copyright 2024 Polyquantique
+# Copyright 2025 Polyquantique
 
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -13,6 +13,14 @@
 # limitations under the License.
 """
 Unitary group Python interface
+
+References
+----------
+    [1] Collins, B. (2003). Moments and cumulants of polynomial random variables on unitarygroups,
+    the Itzykson-Zuber integral, and free probability. International Mathematics Research Notices,
+    2003(17), 953-982.
+    [2] Matsumoto, S. (2013). Weingarten calculus for matrix ensembles associated with compact
+    symmetric spaces. arXiv preprint arXiv:1301.5401.
 """
 
 from math import factorial, prod
@@ -36,12 +44,24 @@
 def representation_dimension(partition: tuple[int], unitary_dimension: Symbol) -> Expr:
     """Returns the dimension of the unitary group U(d) labelled by the input partition
 
-    Args:
-        partition (tuple[int]) : A partition labelling a representation of U(d)
+    Parameters
+    ----------
+        partition (tuple[int]) : a partition labelling a representation of U(d)
         unitary_dimension (Symbol) : dimension d of the unitary matrix U
 
-    Returns:
-        Expr : The dimension of the representation of U(d) labeled by the partition
+    Returns
+    -------
+        Expr : the dimension of the representation of U(d) labeled by the partition
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from haarpy import representation_dimension
+        >>> d = Symbol("d")
+        >>> representation_dimension((2,1,1), 4)
+        15
+        >>> representation_dimension((2,1,1), d)
+        d*(d/2 - 1/2)*(d - 2)*(d + 1)/4
     """
     conjugate_partition = tuple(
         sum(1 for part in partition if i < part) for i in range(partition[0])
@@ -71,16 +91,37 @@ def representation_dimension(partition: tuple[int], unitary_dimension: Symbol) -
 def weingarten_unitary(cycle: Union[Permutation, tuple[int]], unitary_dimension: Symbol) -> Expr:
     """Returns the Weingarten function
 
-    Args:
-        cycle (Permutation, tuple[int]): Permutation from the symmetric group or its cycle-type
-        unitary_dimension (Symbol): Dimension d of the unitary matrix U
+    The function works with both a permutation or a conjugacy class as a partition
 
-    Returns:
-        Expr: The Weingarten function
+    Parameters
+    ----------
+        cycle (Permutation, tuple[int]) : permutation from the symmetric group or its cycle-type
+        unitary_dimension (Symbol) : dimension d of the unitary matrix U
 
-    Raise:
-        TypeError: if unitary_dimension has the wrong type
-        TypeError: if cycle has the wrong type
+    Returns
+    -------
+        Expr : the Weingarten function
+
+    Raise
+    -----
+        TypeError : if unitary_dimension has the wrong type
+        TypeError : if cycle has the wrong type
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from haarpy import weingarten_unitary
+        >>> d = Symbol("d")
+        >>> weingarten_unitary(Permutation(2)(0, 1), 4)
+        Fraction(-1, 180)
+        >>> weingarten_unitary(Permutation(2)(0, 1), d)
+        -1/((d - 2)*(d - 1)*(d + 1)*(d + 2))
+        >>> weingarten_unitary((2, 1), d)
+        -1/((d - 2)*(d - 1)*(d + 1)*(d + 2))
+
+    See Also
+    --------
+        murn_naka_rule, representation_dimension, sympy.utilities.iterables.partitions
     """
     if not isinstance(unitary_dimension, (Expr, int)):
         raise TypeError("unitary_dimension must be an instance of int or sympy.Expr")
@@ -127,16 +168,35 @@ def weingarten_unitary(cycle: Union[Permutation, tuple[int]], unitary_dimension:
 def haar_integral_unitary(sequences: tuple[tuple[int]], unitary_dimension: Symbol) -> Expr:
     """Returns integral over unitary group polynomial sampled at random from the Haar measure
 
-    Args:
-        sequences (tuple[tuple[int]]): Indices of matrix elements
-        unitary_dimension (Symbol): Dimension of the unitary group
-
-    Returns:
-        Expr: Integral under the Haar measure
-
-    Raise:
-        ValueError: If sequences doesn't contain 4 tuples
-        ValueError: If tuples i and j are of different length
+    Parameters
+    ----------
+        sequences (tuple[tuple[int]]) : indices of matrix elements
+        unitary_dimension (Symbol) : dimension of the unitary group
+
+    Returns
+    -------
+        Expr : integral under the Haar measure
+
+    Raise
+    -----
+        ValueError : if sequences do not contain 4 tuples
+        ValueError : if tuples i and j are of different length
+
+    Examples
+    --------
+        >>> from sympy import Symbol
+        >>> from haarpy import haar_integral_unitary
+        >>> d = Symbol("d")
+        >>> seq_i, seq_j = (0, 1, 2), (0, 0, 1)
+        >>> seq_i_prime, seq_j_prime = (0, 1, 2), (0, 1, 0)
+        >>> haar_integral_unitary((seq_i, seq_j, seq_i_prime, seq_j_prime), 5)
+        Fraction(-1, 840)
+        >>> haar_integral_unitary((seq_i, seq_j, seq_i_prime, seq_j_prime), d)
+        -1/(d*(d - 1)*(d + 1)*(d + 2))
+
+    See Also
+    --------
+        stabilizer_coset, weingarten_unitary
     """
     if len(sequences) != 4:
         raise ValueError("Wrong tuple format")
diff --git a/pyproject.toml b/pyproject.toml
index aa4949a..5306055 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,2 +1,6 @@
+[build-system]
+requires = ["setuptools"]
+build-backend = "setuptools.build_meta"
+
 [tool.black]
 line-length = 100