We consider a Hilbert space that consists of two Transmon “qubits” (qudits, actually) with `nq` levels, and a cavity with `nc` levels. Given the dimension of the total Hilbert space, `n = nq*nq*nc`, we want to find `nq` and `nc`, i.e. the tensor structure of the full Hilbert space. We demand that `nq <= nc`.

The decomposition is based on the prime factors of `n`:

```def prime_factors(n):
"""Return the prime decomposition of n"""
i = 2
factors = []
while i * i <= n:
if n % i:
i += 1
else:
n //= i
factors.append(i)
if n > 1:
factors.append(n)
return factors

>>> n = 6*6*20
>>> factors = prime_factors(n)
>>> print factors
[2, 2, 2, 2, 3, 3, 5]
```

Now we put them in two list: the prime factors of the qubit Hilbert space, and the prime factors of the cavity Hilbert space

```qubit_factors = []
cavity_factors = []

for factor in set(factors):
count = factors.count(factor)
while count >= 2:
qubit_factors.append(factor)
count += -2
if count > 0:
cavity_factors.append(factor)

>>> qubit_factors
[2, 2, 3]

>>> cavity_factors

```

We demand that the qubit dimension is not larger than the cavity dimension. In order to assure this, we start moving the smallest prime factors from the qubit to the cavity

```import numpy as np
while np.array(qubit_factors).prod() > np.array(cavity_factors).prod():
factor = qubit_factors.pop(0)
cavity_factors = [factor, factor] + cavity_factors

>>> qubit_factors
[2, 3]

>>> cavity_factors
[2, 2, 5]

>>> nq = np.array(qubit_factors).prod()
>>> nq
6

>>> nc = np.array(cavity_factors).prod()
>> nc
20
```

The total algorithm looks like this (using `prime_factors`, and also allowing to specify the cavity dimension manually):

```def qubit_cavity_decomposition(n, nc=None):
"""
Return (nq, nc) of a decomposition n = nq*nq*nc, where nc >= nq
"""
if nc is not None:
n = n // nc
factors = prime_factors(n)
qubit_factors = []
cavity_factors = []
for factor in set(factors):
count = factors.count(factor)
while count >= 2:
qubit_factors.append(factor)
count += -2
if count > 0:
cavity_factors.append(factor)
# Ensure nc >= nq
if nc is None:
while np.array(qubit_factors).prod() > np.array(cavity_factors).prod():
factor = qubit_factors.pop(0)
cavity_factors = [factor, factor] + cavity_factors
nc = np.array(cavity_factors).prod()
nq = np.array(qubit_factors).prod()
return nq, nc

>>> qubit_cavity_decomposition(5*5*10)
(5, 10)

>>> qubit_cavity_decomposition(5*5*25)
(5, 25)

>>> qubit_cavity_decomposition(10*10*20)
(10, 20)

>>> qubit_cavity_decomposition(8*8*25)
(10, 16)

>>> qubit_cavity_decomposition(8*8*25, nc=25)
(8, 25)

>>> qubit_cavity_decomposition(8*8*30)
(8, 30)
```