def METHOD_NAME(self, f):
    """
    Compute square-free decomposition of the monic ``f`` in ``GF(q)[X]``.
    Notes
    =====
    Uses a modified version of Musser's algorithm for square-free
    decomposition of univariate polynomials over finite fields.
    References
    ==========
    * :cite:`Geddes1992algorithms`, algorithm 8.3
    * :cite:`Musser1971factor`, p.54, algorithm V (modified)
    """
    domain = self.domain
    n, factors, p = 1, [], int(domain.characteristic)
    m = int(domain.order // p)
    while not f.is_ground:
        df = [f.diff(x) for x in self.gens]
        if any(_ for _ in df):
            g = f
            for q in df:
                g = self.gcd(g, q)
            h, f, i = f // g, g, 1
            while h != 1:
                g = self.gcd(f, h)
                h //= g
                if not h.is_ground:
                    factors.append((h, i*n))
                f //= g
                h = g
                i += 1
        n *= p
        g = self.zero
        for monom, coeff in f.items():
            g[tuple(_ // p for _ in monom)] = coeff**m
        f = g
    return factors