-- Example which illustrates the limitations of the classical
-- method of polynomial factorization in Z[x].
-- INTRODUCTION
-- Recall that a polynomial of degree n is completely determined by
-- its values at n+1 distinct points. For simplicity we shall use
-- the points 0,1,2,...,n.
-- If a polynomial G is a factor of F then G(k) divides F(k) for any
-- choice of k; we say that 0 divides 0 in this context.
-- The classical factorization algorithm proceeds as follows.
-- For each degree D starting from 1 up to deg(F)/2 Do
-- For each (D+1)-tuple V[0],...,V[D] satisfying
-- V[0] divides F(0), V[1] divides F(1), etc
-- Interpolate G from the values in V[0],...,V[D].
-- Check whether G divides F.
-- If so, announce the factor G, and replace F by F/G
-- This algorithm is disarmingly simple, and obviously works.
-- The catch lies in the number of (D+1)-tuples which may have to be
-- tried; if F is irreducible, we will try all possible tuples before
-- discovering that it is irreducible.
-- EXERCISE
-- This exercise asks you to estimate (with the help of CoCoA) how many
-- tuples may have to be considered for some quite small polynomials.
-- You are not asked to implement the algorithm.
-----------------------------------------------------------------------------
-- A function which counts how many different factors an integer
-- has; negative factors are counted too.
Define NumFactors(N)
N := Cast(N, INT);
Facs := Comp(SmoothFactor(N, N), 1);
Primes := Set(Facs);
NumFacs := 2;
Foreach P In Primes Do
NumFacs := NumFacs*Count(Facs, P);
EndForeach;
Return NumFacs;
EndDefine; -- NumFactors
PrintLn "4 has ", NumFactors(4); -- there are six factors -4,-2,-1,1,2,4
-----------------------------------------------------------------------------
-- Function to generate random polynomials of degree D.
Define RndPoly(D)
Return Sum([Rand(-999,999)*x^I | I In 0..D]);
EndDefine; -- RndPoly
-- Let F(x) be the product of two cubics with integer coefficients
-- chosen uniformly and randomly in the interval [-999,999].
Use Q[x];
F := RndPoly(3)*RndPoly(3); -- NOT the same as RndPoly(3)^2
-- Now compute the maximum number of tuples we might have to consider.
CountTuples := 0;
For D := 1 To Deg(F)/2 Do
FValues := [Eval(F,[K]) | K In 0..D];
CountTuples := CountTuples + Product([NumFactors(V) | V In FValues]);
EndFor;
PrintLn "Max number of tuples to consider: ", CountTuples
-- The estimate above could be improved. We do not need to consider
-- every tuple since that means we would interpolate each putative
-- factor G twice (once as G,and once as -G). This observation lets
-- us reduce by half the number of tuples to try.
-- Repeat the above count after replacing F by F+120 (which is almost
-- certainly irreducible, so all tuples really would have to be tried).
-- Calculate estimates for the number of tuples to consider for
-- several smallish polynomials.
-- Compare this with the speed of the built-in factorizer of CoCoA:
G := RndPoly(20)*RndPoly(20)*RndPoly(20)*RndPoly(20);
Time NoPrint:=Factor(G); -- should be almost instant
-----------------------------------------------------------------------------
-- Now try factorizing the same polynomial F modulo some small primes;
-- you should avoid any primes which divide the leading coefficient.
-- How often is the modular factorization faithful to the non-modular
-- one? (i.e. the modular factors are two cubics)
Use Z/(101)[x]; -- try factorizing F modulo 101
Factor(BringIn(F)); -- I ought to check that BringIn(F) has full degree.