Polynomial factorization and roots calculator

Integer polynomials in one variable are expressions that include a variable named x, integer numbers and the addition, subtraction and multiplication operators.

They can be expressed in the form: f(x) = f₀ + f₁⁢x + f₂⁢x² + ... + f_n⁢xⁿ.

The number n is the degree of the polynomial. The coefficient f_n is the leading coefficient and the coefficient f₀ is the trailing coefficient. They can be written as deg(f), lc(f) and tc(f) respectively.

Let f(x) and g(x) be two polynomials such that deg(f) ≥ deg(g). We can define addition, subtraction and multiplication as follows:

f(x) + g(x) = h(x) means h_i = f_i + g_i if i ≤ deg(g) and h_i = f_i otherwise.

f(x) − g(x) = h(x) means h_i = f_i − g_i if i ≤ deg(g) and h_i = f_i otherwise.

g(x) − f(x) = h(x) means h_i = g_i − f_i if i ≤ deg(g) and h_i = −f_i otherwise.

f(x) ⁢ g(x) = h(x) means h_i = f₀ ⁢ g_i + f₁ ⁢ g_{i − 1} + ... + f_i ⁢ g₀ if i ≤ deg(g), h_i = f_{i − deg(g)} ⁢ g_deg(g) + f_{i + 1 − deg(g)} ⁢ g_{deg(g) − 1} + ... + f_i ⁢ g₀ otherwise.

The factorization of integer polynomials is a process to find one or more irreducible polynomials whose product is the original polynomial. An irreducible polynomial cannot be expressed as a product of two or more integer polynomials.

For example: x⁴ − 1 = (x² + 1) ⁢ (x + 1) ⁢ (x − 1)

It can be shown that any integer polynomial can be factored in only one way in irreducible polynomials.

Multiplication of polynomials modulo a prime number can be performed in the usual way by multiplying the different terms of the polynomial and then adding the coefficients of the same degree. Finally, the coefficients are reduced modulo that prime.

For example, the product of 3x² + 5x + 1 and 6x² + 4x + 3 modulo 7 is 18x⁴ + (30+12)x³ + (9+20+6)x² + (15+4)x + 3 modulo 7 which equals 4x⁴ + 5x + 3

It can be shown that the polynomials modulo a prime can be factored into the leading coefficient and monic prime polynomials in only one way (monic polynomials have the leading coefficient equal to one.)

It can also be shown that if there are no repeated factors, the polynomial can be factored modulo a power of that prime in only one way.

Use the upper input box to enter the polynomial expression and the lower input box to enter a numerical expression for the modulus, which must be either zero or an integer number greater than 1 of the form prime number raised to an exponent. You can just evaluate the polynomial or evaluate and factor it, by pressing the corresponding button.

Example for integer polynomial factorization: to factor x³⁰ − 1, type xcaret^30minus-1 in the upper input box and 0 in the lower input box. Then press the factor button.

Example for polynomial modulo a prime factorization: copy any of the expressions below in the upper input box and type the number 211 (a prime number) in the lower input box. Then press the button "Factor expression".

• 6x caret^8+x caret^5+3
• 2*left parenthesis(x+6right parenthesis)*left parenthesis(x-5right parenthesis)+xx caret^4+23x

The exponentiation symbol is not present in some mobile devices, so two asterisks ** can by typed as the exponentiation operator. So, equivalent expressions would be:

• 6x**8+x**5+3
• 2*left parenthesis(x+6right parenthesis)*left parenthesis(x-5right parenthesis)+xx**4+23x

You can type a dot (.) and the application will replace it by "x caret^". This saves a lot of time in mobile devices because there is no need to switch between alphabetical and numerical keyboards for simple polynomials.

For the first example you would type:

• 6.8+.5+3

Notice that factoring large degree polynomials will take a lot of time. That's why the applet accepts polynomials of degree up to 1000.

The program has three options for output: with the option Pretty print you can see the exponents as superscripts and the roots have the traditional, printed format; the option TeX shows the output in this format, which is used in many mathematical forums; finally, the option Pari-GP is used when you have to copy the data to another application. The output is compatible with the application Pari-GP, only minor changes are required to copy the data to other math programs.

After the applet finishes factoring, it multiplies the factors in order to validate if they are equal to the input polynomial.

When the modulus box is zero, the application finds the roots of the polynomial factors using radical expressions. The program supports degrees up to 5. By Abel-Ruffini theorem, some polynomials of degree 5 cannot be expressed by radical expressions.

You can use any letter to represent the variable. Multiple variables are not accepted. The program always use the letter x to show the output.

You can also enter expressions that use the following operators and parentheses:

• plus sign+ for addition
• minus sign- for subtraction
• asterisk* for multiplication
• slash/ for division
• percent symbol% for remainder of polynomial long division
• caret or two asterisks^ or ** for exponentiation
• Ans: retrieves the last answer.

Operators and functions for the polynomial (upper box) only:

• Gcd(f,g, ...): Greatest common divisor of these polynomials. Example: GCD(x-1, x^2-1) = x-1.
• Lcm(f,g, ...): Least common multiple of these polynomials. Example: LCM(x+1, x-1, x^2-1) = x^2-1.
• Der(f): Derivative of this polynomial.
• LongDiv(f, g): Long division of these polynomials.
• Random(a,b,c,d): Random polynomial whose degree is selected between the first two parameters and the coefficients between the last two parameters.

Operators and functions for the modulus (lower box) only:

• <, ==, >; <=, >=, != for comparisons. The operators return zero for false and -1 for true.
• AND, OR, XOR, NOT for binary logic. The operations are done in binary (base 2). Positive (negative) numbers are prepended with an infinite number of bits set to zero (one).
• SHL or <<: When b ≥ 0, a SHL b shifts a left the number of bits specified by b. This is equivalent to a × 2^b. Otherwise, a SHL b shifts a right the number of bits specified by −b. This is equivalent to floor(a / 2^−b). Example: 5 SHL 3 = 40.
• SHR or >>: When b ≥ 0, a SHR b shifts a right the number of bits specified by b. This is equivalent to floor(a / 2^b). Otherwise, a SHR b shifts a left the number of bits specified by −b. This is equivalent to a × 2^−b. Example: -19 SHR 2 = -5.
• n!: factorial (n must be greater than or equal to zero). Example: 6! = 6 × 5 × 4 × 3 × 2 = 720.
• n!! ... !: multiple factorial (n must be greater than or equal to zero). It is the product of n times n − k times n − 2k ... (all numbers greater than zero) where k is the number of exclamation marks. Example: 7!! = 7 × 5 × 3 × 1 = 105.
• p#: primorial (product of all primes less or equal than p). Example: 12# = 11 × 7 × 5 × 3 × 2 = 2310.
• B(n): Previous probable prime before n. Example: B(24) = 23.
• F(n): Fibonacci number F_n from the sequence 0, 1, 1, 2, 3, 5, 8, 13, 21, etc. where each element equals the sum of the previous two members of the sequence. Example: F(7) = 13.
• L(n): Lucas number L_n = F_n-1 + F_n+1
• N(n): Next probable prime after n. Example: N(24) = 29.
• P(n): Unrestricted Partition Number (number of decompositions of n into sums of integers without regard to order). Example: P(4) = 5 because the number 4 can be partitioned in 5 different ways: 4 = 3+1 = 2+2 = 2+1+1 = 1+1+1+1.
• Gcd(m,n, ...): Greatest common divisor of these integers. Example: GCD(12, 16) = 4.
• Lcm(m,n, ...): Least common multiple of these integers. Example: LCM(12, 16, 24) = 48.
• Modinv(m,n): inverse of m modulo n, only valid when m and n are coprime, meaning that they do not have common factors. Example: Modinv(3,7) = 5 because 3 × 5 ≡ 1 (mod 7)
• Modpow(m,n,r): finds mⁿ modulo r. Example: Modpow(3, 4, 7) = 4, because 3⁴ ≡ 4 (mod 7).
• Jacobi(m,n): obtains the Jacobi symbol of m and n. When the second argument is prime, the result is zero when m is multiple of n, it is one if there is a solution of x² ≡ m (mod n) and it is equal to −1 when the mentioned congruence has no solution.
• Random(m,n): integer random number between m and n.
• Abs(n): absolute value of n.
• IsPrime(n): returns zero if n is not probable prime, -1 if it is. Example: IsPrime(5) = -1.
• Sqrt(n): Integer part of the square root of the argument.
• Iroot(n,r): Integer r-root of the first argument. Example: Iroot(8, 3) = 2.
• NumDigits(n,r): Number of digits of n in base r. Example: NumDigits(13, 2) = 4 because 13 in binary (base 2) is expressed as 1101.
• SumDigits(n,r): Sum of digits of n in base r. Example: SumDigits(213, 10) = 6 because the sum of the digits expressed in decimal is 2+1+3 = 6.
• RevDigits(n,r): finds the value obtained by writing backwards the digits of n in base r. Example: RevDigits(213, 10) = 312.

