# 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?

edit retag close merge delete

1

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.

( 2021-11-20 23:03:47 +0100 )edit

What should I use instead of Subsets? Maybe Gray code?

( 2021-11-20 23:44:47 +0100 )edit

I wonder if something from the Python itertools module would be faster.

( 2021-11-21 01:17:30 +0100 )edit

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

( 2021-11-22 21:17:40 +0100 )edit

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.

( 2021-11-23 01:44:22 +0100 )edit

Sort by » oldest newest most voted

My best shot so far is iterating over the divisors of the radical of $n$ and saturating prime powers in them:

def unitary_divisors4(n):
fn = factor(n)
D = [(p,k-1) for p,k in fn if k>1]
return sorted( prod(p^k for p,k in D if d%p==0)*d for d in divisors(prod(p for p,_ in fn)) )

more

Possible variations using radical(n) for prod(p for p, _ in factor(n)):

def unitary_divisors5(n):
fn = factor(n)
D = [(p, k - 1) for p, k in fn if k > 1]
return sorted(d*prod(p^k for p, k in D if d % p == 0) for d in divisors(radical(n)))

def unitary_divisors6(n):
D = factor(n // r)  # or list(factor(n // r))
return sorted(d*prod(p^k for p, k in D if d % p == 0) for d in divisors(r))


Computing both radical(n) and factor(n) or factor(n // r) might end up being a waste of time though...

Not sure whether using x.radical(), x.factor(), x.divisors() instead of radical(x), factor(x), divisors(x) would gain anything (assuming only Sage integers, never Python integers, are used).

( 2021-11-21 20:12:59 +0100 )edit

Since we need factorization of $n$ anyway, it's worth to use it for computing the radical (rather than calling function radical, which would internally factor n again) and factorization of n // r.

( 2021-11-22 21:12:02 +0100 )edit
1

Fully agree. A few more things to try...

With fn = factor(n), you can use fn.radical_value() to get the radical (but prod(p for p, _ in fn) works too, of course, and may be faster).

Using fn = list(factor(n)) may win over fn = factor(n) (but then there's no fn.radical_value()).

Using sorted([a for a in b]) may be faster than sorted(a for a in b).

Above all, subsets is a lot faster than Subsets.

One more variation:

def unitary_divisors7(n):
return sorted([prod(p^k for p, k in f)
for f in subsets(factor(n))])

( 2021-11-23 01:23:13 +0100 )edit

This is not a great answer but it's too long to fit in a comment. With N as in the question, it's faster than unitary_divisors1 but slower than unitary_divisors2. Partly taken from the top-ranked answer at https://stackoverflow.com/questions/1..., also appearing as the powerset recipe at https://docs.python.org/3/library/ite....

sage: from itertools import chain, combinations
sage: def unitary_divisors3(n):
....:     s = factor(n)
....:     C = chain.from_iterable(combinations(s, r) for r in range(len(s)+1))
....:     return sorted(prod(p^k for p,k in f) for f in C)
....:

sage: %time L = len(unitary_divisors1(N))
CPU times: user 34.1 s, sys: 242 ms, total: 34.3 s
Wall time: 44.1 s
sage: %time L = len(unitary_divisors2(N))
CPU times: user 1.26 s, sys: 54.2 ms, total: 1.31 s
Wall time: 1.53 s
sage: %time L3 = len(unitary_divisors3(N))
CPU times: user 4.64 s, sys: 74.2 ms, total: 4.71 s
Wall time: 5.88 s


(My system is pretty heavily loaded right now; the timings for the first two were much faster earlier, but these are good for comparisons regardless.)

more

Thanks. This loses a bit to unitary_divisors4 from my answer for squarefree $n$, but a bit better for squareful ones.

( 2021-11-21 03:03:52 +0100 )edit

How fast is Sage at testing whether something is squarefree or similar? You could have several cases for your function depending on properties of n. Certainly in the code I was looking at, factoring n was not the slowest part, so this might make sense.

( 2021-11-21 04:15:02 +0100 )edit

I think testing for square-freeness/-fulness require factorization of a number, and since we need factorization of n anyway, it makes sense to use it as much as possible (rather than calling Sage functions that would internally factor n again).

( 2021-11-22 21:15:31 +0100 )edit