{VERSION 6 0 "IBM INTEL LINUX" "6.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 256 "" 1 24 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 } {PSTYLE "Maple Output" 0 11 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 11 12 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }} {SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT -1 46 "The first procedure in thi s worksheet is the " }{TEXT 257 28 "Extended Euclidean Algorithm" } {TEXT -1 18 " for integers a,b." }}{PARA 0 "" 0 "" {TEXT -1 101 "It co mputes the gcd d. It also computes integers s,t such that d = s*a + t *b. The output is [d,s,t]." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 19 "The procedure uses " }{TEXT 258 18 "integer div ision. " }{TEXT -1 241 "Integer division is the division you learned i n elementary school. Given two integers a,b, if b is not 0 then you c an divide a by b. The result is then the quotient q, and the remainder r, and you have the following relation: a = q*b + r." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 190 "The command in Mapl e for integer division is: iquo, or irem. You can use either one. I f you use irem, then by giving an optional argument 'q' it will also t ell you what the quotient q is." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 465 "my_igcd ex := proc(a,b)\nlocal r,v,q;\n# I assume that a, b are both >= 0\nif \+ a<0 or b<0 then\n error \"I did not implement this case\"\nfi;\nif b =0 then\n [a, 1, 0] # The gcd is d=a. And d = 1*a + 0*b.\nelse\n \+ r := irem(a,b,'q'); # Now a = q*b + r, so r = a-q*b\n v := m y_igcdex(b,r);\n # Now d = v[1] = v[2]*b + v[3]*r\n # So d = \+ v[2]*b + v[3]*(a-q*b)\n # So d = v[3]*a + (v[2]-q*v[3])*b and ou r output is:\n [ v[1], v[3], v[2]-q*v[3]]\n fi\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 74 "The following procedure computes the mult iplicative inverse of a modulo p." }}{PARA 0 "" 0 "" {TEXT -1 58 "In o ther words, it computes s such that s*a is 1 modulo p." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 667 "inverse_mod_p := proc(a, p)\n loc al v;\n # Compute 1/a mod p;\n # Note: In Maple you may type: 1/a m od p\n # and that works. The only purpose of this procedure\n # is t o show you how Maple can compute 1/a mod p.\n #\n # So all you need \+ to do with this procedure is understand\n # how it works. Once you un derstand how 1/a mod p\n # is calculated, then after that it's OK to \+ use it.\n v := my_igcdex(a,p);\n # Now d = igcd(a,p) = 1\n # And s* a + t*p = d.\n # Hence, modulo p we have s*a = d = 1, so s is 1/a mod p.\n # List v contains [d, s, t] so we return v[2] = s = 1/a mod p. \n if v[1]<>1 then\n error \"1 divided by\",a,\"does not exist mo d\",p\n fi;\n v[2];\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 113 "T he following procedure makes a polynomial monic mod p, by multiplying \+ by the inverse of the leading coefficient." }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 260 "make_monic := proc(f,x,p)\n # Make f mod p monic. We could do that like this:\n # f/lcoeff(f,x) mod p;\n # But I wa nted to show how this division 1/lcoeff(f,x) mod p\n # is done. So th at's why I do it like this:\n f * inverse_mod_p(lcoeff(f,x),p) mod p; \nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 129 "The following procedur e divides a by b, where a,b in Q[x].\nThe output is a list [q, r] cont aining the quotient and the remainder." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 937 "my_quo_rem := proc(a,b,x)\n # options trace;\n l ocal q,r,l;\n if b=0 then\n # we don't allow this, so we stop \+ the computation\n # with an error message:\n error \"secon d argument must be nonzero\"\n fi;\n # We want to find q,r such th at: \n # a = q*b + r and want deg(r)= degree(b,x) do\n l := lcoeff(r,x)/lcoeff(b,x) * x^(degre e(r,x)-degree(b,x));\n r := expand(r - l*b);\n q := q + l; \n # It is easy to see that a=q*b+r is still true,\n # and that deg(r) becaome smaller because we cancelled\n # its leadin g coefficient lcoeff(r,x).\n od;\n # Now deg(r) " 0 "" {MPLTEXT 1 0 195 "my_divide := proc(f,g,x)\n local v,r;\n v \+ := my_quo_rem(f,g,x); # v = [q,r]\n r := v[2];\n # Note: The follow ing equals \"true\" if r=0,\n # and equals \"false\" if r is not 0.\n evalb( r = 0 )\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 91 "The fol lowing procedure divides the polynomial a by B, but this is done modul o the prime p." }}{PARA 0 "" 0 "" {TEXT -1 78 "The output is just the \+ remainder (the quotient is not returned in the output)." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 387 "my_rem_modp := proc(a,B,x,p)\n l ocal b,r,l;\n b := B mod p;\n if b=0 then\n error \"second \+ argument must be nonzero\"\n elif lcoeff(b,x) <> 1 then\n err or \"I only programmed the case where b is monic\"\n fi;\n r := a \+ mod p;\n while degree(r,x) >= degree(b,x) do\n l := lcoeff(r,x ) * x^(degree(r,x) - degree(b,x));\n r := Expand(r - l*b) mod p; \n od;\n r\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 93 "The foll owing procedure computes the monic gcd of two polynomials a,b. This is done modulo p." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 983 "my_gcd_ modp := proc(a,b,x,p)\n local bb;\n if b=0 then\n # Make the \+ gcd monic:\n make_monic(a,x,p)\n else\n # Make b monic, th at way my_rem_modp will\n # run faster because it won't need to d o divisions\n # modulo p.\n bb := make_monic(b,x,p);\n \+ # We now replace a by its remainder modulo bb.\n # This way degre es go down (at least if deg(a)>=deg(b)),\n # and we can apply ind uction, which in programming is called\n # recursion:\n my_g cd_modp(bb, my_rem_modp(a,bb,x,p), x, p)\n # Note that I not only replaced a by its remainder and\n # made b monic, I interchanged a,b as well. This way we are\n # sure that next time, the first \+ argument will have higher\n # degree than the second argument, so we know for sure\n # that degrees will go down when we do my_rem _modp, and so\n # the recursion must work because the degrees go \+ down (so\n # that eventually we must hit b=0 and have found the g cd a).\n fi\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 209 "The follo wing is the modular gcd algorithm. It computes the gcd of two polynomi als F1,F2 in Q[x]. It does this by reducing modulo several primes, and combining these primes using the Chinese remainder theorem." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 1207 "my_gcd := proc(F1,F2,x)\nlocal f1,f2,l, p,gp,g,m,g_old;\n f1 := primpart(F1);\n f2 := primpart(F2);\n l := \+ my_igcdex( lcoeff(f1,x), lcoeff(f2,x) )[1];\n # Starting with a large r prime p is better than\n # starting with p=2 because m = product of primes\n # will grow faster this way.\n p := 40009;\n do\n whi le irem(l,p)=0 do\n p := nextprime(p)\n od;\n \+ gp := my_gcd_modp(f1,f2,x,p);\n if not assigned(g) or degree(gp,x )degree(g,x) then\n g_ old := g;\n # Maple's chrem command:\n # g := chrem([g ,l * gp],[m,p]);\n # Lets see if it also works if I use my own code:\n g := my_chrem(g, l*gp, m, p);\n m := m*p;\n \+ g := mods(g,m);\n # Only do the \"divide\"-test if g=g_ old, in other\n # words if g did not change. This way we will \+ only\n # use the \"divide\" algorithm if it is very likely\n \+ # that g is correct, and we may save some CPU time.\n i f g = g_old and (my_divide(f1,g,x) and my_divide(f2,g,x)) then\n \+ return primpart(g)\n fi;\n fi;\n p := nextpr ime(p)\n od\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 59 "Here is an \+ implementation of the Chinese Remainder Theorem:" }}{PARA 0 "" 0 "" {TEXT -1 99 "Note: if you want to see Maple's implementation then type : interface(verboseproc=2); print(chrem);" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 727 "my_chrem := proc(f1,f2,m1,m2)\n local e1, e2;\n \011# Compute idempotents.\n\011# An idempotent is a number e1 that is \n\011# 1 modulo m1 and 0 modulo m2.\n\011# The other idempotent e1 \+ is\n\011# 0 modulo m1 and 1 modulo m2.\n\011#\n\011# I could do this :\n\011# e1 := (1/m1 mod m2) * m2;\n\011#\n\011# But to show that we can also do\n\011# that computation by ourselves I'll do:\n\011#\n \011e1 := inverse_mod_p(m2,m1) * m2;\n\011#\n\011# Now e1 is (1/m2 mod ulo m1) times m2\n\011# which is 1 modulo m1. And because we\n\011# mu ltiplied by m2 it is 0 modulo m2.\n\011#\n\011# We calculate the other idempotent\n\011# as follows:\n\011e2 := 1-e1 ;\n\011#\n\011# Once we have the idempotents then it\n\011# becomes easy:\n\011#\n\011e1*f1 + e2*f2\n\011#\n\011# modulo m1 this is: 1*f1+0*f2 = f1\n\011# modulo m 2 this is: 0*f1+1*f2 = f2\nend;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%) my_chremGf*6&%#f1G%#f2G%#m1G%#m2G6$%#e1G%#e2G6\"F.C%>8$*&-%.inverse_mo d_pG6$9'9&\"\"\"F6F8>8%,&F8F8F1!\"\",&*&F1F89$F8F8*&F:F89%F8F8F.F.F." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 35 "Now lets test my_gcd on an exam ple:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 359 "F:=x^28+14*x^27+37 *x^26-338*x^25-1865*x^24+2710*x^23+32707*x^22+4312*x^21-338333*x^20-27 3080*x^19+2413597*x^18+2718592*x^17-12839085*x^16-15594106*x^15+531856 11*x^14+59131910*x^13-174479807*x^12-149744386*x^11+448712941*x^10+230 741076*x^9-869923283*x^8-123733756*x^7+1178428779*x^6-259180948*x^5-95 8551300*x^4+558328000*x^3+290160000*x^2-345600000*x+86400000:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "my_gcd(F, diff(F,x), x);" }} {PARA 12 "" 1 "" {XPPMATH 20 "6#,L*&\"(+g&>\"\"\"%\"xGF&!\"\"\"'++sF&* &\"'o*)GF&)F'\"#5F&F&*&\"&1V)F&)F'\"#6F&F&*&\"&s'fF&)F'\"#7F&F(*&\"&1v \"F&)F'\"#8F&F(*&\"%swF&)F'\"#9F&F&*&\"%7EF&)F'\"#:F&F&*&\"$:&F&)F'\"# ;F&F(*&\"$S#F&)F'\"#F&F&*$)F'\"# ?F&F&*&\"'+eSF&)F'\"\"#F&F&*&\"(5E0$F&)F'\"\"$F&F&*&\"(LL)=F&)F'\"\"%F &F(*&\"(]7=#F&)F'FKF&F(*&\"(j$*y\"F&)F'\"\"'F&F&*&\"'wj(*F&)F'\"\"(F&F &*&\"'*G3*F&)F'\"\")F&F(*&\"'=4KF&)F'\"\"*F&F(" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 130 "And lets compare this to Maple's gcd. Note: In Mapl e do *NOT* give the variable x in the gcd, it will figure out itself w hat x is:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "gcd(F, diff(F, x));" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#,L*&\"(+g&>\"\"\"%\"xGF&!\"\" \"'++sF&*&\"'o*)GF&)F'\"#5F&F&*&\"&1V)F&)F'\"#6F&F&*&\"&s'fF&)F'\"#7F& F(*&\"&1v\"F&)F'\"#8F&F(*&\"%swF&)F'\"#9F&F&*&\"%7EF&)F'\"#:F&F&*&\"$: &F&)F'\"#;F&F(*&\"$S#F&)F'\"#F&F &*$)F'\"#?F&F&*&\"'+eSF&)F'\"\"#F&F&*&\"(5E0$F&)F'\"\"$F&F&*&\"(LL)=F& )F'\"\"%F&F(*&\"(]7=#F&)F'FKF&F(*&\"(j$*y\"F&)F'\"\"'F&F&*&\"'wj(*F&)F '\"\"(F&F&*&\"'*G3*F&)F'\"\")F&F(*&\"'=4KF&)F'\"\"*F&F(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 18 "OK, it's the same." }}}{EXCHG {PARA 0 "" 0 "" {TEXT 256 11 "Ass ignment:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 114 "Write an algorithm to compute a square-free factorization (just l ike Maple's command sqrfree) using the following:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 36 "differentiation (use Ma ple's diff)" }}{PARA 0 "" 0 "" {TEXT -1 33 "gcd (use the my_gcd pro cedure)" }}{PARA 0 "" 0 "" {TEXT -1 45 "Multiplying polynomials (use \+ Maple's expand)" }}{PARA 0 "" 0 "" {TEXT -1 53 "Dividing polynomials \+ (use the procedure my_quo_rem)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 112 "Write a procedure my_sqrfree. Use my_gcd for gcd's and my_quo_rem for computing a quotient. It should look lik e:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 54 "my_ sqrfree := proc(f,x)\n options trace;\n local ..." }}{PARA 0 "" 0 "" {TEXT -1 6 " ..." }}{PARA 0 "" 0 "" {TEXT -1 4 "end;" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 304 "The line \"option s trace;\" is there only for debugging. It helps you find the bugs bec ause as long as that line is there, all the steps in the algorithm wil l be shown during the computation, which makes it easier to find the e rrors. Once you have corrected all mistakes and tested that it works c orrectly," }}{PARA 0 "" 0 "" {TEXT -1 72 "then you can remove the line \"options trace;\" or place # in front of it." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 77 "The output of your proced ure should be a list:\n [ [f1,e1], [f2,e2], ... ]" }}{PARA 0 "" 0 "" {TEXT -1 10 "such that:" }}{PARA 0 "" 0 "" {TEXT -1 25 " f = f1^e1 * f2^e2 * ..." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 68 "and such that: gcd(f.i, diff(f.i, x) ) = 1 for all f.i i n the list," }}{PARA 0 "" 0 "" {TEXT -1 58 "and such that: gcd(f.i, f. j) = 1 for all i not equal to j." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 63 "So the f.i should have no factors in comm on with the other f.j." }}{PARA 0 "" 0 "" {TEXT -1 92 "And f.i should \+ have no multiple factors, i.e. the gcd of f.i and its derivative shoul d be 1." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 201 "Test your procedure on the polynomial F above. If it doesn't work , then first try very small examples because that makes it easier to f ind the bugs. First try it on F=x or F=x^2 or something like that." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 117 "Unless y ou're an experienced Maple programmer, you should contact me during of fice hours for help with this exercise." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 39 "Here are some hints how to get start ed:" }}{PARA 0 "" 0 "" {TEXT -1 72 "The first step in the exercise sho uld be to compute the gcd of f and f'." }}{PARA 0 "" 0 "" {TEXT -1 156 "If this gcd is 1, then the output should be: [ [f,1] ].\nIf this gcd is not 1, then apply a recursive call to find the squarefree fact orization of this gcd." }}{PARA 0 "" 0 "" {TEXT -1 129 "Suppose this r ecursive call gives you: [[f1,e1], [f2,e2], [f3,e3], ... ]\nNow compu te: f divided by: f1^(e1+1) * f2^(e2+1) * ...." }}{PARA 0 "" 0 "" {TEXT -1 87 "Say that's q. Then your output will be: [[q,1], [f1,e1+1 ], [f2,e2+1], [f3,e3+1], ... ]" }}}}{MARK "0 2 0" 1 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }