-- This example studies whether Chinese remaindering is likely to be
-- beneficial for computing the composition of two polynomials in Z[x].
-- This is a pilot study -- I do not know what result it will give.
-- You will be asked to measure the times for a direct approach and a
-- modular approach, and then to decide which seems better.

-- The aim is to compute f(g(x)) given f,g in Z[x].
-- We shall restrict (initially) to the case of f and g being dense,
-- i.e. most or all of their coefficients are non-zero.  We shall
-- ignore the possibility of being able to do "something clever" if
-- f and/or g has only very few non-zero coefficients (commpared to
-- the degree).

-- We shall use Horner's algorithm for composing the polynomials.
Define HornerCompose(F, G)
  X := Indet(UnivariateIndetIndex(F));
  C := Coefficients(F, X);
  Ans := C[1];
  For I := 2 To Len(C) Do
    Ans := Ans*G + C[I];
  EndFor;
  Return Ans;
EndDefine; -- HornerCompose

-- The test cases should be (initially) some dense random polynomials.
-- Here is a function which generates some random polynomials.
Define RndPoly(D)
  Return Sum([Rand(-999,999)*x^I | I In 0..D]);
EndDefine; -- RndPoly

-- You need to choose the sizes (i.e. degrees) of the polynomials
-- you will use for experimentation.  I suggest choosing degrees
-- which lead to a time of about 10 seconds for the direct computation.
-- Try varying the values of DegF and DegG, but do not increase them
-- too drastically between trials.
DegF := 10;
DegG := 10;
Time FoG := HornerCompose(F, G);

-- Can you determine an empirical formula for the time in terms of
-- DegF and DegG?  Can you improve the formula to account for the
-- coefficient sizes too?


-- Before experimenting with a modular approach, find a bound for
-- the size of the coefficients of F(G(x)) in terms of the coefficients
-- of F and of G.  How does your bound compare with the values actually
-- produced (in some random cases)?
-- Since Z is an integral domain, we know the degree of F(G(x)).
-- Degree and coefficient bound measure the size of the result.

-----------------------------------------------------------------------------
-- Below is some CoCoA code which will estimate how much time would be
-- needed by a modular approach (based on Chinese remaindering).
-- For the moment we shall ignore the reconstruction cost; if the total
-- time for computing the compositions modulo various primes is not much
-- smaller than the time for computing it directly then a CRT modular
-- approach is not going to be useful.
-- [my suspicion is that it will not be useful]

Bound := ???; <--- insert your bound here
M := 1;
P := NextPrime(10000);
-- The loop below tries successive primes until their product exceeds Bound
While M < Bound Do
  M := M*P;
  ModularRing ::= Z/(P)[x];
  Using ModularRing Do
    ModularF := BringIn(F);
    ModularG := BringIn(G);
    Time ModularFoG := HornerCompose(ModularF, ModularG);
  EndUsing;
  P := NextPrime(P);
EndWhile;


-----------------------------------------------------------------------------
-- Harder questions:
-- * Is a CRT approach useful sometimes?  What if DegF is large and DegG is
--   small?  Or vice versa?  What if the coeffs are big (or small)?
-- * Can Hensel lifting be used for computing the composition? [I guess not]
-- * It would be very easy to compute values of F(G(x)) for various values
--   of x.  Is a modular approach based on univariate polynomial
--   reconstruction beneficial?  The answer probably depends criticially
--   on the cost of reconstruction.