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 asunitary_divisors1(N)
. Assembling and iterating over the subsets of size 20 is just going to take a while. Maybe theSubsets
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 theitertools
powerset recipe provides, and it looks like it's faster. I wouldn't be surprised ifSubsets
could be made faster, for example by writing it in Cython, but I don't know if anyone is working on it.