I will try to sketch an algorithm which avoids the computation of the Galois group, since such a computation is very expensive.

The idea is to look at the factorisations of the reductions of $f$ modulo sufficiently many primes $p$. This is the same idea as in the computation of the Galois group, but the point is that we don't need the whole Galois group, only its maximal abelian quotient, so things should be simpler.

We may assume that $f \in \mathbb{Z}[x]$ is monic and irreducible. Let $K = \mathbb{Q}[x]/(f)$ be the number field defined by $f$ and let $L$ be the Galois closure of $K$. We embed both $L$ and $\mathbb{Q}^{\mathrm{ab}}$ into some algebraic closure $\mathbb{Q}^{\mathrm{alg}}$. Let $K^{\mathrm{ab}} = K \cap \mathbb{Q}^{\mathrm{ab}}$ and $L^{\mathrm{ab}} = L \cap \mathbb{Q}^{\mathrm{ab}}$ be the maximal abelian subfields of $K$ and $L$ (sorry for the unusual notation...). Let $f = g_1 g_2 \cdots g_s$ be the factorisation of $f$ over $\mathbb{Q}^{\mathrm{ab}}$. By Galois theory, we know that each $g_i$ has coefficients in $L^{\mathrm{ab}} \cap K = K^{\mathrm{ab}}$.

**Definition.** A *bound* for a finite abelian extension $F$ of $\mathbb{Q}$ is a pair $(N,G)$ where $N \geq 1$ is an integer and $G$ is a subgroup of $(\mathbb{Z}/N\mathbb{Z})^\times$ such that $F \subset \mathbb{Q}(\zeta_N)^G$.

Note that for any abelian $F/\mathbb{Q}$ we always have the trivial bound $(|\Delta_F|,\{1\})$ where $\Delta_F$ is the (absolute) discriminant of $F$. Also, it is not hard to see that among all the bounds for $F$ there is a unique minimal one $(N,G)$. The integer $N$ is the conductor of $F$, and $G$ is the norm group, namely the kernel of $(\mathbb{Z}/N\mathbb{Z})^\times \to \mathrm{Gal}(F/\mathbb{Q})$.

So starting with $f$, we know that $K^{\mathrm{ab}}$ is bounded by $(|\Delta_K|,\{1\})$ (it is possible to simply compute the discriminant of $f$, but this may lead to a worse bound). So we start with $(N,G)=(|\Delta_K|,\{1\})$ and we want to improve the bound. In fact, since $K^{\mathrm{ab}}$ is contained in both $K$ and $\mathbb{Q}(\zeta_N)$, we may take $G$ to be the image of the power map $x \mapsto x^\delta$ in $(\mathbb{Z}/N\mathbb{Z})^\times$, where $\delta= \operatorname{gcd}(\deg(f),\varphi(N))$.

Now we look at the reduction of $f$ modulo some primes, similarly as in Algorithm 4.4.3 in Cohen, *Advanced topics in computational number theory*. Let $p$ be a prime not dividing $N$. Let $\bar{f}$ be the reduction of $f$ mod $p$, and let $\bar{f}=f_1 f_2 \cdots f_r$ be its factorization in $\mathbb{F}_p[x]$. The $f_i$ are pairwise distinct monic irreducible polynomials, say of degree $d_i$. If $K$ were abelian then we would know that $d_1=\ldots=d_r=d$ where $d$ is the order of $p$ in $\mathrm{Gal}(K/\mathbb{Q})$ (seen as a quotient of $(\mathbb{Z}/N\mathbb{Z})^\times$). We are not in this case, but still know by multiplicativity of the inertia degree that the order $d$ of $p$ in $\mathrm{Gal}(K^{\mathrm{ab}}/\mathbb{Q})$ divides each $d_i$. So we adjoin to the group $G$ the element $p^h$ where $h$ is the gcd of $d_1,\ldots,d_r$.

Doing this for enough primes $p$, we should get a reasonable bound $(N,G)$. Now we can replace $(N,G)$ by a minimal equivalent bound $(N',G')$ where $N'$ divides $N$. (Precisely, $N'$ is the minimal divisor of $N$ such that $G$ contains the kernel of reduction mod $N'$, and $G'$ is the image of $G$ mod $N'$.) At this point it remains to factor $f$ over $K'=\mathbb{Q}(\zeta_N)^G$. By the theory of Gauss sums, we know a primitive element $\alpha$ in $K'$ (written in terms of $\zeta_N$) and its minimal polynomial, so we simply factor $f$ over the number field $\mathbb{Q}(\alpha)$.

A natural question is how many primes to use. In fact, if we are only interested in having an algorithm, this doesn't matter: we could simply have used the trivial bound. To get a more efficient algorithm, we may use some primes but it's not clear to me, for example, how many primes are needed to ensure that we get $K'=K^{\mathrm{ab}}$. In general, it is known that the norm group of an abelian extension is generated by the prime powers $p^d$ as above when $p$ runs through the primes less than an explicit bound $B$, see Algorithm 4.4.5 in Cohen. But the bound $B$ is a priori exponential in the size of the discriminant, while there is a much better polynomial bound assuming GRH, so some choices need to be made.

Another difficulty is that the gcd $h$ used above need not equal $d$. Maybe the Cebotarev theorem (an effective version!) could be used to prove this doesn't matter, but this may be a non-trivial task and I haven't thought about it. So I would make the following compromise: simply use a certain amount of primes without caring whether $K'=K^{\mathrm{ab}}$ or not. At the end, if really needed, one could determine over which minimal extension the factorisation of $f$ actually takes place. The "right" number of primes will depend on the parameters but could be determined experimentally or by finding some heuristics.