{VERSION 3 0 "SGI MIPS UNIX" "3.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }} {SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT -1 527 "Greatest Common Divisor. \+ GCD. In Maple: gcd for polynomials, and igcd for integers.\nIf a and b are two nonnegative integers then the gcd g is such that the ideal (a ,b) = Z*a+Z*b is generated by g, so (g)=(a,b). Since the ideal (a,b) e quals the ideal (a-b,b) and also (a,b-a) we see that by repeated subst ractions we can make the numbers a,b smaller until one of them vanishe s. Then we've found generators a', b' of this ideal, with one of these two equal to zero. So we've found just 1 generator of the ideal, whic h is the gcd." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 170 "GCD := pr oc(a::nonnegint, b::nonnegint)\n if a=0 then\n b\n elif b=0 \+ then\n a\n elif a>b then\n procname(a-b,b)\n else\n \+ procname(a,b-a)\n fi\nend;\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "GCD( 12, 18);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 120 "Now the same algorithm, but a little bit faster, based on the fac t that the ideal (a,b) equals the ideal (b, irem(a,b))." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 434 "GCD := proc(a::nonnegint, b::nonne gint)\n if b=0 then\n a\n else\n # Now b<>0 so I can rep lace a by irem(a,b)\n # without changing the ideal (a,b) = Z*a+Z* b.\n # I'll also interchange a and b, this way from\n # this point on, b will always be smaller than a,\n # which causes that the irem in the next step will\n # always make a smaller, so the algorithm will end.\n procname(b,irem(a,b))\n fi\nend;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "GCD(100!, 12345);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 156 "Now lets do the same for polynomi als. If f,g in Q[x] then we will compute the gcd d. So (f,g) = (d) whe re (f1..fn) stands for the ideal Q[x]*f1+...+Q[x]*fn." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "GCD:=proc(f,g,x::name)\n if g=0 t hen\n f\n else\n procname(g, rem(args), x)\n fi\nend;" } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "f:=expand( mul(x-i,i=-2..6));" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 30 "g:=expand( mul(x-i,i=4..10) );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "GCD(f,g,x);" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 13 "content(%,x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 95 "The content C of a polynomial F in Q[x] is the largest ra tional number C such that F/C in Z[x]." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "%% / %;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 207 "The p olynomial F/C is called the primitive part of F. It is called primpart in Maple. We can use primpart to get rid of the fractions, and we can also use it to remove common integer factors of a polynomial." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "GCD(f,g,x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "primpart(%);" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 27 "primpart( 15*x+27*x^2 , x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "f:=randpoly(x,degree=3);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "f,g := seq(expand(f*randpoly(x,degr ee=15)) , i=1..2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "GCD(f ,g,x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "primpart(%);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 314 "As you can see, the output that o ur GCD produces is not very nice, although in principle it is correct \+ because it is a generator of the Q[x]-ideal (f,g) and so we can call i t the gcd of f and g. Note that this gcd is only determined up to a no n-zero factor in Q. We could get rid of these large numbers as follows :" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 108 "GCD:=proc(f,g,x::name )\n if g=0 then\n primpart(f,x)\n else\n procname(g, rem (args), x)\n fi\nend;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 " GCD(f,g,x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 217 "Now it looks nice r, but it's still not very good. This algorithm will be very slow on l arger inputs because it internally computes with very large rational n umbers. Only and the end they get reduced to small integers." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "f:=randpoly(x,degree=3);" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "f,g := seq(expand(f*randpol y(x,degree=40)) , i=1..2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "GCD(f,g,x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "time();" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 106 "That already took quite a whil e. For higher degrees it will get much worse. Lets see why it takes so long:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 136 "GCD:=proc(f,g,x: :name)\n lprint( length( [args] ));\n if g=0 then\n primpart( f)\n else\n procname(g, rem(args), x)\n fi\nend;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "GCD(f,g,x);" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 18 "length( [f,g,x] );" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 154 "So the sum of the number of digits of the two polynomial s went up to about 70000, even though in the input its only 400 someth ing, and in the output it's:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "length(%%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 236 "So our GCD \+ algorithm is very slow, in fact, it's exponentially slow. Try to run s ome examples for degree 50, 100, 150, 200,... and see what happens to \+ the number of digits. Lets also just look at the largest rational numb er that occurs." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 207 "GCD:=pr oc(f,g,x::name)\n local i;\n lprint( [length( [args] ),\n max( seq(length(i),i=[coeffs(f,x),coeffs(g,x)]))] );\n if g=0 then\n \+ primpart(f)\n else\n procname(g, rem(args), x)\n fi\nend; " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "GCD(f,g,x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 349 "We shouldn't use \"options trace\" in th is procedure, the numbers are so large that when you'd print them on t he screen it would be a big mess. Some of those coefficients had 5000 \+ digits (in numerator and denominator combined). Some reduction during \+ the computation will help keeping coefficients smaller, and so that sh ould make the algorithm faster." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 219 "GCD:=proc(f,g,x::name)\n local i;\n lprint( [length( [arg s] ),\n max(seq(length(i),i=[coeffs(f,x),coeffs(g,x)]))] );\n \+ if g=0 then\n primpart(f)\n else\n procname(g, primpart(re m(args),x), x)\n fi\nend;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "GCD(f,g,x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 286 "Now we're loo king at close to 150 digits, which is better, but still not very good \+ when you consider that f, g and the gcd all had small coefficients. Ho wever, just adding this primpart did already speed up the code. The fo llowing example would take *much* longer without that primpart." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "f:=randpoly(x,degree=4);" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "f,g := seq(expand(f*randpol y(x,degree=50,terms=10)) , i=1..2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "GCD(f,g,x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 263 " Now lets see how long it takes in Maple. Note that in maple's gcd comm and you should *NOT* give the third variable x, that leads to strange \+ errors. Maple's gcd is in fact for the ring Q[several variables] and i t will figure out by itself what those variables are." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "gcd(f,g);" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 149 "So this goes pretty much instantly. Why is this so muc h faster? Optimized and compiled code or something like that? If you t hink that, then try this:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "f:=randpoly(x,degree=250,dense):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "f,g := seq(expand(f*randpoly(x,degree=250)) , i=1..2) :" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 482 "If we applied our algorithm GCD on these f and g, w hich have degree 500, then we would make the computer work with intege rs that have many million digits. Even written in machine language or \+ C and highly optimized, that would still be an extremely lenghty compu tation. And besides, if you have say 128 Megabyte, you can not store i ntegers of more than 128*10^6*log(256)/log(10) digits so even if the C PU was infinitely fast the computation would still fail, you'd run out of memory." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "d:=gcd(f,g): " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 288 "As we see, the computation s till takes less than a second. If we had computed over the ring Z[x], \+ which our algorithm GCD does, it would have been impossible to do this , to handle the many-millions of digits in such a short time. So we mu st conclude that Maple uses a different algorithm." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 110 "GCD_mod_p:=proc(f,g,x,p)\n if g=0 then \n f\n else\n procname(g, Rem(f,g,x) mod p, x, p)\n fi\n end;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "d13:=GCD_mod_p(f,g, x,13);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 575 "As we see, this algorithm GCD_mod_p is *much* \+ faster than GCD, while the only difference is that we reduce coefficie nts modulo the prime number p. This makes a big difference because by \+ reducing modulo p our coefficients can never grow large, and whether w e use a simple multiplication algorithm or Karatsuba or even more adva nced algorithms for integer arithmetic, it's obvious that computing wi th just a few digits is always much faster than with very large intege rs. We will see later in this course how gcd's can be computed by comp uting a few gcd's modulo prime numbers." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 689 "Let's first look a little closer at this remainder procedure. How to compute the quotient and remainder o f polynomials. Basically, the idea is exactly the same as the integer \+ division. When you do an integer division of positive integers a and b , then as long as you can substract a multiple of b from a you do so. \+ Then what you end up with is the remainder of a modulo b. And, if you \+ kept track of which multiple of b you substracted from a, then you hav e the quotient too. We will see here that for polynomials it works exa ctly the same. If F,G are polynomials and I want to divide F by G, the n as long as the degree of F is >= the degree of G, I can substract so me multiple of G from F." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "F:=3*x^5+2*x^2+1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "G:=2* x^2+x+1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "a := lcoeff(F,x )/lcoeff(G,x) * x^(degree(F,x)-degree(G,x));" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 133 "Now we see that we can substract a multiple of G from \+ F in such a way that the degree of G will get smaller, namely we subst ract a*G." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "F:=expand(F-a* G);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "a := lcoeff(F,x)/lco eff(G,x) * x^(degree(F,x)-degree(G,x));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "F:=expand(F-a*G);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "a := lcoeff(F,x)/lcoeff(G,x) * x^(degree(F,x)-degree( G,x));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "F:=expand(F-a*G); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "a := lcoeff(F,x)/lcoeff (G,x) * x^(degree(F,x)-degree(G,x));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "F:=expand(F-a*G);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 254 "Now the degree of F is smaller than the degree of G. So we can no further reduce F modulo G, we have found the remainder of the divi sion of F by G. Of course, the quotient is just the sum of all the a's we computed. This leads to the following algorithm." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 265 "quo_rem:=proc(f,g,x)\n local F,q,a;\n \+ F:=f;\n q:=0;\n while degree(F,x) >= degree(g,x) do\n a := l coeff(F,x)/lcoeff(g,x) * x^(degree(F,x)-degree(g,x));\n q := q+a; \n F:=expand(F-a*g)\n od;\n # Now F is the remainder and q is \+ the quotient.\n q,F\nend;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "f:=x^12+x^4+1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "g:=2* x^3+3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "q, r :=quo_rem(f, g,x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "q*g+r;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "expand(%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 51 "For the algorithm to be correct we need two things:" } }{PARA 0 "" 0 "" {TEXT -1 145 "1) degree(r,x) < degree(g,x).\n2) f = \+ q*g + r.\nBoth hold, so the algorithm is correct. Note that q, r are \+ uniquely defined by conditions 1), 2)." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "f:=rand poly(x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "c:=12;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "g:=x-c;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "subs(x=c,f);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "rem(f,g,x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "subs(x=c, f^3+x^3*f+3 );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "rem( f^3+x^3*f+3 ,g,x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 273 "We see that taking the remainder modulo x-c is the same \+ thing as substituting x=c.\nNow how about taking the remainder modulo \+ a polynomial of degree > 1?\nLet's take the polynomial g:=x^2+1\nIn Ma ple the complex number i is denoted by I. So I is one of the two roots of x^2+1." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2 "f;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "subs(x=I,f);" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 9 "g:=x^2+1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "rem(f,g,x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 78 "N ow another example. Lets take the number 2^(1/3). This is a root of g: =x^3-2." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "f:=randpoly(x,de gree=10);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "c:=2^(1/3);" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "g:=x^3-2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "subs(x=c,f);" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 11 "rem(f,g,x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "subs(x=c,%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 110 "The conc lusion is the following.\nIf f and g are polynomials, and r:=rem(f,g,x ), and if c is a root of g, then:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 29 "subs(x=c,f) = subs(x=c, r)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 106 "The explanation \+ is easy,\n f = q*g+r, subs(x=c,f) = subs(x=c,q)*subs(x=c,g) + subs( x=c,r) = subs(x=c,r)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 358 "So we see that computing in the ring Q[ 2^(1/3) ], th e set of all sums of rational multiples of powers of 2^(1/3), is compl etely equivalent to computing in the ring Q[x] modulo the polynomial x ^3-2. Computing modulo x^3-2 simply means that whenever you see a poly nomial with degree >= degree(x^3-2,x) then you just replace it with th e remainder modulo x^3-2." }}{PARA 0 "" 0 "" {TEXT -1 200 "In fact thi s proves that the ring Q[2^(1/3)] (the rational numbers extended with \+ the number 2^(1/3))\nis isomorphic to the ring Q[x]/(x^3-2) (polynomia ls in the variable x modulo the polynomial x^3-2)." }}{PARA 0 "" 0 "" {TEXT -1 88 "Since x^3-2 is irreducible we know from algebra that thes e two rings are in fact fields." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 668 "Let m be a positive integer. Previously \+ we have already in the ring Z/(m), the integers modulo the number m. C omputing modulo m meant that if a number was negative, we add a multip le of m to make it non-negative, and if the number is >=m we substract multiples of m to make it smaller than m. This is done by the command mod. If a is an integer then \"a mod m\" will produce the remainder o f a modulo m, and the result is in the range 0..m-1. So every element \+ of Z/(m) can be represented uniquely by a number in this range. So the ring Z/(m) contains m elements, because in Z/(m) we will only conside r integers to be different when their remainder modulo m is different. " }}{PARA 0 "" 0 "" {TEXT -1 392 "From algebra we know that Z/(m) is a field (meaning that we can divide by non-zero elements) if and only i f m is irreducible, so m has to be a prime number p. When m is a prime number p we will denote this field by Fp. So F5 contains five element s that we can represent by 0,1,2,3,4. The product 2*3 for example will be 1 because 2*3 reduces to 1 modulo 5. So 1/2 will be 3 and 1/3 will be 2." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}}}{MARK "99" 0 }{VIEWOPTS 1 1 0 1 1 1803 }