{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 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 }{CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 256 "" 1 24 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 268 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 272 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 273 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 274 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 275 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 276 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 277 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 278 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 279 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 280 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 281 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 282 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 283 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 284 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 285 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 286 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 287 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 288 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 289 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 290 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 291 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 292 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 293 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 294 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 295 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 296 "" 0 1 0 0 0 0 1 1 0 0 0 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 }{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }1 0 0 0 8 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 2" 3 4 1 {CSTYLE "" -1 -1 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 } 0 0 0 -1 8 2 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 }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 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }} {SECT 0 {SECT 1 {PARA 3 "" 0 "" {TEXT -1 0 "" }{TEXT 256 55 "Factoring in Q[x] with Zassenhaus' algorithm, overview." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 39 "Let Q be the field of rat ional numbers." }}{PARA 0 "" 0 "" {TEXT -1 30 "Let Z be the ring of in tegers." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 52 "Let Fp = Z/pZ be the finite field with p elements." }}{PARA 0 " " 0 "" {TEXT -1 39 "Let Z_p be the ring of p-adic integers." }}{PARA 0 "" 0 "" {TEXT -1 39 "Let Q_p be the field of p-adic numbers." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 55 "Then Z is a subring of Z_p, and Q is a subfield of Q_p." }}{PARA 0 "" 0 "" {TEXT -1 99 "Note that if a is in Q, then a is in Z_p if and only if t he denominator of a is not divisible by p." }}{PARA 0 "" 0 "" {TEXT -1 45 "We have Fp is isomorphic to Z_p modulo p Z_p." }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 46 "If f in Q[x] then we wi ll factor f as follows:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 45 "*) First do a square-free factorization of f." }} {PARA 0 "" 0 "" {TEXT -1 46 " Then we have f = f1^1 * f2^2 * f3^3 \+ * ..." }}{PARA 0 "" 0 "" {TEXT -1 74 " and each of these factors \+ f.i is squarefree, so gcd(f.i, f.i') = 1." }}{PARA 0 "" 0 "" {TEXT -1 160 " Another way to say that is that discriminant(f.i, x) <> \+ 0.\n The discriminant of a polynomial is 0 if and only if f.i has roots with multiplicity > 1." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 52 "*) After that, we have a partial factoriz ation of f." }}{PARA 0 "" 0 "" {TEXT -1 51 " But the f.i's are not necessarily irreducible." }}{PARA 0 "" 0 "" {TEXT -1 89 " The sqr free algorithm can not factor f.i because it's main step is g:=gcd(f.i , f.i')" }}{PARA 0 "" 0 "" {TEXT -1 47 " which only gives us a tri vial factor g:=1." }}{PARA 0 "" 0 "" {TEXT -1 62 " So we must fact or f.i but we know that f.i is squarefree." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 77 "*) So we need an algorithm for factoring a square-free polynomial f in Q[x]." }}{PARA 0 "" 0 "" {TEXT -1 82 " First we divide f by its integer-content. The integ er-content is the largest" }}{PARA 0 "" 0 "" {TEXT -1 46 " ration al number c such that f/c in Z[x]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 41 "*) Then we may assume that f is in Z[x]. " }}{PARA 0 "" 0 "" {TEXT -1 31 " Then f is also in Z_p[x]." }} {PARA 0 "" 0 "" {TEXT -1 14 " Denote " }{TEXT 257 5 "f = " } {TEXT -1 31 "f mod p = f + (p) in Fp[x]." }}{PARA 0 "" 0 "" {TEXT -1 13 " Denote " }{TEXT 259 2 "l " }{TEXT -1 14 "= lcoeff(f,x)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 10 " N ow " }{TEXT 260 2 "l " }{TEXT -1 35 " is nonzero mod p <====> degre e(" }{TEXT 258 2 "f " }{TEXT -1 67 ",x) = degree(f,x).\n Assume f rom now on that p does not divide " }{TEXT 261 3 "l ." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 83 " Furthermore: di scriminant(f,x) is not zero (because f is squarefree, and the" }} {PARA 0 "" 0 "" {TEXT -1 86 " discriminant is defined in such a w ay that it is zero if and only if f has roots" }}{PARA 0 "" 0 "" {TEXT -1 29 " with multiplicity > 1)." }}{PARA 0 "" 0 "" {TEXT -1 19 " Then we have:" }}{PARA 0 "" 0 "" {TEXT -1 51 " \+ discriminant is nonzero mod p <====> " }{TEXT 262 2 "f " }{TEXT -1 23 "in Fp[x] is squarefree." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 37 " In the algorithm we'll test if " }{TEXT 263 2 "l " }{TEXT -1 45 "or the discriminant is zero mod p, and if so, " }}{PARA 0 "" 0 "" {TEXT -1 30 " we take another prime p." }} {PARA 0 "" 0 "" {TEXT -1 32 " After that, we may assume:" }} {PARA 0 "" 0 "" {TEXT -1 12 " " }{TEXT 264 3 "f " }{TEXT -1 43 "is squarefree and has the same degree as f." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 76 "*) According to Hensel' s lemma, we can then obtain all irreducible factors" }}{PARA 0 "" 0 " " {TEXT -1 43 " of f in Z_p[x] in the following way:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 20 " *) Fact or " }{TEXT 265 2 "f " }{TEXT -1 9 "in Fp[x]." }}{PARA 0 "" 0 "" {TEXT -1 85 " This gives us all irreducible factors of f in Z_p[x] up to accuracy 1." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 25 " *) Hensel lift." }}{PARA 0 "" 0 "" {TEXT -1 77 " This way we can increase the accuracy of o ur approximations of" }}{PARA 0 "" 0 "" {TEXT -1 43 " t he factors of f in Z_p[x]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 81 " Note that since the elements of Q_p or Z_p are infinite power series in p," }}{PARA 0 "" 0 "" {TEXT -1 80 " \+ we can not calculate the factors of f in Z_p[x] exactly, we can only find" }}{PARA 0 "" 0 "" {TEXT -1 39 " them up to some finite ac curacy." }}{PARA 0 "" 0 "" {TEXT -1 79 " The same is true for re al numbers, if r is a real number, and we have no" }}{PARA 0 "" 0 "" {TEXT -1 81 " additional information about r, then no matter how many digits we compute," }}{PARA 0 "" 0 "" {TEXT -1 54 " it's n ever enough to find the exact value of r." }}{PARA 0 "" 0 "" {TEXT -1 40 " However, if we know some integer " }{TEXT 266 2 "l " } {TEXT -1 17 "and we know that " }{TEXT 267 2 "l " }{TEXT -1 7 "* r is: " }}{PARA 0 "" 0 "" {TEXT -1 57 " (**) an integer with no \+ more than 100 digits" }}{PARA 0 "" 0 "" {TEXT -1 80 " and if we \+ can calculate r with 101 digits accuracy, then we can calculate" }} {PARA 0 "" 0 "" {TEXT -1 75 " the exact value of r. So whenever r is some real number and we know" }}{PARA 0 "" 0 "" {TEXT -1 80 " \+ nothing else, then no accuracy is high enough to prove that we hav e found" }}{PARA 0 "" 0 "" {TEXT -1 77 " the exact value of r. H owever, if we know that r is a rational number," }}{PARA 0 "" 0 "" {TEXT -1 28 " and we know a number " }{TEXT 268 1 "l" }{TEXT -1 46 " that is divisible by denom(r), and we know a" }}{PARA 0 "" 0 "" {TEXT -1 51 " bound for the absolute value of the integer " } {TEXT 269 2 "l " }{TEXT -1 29 "* r, then as soon as we know" }}{PARA 0 "" 0 "" {TEXT -1 81 " r up to some sufficiently high accuracy, we can figure out the exact value" }}{PARA 0 "" 0 "" {TEXT -1 78 " \+ of r. That's why \"finite accuracy\" + \"inside information like \+ in (**)\"" }}{PARA 0 "" 0 "" {TEXT -1 68 " results in \"infinite accuracy\", i.e. the precise value for r." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 71 " Now note that the same a rgument is true when we use not the real" }}{PARA 0 "" 0 "" {TEXT -1 37 " numbers but the p-adic number." }}{PARA 0 "" 0 "" {TEXT -1 81 " In the real numbers, if we had a floating point approximat ion of r with n" }}{PARA 0 "" 0 "" {TEXT -1 50 " digits accurate , and if we knew an integer " }{TEXT 270 2 "l " }{TEXT -1 22 "for whic h we knew that" }}{PARA 0 "" 0 "" {TEXT -1 7 " " }{TEXT 271 2 "l " }{TEXT -1 55 "* r was some integer with no more than 100 digits, th en" }}{PARA 0 "" 0 "" {TEXT -1 57 " if we knew r up to n=80 digi ts, we can not find r." }}{PARA 0 "" 0 "" {TEXT -1 56 " But if w e knew r up to n=110 digits, we can find " }{TEXT 272 2 "l " }{TEXT -1 17 "* r by computing:" }}{PARA 0 "" 0 "" {TEXT -1 14 " round( " }{TEXT 273 2 "l " }{TEXT -1 6 "* r )." }}{PARA 0 "" 0 "" {TEXT -1 78 " Here the Maple round command determines the nearest integer . Obviously," }}{PARA 0 "" 0 "" {TEXT -1 80 " if we have an 100- digit integer for which we have an 110-digits accurate," }}{PARA 0 "" 0 "" {TEXT -1 78 " floating point approximation, then we can mak e that approximation exact" }}{PARA 0 "" 0 "" {TEXT -1 46 " by r ounding off to the nearest integer." }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 0 "" 0 "" {TEXT -1 82 " With p-adic approximations it wor ks very similar, but we have to use \"mods\"" }}{PARA 0 "" 0 "" {TEXT -1 24 " instead of round." }}{PARA 0 "" 0 "" {TEXT -1 11 " \+ If " }{TEXT 274 2 "l " }{TEXT -1 77 "* r is an integer whose absolut e value is less than B, then take an integer n" }}{PARA 0 "" 0 "" {TEXT -1 47 " such that p^n > 2*B + 1. Then, compute " }{TEXT 275 2 "l " }{TEXT -1 37 "* r up to p-adic accuracy n, that is," }} {PARA 0 "" 0 "" {TEXT -1 15 " compute " }{TEXT 276 2 "l " } {TEXT -1 67 "* r mod p^n. Now instead of \"round\" as for floating po int numbers," }}{PARA 0 "" 0 "" {TEXT -1 43 " we have to use \"m ods\". Take: mods(" }{TEXT 277 2 "l " }{TEXT -1 16 "* r, p^n). If \+ " }{TEXT 278 2 "l " }{TEXT -1 27 "* r really was an integer" }} {PARA 0 "" 0 "" {TEXT -1 80 " with absolute value < B, then this command results in the exact value of " }{TEXT 279 2 "l " }{TEXT -1 4 "* r." }}{PARA 0 "" 0 "" {TEXT -1 48 " We can then easily find r by dividing by " }{TEXT 282 2 "l." }}{PARA 0 "" 0 "" {TEXT -1 16 " \+ But if " }{TEXT 280 2 "l " }{TEXT -1 72 "* r was not an integ er (for example if r was not a rational number), or" }}{PARA 0 "" 0 " " {TEXT -1 86 " if it was an integer whose absolute value is bi gger than B, then the result of" }}{PARA 0 "" 0 "" {TEXT -1 13 " \+ mods(" }{TEXT 281 2 "l " }{TEXT -1 70 "* r, p^n) is in essence garba ge. So the bound B is very important for" }}{PARA 0 "" 0 "" {TEXT -1 54 " turning an approximation into an exact result." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 82 "*) We will use \+ factoring in Fp[x] and Hensel lifting not for f but to factor f/" } {TEXT 283 2 "l." }}{PARA 0 "" 0 "" {TEXT -1 13 " Now f/" }{TEXT 284 2 "l " }{TEXT -1 57 "is monic and that is convenient. Since p does not divide " }{TEXT 285 3 "l, " }{TEXT -1 7 "we have" }}{PARA 0 "" 0 "" {TEXT -1 14 " that 1/" }{TEXT 286 3 "l " }{TEXT -1 27 "is a \+ p-adic integer, so f/" }{TEXT 287 3 "l " }{TEXT -1 30 "is still an e lement of Z_p[x]." }}{PARA 0 "" 0 "" {TEXT -1 43 " We find the f ollowing factorization:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 13 " f/" }{TEXT 288 7 "l = " }{TEXT -1 17 " f1 * f2 * .. * fn" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 70 " where f1, f2,..,fn are all monic irreducible facto rs in Z_p[x]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 28 "*) All monic factors of f/" }{TEXT 289 2 "l " }{TEXT -1 40 "in Q_p[x] are in Z_p[x] (Gauss' lemma)." }}{PARA 0 "" 0 "" {TEXT -1 64 " There are precisely 2^n such factors. The smallest is g=1" }}{PARA 0 "" 0 "" {TEXT -1 32 " and the largest is g = \+ f/" }{TEXT 290 2 "l." }}{PARA 0 "" 0 "" {TEXT -1 73 " To find al l of these factors, take all subsets of \{f1, f2,.., fn\}." }}{PARA 0 "" 0 "" {TEXT -1 43 " For each subset S, take the product." }} {PARA 0 "" 0 "" {TEXT -1 53 " That way, we find all monic factor s in Z_p[x]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 61 " If there are k irreducible monic factors in Q[x], say:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 12 " \+ f/" }{TEXT 291 3 "l " }{TEXT -1 41 "= g1 * g2 * .. * gk where g .i in Q[x]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 70 " then of these 2^n monic factors in Z_p[x] there are pr ecisely" }}{PARA 0 "" 0 "" {TEXT -1 71 " 2^k that are in Q[x]. \+ This is because Q is contained in Q_p and" }}{PARA 0 "" 0 "" {TEXT -1 133 " so g1 .. gk are in Q_p[x], and hence in Z_p[x] by Gauss' \+ lemma,\n and so each g.i must be a product of some of the f.j. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 64 "*) \+ We are now facing a problem: We can only compute the f.i up" }}{PARA 0 "" 0 "" {TEXT -1 73 " to finite accuracy. Given a subset of \{ f1,..,fn\} we have to decide" }}{PARA 0 "" 0 "" {TEXT -1 70 " if the product of this subset is in Q[x] or not. But how can we" }} {PARA 0 "" 0 "" {TEXT -1 63 " decide that if we are only working with finite accuracy?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 46 " Zassenhaus solved this problem in 1969:" }} {PARA 0 "" 0 "" {TEXT -1 65 " Let g be the product of that s ubset. We then now g not" }}{PARA 0 "" 0 "" {TEXT -1 55 " pr ecisely, but only up to a finite accuracy." }}{PARA 0 "" 0 "" {TEXT -1 70 " Now g is a factor of f. If g were an element of Q[x ], then" }}{PARA 0 "" 0 "" {TEXT -1 40 " by Gauss' lemma we \+ have that " }{TEXT 292 2 "l " }{TEXT -1 26 "*g is in Z[x]. Then by the " }}{PARA 0 "" 0 "" {TEXT -1 63 " Landau-Mignotte bound, we \+ can find a bound B for the" }}{PARA 0 "" 0 "" {TEXT -1 49 " \+ absolute value of the coefficients of " }{TEXT 293 2 "l " }{TEXT -1 3 "*g." }}{PARA 0 "" 0 "" {TEXT -1 66 " Now, if the p-power is greater than 2*B+1, then we have" }}{PARA 0 "" 0 "" {TEXT -1 38 " \+ enough accuracy to recover " }{TEXT 294 2 "l " }{TEXT -1 26 "*g with infinite accuracy." }}{PARA 0 "" 0 "" {TEXT -1 36 " Th en we can divide it by " }{TEXT 295 2 "l " }{TEXT -1 34 "and we find g . If we want to avoid" }}{PARA 0 "" 0 "" {TEXT -1 58 " frac tions then we could also not divide it by " }{TEXT 296 2 "l " }{TEXT -1 20 "but divide it by its" }}{PARA 0 "" 0 "" {TEXT -1 28 " \+ integer-content." }}{PARA 0 "" 0 "" {TEXT -1 43 " Then we have two possibilities:" }}{PARA 0 "" 0 "" {TEXT -1 58 " \+ *) Either we found an exact factor in Q[x]." }}{PARA 0 "" 0 "" {TEXT -1 74 " *) Or: g wasn't in Q[x] at all, in which c ase the result of" }}{PARA 0 "" 0 "" {TEXT -1 47 " \+ the computation is garbage." }}{PARA 0 "" 0 "" {TEXT -1 74 " \+ To test which of these two is true, all we have to do is check" }} {PARA 0 "" 0 "" {TEXT -1 76 " if the factor that we found i s a true exact factor of f in Q[x]." }}{PARA 0 "" 0 "" {TEXT -1 69 " \+ We can do this by the \"divide\" command (or we can use our " }}{PARA 0 "" 0 "" {TEXT -1 68 " own implementation, we ca n copy the quo_rem command from" }}{PARA 0 "" 0 "" {TEXT -1 33 " \+ an earlier worksheet)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 74 " For the 2^n subsets of \{f1,..fn\} there \+ will be 2^k subsets where the" }}{PARA 0 "" 0 "" {TEXT -1 72 " re sult is not garbage, in which case divide(f, factor_we_found)" }} {PARA 0 "" 0 "" {TEXT -1 30 " gives the output \"true\"." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 64 "*) Note \+ that we do not want to find all 2^k monic factors of f," }}{PARA 0 "" 0 "" {TEXT -1 53 " (which correspond to all subsets of \{g1,..,gk \})" }}{PARA 0 "" 0 "" {TEXT -1 34 " we only want to find g1..gk. " }}{PARA 0 "" 0 "" {TEXT -1 63 " So we want to avoid finding big ger factors than g1 .. gk." }}{PARA 0 "" 0 "" {TEXT -1 45 " This \+ is easily accomplished as follows:" }}{PARA 0 "" 0 "" {TEXT -1 69 " \+ *) First try all subsets of \{f1,..,fn\} with only 1 element." }}{PARA 0 "" 0 "" {TEXT -1 53 " *) Then those subsets that ha ve 2 elements." }}{PARA 0 "" 0 "" {TEXT -1 46 " *) Then 3,\n \+ *) then 4, etc." }}{PARA 0 "" 0 "" {TEXT -1 47 " What we also have to do is the following:" }}{PARA 0 "" 0 "" {TEXT -1 81 " \+ *) As soon as some subset S of \{f1,..,fn\} gives us an exact t rue factor" }}{PARA 0 "" 0 "" {TEXT -1 75 " g in Q[x] of f, then we remove the f.i in S from \{f1,..,fn\}." }}{PARA 0 "" 0 "" {TEXT -1 152 " If we set up the algorithm in that way, then we'l l only find irredudible\n factors g1 .. gk in Q[x] and we will n ot find products of g1..gk." }}{PARA 0 "" 0 "" {TEXT -1 7 " " }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 428 "my_factor : = proc(f)\n local i,x;\n\n # Figure out all names that occur in f.\n x := indets(f);\n if x=\{\} then\n # f has no variables, i.e. f \+ is a constant,\n # so its factorization in Q[x] is trivial:\n RE TURN(f)\n fi;\n for i in x do if not type(i,'name') then\n ERROR (f, \"should be a polynomial\")\n fi od;\n\n if nops(x)=1 then\n \+ factor_univariate_case(f, x[1])\n else\n factor_multivariate_cas e(f, x)\n fi:\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 194 "factor_multivariate_case : = proc()\n ERROR(\"more than one variable not yet implemented\")\n # We will do this later. We have already seen\n # partially h ow to do it for 2 variables.\nend:" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 23 "Hensel lift algorithms." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {EXCHG {PARA 0 "" 0 "" {TEXT -1 331 "See section 6.6 in the book for t he bound used in the following algorithm \"required_accuracy\". Note \+ that we multiply the bound by l because we're multiplying the product \+ g in the procedure factor_sqrfree by l. That way the rational number c oefficients become integer coefficients by Gauss' lemma, which is sect ion 6.2 in the book." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 144 "re quired_accuracy := proc(f,l,x,p)\n local B,n;\n n := degree(f,x);\n \+ B := l * 2^n*sqrt(n+1)*maxnorm(f);\n ceil( evalf(log(2*B+1)/log(p)) \+ )\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 182 "The following algorit hm is there so we can choose which Hensel-algorithm we'll use. That wa y we can compare them more easily, and we can then test experimentally which one is faster." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 288 "H enselLift := proc()\n # To change the algorithm that is used,\n # ju st uncomment (remove #) from one line\n # and comment (place #) the o ther lines:\n\n # HenselLift_linalg(args)\n # HenselLift_linalg_bigg ersteps(args)\n # HenselLift_gcdex(args)\n HenselLift_gcdex_biggerst eps(args)\n\nend:" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 54 "Hensel lift ing using the extended Euclidean Algorithm." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 944 "Hensel_lift_two_factors := proc(f,g,h,x,N,p)\n # I nput: a positive integer n, and\n # monic polynomials f,g,h in Z[x] \+ such that:\n # h is f*g mod p.\n\n # Output: monic polynomials F, \+ G such that:\n # h is F * G mod p^N.\n local s,t,H,f0,g0,R,q,F,G,n ;\n \n # Compute s and t:\n if Gcdex(f,g,x,'s','t') mod p<>1 the n\n ERROR( \"f and g should have Gcd 1 modulo p\")\n fi;\n F, G := f, g; # Now the accuracy is 1.\n H := h mod p^N;\n for n fro m 1 to N-1 do\n # Compute f0, g0 such that (F+f0*p^n)*(G+g0*p^n ) is\n # H modulo p^(n+1).\n # Take f0, g0 = R*t, R*s as explained above.\n R:=( Expand(H-F*G) mod p^(n+1))/p^n;\n \+ f0, g0 := R*t, R*s;\n # Now reduce f0, g0 to meet the degree conditions,\n # while at the same time keeping f0*g+g0*f the\n # same, namely R.\n f0:=Rem(f0,f,x,'q') mod p;\n \+ g0:=Expand(g0+q*g) mod p;\n F, G := F+f0*p^n, G+g0*p^n\n od; \n F, G\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1005 "HenselL ift_gcdex:=proc(L::list,h,x::name,N::posint,p::posint)\n # Recursive algorithm for Hensel lifting when there are\n # more than 2 factors . Note: This is perhaps not the most\n # efficient way to do it but \+ it's fine for us. For more\n # details see chapter 15 in the book.\n #\n # Input: a list L=[f1,f2,...fk] of polynomials such that\n \+ # f1*f2*...*fk is h modulo p, and all f.i have Gcd 1.\n #\n # Outp ut: a list [F1,F2,..Fk] such that\n # F1*F2*..*Fk is h modulo p^n.\n local l,L1,L2,h1,h2;\n l:=nops(L);\n if l=0 then\n ERROR( \"wrong input\")\n elif l=1 then\n [h]\n elif l=2 then\n \+ [Hensel_lift_two_factors(op(L),h,x,N,p)]\n else\n L1, L2 := L[1..iquo(l,2)], L[iquo(l,2)+1..-1];\n # Take the products.\n \+ h1, h2 := convert(L1,`*`), convert(L2,`*`);\n # Lift the pr oducts:\n h1, h2 := Hensel_lift_two_factors(h1,h2,h,x,N,p);\n \+ # Lift L1 and L2 and return the answer\n [op(procname(L1,h1, x,N,p)), op(procname(L2,h2,x,N,p))]\n fi\nend:" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 43 "Hensel lifting by solving linear equations." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1629 "HenselLift_linalg:=proc(L, f,y,n,p)\n # Inpute f in Z[y] monic in y.\n # L a list of factors \+ of f mod p^1\n # (so the elements in L have accuracy 1).\n # The e lements in L must have gcd 1 in Fp[y].\n #\n # n a positive intege r.\n #\n # Note: to make it easy I only implemented the case\n # where f is monic, and where the elements of L\n # are also monic (a s polynomials in y).\n #\n local i,j,k,c,Li,pol,cpol;\n\n Li:=L; # accuracy right now is i=1.\n\n for i from 1 to n-1 do\n \n \+ Li := [seq(\n Li[k] + p^i * add(c[j,k]*y^j,\n \+ j=0..degree(L[k],y)-1)\n ,k=1..nops(L))];\n \n # Now Li mod p^i did not change because we added\n # some thing that is 0 mod p^i.\n #\n # Li mod p^(i+1) did change. \+ We did not know any of\n # the coefficients of p^i*y^j in any of \+ the\n # factors. So we used an unknown c[j,k] for the\n # th e coefficient of p^i*y^j in the k'th factor. \n #\n # Now le ts compute equations for the c[j,k].\n # We know that the product convert(Li,`*`) is the\n # same as f up to accuracy i, so the fo llowing is\n # a polynomial:\n pol:=normal((convert(Li,`*`)- f)/p^i);\n # We want to choose values for the unknowns c[j,k]\n \+ # in such a way that Li is OK with accuracy i+1,\n # which me ans that pol is 0 mod x^1.\n pol:=Expand(pol) mod p;\n # Thi s pol should be the zero polynomial,\n # which means that the fol lowing coeffs\n # should be zero:\n cpol:=\{coeffs(pol,y)\}; \n # Solve and increase the accuracy of Li to i+1:\n Li:=sub s(solve(cpol) mod p,Li)\n od;\n Li\nend:" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 0 "" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 90 "Hensel liftin g by solving linear equations, increasing the accuracy more than 1 at \+ a time." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1627 "HenselLift_linalg_biggersteps:=proc(L,f,y,n,p)\n \+ # options trace;\n local i,j,k,c,Li,pol,cpol, inc_i;\n\n i := 1;\n Li:=L; # accuracy right now is i=1.\n\n while i < n do\n \n \+ Li := [seq(\n Li[k] + p^i * add(c[j,k]*y^j,\n \+ j=0..degree(L[k],y)-1)\n ,k=1..nops(L))];\n \n pol:=normal((convert(Li,`*`)-f)/p^i);\n # We want to choo se values for the unknowns c[j,k]\n # in such a way that Li is OK with accuracy i + inc_i\n # where inc_i > 0.\n # Since we'r e only going to accuracy n, we see that\n # there is no reason fo r inc_i to be bigger than n-i.\n # So inc_i <= n-i.\n\n # Ho wever, we want that the equations for the\n # unknowns c[j,k] to \+ be linear: If they are not\n # linear then they are difficult to \+ solve.\n\n # Now \"pol\" is not linear in the c[j,k], it contains \n # products of the unknowns c[j,k]. However, these\n # pro ducts have p-powers at least p^(2*i) in front\n # of them. So if \+ we mod out by p^(2*i), then \"pol\"\n # becomes linear in c[j,k]. \n # This puts another condition on inc_i. We must have:\n # inc_i <= i. Otherwise in the equations \"cpol\"\n # below will n ot be linear.\n\n # So we take something that's <=i and <= n-i. \n inc_i := min(i, n-i);\n \n pol:=Expand(pol) mod p^in c_i;\n\n # This pol should be the zero polynomial,\n # which means that the following coeffs\n # should be zero:\n cpol: =\{coeffs(pol,y)\};\n # Solve and increase the accuracy of Li to \+ i+1:\n Li:=subs( solve(cpol) mod p^inc_i, Li);\n i := i + in c_i\n od;\n Li\nend:" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 74 "Hensel lifting, extended Euclidean Algori thm, more than 1 step at at time." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 574 "Hensel_lift_two_factors_b := proc(f,g,h,x,N,p)\n local s,t,H,f0 ,g0,R,q,F,G,n, inc_n;\n \n # Compute s and t:\n if Gcdex(f,g,x,' s','t') mod p^floor(N/2)<>1 then\n ERROR( \"f and g should have G cd 1 modulo p\")\n fi;\n F, G := f, g; # Now the accuracy is 1.\n \+ H := h mod p^N;\n n:=1;\n while n < N do\n inc_n := min(n , N-n);\n R:=( Expand(H-F*G) mod p^(n+inc_n))/p^n;\n f0, g0 := R*t, R*s;\n f0:=Rem(f0,f,x,'q') mod p^inc_n;\n g0 :=Expand(g0+q*g) mod p^inc_n;\n F, G := F+f0*p^n, G+g0*p^n;\n \+ n := n + inc_n\n od;\n F, G\nend:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 491 "HenselLift_gcdex_biggersteps:=proc(L::list,h,x::na me,N::posint,p::posint)\n local l,L1,L2,h1,h2;\n l:=nops(L);\n i f l=0 then\n ERROR(\"wrong input\")\n elif l=1 then\n [h ]\n elif l=2 then\n [Hensel_lift_two_factors(op(L),h,x,N,p)]\n else\n L1, L2 := L[1..iquo(l,2)], L[iquo(l,2)+1..-1];\n \+ h1, h2 := convert(L1,`*`), convert(L2,`*`);\n h1, h2 := Hensel_ lift_two_factors(h1,h2,h,x,N,p);\n [op(procname(L1,h1,x,N,p)), o p(procname(L2,h2,x,N,p))]\n fi\nend:" }}}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 71 "Compute all monic irreduc ible factors in Z_p[x] for a suitable prime p." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1537 "monic_irreducible_p_adic_factors:=proc(f,x)\n l ocal p,l,fm,n,v,ppower,i;\n l := lcoeff(f,x);\n\n p := take_the_firs t_suitable_prime(f,l,x);\n \n # Now we factor f mod p, but lets first make it monic:\n # We know that Maple can calculate \"1/l mod p\" by \n # computing igcdex(l,p,'s','t'), and we know how\n # that algorit hm works. So we know how Maple can\n # compute the following: first \+ \"1/l mod p\" (which is\n # the number s in igcdex(l,p,'s','t')) and \+ then\n # multiply that by f, and then take the remainder\n # of the \+ coefficients mod p.\n fm := f/l mod p; \n\n # Now we have a monic p olynomial fm in Fp[x] where Fp\n # is a finite field with p elements. Furthermore, fm\n # has no multiple factors, all factors appear with \n # multiplicity 1.\n # We have not yet studied how to factor in Fp [x], so\n # for the moment we'll use Maple to do that for us:\n v := Factors(fm) mod p;\n\n # We have: v = [lc , [ [f1, m1], [f2,m2], ... ] ]\n # but we know that this leading coefficient is lc=1\n # and a ll multiplicities m1, m2,.. are also 1 because\n # we did not take a \+ prime p where the discriminant\n # Discrim(f,x) is zero mod p.\n\n v := [seq(i[1],i=v[2])]; # now v = [f1,f2,...]\n\n n := required_accur acy(f,l,x,p);\n ppower := p^n;\n v := HenselLift(v, f/l, x, n, p);\n v := \{op(v)\}; # Turn v into a set.\n # We've factored f/l in Z_p[ x] where Z_p = the ring\n # of p-adic integers. We've determined the \+ factors\n # up to accuracy n, i.e. we've computed them\n # modulo pp ower.\n\n # Now return v and the modulus ppower:\n v, ppower\nend:\n " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 129 "In the algorithm below, If l mod p is 0 then the degree of f mod p\n is smaller than the degree of f. Then we take the next prime." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 71 "If Discrim(f,x) mod p is 0 then f is not \+ square-free mod p. If we would" }}{PARA 0 "" 0 "" {TEXT -1 73 "factor \+ f mod p there would be factors of multiplicity > 1. So some factor" }} {PARA 0 "" 0 "" {TEXT -1 69 "would appear more than once in the factor ization. That's a problem in" }}{PARA 0 "" 0 "" {TEXT -1 71 "Hensel li fting because we can only Hensel lift factors that have gcd 1," }} {PARA 0 "" 0 "" {TEXT -1 80 "so the irreducible factors must all be di stinct (they must have multiplicity 1)." }}{PARA 0 "" 0 "" {TEXT -1 71 "Now the discriminant is something that is zero if and only if ther e are" }}{PARA 0 "" 0 "" {TEXT -1 80 "factors of multiplicity > 1, whi ch is true if and only if f and f' have a common" }}{PARA 0 "" 0 "" {TEXT -1 88 "factor, which is true if and only if the resultant of f a nd f' is zero (the discriminant" }}{PARA 0 "" 0 "" {TEXT -1 53 "of f i s some number times the resultant of f and f')." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 156 "take_the_f irst_suitable_prime := proc(f,l,x)\n local p;\n p := 2;\n while l m od p = 0 or Discrim(f,x) mod p = 0 do\n p := nextprime(p)\n od; \n p\nend:\n " }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 1 {PARA 3 " " 0 "" {TEXT -1 82 "Factoring univariate polynomials, reducing the pro blem to square-free polynomials." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 994 "factor_univariate_case := proc(f,x)\n local v,result,i,w;\n\n # Determine a square-free factorization. We know by\n # now of course \+ how that works. It's mainly based on\n # gcd(f, f') and this computat ion is fast because we\n # use modular arithmetic which prevents exce ssive\n # growth of coefficients, and if the coefficients don't\n # \+ grow then the Euclidean algorithm will be fast.\n v := sqrfree(f,x); \n\n result := v[1]; # v[1] is a constant, v[2] is\n \+ # a list of factors and multiplicities.\n\n for i in v[2] do\n \+ # Now i[1] is a square-free polynomial. So\n # gcd( i[1], i[1]' ) = 1 so the sqrfree\n # algorithm can not factor i[1] any further .\n # To factor a squarefree polynomial i[1] we'll\n # use t he following algorithm:\n w := factor_sqrfree(i[1],x);\n\n # Now we raise all of the factors we found\n # by the i[2]-th powe r because i[1] appeared\n # in f with multiplicity i[2].\n r esult := result * w^i[2]\n od;\n result\nend:\n" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 45 "Factoring univar iate square-free polynomials." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1482 "factor_sqrfree := proc(f,x)\n # Input: a squarefree polynomial \+ f in Q[x]\n # Output: The factorization.\n # options trace;\n local c,ff,v,ppower,result,n,l,su,S,g;\n if degree(f,x) < 2 then\n \+ # A reducible polynomial must have degree\n # at least 2, so we 're done:\n RETURN( f )\n fi;\n\n c := icontent(f);\n ff := \+ f/c;\n # Since we divided by the integer-content c, we\n # now have \+ a polynomial with integer coefficients.\n\n v, ppower := monic_irredu cible_p_adic_factors(ff,x);\n \n result := c; # That's the constant \+ that we divided\n # f by to obtain ff.\n n := 0;\n whi le v <> \{\} do\n l := lcoeff(ff,x);\n n := n+1;\n su : = [subsets(v,n)];\n # Go through all subsets of v with n elements :\n for S in su do\n g := convert(S,`*`);\n g := mods(l * Expand(g),ppower);\n # Now divide g by the gcd \+ of its coefficients:\n g := g/icontent(g);\n # We \+ know that modulo ppower, g divides ff.\n # So g is a factor \+ of ff, at least up to some\n # finite accuracy. We will now \+ test if g\n # divides ff up to infinite accuracy, that is,\n # test if q = ff/g is a polynomial or not: \n if \+ divide(ff,g,'q') then\n result := result * g;\n \+ ff := q; # ff = ff/g\n v := v minus S\n fi \n od\n od;\n if not member(ff, \{1,-1\}) then\n ERROR(\"Al l factors should have been removed by now\")\n fi;\n ff*result\nend: \n\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 75 "Compute the sequence S1,S 2,.. of all subsets with n elements of some set v." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 263 "subsets:=proc(v::set, n::posint)\n local \+ i,e,T;\n if n=1 then\n seq( \{e\}, e=v);\n elif n>nops(v) then\n NULL\n elif n=nops(v) then\n v\n else\n e := v[1];\n \+ T := v minus \{e\};\n subsets(T,n),\n seq(\{e\} union i, i =[subsets(T, n-1)])\n fi\nend:" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 9 "Examples." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 219 "Note: If you want to see how an algorithm, for example factor_ sqrfree, works then type: trace(factor_sqrfree) and try an example. No te: it's best to do that with small examples, then it's easier to see \+ what is going on." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "my_factor(x^2-1,x);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#*&,&%\"xG\"\"\"F&F&F&,&F%F&!\"\"F&F&" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "my_factor(x^12+1);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#*&,&*$)%\"xG\"\"%\"\"\"\"\"\"F*F*F*,(*$)F'\"\")F)F*F%! \"\"F*F*F*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "my_factor(16* x^2-1);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#*&,&%\"xG\"\"%!\"\"\"\"\"F( ,&F%F&F(F(F(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "my_factor(x ^30-2*x^15+1);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#**),,*$)%\"xG\"\"%\" \"\"\"\"\"*$)F(\"\"$F*F+*$)F(\"\"#F*F+F(F+F+F+F1F*),&F(F+!\"\"F+F1F*), (F/F+F(F+F+F+F1F*),0*$)F(\"\")F*F+*$)F(\"\"(F*F4*$)F(\"\"&F*F+F&F4F,F+ F(F4F+F+F1F*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "expand(%); " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#,(*$)%\"xG\"#I\"\"\"\"\"\"*$)F&\"# :F(!\"#F)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "T := time(): my_factor(x^52-1); computation_time = time()-T;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#*.,&*$)%\"xG\"\"#\"\"\"\"\"\"F*F*F*,&F'F*F*F*F*,&F'F*! \"\"F*F*,F*F5F*F8F*F%F*F'F*F*F*F*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%1computation_timeG$\"%!4&!\"$" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 25 "Rema rks and improvements." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 62 "There are several ways in which the algorithm can be i mproved:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 55 "First of all: At the bottom of factor_sqrfree you see:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 137 " while v <> \+ \{\} do\n n := n+1;\n su := [subsets(v,n)];\n # Go thro ugh all subsets of v with n elements:\n for S in su do\n" }} {PARA 0 "" 0 "" {TEXT -1 78 "Here v is the list of monic factors in Z_ p[x] that have not yet been used (the" }}{PARA 0 "" 0 "" {TEXT -1 86 " set of all monic irreducible factors in Z_p[x] of the polynomial ff in the algorithm)." }}{PARA 0 "" 0 "" {TEXT -1 72 "The improvement we ca n make is the following: As soon as n > nops(v)/2" }}{PARA 0 "" 0 " " {TEXT -1 82 "we can stop. After all, if there would be a factor g of ff in Q[x] consisting of n" }}{PARA 0 "" 0 "" {TEXT -1 69 "factors in Z_p[x], then there would also be another one (namely ff/g)" }}{PARA 0 "" 0 "" {TEXT -1 81 "that only consists of the other nops(v)-n facto rs. So if nops(v) - n < n then we" }}{PARA 0 "" 0 "" {TEXT -1 78 "can not possibly find something with n; because if that were so we should have" }}{PARA 0 "" 0 "" {TEXT -1 70 "already found something earlier \+ when we were trying g's that consisted" }}{PARA 0 "" 0 "" {TEXT -1 36 "of only nops(v)-n factors in Z_p[x]." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 70 "In other words: as soon as n > nops(v )/2 we know that whatever is left" }}{PARA 0 "" 0 "" {TEXT -1 72 "(tha t's called ff in the algorithm) must be irreducible, so we can stop." }}{PARA 0 "" 0 "" {TEXT -1 73 "In essence, the idea is the following: \+ whenever a subset S of v has been" }}{PARA 0 "" 0 "" {TEXT -1 85 "tri ed and did not lead to a factor in Q[x], then it's pointless to try th e complement" }}{PARA 0 "" 0 "" {TEXT -1 69 "of that set S. This way w e can save at most half the number of cases." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 73 "However, there are other improv ements possible that are more significant." }}{PARA 0 "" 0 "" {TEXT -1 81 "One of those is the following: Don't try just one prime p, but \+ try several primes" }}{PARA 0 "" 0 "" {TEXT -1 82 "p, and then take th e one that leads to the fewest number of p-adic factors. Trying" }} {PARA 0 "" 0 "" {TEXT -1 78 "a couple of p's can decrease nops(v), and because the cost of the algorithm is" }}{PARA 0 "" 0 "" {TEXT -1 79 " exponential in nops(v), even a small decrease in nops(v) can make a su bstantial" }}{PARA 0 "" 0 "" {TEXT -1 11 "difference." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 82 "Another improvement: w hen S is a subset of the set v of monic irreducible factors," }}{PARA 0 "" 0 "" {TEXT -1 60 "then it is relatively expensive to compute the \+ product of S." }}{PARA 0 "" 0 "" {TEXT -1 78 "The idea to speed it up \+ is the following: Just compute one coefficient of that" }}{PARA 0 "" 0 "" {TEXT -1 84 "product, see if it satisfies a certain bound, and on ly if that is the case, then you" }}{PARA 0 "" 0 "" {TEXT -1 76 "compu te the product g of S. That way, whenever g is not a factor in Q[x], i n" }}{PARA 0 "" 0 "" {TEXT -1 72 "most cases you can avoid computing g . It is possible design an extremely" }}{PARA 0 "" 0 "" {TEXT -1 93 "f ast test in this way that skips almost all of the cases that don't lea d to a factor in Q[x]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 76 "When this test is implemented, then in most factoriz ation problems, the most" }}{PARA 0 "" 0 "" {TEXT -1 50 "time consumin g portion is the Hensel lifting part." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 80 "However, when nops(v) is very large, \+ then the number of subsets is exponentially" }}{PARA 0 "" 0 "" {TEXT -1 76 "large. Even though with a clever test we can discard each subse ts S within a" }}{PARA 0 "" 0 "" {TEXT -1 78 "very short time, the fac t that there are so many subsets S makes the algorithm" }}{PARA 0 "" 0 "" {TEXT -1 57 "exponentially slow, because the computation time wil l be:" }}{PARA 0 "" 0 "" {TEXT -1 39 " \" time to do 1 test\" * \"nu mber of S\"" }}{PARA 0 "" 0 "" {TEXT -1 63 "which is: \"a very small \+ constant\" * \"an exponential function\"." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 79 "Zassenhaus algorithm will be fast \+ in practise because usually nops(v) is not so" }}{PARA 0 "" 0 "" {TEXT -1 85 "large. But at the same time it can be exponentially slow, it is possible to construct" }}{PARA 0 "" 0 "" {TEXT -1 71 "examples \+ where nops(v) is large. And a \"small constant\" (small when you" }} {PARA 0 "" 0 "" {TEXT -1 83 "implement a clever test) times an exponen tial function, that's still an exponential" }}{PARA 0 "" 0 "" {TEXT -1 79 "function, so when nops(v) is large enough then sooner or later \+ that has to slow" }}{PARA 0 "" 0 "" {TEXT -1 75 "down the computation. Even though Zassenhaus algorithm dates back to 1969," }}{PARA 0 "" 0 "" {TEXT -1 77 "this combinatorial problem (which elements of v to c ombine) has only recently" }}{PARA 0 "" 0 "" {TEXT -1 12 "been solved. " }}{PARA 0 "" 0 "" {TEXT -1 78 "An example where nops(v) is large is \+ mentioned on page 439 in the book. In the" }}{PARA 0 "" 0 "" {TEXT -1 86 "last sentence on that page they mention a polynomial that took 2 d ays to factor in NTL" }}{PARA 0 "" 0 "" {TEXT -1 85 "which has the fas test factoring implementation. But the book is from 1999, before the" }}{PARA 0 "" 0 "" {TEXT -1 85 "combinatorial problem was solved. Right now, on the same computer, that factorization" }}{PARA 0 "" 0 "" {TEXT -1 68 "takes only 1 second. We'll see later in this course how t hat works. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 75 "You may notice that Maple's \"factor\" command is faster than t he \"my_factor\"" }}{PARA 0 "" 0 "" {TEXT -1 86 "command in this works heet. The reason for that is the following: In our algorithm, the" }} {PARA 0 "" 0 "" {TEXT -1 87 "datastructure we use is: Maple's polynomi als with integer coefficients. Now in computer" }}{PARA 0 "" 0 "" {TEXT -1 84 "algebra, integers are fairly complicated expressions beca use they can have arbitrary" }}{PARA 0 "" 0 "" {TEXT -1 81 "length. Ho wever, as soon as you work modulo m (where m is a power of p), then th e" }}{PARA 0 "" 0 "" {TEXT -1 51 "expressions (the integers mod m) hav e fixed length." }}{PARA 0 "" 0 "" {TEXT -1 84 "One can obtain an effi ciency advantage from that. What Maple did is implement a fast" }} {PARA 0 "" 0 "" {TEXT -1 91 "datastructure for polynomials with coeffi cients mod m. Within this datastructure everything" }}{PARA 0 "" 0 "" {TEXT -1 93 "runs faster, a multiplication of polynomials is perhaps 5 0 times faster in that datastructure" }}{PARA 0 "" 0 "" {TEXT -1 90 "t han in the datastructure we are using. So if we did exactly the same c omputational steps," }}{PARA 0 "" 0 "" {TEXT -1 91 "then we would be a round 50 times slower (that number 50 may vary in different examples, \+ I'm" }}{PARA 0 "" 0 "" {TEXT -1 37 "basically just guessing that numbe r)." }}{PARA 0 "" 0 "" {TEXT -1 85 "So in essence, Maple's factor comm and does the same (perhaps with some minor changes)" }}{PARA 0 "" 0 " " {TEXT -1 90 "as what we're doing in this worksheet, however, each op eration + or * goes 50 times faster" }}{PARA 0 "" 0 "" {TEXT -1 85 "du e to the fact that we're using Maple's \"general integer polynomials\" and Maple uses" }}{PARA 0 "" 0 "" {TEXT -1 91 "\"specific univariate \+ modular polynomials\". So Maple's factor perform the same computationa l" }}{PARA 0 "" 0 "" {TEXT -1 95 "steps as we do, but each such step + or * is much faster, so the overall computation is faster." }}{PARA 0 "" 0 "" {TEXT -1 89 "Note that if you program in C instead of Maple, then you can make it even faster. You can" }}{PARA 0 "" 0 "" {TEXT -1 88 "download software from the web that runs yet another 50 times f aster than Maple's factor" }}{PARA 0 "" 0 "" {TEXT -1 91 "command. Mak ing the software run so fast is hard. If you want to get some idea abo ut what's" }}{PARA 0 "" 0 "" {TEXT -1 93 "involved in optimizing softw are you should read section 9.7 in the book. In this course we're" }} {PARA 0 "" 0 "" {TEXT -1 84 "more concerned with the mathematical aspe cts of the algorithms and not so much as to" }}{PARA 0 "" 0 "" {TEXT -1 37 "how to fully optimize the algorithms." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 84 "Improvement 1, try several primes and choose one where th ere are few p-adic factors." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 233 "Now, instead of just taking some prime, we try 5 primes, and take the one where we have the smallest number of p-adi c factors. That will lead to fewer p-adic factors, which greatly reduc es the number of combinations we need to check." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 261 "monic_irred ucible_p_adic_factors:=proc(f,x)\n local p,l,n,v,ppower;\n l := lcoe ff(f,x);\n p,v := try_some_primes(f,x);\n n := required_accuracy(f,l ,x,p);\n ppower := p^n;\n v := HenselLift(v, f/l, x, n, p);\n v := \+ \{op(v)\}; # Turn v into a set.\n v, ppower\nend:\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 210 "In this algorithm, instead of taking the first usable prime, we take the first 5 usable primes and then we choose th e best one, the one where we have the fewest irreducible factors in Fp [x] and hence in Z_p[x]." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 604 "try_some_primes := proc(f,x)\nlocal p,l,number_factors,v,i,j,best _prime_so_far,factors_so_far; \n p := 2;\n l := lcoeff(f,x);\n num ber_factors := infinity;\n for i to 5 do\n while l mod p = 0 or Di scrim(f,x) mod p = 0 do\n p := nextprime(p)\n od;\n v := \+ Factors(f/l) mod p;\n # lprint( nops(v[2]) );\n if nops(v[2]) < \+ number_factors then\n # This prime is better than what we had b efore:\n best_prime_so_far := p;\n number_factors := nop s(v[2]);\n factors_so_far := [seq(j[1],j=v[2])];\n fi;\n \+ p := nextprime(p);\n od;\n best_prime_so_far, factors_so_far\nend: " }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 9 "Examples." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 175 "In the following \+ example we see that improvement 1 gives a noticable speedup (about a f actor 2 or 3 on my computer). So trying a few primes instead of one pr ime is worthwhile." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "T := \+ time(): my_factor(x^52-1): computation_time = time()-T;" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#/%1computation_timeG$\"%:8!\"$" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 30 "Now lets try a harder example:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 287 "f:=x^54-12*x^51+273*x^48+55172*x^4 5+636717*x^42-69118164*x^39+626174940*x^36+9955088304*x^33+16113952082 4*x^30+82525678888*x^27+147157350765*x^24-7164210121644*x^21-139371164 83437*x^18-39511851689016*x^15-17611233093267*x^12-12897223914040*x^9- 77448210738*x^6-47689721112*x^3-5649262541:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "T := time(): my_factor(f); computation_time = time()-T;" }}{PARA 12 "" 1 " " {XPPMATH 20 "6#*&,6*$)%\"xG\"#F\"\"\"\"\"\"*$)F'\"#:F)\"''[u\"*$)F' \"#=F)\"&3k\"*$)F'\"#7F)!'*)\\R*$)F'\"\"$F)!(sWx#*$)F'\"\"'F)!(,*zn*$) F'\"\"*F)!':'G%!(`8d#F**$)F'\"#CF)!#L*$)F'\"#@F)\"%D#F*F*" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#/%1computation_timeG$\"&pg\"!\"$" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 " " {TEXT -1 68 "Improvement 2a, test one coefficient of g before comput ing all of g." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 100 "The \"rootbound\" algorithm calculates a simple upper \+ bound for the absolute values of the roots of f." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "readlib(rootbound): # not necessary in Maple \+ 6, 7, .." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 373 "In the following algorithm, when we have a set S of fact ors in Z_p[x], before we compute all coefficients of the product g, we just compute one coefficient of g, and then compare that coefficient \+ to a bound. If the coefficient is larger than the bound, then \"test_s ome_coefficients\" returns \"true\" and then we can skip S without cal culating the remaining coefficients of g." }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 948 "factor_sqrfree := proc(f,x)\n local c,ff,v,ppower ,result,n,l,su,S,g,B,tc;\n if degree(f,x) < 2 then\n RETURN( \+ f )\n fi;\n\n c := icontent(f);\n ff := f/c;\n v, ppower := monic_ irreducible_p_adic_factors(f,x);\n result := c; \n n := 0;\n B := r ootbound(f,x);\n while v <> \{\} do\n l := lcoeff(ff,x);\n \+ tc := coeff(ff,x,0);\n n := n+1;\n su := [subsets(v,n)];\n \+ # Go through all subsets of v with n elements:\n for S in su \+ do\n if n>1 and test_some_coefficients(S, x, ppower, l, tc, \+ B) then\n next\n fi;\n g := convert (S,`*`);\n g := mods(l * Expand(g),ppower);\n g := g/icontent(g);\n if divide(ff,g,'q') then\n re sult := result * g;\n ff := q; # ff = ff/g\n \+ v := v minus S\n fi\n od\n od;\n if not member(ff, \+ \{1,-1\}) then\n ERROR(\"All factors should have been removed by n ow\")\n fi;\n ff*result\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 345 "Here we'll implement the \"d-1 test\". That means that if g has d egree d, then we calculate coeff(g,x,d-1). Now if g is monic, then tha t coefficient is (-1) times the sum of the roots. So it is bounded by: degree * rootbound. If g is not monic, but has leading coefficient l , then abs( coeff(g,x,d-1) ) must be no greater than l*degree*rootboun d." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 352 "test_some_coefficien ts := proc(S, x, ppower, l, tc, RootBound)\n local bd, l_times_coeff, i;\n # Bound for the second highest coefficient is:\n bd := RootBoun d * add(degree(i,x),i=S);\n # l times the actual second highest coeff icient is:\n l_times_coeff := mods(l *\n add(coeff(i,x,degree(i,x )-1), i=S), ppower);\n evalb( l_times_coeff > l * bd)\nend:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 " " {TEXT -1 23 "Try last example again." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 34 "Now try the same polynomial again:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "T := time(): my_factor(f): computation_time = time()- T;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%1computation_timeG$\"&r<\"!\"$ " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 48 "This example runs faster now \+ than it did before." }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 43 "Improvement 2b, test two coefficients of \+ g." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {EXCHG {PARA 0 "" 0 "" {TEXT -1 386 "Now do the \"d-1 test\" as well a s the \"constant term test\". If g is the product of the elements of S , then coeff(g,x,0) is easily calculated: just take the product of the constant coefficients of the elements of S. Now if g is primitive (th e gcd of the coeffs is 1) then coeff(g,x,0) must divide coeff(f,x,0) i f g divides f. So abs(coeff(g,x,0)) can not be bigger than abs(coeff(f ,x,0))." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 567 "test_some_coeff icients := proc(S, x, ppower, l, tc, RootBound)\n local bd, l_times_c oeff,i;\n # Bound for the second highest coefficient is:\n bd := Roo tBound * add(degree(i,x),i=S);\n # l times the actual second highest \+ coefficient is:\n l_times_coeff := mods(l *\n add(coeff(i,x,degre e(i,x)-1), i=S), ppower);\n if l_times_coeff > l * bd then\n RETU RN(true)\n fi;\n i := mods(l * mul(coeff(i,x,0), i=S),ppower):\n ev alb(i <> 0 and irem(l*tc, i) <> 0)\n # So if lc*i is not divisible by i (when irem<>0)\n # then we return \"true\" which means: skip this \+ S.\nend:" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 23 "Try last example ag ain." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "T := time(): my_factor(f): computation_time = time()-T;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%1co mputation_timeG$\"%gj!\"$" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 101 "If \+ we have a subset S whose product is g, and if the degree of g is d, th en we can quickly calculate:" }}{PARA 0 "" 0 "" {TEXT -1 36 " coeff (g,x,d-1) and coeff(g,x,0)." }}{PARA 0 "" 0 "" {TEXT -1 112 "We can co mpare those numbers to some a-priori bounds, and if they are bigger th an the bound, then we can skip S." }}{PARA 0 "" 0 "" {TEXT -1 111 "In \+ that case we save something, namely we save the computation of all oth er coefficients of the product g of S." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}}{MARK "1" 0 }{VIEWOPTS 1 1 0 1 1 1803 }