You can use the prefix 0x for hexadecimal numbers, for example 0x38 is equal to 56.

Division of polynomials modulo a prime number

The polynomial division requires multiple modular divisions where the divisor is the leading coefficient of the divisor polynomial.

To compute the modular division a / b (mod p), first the modular multiplicative inverse c is found. This number has the property b⁢c ≡ 1 (mod p) and it can be found using the extended Euclidean algorithm as follows:

function modInv(value, modulus)
{
  V1 ← 1; V3 ← value;
  U1 ← 0; U3 ← modulus;
  while V3 ≠ 0
  {
    QQ ← U3 / V3;
    T1 ← U1 − V1 * QQ;
    T3 ← U3 − V3 * QQ;
    U1 ← V1; U3 ← V3;
    V1 ← T1; V3 ← T3;
  }
  return U1;
}

Then the division is equal to the product a⁢c (mod p).

To minimize the number of modular divisions (which are expensive), we can multiply all coefficients of both polynomials (dividend and divisor) by the multiplicative inverse of the leading coefficient of the divisor polynomial. In this way, we will divide by a monic polynomial, that is, a polynomial whose leading coefficient equals 1. This will not change the quotient, but the remainder will need to be multiplied by the leading coefficient of the divisor polynomial.

Example: perform the division (3⁢x³ + 7⁢x² + 5⁢x + 6) / (4⁢x² + 3⁢x + 10) (mod 11):

First we have to find the multiplicative inverse of 4 (mod 11), which equals 3, because 4 × 3 = 12 ≡ 1 (mod 11). Multiplying all coefficients by 3 we get:

(9⁢x³ + 10⁢x² + 4⁢x + 7) / (x² + 9⁢x + 8) (mod 11)

Dividing the leading coefficient of the dividend polynomial over the leading coefficient of the divisor polynomial: 9⁢x³ / x² ≡ 9⁢x (mod 11).

Now we subtract the product of what we have just found by the divisor polynomial from the dividend polynomial. We obtain:

9⁢x³ + 10⁢x² + 4⁢x + 7 - 9⁢x(x² + 9⁢x + 8) ≡ 6⁢x² + 9⁢x + 7 (mod 11). Notice that 10 − 9 × 9 ≡ 6 (mod 11) and 4 − 9 × 8 (mod 11) ≡ 9 (mod 11).

Dividing the leading coefficient of the remainder polynomial we just found over the leading coefficient of the divisor polynomial: 6⁢x² / x² ≡ 6 (mod 11).

Now we subtract the product of what we have just found by the divisor polynomial from the remainder polynomial. We obtain:

