Add function to get unitary divisors

[jm58660]

Add a function to get unitary divisors of an integer / of an element of UFD. For definition see https://en.wikipedia.org/wiki/Unitary_divisor. Trivial way to make this function is

def unitary_divisors(n):
    return [prod([p[0]^p[1] for p in f]) for f in Set(factor(n)).subsets()]

Before proceeding we should discuss about what to put in the global namespace. For that discussion see https://groups.google.com/forum/#!topic/sage-devel/DHMCJydPqJU.

CC: @bgrenet

Component: number theory

Issue created by migration from https://trac.sagemath.org/ticket/20369