-- 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.