{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 56 "Let f be an element of Q[x ], and p be a rational number." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 16 "As we have seen:" }}{PARA 0 "" 0 "" {TEXT -1 17 " rem(f, x-p, x);" }}{PARA 0 "" 0 "" {TEXT -1 56 "gives y ou the value of f at the point p. More generally," }}{PARA 0 "" 0 "" {TEXT -1 20 " rem(f, (x-p)^n, x);" }}{PARA 0 "" 0 "" {TEXT -1 267 "giv es you an approximation of f in Q[[ x-p ]] of accuracy n, which means \+ the taylor series of f at x=p of order n (i.e. you give the coefficien ts of (x-p)^i for i=0..n-1, or another way to look at is is that you \+ give the derivatives of f at x=p of order 0,1,..,n-1)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 51 "Let p and q be two d ifferent rational numbers. Now:" }}{PARA 0 "" 0 "" {TEXT -1 16 " rem(f , x-p, x);" }}{PARA 0 "" 0 "" {TEXT -1 22 "is the value at p, and" }} {PARA 0 "" 0 "" {TEXT -1 14 " rem(f,x-q,x);" }}{PARA 0 "" 0 "" {TEXT -1 217 "is the value at q. So these two remainders tell you no more an d no less about f than just the values at p and q. We have seen that i f you have a polynomial m, and p and q are roots of m, then the follow ing polynomial:" }}{PARA 0 "" 0 "" {TEXT -1 12 " rem(f,m,x);" }}{PARA 0 "" 0 "" {TEXT -1 78 "has the same values at x=p,q as f. Therefore, i f we take m:=(x-p)*(x-q), then:" }}{PARA 0 "" 0 "" {TEXT -1 12 " rem(f ,m,x);" }}{PARA 0 "" 0 "" {TEXT -1 86 "contains exactly the same infor mation (namely the values of f at p and q) as the pair:" }}{PARA 0 "" 0 "" {TEXT -1 28 " rem(f,x-p,x), rem(f,x-q,x)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 74 "More generally, if p1, p2 and p3 are three different rational numbers, and" }}{PARA 0 "" 0 "" {TEXT -1 15 " m1:=(x-p1)^n1;" }}{PARA 0 "" 0 "" {TEXT -1 15 " m2:=(x-p 2)^n2;" }}{PARA 0 "" 0 "" {TEXT -1 15 " m3:=(x-p3)^n3;" }}{PARA 0 "" 0 "" {TEXT -1 17 " r1:=rem(f,m1,x);" }}{PARA 0 "" 0 "" {TEXT -1 17 " r 2:=rem(f,m2,x);" }}{PARA 0 "" 0 "" {TEXT -1 17 " r3:=rem(f,m3,x);" }} {PARA 0 "" 0 "" {TEXT -1 212 "then r1 gives you the value of f and the value of the derivatives of f (up to order n1 - 1) at the point p1. L ikewise for r2 and r3. But all of this information is also still conta ined in the following remainder:" }}{PARA 0 "" 0 "" {TEXT -1 13 " m:=m 1*m2*m3;" }}{PARA 0 "" 0 "" {TEXT -1 15 " r:=rem(f,m,x);" }}{PARA 0 " " 0 "" {TEXT -1 60 "because r1=rem(r,m1,x), r2=rem(r,m2,x) and r3=re m(r,m3,x);" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 307 "In a way, you can think of r as some kind of interpolation of \+ the information given by r1, r2 and r3. When n1=n2=n3=1 then it really is an interpolation, because then r1, r2 and r3 are just the values o f f at p1, p2, p3, and r is the polynomial that you get by interpolati on of these values at these points." }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 0 "" 0 "" {TEXT -1 139 "We see that going from r to r1, r2, r3 \+ is done simply by taking the remainder. In the special case n1=n2=n3 t his is the same as evaluation." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 437 "The inverse process, going from r1, r2, \+ r3 to r is called \"Chinese Remainder\", We see that in the special ca se n1=n2=n3 the Chinese Remainder is the same thing as interpolation. \+ We have already seen the usefullness of the pair \"evaluation\" and \" interpolation\" in the fast multiplication of polynomials by FFT. Now \+ Chinese Remainder is simply a generalization of this same principle, a nd we'll see that there are many useful applications." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 161 "Let R be a ring, and \+ let m1, m2, ..., mk in R, all having gcd 1. (e.g. R = Z and m1=6, m2 =7, m3=25. Or R=Q[x] and m1=(x-3), m2=(x^2+2)^2, m3=(x-5)^3*(x-6)^2). " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 25 "Let m =m1*m2*...*mk. Then:" }}{PARA 0 "" 0 "" {TEXT -1 26 "CHINESE REMAINDER THEOREM:" }}{PARA 0 "" 0 "" {TEXT -1 76 "The ring R/(m) is isomorphic to the ring S = R/(m1) * R/(m2) * ... * R/(mk)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 195 "We know how we can compute in R/(m) when R is a ring like Z or Q[x], we compute just like in R, but we always reduce elements of R t o a unique normal form modulo m, namely the remainder modulo m." }} {PARA 0 "" 0 "" {TEXT -1 72 "Likewise we can compute in each of the ri ngs R/(m1), R/(m2), ... R/(mk)." }}{PARA 0 "" 0 "" {TEXT -1 86 "An ele ment of the ring S, which is a product of the rings R/(m1), R/(m2),.. . , R/(mk)" }}{PARA 0 "" 0 "" {TEXT -1 110 "is represented by a list ( r1, r2, ..., rk) where r1 is an elements of R/(m1), r2 is an element o f R/(m2), etc." }}{PARA 0 "" 0 "" {TEXT -1 126 "If you have two elemen ts of S, say (r1, ..,rk) and (s1, ..,sk) then the sum and product in S are defined in the following way:" }}{PARA 0 "" 0 "" {TEXT -1 47 "(r1 ,..,rk) + (s1,...,sk) = (r1+s1, ... , rk+sk)" }}{PARA 0 "" 0 "" {TEXT -1 3 "and" }}{PARA 0 "" 0 "" {TEXT -1 46 "(r1, ..,rk)*(s1,..,sk) = (r1 *s1, ... , rk*sk)." }}{PARA 0 "" 0 "" {TEXT -1 121 "So one addition (o r multiplication) in S is the same as one addition (or multiplication ) in each of the R/(m.i), i=1..k." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 137 "The isomorphism from R/(m) to the ring S is given in the following way: if r in R/(m) then the image of r in S under this isomorphism is:" }}{PARA 0 "" 0 "" {TEXT -1 41 "(r mod m1, r mod m2, ..., r mod mk) in S." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 299 "The Chinese Remainder Theorem says that \+ this map from R/(m) to S is an isomorphism, so in particular there is \+ an inverse of this map. This map, for the case that the ring R is the \+ ring of integers, is given in Maple by the command chrem. Lets try som e examples to see if it really is an isomorphism." }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 7 "m1:=80;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "m2:=101;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 " m3:=103;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "igcd(m1,m2);" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "igcd(m1,m3);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "igcd(m2,m3);" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 12 "m:=m1*m2*m3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "ms:=[m1,m2,m3];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "r:=rand(0..10^5)() mod m;" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 24 "P:=randpoly(x,degree=5);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "P:=unapply(P,x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 246 "Suppose we want to do some computation in Z/(m). For example, \+ we want to compute P(r) in Z/(m). We can do that directly, but we can \+ also do it in the ring S = Z/(m1) * Z/(m2) * Z/(m3) because that ring \+ is isomorphic. And the results should match." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "P(r) mod m;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "s:=[seq(r mod i,i=ms)];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "r=chrem(s,ms);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "123*s;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "chrem(%,ms );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "123*r mod m;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "Ps:=[seq(P(s[i]) mod ms[i],i =1..3)];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "chrem(Ps,ms);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "P(r) mod m;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 207 "A typical application of Chinese Remaindering is to speed up a lgorithms. Lets suppose you have a very large matrix A with integer en tries and you want to compute the determinant of A. It is easy to see \+ that:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 70 " Hadamard bound for the determinant (see also the help page ?hadamard): " }}{PARA 0 "" 0 "" {TEXT -1 82 "abs(det(A)) <= length of the vector c ol(A,1) * length of the vector col(A,2) * ..." }}{PARA 0 "" 0 "" {TEXT -1 178 "where the length of the i'th column equals sqrt(sum of s quares of entries of i'th column). The bound is sharp for orthogonal m atrices, when all vectors have angles of 90 degrees." }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 103 "Lets say this bound is B. If you compute d is det(A) modulo m, and you know that abs(det(A)) <= B, then" }}{PARA 0 "" 0 "" {TEXT -1 18 "d := det(A) mod m;" }} {PARA 0 "" 0 "" {TEXT -1 64 "lets you determine det(A) precisely, prov ided that m >= 2*B + 1." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 391 "If we compute det(A ) directly by elimination over Z one can get exponential coefficients \+ growth and the computation for a 100 by 100 matrix can take \"forever \". If you do the same computation over Z/(m) it will go much faster, \+ you're avoiding the many-million digit integers that way. But now supp ose that 2*B+1 is about 10^1000. So you're still computing with intege rs that have 1000 digits." }}{PARA 0 "" 0 "" {TEXT -1 288 "The compute r can work very fast with integers that are small enough so that the p rocessor can multiply them in 1 single machine command. Lets assume th at the computer is very fast when you're computing modulo a number m.i that are smaller than 10^11. Then what you do is you take numbers:" } }{PARA 0 "" 0 "" {TEXT -1 17 "m1, m2, ..., m100" }}{PARA 0 "" 0 "" {TEXT -1 94 "that all have gcd 1 (the easiest is just to take prime nu mbers < 10^11) such that the product:" }}{PARA 0 "" 0 "" {TEXT -1 17 " m = m1*m2*..*m100" }}{PARA 0 "" 0 "" {TEXT -1 18 "is at least 2*B+1." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 257 "Then y ou compute det(A) first in Z/(m1), then in Z/(m2) etc. (note that if \+ your computer has more than 1 processor then you can do some of these \+ these determinants modulo m.i in parallel, i.e. you compute several at the same time, on different processors)." }}{PARA 0 "" 0 "" {TEXT -1 64 "When r1,r2,.. are the determinants modulo m1, m2,.. then you do:" }}{PARA 0 "" 0 "" {TEXT -1 34 "d:=chrem([r1,r2,..], [m1,m2,...]);" }} {PARA 0 "" 0 "" {TEXT -1 37 "and you get the determinant modulo m." }} {PARA 0 "" 0 "" {TEXT -1 56 "Then you do mods(d,m) and you have the de terminant in Z." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 16 "----------------" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 77 "We will b e looking more closely at another application, namely the following:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 241 "Let f, g in Z[x] be polynomials with integer coefficients, and let d in Z[x] be the gcd taken in the ring Z[x] (so not Q[x] as we had before), whi ch means that f=f'*d and g=g'*d for some f', g' in Z[x] (not just f', \+ g' in Q[x]). Furthermore:" }}{PARA 0 "" 0 "" {TEXT -1 51 " content(d,x ) = igcd( content(f,x), content(g,x) );" }}{PARA 0 "" 0 "" {TEXT -1 293 "There exists a bound for the size of the coefficients of any fact or of f in Z[x], called the Landau-Mignotte bound. So the coefficients of d have absolute value <= B where B is the Landau-Mignotte bound. T herefore we can determine the gcd d with the mods command when we have d modulo m, where" }}{PARA 0 "" 0 "" {TEXT -1 11 "m >= 2*B+1." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 122 "We will \+ take prime numbers p1, p2, ... such that the product is m >= 2*B+1, an d we will compute the gcd modulo p1, p2, ..." }}{PARA 0 "" 0 "" {TEXT -1 69 "As we have seen, computing the gcd modulo p can be done very qu ickly." }}{PARA 0 "" 0 "" {TEXT -1 448 "Since computing with 5 or 10 d igits costs about the same time as computing with 1 or 2 digits (this \+ depends on the processor, if you had a processor that could only do ve ry small integers in 1 operation then one would notice a difference be tween 10 digits and 2 digits), it will be a good idea to use prime num bers with not too few digits. Using very small prime numbers would mea n that we need more primes in order for the product m to be >=2*B+1." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 32 "There a re a few problems though:" }}{PARA 0 "" 0 "" {TEXT -1 16 "What we need is:" }}{PARA 0 "" 0 "" {TEXT -1 7 "d mod p" }}{PARA 0 "" 0 "" {TEXT -1 54 "for sufficiently many (i.e. product >=2*B+1) primes p." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 18 "However, \+ if we do:" }}{PARA 0 "" 0 "" {TEXT -1 21 "dp := Gcd(f,g) mod p;" }} {PARA 0 "" 0 "" {TEXT -1 77 "then there are three possible reasons why dp could be different from d mod p." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 74 "Bad case 1): p divides lcoeff(d,x). Man y things can go wrong in this case." }}{PARA 0 "" 0 "" {TEXT -1 145 "B ad case 2): degree(dp,x) > degree(d,x). Cases 1) and 2) are called \"b ad primes\". There are only finitely many bad primes as we will prove \+ later." }}{PARA 0 "" 0 "" {TEXT -1 81 "Bad case 3): d.i equals some c onstant c in Fp.i times the polynomial d mod p.i." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 117 "Note: Case 1) can of cou rse easily be avoided, by avoiding all primes that divide dl=gcd( lcoe ff(f,x), lcoeff(g,x) )." }}{PARA 0 "" 0 "" {TEXT -1 341 "When we exclu de case 1), then it must be so that \"d mod p\" is a divisor in Fp[x] \+ of degree = degree(d) of both polynomials f and g in Zp[x]. Hence, Gcd (f,g) mod p; will have degree at least degree(d), and furthermore Gcd( f,g) mod p; will be some multiple of \"d mod p\". Therefore, we must h ave either case 2) or case 3) when we avoid case 1)." }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 237 "We know that lcoeff(d, x) divides dl:=igcd( lcoeff(f,x), lcoeff(g,x). So, if dl=1 then (we ma y take lcoeff(d,x)>0) we know that lcoeff(d,x)=1. In this case it is e asy to avoid case 3), by always making the result of Gcd(f,g) mod p; m onic." }}{PARA 0 "" 0 "" {TEXT -1 819 "If dl>1 then it is slightly har der. What we do then is instead of making Gcd(f,g) mod p; monic, we mu ltiply it by an element in Fp such that the result will have lcoeff = \+ dl mod p. Then we apply Chinese remainder on all of these polynomials \+ with lcoeff=dl mod p, so the result will be a polynomial in Z[x]/(p1*p 2*..) with lcoeff=dl. Then we apply mods(.., m) to get the polynomial \+ in Z[x], and then we divide it by it's content, and multiply it with c ontent(d,x), (note that we already know a priori what content(d,x) is, see above). This way the case dl>1 can be handled too. Notice that in the Hensel lifting worksheets the same problem could occur, that is w hy we always assumed the polynomial f we wanted to factor to be monic. If f is not monic, then that can be handled in pretty much the same w ay as we did here." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 244 "From now we will assume for simplicity that dl=1, so cas e 3) can be avoided by making Gcd(f,g) mod p; monic, and case 1) is au tomatically avoided when dl=1. How about the \"bad primes\" for case t wo? Are those easy to recognize? The answer is no:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "f:=x^5+23*x+12;" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 18 "g:=x^6+15*x^3+2*x;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "ifactor( resultant(f,g,x) );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "p:=46877;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "f:=expand( f*(x^2-2) );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "g:=expand( g*(x^2-2) );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "gcd(f,g);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "Gcd(f,g ) mod p;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "p:=rand(0..10^3 )();" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "p:=nextprime(p);" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "Gcd(f,g) mod p;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "mods(%,p);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 73 "In the following example you can see more easil y what the bad primes are:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "f:=(x-1)*(x+1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "g:=( x-6)*(x+5);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "sol_f, sol_g := [solve(f)] , [solve(g)];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "R:=mul(mul( i-j ,i=sol_f), j=sol_g);" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 20 "R:=resultant(f,g,x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "ifactor(R);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "[f,g] mod 2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "[f,g] mod 3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "[f,g] mod 5;" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "[f,g] mod 7;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 483 "We see that in the example f,g above, we have gcd(f,g)=1, but \"Gcd(f,g) mod p\" is not 1 for all primes divid ing the resultant R. Here the resultant is two monic polynomials f,g i s defined as the product of all i-j for all roots i of f and roots j o f g. We see that when the resultant is 0, then one of the i-j=0, and t hen f,g have a common factor x-i=x-j. Therefore, when the gcd of two m onic polynomials f and g is 1, the bad primes are exactly the primes t hat divide the resultant." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 390 "Computing a resultant is about the same cost as c omputing a gcd, so we should not compute a resultant in our modular_gc d algorithm, because that would roughly slow it down a factor 2. Besid es, if the gcd is not 1 then the resultant will be equal to 0 because \+ the resultant contains all factors i-j where i=[roots of f], j=[roots \+ of g], and so when the gcd<>1 then one of those factors is 0." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 97 "Primes th at divide igcd(lcoeff(f,x), lcoeff(g,x)) can be very bad. Look at the \+ following example:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "f:=(3 *x-1)*(x-1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "g:=(3*x-1)* (x+2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "[f,g] mod 3;" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 715 "We see that when p divides lcoeff (f,x) and lcoeff(g,x), then Gcd(f,g) mod p can still have degree equal to degree gcd(f,g). But nevertheless, the prime p=3 in the above exam ple is bad for two reasons. The first one is that gcd(f,g) contains th e factor 3*x-1 which is missing from \"Gcd(f,g) mod p\", and the secon d reason p is bad is that \"Gcd(f,g) mod p\" has a factor (x-1)=(x+2) which does not appear in gcd(f,g). If you would use one of such prime in the modular_gcd algorithm, then when you do m:=mods(d,m) you are b ound to find the wrong answer, because for the bad prime p we have tha t gcd(f,g) reduced modulo p is not the same as Gcd(f,g) mod p, and so \+ we could never reconstruct gcd(f,g) from Gcd(f,g) mod p." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "gcd(f,g);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "% mod 3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "Gcd(f,g) mod 3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1403 "modular_gcd := proc(f,g,x, B,initial_p)\n # Input f,g in Z[x]\n # bound B for absolute valu es of the coeffs of the gcd\n # Assumption: igcd( lcoeff(f,x), lcoe ff(g,x) ) = 1\n local p,d,m,dp;\n if nargs=5 then\n p:=i nitial_p\n else\n # Note: it would be faster to start with \+ larger values of p,\n # because m would go up quicker when we \+ have higher p's.\n p:=2;\n fi;\n d:=Gcd(f,g) mod p;\n \+ m:=p;\n while m<2*B+1 do\n p:=nextprime(p);\n dp:=Gcd( f,g) mod p;\n if degree(dp,x) > degree(d,x) then\n # \+ The prime p is bad, skip p:\n next\n elif degree(dp,x ) < degree(d,x) then\n # All primes up to this point were b ad.\n RETURN( procname(f,g,x,B,p) )\n fi;\n # T he following command applies Chinese Remaindering\n # on all coe fficients of d and dp:\n d:=chrem([d,dp],[m,p]);\n m:=m*p \n od;\n # Necessary in case the gcd has negative coefficients: \n d:=mods(d,m);\n # Note that at this point one could imagine t hat d is not the\n # gcd because degree(d,x)>degree(gcd). In that c ase:\n # d doesn't divide f (i.e. f/d not in Z[x]) or d doesn't div ide g.\n if divide(f,d) and divide(g,d) then\n # Now it is c ertain that d is the gcd\n d\n else\n # All primes so far were bad, so use higher primes:\n procname(f,g,x,B,nextpri me(p))\n fi\nend;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "tra ce(modular_gcd); # To see what happens." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "B:=100;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "f:=expand( (x^2+x+1)*(x -15) );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "g:=expand( (x^2+ x+1)*(x+15) );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "modular_g cd(f,g,x,B); # we see that 2,3 and 5 are bad primes." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "f:=x^2-1; # = (x-1)*(x+1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 154 "g:=x^2+x-30; # = (x-6)*(x+5);\n# this is the same \+ example as before, recall that the\n# primes 2,3,5, and 7 are the prim es that\n# divide resultant(f,g,x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "B:=1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "mo dular_gcd(f,g,x,B); # we see that 2,3,5,7 are bad primes." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "N:=100; g:=randpoly(x,degree=N,den se):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "f1:=expand(randpoly (x,degree=N,dense)*g):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "f 2:=expand(randpoly(x,degree=N,dense)*g):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "B:=10^10;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 135 "Now modular_gcd won't work, because f1 and f2 and the gcd g dont' have to be monic. The following will also work in the non-monic case:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1300 "modular_gcd := proc(f,g,x, B,initial_p)\n # Input f,g in Z[x]\n # bound B for absolute valu es of the coeffs of the gcd\n local p,d,m,dp,dl;\n if nargs=5 th en\n p:=initial_p\n else\n p:=2;\n fi;\n d:=G cd(f,g) mod p;\n m:=p;\n while m<2*B+1 do\n p:=nextprime(p );\n dp:=Gcd(f,g) mod p;\n if degree(dp,x) > degree(d,x) t hen\n # The prime p is bad, skip p:\n next\n \+ elif degree(dp,x) < degree(d,x) then\n # All primes up t o this point were bad.\n RETURN( procname(f,g,x,B,p) )\n \+ fi;\n d:=chrem([d,dp],[m,p]);\n m:=m*p\n od;\n # So far d is monic. But the lcoeff of d need not be monic,\n # in g eneral it divides dl:\n dl:=igcd(lcoeff(f,x),lcoeff(g,x));\n d:= d*dl;\n d:=mods(d,m);\n # It could be that multiplying d by dl w as too much, we\n # should have multiplied d by something that give s us the\n # correct content. The content of d should be the follow ing:\n dl:=igcd(content(f,x),content(g,x));\n d:=primpart(d)*dl; \n # Now content(d)=dl=content of the gcd.\n if divide(f,d) and \+ divide(g,d) then\n # Now it is certain that d is the gcd\n \+ d\n else\n # All primes so far were bad, use higher prime s.\n procname(f,g,x,B,nextprime(p))\n fi\nend;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 85 "modular_gcd(f1,f2,x,B,nextprime(10^ 4)); # Since (10^4)^3 > 2*B+1 we'll need 3 primes." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "gcd(f1,f2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "77" 0 }{VIEWOPTS 1 1 0 1 1 1803 }