(6⁢x² + 9⁢x + 7) − 6⁢(x² + 9⁢x + 8) ≡ 10x + 3 (mod 11). Notice that 9 − 6 × 9 ≡ 10 (mod 11) and 7 − 6 × 8 ≡ 3 (mod 11).

The quotient polynomial equals 9⁢x + 6 and the remainder polynomial equals 4⁢(10x + 3) ≡ 7 x + 1 (mod 11).

Division of integer polynomials

The division proceeds in the same way as in the previous subsection, with the difference that there is no multiplicative inverse of the leading coefficient of the divisor polynomial. So, for each iteration of the algorithm, a division is required. If the divisor polynomial is not monic, sometimes the division cannot be performed. This occurs when the leading coefficient of the remainder is not multiple of the leading coefficient of the divisor.

The greatest common divisor of two polynomials is the polynomial of the highest possible degree, that divides both polynomials.

For example: gcd(x² + 6⁢x + 5, 2⁢x² + 13⁢x + 15) = x + 5

GCD of polynomials modulo a prime number

The greatest common divisor can be found using the Euclidean algorithm as follows:

function gcdm(poly1, poly2, p)
{
  a ← poly1;
  b ← poly2;
  while b ≠ 0
    t ← b;
    b ← remainder(a, b) (mod p);
    a ← t;
  return a;
}


GCD of integer polynomials

While the same algorithm could be used for finding the gcd of two integer polynomials, the coefficients of the intermediate polynomials increase very quickly, so this is not useful.

There are two methods to find the gcd: the subresultant algorithm and the modular algorithm. The latter, invented by William Brown, uses gcd of polynomials modulo a prime, so this is better for our purposes.

function gcdp(poly1, poly2)
{
  c1 ← cont(poly1);
  c2 ← cont(poly2);
  c ← gcd(c1, c2);
  a ← poly1 / c1;
  b ← poly2 / c2;
  h ← gcd(lc(poly1), lc(poly2));
  d ← min(deg(poly1), deg(poly2));
  m ← 1;
  gm ← 0;
  repeat
  {
    do
    {
      do
      {
        p ← nextPrime();
      } while remainder(m*h, p) = 0;
      r ← poly1 (mod p);
      s ← poly2 (mod p);
      gp ← gcdm(r, s, p);
      if deg(gp) = 0
      {
        output c; stop;
      }
    } while deg(gp) > d;
    gp ← (h mod p)/lc(gp) * gp;
    if deg(gp) < d
    {
      m ← 1;
      gm ← 0;
      d ← deg(gp);
    }
    h ← CRA([gp, p], [gm, m]);
    if h = gm
    {
      h ← h / cont(h);
      if remainder(a, h) = 0 and remainder(b, h) = 0
      {
        output c*h; stop;
      }
    }
    m ← p * m;
    gm ← h;
  }
}

The content of a polynomial is the greatest common divisor of all coefficients of that polynomial. The sign of the content matches the sign of the leading coefficient. It is expressed as cont(f) in the previous algorithm.

The coefficients of the result of the Chinese Remainder Algorithm (function CRA) computed modulo mp, must be in the range −mp/2 to mp/2.

The idea behind this algorithm is to compute several gcd of the input polynomials modulo different primes. We can see that their degrees are greater or equal than the degree of the polynomial gcd.

There are three cases:

• The degree of the modular gcd is greater than the degree of previous modular gcd: the new result is discarded because it has incorrect degree.
• The degree of the modular gcd is less than the degree of previous modular gcd: the old result is discarded because it has incorrect degree. The new result is used instead.
• Both degrees are equal: the coefficients of both gcds are merged using the Chinese Remainder Theorem. This algorithm computes g (mod mp) given g_m (mod m) and g_p (mod p).

The algorithm continues until the computed polynomial divides both input polynomials.

The factoring algorithm can be divided in four steps: square-free factorization, distinct degree factorization, equal degree factorization and factor lift. All steps require monic polynomials, so before starting, we will have to multiply all coefficients by the modular multiplicative inverse of the leading coefficient of the polynomial.

