Algorithm to compute nth root (radical) $\sqrt[n]{p(X)}$ of polynomial

Question

I am making a symbolic computation library which supports symbolic polynomials (both univariate and multivariate) and among other things I would like to support (possibly truncated) nth root (radical) computation of a (univariate or multivariate) polynomial since this computation is needed for some other things.

So lets say we have a polynomial $p(X) = \sum_{i}{a_{i}X^{i}}$ and I need an algorithm for $\sqrt[n]{p(X)}$.

I have found this post on computing the square root of a polynomial but I dont quite get it but also I would like an algorithm for nth root.

Library can compute all primitive operations of polynomials (both univariate and multivariate) on a general ring of coeffcients (ie add/sub/mul/div/mod/pow/derivative). (For reference library is on github)

Please provide a step-by-step algorithm in your answer or a link to such algorithm so I can test it.

PS. I tried using a variant of Newton's method for computing nth root and adapted it to polynomial operations but result is completely off, not even close.

Note: If polynomial is not a perfect nth power then a truncated power series can be computed up to some limit (eg truncated Taylor series) as an approximation. For example the square root algorithm on PlanetMath computes the Taylor series if poly is not perfect power of 2.

if poly is not a perfect nth power or power series approximation cannot be computed is fine for me if algorithm throws an error.

Thank you.

Answer 1

I managed to decipher the square root algorithm referenced in the question and generalise it to n-th root algorithm.

Algorithm: Compute nth root/radical of polynomial $p(X) = \sum_{i}a_iX^i$

Preliminary: Arrange the polynomial from lowest degree term up to highest degree term or if multivariate do same according to some monomial ordering (eg LEX). $LT(p(X))$ refers to the leading term according to monomial order and $TT(p(X))$ refers to the tail / last term according to monomial order. Algorithm works both if LT is used instead of TT below but with TT power series approximation (if $p(X)$ is not a perfect nth power) is easier to compute up to any desired accuracy. $maxdeg(p(X))$ refers to maximum power that exists in $p(X)$ in any variable (if multivariate). $maxterms$ is a user-defined limit and refers to maximum amount of terms to be computed if a power series approximation is the result. Eg up to $6$ terms of Taylor series expansion (if possible).

• Check: If $p(X)=0$ or $p(X)=1$ return $p(X)$
• Init: $i \leftarrow 0, nterms \leftarrow 0, \text{ } r_{i}(X) \leftarrow \sqrt[n]{TT(p(X))}$
• Step 1: if $nterms \geq maxterms$ return $r_{i}(X)$ else goto to Step 2
• Step 2: $d_{i}(X) \leftarrow p(X)-r_{i}(X)^{n}$
• step 3: if $d_{i}(X)=0$ return $r_{i}(X)$ else goto to Step 4
• Step 4: $q_{i}(X) \leftarrow TT(d_{i}(X)) \text{ }/\text{ } TT(nr_{i}(X)^{n-1})$
• Step 5: if $q_{i}(X)=0$ return $r_{i}(X)$ else goto to Step 6
• Step 6: $r_{i+1}(X) \leftarrow r_{i}(X) + q_{i}(X)$
• Step 7: if $n \cdot maxdeg(r_{i+1}(X)) > maxdeg(p(X))$ this means that $p(X)$ is not a perfect nth power and that a power series approximation is being computed. $nterms \leftarrow nterms+1$
• Step 8: set $i \leftarrow i+1$ and goto Step 1

Optionally one can normalise returned $r(X)$ to have positive leading term coefficient if $n$ is a multiple of $2$, since in that case both $r(X)$ and $-r(X)$ are roots.

My library is updated on github to use the above algorithm for computing nth radicals of polynomials.

Above algorithm correctly passes the example tests:

$p \in Q[x]$

$\sqrt{x^2} = x$

$\sqrt{(x^2)^2} = x^2$

$\sqrt[5]{(x+1)^5} = x+1$

$\sqrt{9x^4+6x^3-11x^2-4x+4} = -3x^2-x+2$

$\sqrt{x+1}=\frac{7}{256}x^5-\frac{5}{128}x^4+\frac{1}{16}x^3-\frac{1}{8}x^2+\frac{1}{2}x+1$ (truncated Taylor series up to maxterms)

$p \in Q[x, y]$

$\sqrt{4x^2-12xy+9y^2} = -2x+3y$

Answer 2

Suggestion:

Write $p(x)$ as $p_0(x)(1+x\,q(x))$ and compute

$$\sqrt[n]{p(x)}=\sqrt[n]{p_0(x)}\sum_{k=0}^d\binom{1/n}kx^k(q(x))^k.$$

You don't need to fully compute the powers of $q$, you can stop at the highest power that does not exceed $d-k$.