-- Example of Ben-Or and Tiwari's sparse interpolation algorithm.
-- You might like to look at the motivating example in MultivariateHensel.coc
-- before going through this exercise.

-----------------------------------------------------------------------------
-- I suggest you run through the example below fairly quickly, and then
-- consider some of the variations and questions at the end.  Keep your
-- notes about BenOr-Tiwari's method handy, to help you understand what
-- follows.

-- Here is the polynomial we shall reconstruct; of course, we must
-- pretend we do not know its representation.
Use Q[x,y,z];
F := 1+2*x^10+3*y^100+4*z^1000;

-- T is the number of terms in F; we suppose we know the exact value.
T := Len(F);

-----------------------------------------------------------------------------
-- Now we start the algorithm.
-- Evaluate F at 2*T+1 specially chosen points -- we use F only here.
V := [Eval(F,[2^I,3^I,5^I]) | I In 0..(2*T)];

-- Use LinKer to find the coefficients of the recurrence relation.
M1 := Mat([[V[I+J-1] | I In 1..(T+1)] | J In 1..T]);
K := LinKer(M1);
If Len(K) <> 1 Then Error("Estimated number of terms is wrong."); EndIf;

-- The coefficients of the recurrence relation are in K;
-- build an auxiliary polynomial with these coefficients and factorize it.
AuxPoly := Sum([K[1][I]*x^(I-1) | I In 1..Len(K[1])]);
FacPows := Factor(AuxPoly);

-- The values of the power products of f at the point (2,3,5) are
-- the roots of the (linear) factors...
PPImages := [-LC(Subst(FacPow[1],x,0)) | FacPow In FacPows And Deg(FacPow[1]) = 1];

-- Determine the coefficients of the polynomial by solving a linear system.
M2 := Mat([[PPImages[I]^(J-1) | I In 1..T] | J In 1..T]);
RHS := [V[I] | I In 1..T];
Coeffs := LinSol(M2, RHS);

-- We have the complete answer now, so print it out.
PrintLn "The polynomial is the sum of:";
For I := 1 To T Do
  Xexp := FactorMultiplicity(Num(PPImages[I]), 2);
  Yexp := FactorMultiplicity(Num(PPImages[I]), 3);
  Zexp := FactorMultiplicity(Num(PPImages[I]), 5);
  PrintLn Coeffs[I], " * ", "x^",Xexp, " * y^",Yexp, " * z^",Zexp;
EndFor;

-----------------------------------------------------------------------------
-- Notice how the method is (largely) unaffected by the degree of the
-- polynomial to be reconstructed.  A direct application of multivariate
-- Hensel lifting would be utterly impossible.  Even Zippel's method
-- would have to work hard.

-- Recently W-S Lee has combined Zippel's method with BenOr-Tiwari's.
-- The idea is to use BenOr-Tiwari for sparse univariate reconstruction
-- inside Zippel's method.  The combined method works very well for
-- extremely sparse polynomials of high degree.  Lee has also studied
-- applying sparse multivariate reconstruction to the case where the
-- values of the polynomial can be found only approximately (e.g. floating
-- point numbers).


-- The weak point of BenOr-Tiwari's method is the size of the numbers it
-- generates (primarily because the later evaluation points have large
-- coordinates).  Have a look at the size of the later values in the
-- list V -- you might find the function FloatStr helpful here.
-- Have a look at the coefficients in AuxPoly, the one whose roots give
-- us the power products in the polynomial being reconstructed.
-- Finally have a look at the entries in the matrix M2 from the linear
-- system we solved to find the coefficients of our polynomial.

-- Try running the method to reconstruct some other sparse multivariate
-- polynomials (but limit yourself to a maximum of 10 terms, otherwise
-- the computations may become lengthy).  You will need to modify the
-- code if you want to change the number of indeterminates.


-- Some questions.
-- (1) It would be nice if BenOr-Tiwari could use smaller evaluation
--     points.  Find a family of univariate binomials (i.e. just 2 terms)
--     which shows that there is no (sparse) reconstruction algorithm
--     which uses only the five evaluation points -2,-1,0,1,2.
-- (2) Let zeta be an r-th root of unity (i.e. zeta^r-1 = 0).
--     Find a family of univariate binomials (i.e. just 2 terms)
--     which shows that there is no (sparse) reconstruction algorithm
--     which uses as its five evaluation points zeta^e1, zeta^e2,... zeta^e5.
--     (where e1,...,e5 are integers).

-- HARD QUESTIONS (I don't know how to answer them).
-- The bane of BenOr-Tiwari's method is the size of numbers it uses;
-- this limits drastically its practical usefulness.  Can you modify the
-- method so that it uses smaller evaluation points?  Maybe there is a
-- way for sparse univariate polynomials?  What about the bivariate case?
-- What about binomials?

-- BenOr-Tiwari's method obviously exhibits "intermediate expression swell"
-- so appears to be an ideal candidate for being "modularized"... except
-- for one aspect: FactorMultiplicity will not work in a modular context.
-- W-S Lee has found an interesting way of adapting BenOr-Tiwari to work
-- for polynomials over Z/(P).  Can you find a way?
-- NB we must assume that the degree of the polynomial is less than P
-- because x^P-x evaluates to zero at every point in Z/(P).