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