Square-free factorization

The next steps do not work if there are repeated factors, so the first step is to factor the polynomial in such a way that there are no repeated factors.

The derivative of the polynomial f(x) = f₀ + f₁⁢x + f₂⁢x² + ... + f_n⁢xⁿ is f'(x) = f₁ + 2⁢f₂⁢x + ... + n⁢f_n⁢x^{n − 1}

The recursive algorithm is:

function SFF(f)
{
  i ← 1; g ← f';
  if g ≠ 0
  {
    c ← gcd(f, g);
    w ← f / c;
    while w ≠ 1
    {
      y ← gcd(w, c); z ← w / y;
      Output(zi); i ← i + 1;
      w ← y; c ← c / y;
    }
    if c ≠ 1
    {
      c ← c1/p;
      Output(SFF(c)p);
    }
  }
  else
  {
    f ← f1/p;
    Output(SFF(f)p);
  }
}

For each polynomial on the output of this algorithm we have to perform the next steps.

Distinct degree factorization

This is a method that exploits the fact that the product of all monic irreducible polynomials of degrees that divide d (mod p) equals x^e − x where e = p^d.

function DDF(f, p)
{
  i ← 1;
  h ← f;
  j ← x;
  q ← 0;
  while deg(h) ≥ 2⁢i
  {
    j ← jp (mod h);
    g ← gcdm(h, j − x);
    if g ≠ 1
    {
      Output(g, i);
      q ← 1;
      h ← h / g;
    }
  }
  if h ≠ 1
  {
    Output(h, deg(h));
    q ← 1;
  }
  if q = 0
  {
    Output(f, 1);
  }
}

The pairs that form the output of this function are the input for the next step. The first element of the pair is a factor of f which is the product of all factors of the degree specified in the second element of the pair.

Equal degree factorization

This is a probabilistic method due to David Cantor and Hans Zassenhaus that finds all factors of the same degree from the output of the previous algorithm:

function EDF(f, d, p)
{
  n ← deg(f);
  set list of factors to f;
  while size(list of factors) < n/d
  {
    h ← random polynomial with deg(h) < n;
    e ← (qd − 1) / 2;
    g ← he − 1 (mod f);
    for each element u of list of factors
    {
      if deg(u) > d
      {
        j ← gcdm(g, u);
        if j ≠ 1 and j ≠ u
        {
          discard u from list of factors;
          add j and u/j to list of factors;
        }
      }
    }
  }
  Output list of factors
}

Polynomial lift

Once all irreducible factors of the polynomial modulo p are found, we can find the factor of the polynomial modulo pⁿ if there are no repeated factors. The following algorithm can lift from modulo m to m², so we can execute it several times until the complete lifting is done.

On input: f = g*h (mod m), s*g + t*h = 1 (mod m)

On output: f = g'*h' (mod m²), s'*g' + t'*h' = 1 (mod m²) with deg(s') < deg(h'), deg(t') < deg(g')

All calculations below are performed mod m².

  e ← f - g*h
  Compute q, r such that s*e = q*h + r
  g' ← g + t*e + q*g
  h' ← h + r

  e ← s*g' + t*h' − 1
  Compute q, r such that s*e = q*h + r
  s' ← s − r
  t' ← t − t*e − q*g'

The initial values of s and t are obtained from g and h using Extended Euclidean algorithm.

We can use the output of the previous section to factor integer polynomials.

First of all, we need to divide the polynomial by its content (the greatest common divisor of all coefficients with the sign of the leading coefficient). The result is the principal part, named pp(f).

The methods shown below do not work if the factors are repeated. To separate them, we can use Yun's algorithm as follows:

Let f = a₁ ⁢ a₂² ⁢ a₃³ ⁢ ... ⁢ a_nⁿ

