-- This exercise shows the simple reconstruction method from Z to Z[x].
-- In fact we shall also use it to reconstruct from Z[y] to Z[x,y].

-----------------------------------------------------------------------------
-- We start with the definitions of two functions which effect the
-- reconstruction.  Then we verify correct working on a simple example.

-- This function accepts two integers Xval and Yval.
-- We suppose that Yval = f(Xval) for some univariate
-- polynomial whose coefficients have magnitude at most
-- Xval/2.  The function computes the unique polynomial
-- with coefficients less than Xval/2 whose value at
-- x=Xval is Yval.
Define ReconstructUniv(Xval, Yval)
  -- Student should write this function.
  F := 0;
  Xpower := 1;
  While Yval <> 0 Do
    C := Mod(Yval, Xval);
    Yval := (Yval-C)/Xval;
    If C > Xval/2 Then
      C := C-Xval;
      Yval := Yval+1;
    EndIf;
    F := F+Xpower*C;
    Xpower := Xpower*x;
  EndWhile;
  Return F;
EndDefine; -- ReconstructUniv

-- This function reconstructs a bivariate polynomial from
-- its univariate image under x |-> Xval (where Xval is a
-- positive integer).  The coefficients of the bivariate
-- polynomial have absolute value less than Xval/2.
Define Univ2Biv(F, Xval)
  Ans := 0;
  Foreach CPP In Monomials(F) Do
    Ans := Ans + ReconstructUniv(Xval, LC(CPP))*LPP(CPP);
  EndForeach;
  Return Ans;
EndDefine; -- Univ2Biv

-- Check that ReconstructUniv gives the right answer in
-- an easy case.  The result here should be 3x^3-x+7
ReconstructUniv(100,2999907);


-----------------------------------------------------------------------------
-- Now we apply the reconstruction routines to enable us to factorize
-- bivariate polynomials delegating to CoCoA the task of factorizing
-- univariate polyomials.

-- Create two random dense bivariate polynomials F, G.
-- Put their product in H, the ask CoCoA to factorize H.
Skel := (x^21-1)/(x-1)*(y^21-1)/(y-1);
F := Randomized(Skel);
G := Randomized(Skel);
H := F*G;
Time Facs := Factor(H); -- Note how long this took.

-- Now repeat the factorization but using a different way.
-- For simplicity we assume that no factor of H has a coefficient
-- greater than 0.5*10^25; ideally we should use a proper bound.
Xvalue := 10^25;

-- HUniv is the univariate image of H under x |-> Xvalue (Xvalue=10^25).
-- Factorize HUniv, and reconstruct bivariate factors from
-- the univariate ones using ReconstructUniv.
Time HUniv := Subst(H, x, Xvalue);
Time FacsUniv := Factor(HUniv);
Time Facs2 := [Univ2Biv(FacPow[1], Xvalue) | FacPow In FacsUniv];
-- Verify the answer.
Product(Facs2) = H;

-- Now repeat the above for several different pairs of random
-- bivariate polynomials.  Does this approach work?
-- How does it compare with CoCoA's built-in algorithm for speed?
-- How does varying the degree in x affect the relative speeds?
-- How does varying the degree in y affect the relative speeds?
-- What happens if the factors are quite sparse? (i.e. most of the
-- coeffs are zero)

-----------------------------------------------------------------------------
-- Now try starting from these two bivariate polynomials.
F := 2*y+x^2+x+2;
G := F+2;

-- The approach in this exercise gives an obviously wrong answer: why?
-- How can the method be modified so that it always gives the right result?