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