a[0] = gcd(f, der(f));
b[1] = f / a[0];
c[1] = der(f) / a[0];
d[1] = c[1] - der(f[1]);
i = 1;
repeat
  a[i] = gcd(b[i], d[i]);
  b[i+1] = b[i] / a[i];
  c[i+1] = d[i] / a[i];
  i = i + 1;
  d[i] = c[i] - der(b[i]);
until b = 1;
output a[1], ..., a[i-1];

At this moment we know that there are no repeated factors. We have to find a small prime p for which the polynomial has no repeated factors. This can be easily found by checking that gcd(f, f') ≡ 1 (mod p).

We need to find a bound of the coefficient of factors. The Knuth-Cohen bound for coefficients of polynomial factors can be computed as follows:

If polynomial G divides F we have for all j:

|G_j| ≤ binomial(n − 1, j)*(Σ_i |F_i|²)^1/2 + binomial(n − 1, j −1) * |F_m|

where m is the degree of F and n is the degree of G. The maximum degree to be considered is n = ceil(m/2).

Once the bound B is found, we have to compute the least value of e such that 2⁢B < p^e.

Now we factor the polynomial mod p^e which has r different factors. If r > 10, we can try a different value of p, so we minimize the number of factors found. The application tries up to 5 different primes.

We will now combine the factors found modulo p^e to get integer polynomial factors. There are 2^r possible combinations, but most of them can be easily discarded.

Let the combination of factors found are f₀, f₁, ..., f_s. Compute a ≡ lc(f) ⁢ tc(f₀)⁢ tc(f₁) ⁢ ... ⁢tc(f_s) (mod p^e) and b ≡ lc(f) ⁢ tc(f) (mod p^e). In both cases the products must be in the range −p^e/2 to p^e/2. If a does not divide b, we can discard that combination.

For the very few combinations that remain, compute a ≡ lc (f) ⁢ f₀ ⁢ f₁ ⁢ ... ⁢ f_s (mod p^e). Again, this time the coefficients of this polynomial must be in the range −p^e/2 to p^e/2. Compute b = gcdp(a, lc(f) ⁢ f). If the degree of b is zero, discard this combination. Otherwise, the polynomial b is a proper factor of f.

Even though the steps are very fast, for some polynomials the number of combinations is very huge. So, the program uses Van Hoeij algorithm which uses LLL algorithm to determine quickly the combinations to use.

The program finds the roots from the irreducible polynomials obtained in the previous section.

For polynomials of degrees up to 4, the standard formulas are used. No cube roots of complex numbers are used. In these cases (named casus irreducibilis), the program uses trigonometric functions.

For polynomials of degree 5, the application uses the results of D. S. Dummit's paper Solving Solvable Quintic, Mathematics of Computation volume 57, number 195, July 1991, pp. 387-401.

The program can determine whether an irreducible polynomial is cyclotomic, i. e., a divisor of xⁿ − 1. In this case the program shows the complex roots of 1:

cos cosine of m times pi πn + i sin times sine of m times pi πn

where m and n are coprime.

If n = 2^q × 3^r × 5^s × 17^t where r < 2, s < 2 and t < 2, the program can find the roots as radical expressions.

For irreducible polynomials of the form a xⁿ + b, the roots can be computed as the product of a real number by a complex root of 1, so the method used in the previous paragraph is used.

For a x²ⁿ + b xⁿ + c there are two cases. If the quadratic polynomial (discarding the variable n) has real roots, again we can use the methods of the previous paragraph. Otherwise, only trigonometric functions are shown.

For almost all polynomials with degree 5 or greater, the roots cannot be expressed as radical expressions. To detect these cases with high probability, the program uses the results of two papers about Galois theory:

• J.H. Davenport, G.C. Smith, Fast recognition of alternating and symmetric Galois groups, Journal of Pure and Applied Algebra 153 (2000) 17-25.
• K. Conrad, Recognizing Galois groups S_n and A_n

You can download the source of the current program and the old sum polynomial factorization applet from GitHub. Notice that the source code is in C language and you need the Emscripten environment in order to generate JavaScript.