-- Exercise about Euclid's Algorithm; it illustrates why a direct approach
-- does not work well.

-- Function which computes the GCD of univariate polynomials F and G
-- using Euclid's algorithm
Define EA(F, G)
  -- Note that NR(F,[G]) computes the remainder of F modulo G
  While G <> 0 Do
    Swap := NR(F,[G]);
    F := G;
    G := Swap;
  EndWhile;
  Return F;
EndDefine; -- EA


-- Compare the direct implementation of Euclid's algorithm
-- with CoCoA's built-in GCD.  We compute the GCD of two
-- random polynomials of degree D.
D := 100;
T := Randomized((x+1)^D); -- random poly of degree D
F := T*Randomized((x+1)^D); -- random poly of degree D
G := T*Randomized((x+1)^D); -- random poly of degree D

-- First we compute with the built-in function.
Time H := GCD(F,G);     -- should be almost instant
PrintLn "Result is ",H; -- probably just 1


-- Now try the same computation using a naive implementation.
Time K := EA(F,G);
-- It probably took a lot longer than the built-in GCD function
-- Result is essentially a rational number
PrintLn "Numerator of result is about ",FloatStr(Num(LC(K)));
PrintLn "Denominator of result is about ",FloatStr(Den(LC(K)));

-- Before reading on, think about why EA is slow.
-- Can you think of a way of making it a bit faster (without
-- turning to modular methods)?

-----------------------------------------------------------------------------
--  Here we have a smarter implementation of Euclid's Algorithm.
--  It assumes that we start with polynomials with integer coeffs,
--  and is written to avoid completely all rational numbers.

-- IRem computes the remainder of N*F modulo G where N is an
-- integer chosen to ensure that the result has integer coeffs
-- (N may not be the smallest possible choice).
Define IRem(F, G)
  X := Indet(UnivariateIndetIndex(F));
  While F <> 0 And Deg(F) >= Deg(G) Do
    H := GCD(Num(LC(F)), Num(LC(G)));
    F := (LC(G)/H)*F - (LC(F)/H)*G*X^(Deg(F)-Deg(G));
  EndWhile;
  Return F;
EndDefine; -- IRem

Define Content(F)
  Return GCD([Cast(N, INT) | N In Coefficients(F)]);
End;

-- Same as EA except we use IRem to compute the remainder.
Define IEA(F, G)
  While G <> 0 Do
    Swap := IRem(F,G);
    If Swap <> 0 Then 
      Swap := Swap/Content(Swap);
    EndIf;
    F := G;
    G := Swap;      
  EndWhile;
  Return F;
EndDefine; -- EA


-- Let us see how the smart implementation performs.
Time L := IEA(F,G); -- should be fairly quick
PrintLn "Result is a large integer; value is about ", FloatStr(L);

-- Can you think why IEA is so much faster than EA?
-- Bear in mind that the integer coeffs appearing during the
-- computation in IEA are about as big as the rationals appearing
-- during the run of EA; here the size of an integer is it number
-- of digits, and the size of a rational is the sum of the sizes
-- of its numerator and its denominator.


-----------------------------------------------------------------------------
-- How does the built-in GCD function manage to be so fast?
-- In part it is a trick: if the GCD is 1 then the built-in
-- function is especially fast.  Try the following.

P := NextPrime(10000);
ModularRing ::= Z/(P)[x,y,z];
Using ModularRing Do
  ModularF := BringIn(F);
  ModularG := BringIn(G);
  ModularH := GCD(ModularF, ModularG); -- This should be fast because
                                       -- coeffs cannot become enormous.
  PrintLn "GCD modulo ",P," is ",ModularH;
EndUsing;