Ask Your Question

efficient generation of unitary divisors

asked 2021-11-20 21:48:51 +0200

Max Alekseyev gravatar image

updated 2021-11-20 21:56:41 +0200

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

sage: %time len(unitary_divisors2(N))                                                                                                                                                                              
CPU times: user 672 ms, sys: 16 ms, total: 688 ms
Wall time: 688 ms

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 flag offensive close merge delete



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.

John Palmieri gravatar imageJohn Palmieri ( 2021-11-20 23:03:47 +0200 )edit

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

Max Alekseyev gravatar imageMax Alekseyev ( 2021-11-20 23:44:47 +0200 )edit

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

John Palmieri gravatar imageJohn Palmieri ( 2021-11-21 01:17:30 +0200 )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.

Max Alekseyev gravatar imageMax Alekseyev ( 2021-11-22 21:17:40 +0200 )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.

John Palmieri gravatar imageJohn Palmieri ( 2021-11-23 01:44:22 +0200 )edit

2 Answers

Sort by ยป oldest newest most voted

answered 2021-11-21 03:00:40 +0200

Max Alekseyev gravatar image

updated 2021-11-21 03:08:08 +0200

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)) )
edit flag offensive delete link 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):
    r = radical(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).

slelievre gravatar imageslelievre ( 2021-11-21 20:12:59 +0200 )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.

Max Alekseyev gravatar imageMax Alekseyev ( 2021-11-22 21:12:02 +0200 )edit

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))])
slelievre gravatar imageslelievre ( 2021-11-23 01:23:13 +0200 )edit

answered 2021-11-21 01:25:52 +0200

updated 2021-11-21 01:27:49 +0200

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, also appearing as the powerset recipe at

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.)

edit flag offensive delete link more


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

Max Alekseyev gravatar imageMax Alekseyev ( 2021-11-21 03:03:52 +0200 )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.

John Palmieri gravatar imageJohn Palmieri ( 2021-11-21 04:15:02 +0200 )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).

Max Alekseyev gravatar imageMax Alekseyev ( 2021-11-22 21:15:31 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower


Asked: 2021-11-20 21:48:51 +0200

Seen: 255 times

Last updated: Nov 21 '21