I just found a really cool polynomial decomposition algorithm which works beautifully for polys over GF(p) whose degree is <p, i.e., particularly those which are permutation polynomials. I believe this algorithm is due to Ms. Susan Landau (now @ WPI), although I can't be 100% sure about this w/o asking her. Decompose monic f(x) = g(h(x)), where deg(f)=n=r*s=deg(g)*deg(h), where g(x) and h(x) are also monic. Compute fbar(x) = reverse(f(x)) = x^n*f(1/x). Compute hbar(x) = r'th root of fbar(x) to precision (degree) s=n/r. This is well defined if r<p because hbar(0)=fbar(0)=1, as f(x) is monic, so fbar(0)=1. Use Newton iteration on power series approx. to r'th root of fbar(x). Compute h(x) = reverse(hbar(x)) = x^s*hbar(1/x). Express f(x) in base-h(x) notation using classical number base conversion algorithm. Each of the "digits" is a polynomial of degree <deg(h)=s. If each of the digits is actually a *constant* -- i.e., free of the variable x, then these "digits" are the coefficients of a polynomial g(x), and we have shown that f(x)=g(h(x)). Note that decomposition is essentially *cheaper* than factoring. Who knew? Cool ! (Note that although the paper looks to be behind a paywall, I was still able to download it, and can email you a copy if necessary.) J. Symbolic Computation (1990) 9, 281-299 Functional Decomposition of Polynomials: the Tame Case JOACHIM VON ZUR GATHEN Department of Computer Science, University of Toronto Toronto, Ontario MSS 1A4, Canada Using fast polynomial arithmetic as quoted above, we can implement the algorithm of Kozen & Landau (1989) very efficiently. For a polynomial f = x^n+a_n-1 x^(n-1) +...+ a_0 in F[x], we denote by fbar = a_0 x^n +...+a_n-1 x + 1 = x^n f(1/x) the reversal of f. Algorithm: Univariate decomposition. Input: f in F[x], monic of degree n = r s, and r in N, not divisible by char(F). Output: The unique normal decomposition (g,h) of f with deg g = r, if such a decomposition exists, and "no decomposition" otherwise. 1. Let fbar be the reversal of f, and compute hbar in F[x] of degree less than s with hbar^r = fbar mod x^s and h(0) = 1. Let h = x^s hbar(1/x) in F[x]. 2. Compute b_0,..., b_r in F[x] as in "Taylor expansion". 3. If b_0,..., b_r in F, set g = sum(0<=i<=r, b_i x^i) in F[x] and return (g,h). Otherwise return "no decomposition".