This method will return the maximum number of edge disjoint Hamiltonian cycles. This is certainly not the best possible formulation.

```
def disjoint_hamiltonian_cycles(g, solver=None, verbose=0,
*, integrality_tolerance=1e-3):
"""
"""
from sage.numerical.mip import MixedIntegerLinearProgram
from sage.numerical.mip import MIPSolverException
from sage.graphs.graph import Graph
UB = min(g.degree()) // 2
if not UB:
raise EmptySetError("the given graph is not Hamiltonian")
eps = 1 / (2*Integer(g.order()))
x = next(g.vertex_iterator())
while True:
# Search for UB edge disjoint Hamiltonian cycles
p = MixedIntegerLinearProgram(solver=solver)
f = p.new_variable(binary=True)
r = p.new_variable(nonnegative=True)
# An edge belongs to a unique cycle
for e in g.edge_iterator(labels=False):
p.add_constraint(p.sum(f[frozenset(e), i] for i in range(UB)) <= 1)
# All the vertices of cycle i have degree 2
for i in range(UB):
for v in g:
p.add_constraint(p.sum(f[frozenset((u, v)), i] for u in g.neighbor_iterator(v)), min=2, max=2)
# Connectivity constraints based on the maximum average degree
for i in range(UB):
# r is greater than f
for u, v in g.edge_iterator(labels=None):
p.add_constraint(r[u, v, i] + r[v, u, i] - f[frozenset((u, v)), i], min=0)
# no cycle which does not contain x
for v in g:
if v != x:
p.add_constraint(p.sum(r[u, v, i] for u in g.neighbor_iterator(v)), max=1 - eps)
try:
p.solve(log=verbose)
f_val = p.get_values(f, convert=bool, tolerance=integrality_tolerance)
tsp = [Graph() for i in range(UB)]
for i in range(UB):
tsp[i].add_vertices(g.vertex_iterator())
tsp[i].set_pos(g.get_pos())
for (e, i), b in f_val.items():
if b:
tsp[i].add_edge(e)
return tsp
except MIPSolverException:
UB -= 1
if not UB:
raise EmptySetError("the given graph is not Hamiltonian")
```

You get:

```
sage: G = graphs.CompleteGraph(5)
sage: T = disjoint_hamiltonian_cycles(G)
sage: T
[Graph on 5 vertices, Graph on 5 vertices]
sage: [t.edges(labels=False, sort=True) for t in T]
[[(0, 1), (0, 3), (1, 4), (2, 3), (2, 4)],
[(0, 2), (0, 4), (1, 2), (1, 3), (3, 4)]]
sage: G = graphs.CompleteGraph(7)
sage: T = disjoint_hamiltonian_cycles(G)
sage: [t.edges(labels=False, sort=True) for t in T]
[[(0, 1), (0, 4), (1, 6), (2, 3), (2, 5), (3, 6), (4, 5)],
[(0, 2), (0, 5), (1, 3), (1, 4), (2, 6), (3, 5), (4, 6)],
[(0, 3), (0, 6), (1, 2), (1, 5), (2, 4), (3, 4), (5, 6)]]
```

We don't have such method in Sagemath. This is a hard optimisation problem. You can use integer linear programming to model and solve it.

Also, your question is not entirely clear. You want the maximal number of pairwise disjoint hamiltonian cycles? I doubt that starting with a random hamiltonian cycle (such as given by

`G.hamiltonian_cycle()`

), deleting it and repeating the procedure would give this maximum. If you want a maximal collection of pairwise disjoint hamiltonian cycles (ie such that what remains is not hamiltonian) your suggestion would certainly works.@vdelecroix yes, I want the maximum number of disjoint hamiltonian cycles in a graph. Yes, my method would work for maximal, but I actually wanted maximum of such cycles. edited the question

@David Coudert thanks, but could you just slightly elaborate?