# efficient generation of unitary divisors

A divisor $d\mid n$ is called unitary if it is coprime to its cofactor, i.e. $\gcd(d,\frac{n}d)=1$. Since for each prime power $p^k\| n$, we either have $p^k\|d$ or $p\nmid d$, which lead to the following seemingly-efficient code to generate all unitary divisors:

```
def unitary_divisors1(n):
return sorted( prod(p^k for p,k in f) for f in Subsets(factor(n)) )
```

For comparison, the naive approach by filtering the divisors of $n$ would be:

```
def unitary_divisors2(n):
return [d for d in divisors(n) if gcd(d,n//d)==1]
```

From theoretical perspective, `unitary_divisors1(n)`

is much more efficient than `unitary_divisors2(n)`

for powerful/squareful numbers $n$ and should be comparable for square-free numbers $n$. However, in practice `unitary_divisors1(n)`

loses badly to `unitary_divisors2(n)`

when $n$ is square-free (in which case every divisor of $n$ is unitary):

```
sage: N = prod(nth_prime(i) for i in (1..20))
sage: %time len(unitary_divisors1(N))
CPU times: user 26.8 s, sys: 64 ms, total: 26.9 s
Wall time: 26.9 s
1048576
sage: %time len(unitary_divisors2(N))
CPU times: user 672 ms, sys: 16 ms, total: 688 ms
Wall time: 688 ms
1048576
```

In this example with `N`

being the product of first 20 primes, `unitary_divisors1(N)`

is over 40 times slower than `unitary_divisors2(N)`

.

Is there a way to generate unitary divisors that works with efficiency close to `divisors(n)`

when $n$ is squarefree or almost squarefree, and takes advantage of the known structure of unitary divisors?

Essentially all of the time in

`unitary_divisors1`

is spent in`[... for f in Subsets(factor(n))]`

: on my machine,`[f for f in Subsets(factor(N))]`

took about the same amount of time as`unitary_divisors1(N)`

. Assembling and iterating over the subsets of size 20 is just going to take a while. Maybe the`Subsets`

iterator could be made more efficient.What should I use instead of Subsets? Maybe Gray code?

I wonder if something from the Python

`itertools`

module would be faster.Can

`Subsets`

construction be improved? I'd prefer to rely on its efficient implementation rather than trying to optimize simple stuff here and there.I don't know the

`Subsets`

code. It looks like all you need from that is to iterate over all subsets, which the`itertools`

powerset recipe provides, and it looks like it's faster. I wouldn't be surprised if`Subsets`

could be made faster, for example by writing it in Cython, but I don't know if anyone is working on it.