{VERSION 6 0 "SUN SPARC SOLARIS" "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 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 268 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 272 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 273 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 274 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 275 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 276 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 277 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 278 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 279 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 280 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 281 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 282 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 283 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 284 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 285 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 286 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 287 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 288 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 289 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 290 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 291 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 292 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 293 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 294 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 295 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 296 "" 0 1 0 0 0 0 1 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 "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }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 1 }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 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 {SECT 0 {PARA 3 "" 0 "" {TEXT -1 0 "" }{TEXT 256 55 "Factoring in Q[x] with Zassenhaus' algorithm, overview." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 39 "Let Q be the field of ra tional numbers." }}{PARA 0 "" 0 "" {TEXT -1 30 "Let Z be the ring of i ntegers." }}{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 66 " Let g be the product of that s ubset. We then know g not" }}{PARA 0 "" 0 "" {TEXT -1 55 " p recisely, 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 426 "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 f actor_univariate_case(f, x[1])\n else\n factor_multivariate_case( f, x)\n fi:\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 125 "fact or_multivariate_case := proc()\n error \"more than one variable no t yet implemented\"\n # We will do this later.\nend:" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 23 "Hensel lift algorithms." }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 126 "For monic polyn omials, the bound on the size of the coefficients of a factor of f is: 2^n times the biggest coefficient of f." }}{PARA 0 "" 0 "" {TEXT -1 201 "Note that we multiply the bound by l because we're multiplying th e product g in the procedure factor_sqrfree by l. That way the rationa l number coefficients become integer coefficients by Gauss' lemma." }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 134 "required_accuracy := proc( f,l,x,p)\n local B,n;\n n := degree(f,x);\n B := l * 2^n*maxnorm(f) ;\n ceil( evalf(log(2*B+1)/log(p)) )\nend:" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 182 "The following algorithm is there so we can choose whic h Hensel-algorithm we'll use. That way we can compare them more easily , and we can then test experimentally which one is faster." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 287 "HenselLift := proc()\n # To chang e the algorithm that is used,\n # just uncomment (remove #) from one \+ line\n # and comment (place #) the other lines:\n\n # HenselLift_lin alg(args)\n # HenselLift_linalg_biggersteps(args)\n # HenselLift_gcd ex(args)\n HenselLift_gcdex_biggersteps(args)\nend:" }}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 54 "Hensel lifting using the extended Euclide an Algorithm." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 942 "Hensel_lif t_two_factors := proc(f,g,h,x,N,p)\n # Input: a positive integer n, \+ and\n # monic polynomials f,g,h in Z[x] such that:\n # h is f*g mo d 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 then\n error \"f and g sho uld 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 from 1 to N-1 do\n # Compu te 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 954 "HenselLift_gcdex:=proc(L::list,h,x ::name,N::posint,p::posint)\n # Recursive algorithm for Hensel lifti ng 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.\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 # Output: a li st [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 [Hense l_lift_two_factors(op(L),h,x,N,p)]\n else\n L1, L2 := L[1..iqu o(l,2)], L[iquo(l,2)+1..-1];\n # Take the products.\n h1, \+ h2 := convert(L1,`*`), convert(L2,`*`);\n # Lift the products:\n h1, h2 := Hensel_lift_two_factors(h1,h2,h,x,N,p);\n # Lif t 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 0 {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 mo d p^1\n # (so the elements in L have accuracy 1).\n # The elements in L must have gcd 1 in Fp[y].\n #\n # n a positive integer.\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 (as polyn omials in y).\n #\n local i,j,k,c,Li,pol,cpol;\n\n Li:=L; # accu racy 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 # something th at is 0 mod p^i.\n #\n # Li mod p^(i+1) did change. We did n ot 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 # the coeffi cient of p^i*y^j in the k'th factor. \n #\n # Now lets compu te 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 following \+ 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 means that pol is 0 mod x^1.\n pol:=Expand(pol) mod p;\n # This pol sh ould be the zero polynomial,\n # which means that the following c oeffs\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,Li)\n od;\n Li\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 90 "Hensel lifting by sol ving linear equations, increasing the accuracy more than 1 at a time. " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1627 "HenselLift_linalg_bigg ersteps:=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 choose 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're only going to accuracy n, we see that \n # there is no reason for inc_i to be bigger than n-i.\n # So inc_i <= n-i.\n\n # However, we want that the equations for t he\n # unknowns c[j,k] to be linear: If they are not\n # lin ear then they are difficult to solve.\n\n # Now \"pol\" is not li near in the c[j,k], it contains\n # products of the unknowns c[j, k]. However, these\n # products 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 o n inc_i. We must have:\n # inc_i <= i. Otherwise in the equations \"cpol\"\n # below will not be linear.\n\n # So we take som ething that's <=i and <= n-i.\n inc_i := min(i, n-i);\n \n \+ pol:=Expand(pol) mod p^inc_i;\n\n # This pol should be the z ero polynomial,\n # which means that the following coeffs\n \+ # should be zero:\n cpol:=\{coeffs(pol,y)\};\n # Solve and i ncrease the accuracy of Li to i+1:\n Li:=subs( solve(cpol) mod p^ inc_i, Li);\n i := i + inc_i\n od;\n Li\nend:" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 74 "Hensel lift ing, extended Euclidean Algorithm, more than 1 step at at time." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 572 "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 Gcd 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 wh ile 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 490 "HenselLift_gcde x_biggersteps:=proc(L::list,h,x::name,N::posint,p::posint)\n local l ,L1,L2,h1,h2;\n l:=nops(L);\n if l=0 then\n error \"wrong in put\"\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)), op(procname(L2,h2,x,N,p))]\n fi\nen d:" }}}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 71 "Compute all monic irreducible factors in Z_p[x] for a sui table prime p." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1537 "monic_ir reducible_p_adic_factors:=proc(f,x)\n local p,l,fm,n,v,ppower,i;\n l := lcoeff(f,x);\n\n p := take_the_first_suitable_prime(f,l,x);\n \n \+ # Now we factor f mod p, but lets first make it monic:\n # We know t hat Maple can calculate \"1/l mod p\" by\n # computing igcdex(l,p,'s' ,'t'), and we know how\n # that algorithm works. So we know how Maple can\n # compute the following: first \"1/l mod p\" (which is\n # th e number s in igcdex(l,p,'s','t')) and then\n # multiply that by f, a nd then take the remainder\n # of the coefficients mod p.\n fm := f/ l mod p; \n\n # Now we have a monic polynomial fm in Fp[x] where Fp \n # is a finite field with p elements. Furthermore, fm\n # has no m ultiple 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 thi s leading coefficient is lc=1\n # and all multiplicities m1, m2,.. ar e 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_accuracy(f,l,x,p);\n ppower := p^n ;\n v := HenselLift(v, f/l, x, n, p);\n v := \{op(v)\}; # Turn v int o a set.\n # We've factored f/l in Z_p[x] where Z_p = the ring\n # o f p-adic integers. We've determined the factors\n # up to accuracy n, i.e. we've computed them\n # modulo ppower.\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 pr ime." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 71 "I f 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 "w ould appear more than once in the factorization. That's a problem in" }}{PARA 0 "" 0 "" {TEXT -1 71 "Hensel lifting because we can only Hens el lift factors that have gcd 1," }}{PARA 0 "" 0 "" {TEXT -1 80 "so th e irreducible factors must all be distinct (they must have multiplicit y 1)." }}{PARA 0 "" 0 "" {TEXT -1 71 "Now the discriminant is somethin g that is zero if and only if there are" }}{PARA 0 "" 0 "" {TEXT -1 80 "factors of multiplicity > 1, which is true if and only if f and f' have a common" }}{PARA 0 "" 0 "" {TEXT -1 88 "factor, which is true i f and only if the resultant of f and f' is zero (the discriminant" }} {PARA 0 "" 0 "" {TEXT -1 53 "of f is some number times the resultant o f f and f')." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 153 "take_the_first_suitable_prime := proc(f,l,x)\n \+ local p;\n p := 2;\n while l mod p = 0 or Discrim(f,x) mod p = 0 do \n p := nextprime(p)\n od;\n p\nend:" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 82 "Factoring univariate p olynomials, reducing the problem to square-free polynomials." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 994 "factor_univariate_case := p roc(f,x)\n local v,result,i,w;\n\n # Determine a square-free factori zation. We know by\n # now of course how that works. It's mainly base d on\n # gcd(f, f') and this computation is fast because we\n # use \+ modular arithmetic which prevents excessive\n # growth of coefficient s, and if the coefficients don't\n # grow then the Euclidean algorith m 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 multipli cities.\n\n for i in v[2] do\n # Now i[1] is a square-free polyn omial. So\n # gcd( i[1], i[1]' ) = 1 so the sqrfree\n # algo rithm can not factor i[1] any further.\n # To factor a squarefree polynomial i[1] we'll\n # use the 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 power because i[1] appeared\n # \+ in f with multiplicity i[2].\n result := result * w^i[2]\n od;\n result\nend:\n" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 45 "Factoring univariate square-free polynomials." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1477 "factor_sqrfree := proc(f,x )\n # Input: a squarefree polynomial f in Q[x]\n # Output: The facto rization.\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 ha ve degree\n # at least 2, so we're done:\n return f\n f i;\n\n c := icontent(f);\n ff := f/c;\n # Since we divided by the i nteger-content c, we\n # now have a polynomial with integer coefficie nts.\n\n v, ppower := monic_irreducible_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 while v <> \{\} do\n l := lcoeff( ff,x);\n n := n+1;\n su := [subsets(v,n)];\n # Go throu gh 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(f f, \{1,-1\}) then\n error \"All factors should have been removed b y now\"\n fi;\n ff*result\nend:\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 75 "Compute the sequence S1,S2,.. of all subsets with n elements of some set v." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 263 "subsets:=p roc(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 0 {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. Note: it's best to do that with sm all 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_fact or(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)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "my_factor(16*x^2-1);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#*&,&*&\"\"%\"\"\"%\"xGF'F'F'!\"\"F',&*&F&F'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**$)F(\"\"$F*F*F&F*F(F*F*F*F)F*),&F(F*F*!\"\"F)F*),0*$)F(\" \")F*F**$)F(\"\"(F*F5*$)F(\"\"&F*F*F-F5F0F*F(F5F*F*F)F*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "expand(%);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#,(*$)%\"xG\"#I\"\"\"F(*&\"\"#F()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), " 0 "" {MPLTEXT 1 0 0 " " }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 25 "Remarks and improvements." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 62 "There a re several ways in which the algorithm can be improved:" }}{PARA 0 "" 0 "" {TEXT -1 55 "First of all: At the bottom of factor_sqrfree you s ee:" }}{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 through all subsets of v with n elements:\n for S in su d o\n" }}{PARA 0 "" 0 "" {TEXT -1 78 "Here v is the list of monic factor s 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 polynomia l ff in the algorithm)." }}{PARA 0 "" 0 "" {TEXT -1 72 "The improvemen t we can 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 facto r g of ff in Q[x] consisting of n" }}{PARA 0 "" 0 "" {TEXT -1 69 "fact ors 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 factors. 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 e arlier 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 a s n > nops(v)/2 we know that whatever is left" }}{PARA 0 "" 0 "" {TEXT -1 72 "(that's called ff in the algorithm) must be irreducible, \+ so we can stop." }}{PARA 0 "" 0 "" {TEXT -1 73 "In essence, the idea i s the following: whenever a subset S of v has been" }}{PARA 0 "" 0 " " {TEXT -1 85 "tried and did not lead to a factor in Q[x], then it's p ointless to try the complement" }}{PARA 0 "" 0 "" {TEXT -1 69 "of that set S. This way we 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 improvements 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 the one that leads to the fewest number of p-adic factors. Trying" }}{PARA 0 "" 0 "" {TEXT -1 78 "a couple of p's can d ecrease 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 n ops(v) can make a substantial" }}{PARA 0 "" 0 "" {TEXT -1 11 "differen ce." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 82 "An other improvement: when S is a subset of the set v of monic irreducibl e factors," }}{PARA 0 "" 0 "" {TEXT -1 60 "then it is relatively expen sive 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 ce rtain bound, and only if that is the case, then you" }}{PARA 0 "" 0 " " {TEXT -1 76 "compute the product g of S. That way, whenever g is not a factor in Q[x], in" }}{PARA 0 "" 0 "" {TEXT -1 72 "most cases you c an avoid computing g. It is possible design an extremely" }}{PARA 0 " " 0 "" {TEXT -1 93 "fast test in this way that skips almost all of the cases that don't lead 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 factorization problems, the most" }}{PARA 0 "" 0 "" {TEXT -1 50 "time consuming 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 exponentiall y" }}{PARA 0 "" 0 "" {TEXT -1 76 "large. Even though with a clever tes t we can discard each subsets S within a" }}{PARA 0 "" 0 "" {TEXT -1 78 "very short time, the fact that there are so many subsets S makes t he algorithm" }}{PARA 0 "" 0 "" {TEXT -1 57 "exponentially slow, becau se the computation time will be:" }}{PARA 0 "" 0 "" {TEXT -1 39 " \" time to do 1 test\" * \"number 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 "Zassenh aus 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 te st) times an exponential 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 (whi ch elements of v to combine) has only recently" }}{PARA 0 "" 0 "" {TEXT -1 12 "been solved." }}{PARA 0 "" 0 "" {TEXT -1 78 "There are ex amples where nops(v) is large, and such examples used to run slow." }} {PARA 0 "" 0 "" {TEXT -1 82 "However, nowadays there is a fast alterna tive for the case where nops(v) is large," }}{PARA 0 "" 0 "" {TEXT -1 36 "ask me if you're interested in that." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 75 "You may notice that Maple's \"fact or\" command is faster than the \"my_factor\"" }}{PARA 0 "" 0 "" {TEXT -1 86 "command in this worksheet. The reason for that is the fol lowing: In our algorithm, the" }}{PARA 0 "" 0 "" {TEXT -1 87 "datastru cture we use is: Maple's polynomials with integer coefficients. Now in computer" }}{PARA 0 "" 0 "" {TEXT -1 84 "algebra, integers are fairly complicated expressions because they can have arbitrary" }}{PARA 0 " " 0 "" {TEXT -1 81 "length. However, as soon as you work modulo m (whe re m is a power of p), then the" }}{PARA 0 "" 0 "" {TEXT -1 51 "expres sions (the integers mod m) have fixed length." }}{PARA 0 "" 0 "" {TEXT -1 84 "One can obtain an efficiency advantage from that. What Ma ple did is implement a fast" }}{PARA 0 "" 0 "" {TEXT -1 91 "datastruct ure for polynomials with coefficients mod m. Within this datastructure everything" }}{PARA 0 "" 0 "" {TEXT -1 93 "runs faster, a multiplicat ion of polynomials is perhaps 50 times faster in that datastructure" } }{PARA 0 "" 0 "" {TEXT -1 90 "than in the datastructure we are using. \+ So if we did exactly the same computational steps," }}{PARA 0 "" 0 "" {TEXT -1 87 "then we would be around 50 times slower (that number 50 m ay vary in different examples," }}{PARA 0 "" 0 "" {TEXT -1 41 "I'm bas ically just guessing that number) " }}{PARA 0 "" 0 "" {TEXT -1 85 "So \+ in essence, Maple's factor command does the same (perhaps with some mi nor changes)" }}{PARA 0 "" 0 "" {TEXT -1 90 "as what we're doing in th is worksheet, however, each operation + or * goes 50 times faster" }} {PARA 0 "" 0 "" {TEXT -1 85 "due to the fact that we're using Maple's \+ \"general integer polynomials\" and Maple uses" }}{PARA 0 "" 0 "" {TEXT -1 92 "\"specific univariate modular polynomials\". So Maple's f actor performs the same computational" }}{PARA 0 "" 0 "" {TEXT -1 95 " steps as we do, but each such step + or * is much faster, so the overa ll 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 faster than Maple's factor" }}{PARA 0 "" 0 "" {TEXT -1 91 "command. Making the software run so fast is hard. If you want to get some idea about what's" }}{PARA 0 "" 0 "" {TEXT -1 87 "involved in optimizing software you should read section 9.7 in \+ the book Modern Computer" }}{PARA 0 "" 0 "" {TEXT -1 83 "Algebra by Ge rhard and von zur Gathen. In this course we're more concerned with the " }}{PARA 0 "" 0 "" {TEXT -1 98 "mathematical aspects of the algorithm s and not so much as to how to fully optimize the algorithms." }}} {SECT 0 {PARA 3 "" 0 "" {TEXT -1 84 "Improvement 1, try several primes and choose one where there are few p-adic factors." }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 233 "Now, instead of just t aking some prime, we try 5 primes, and take the one where we have the \+ smallest number of p-adic factors. That will lead to fewer p-adic fact ors, which greatly reduces the number of combinations we need to check ." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 261 "monic_irreducible_p_adic_factors:=proc(f,x)\n local p,l,n,v,ppower;\n l := lcoeff(f,x);\n p,v := try_some_primes(f,x); \n n := required_accuracy(f,l,x,p);\n ppower := p^n;\n v := HenselL ift(v, f/l, x, n, p);\n v := \{op(v)\}; # Turn v into a set.\n v, pp ower\nend:\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 210 "In this algorith m, instead of taking the first usable prime, we take the first 5 usabl e primes and then we choose the best one, the one where we have the fe west 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)\nloca l p,l,number_factors,v,i,j,best_prime_so_far,factors_so_far; \n p := 2;\n l := lcoeff(f,x);\n number_factors := infinity;\n for i to 5 \+ do\n while l mod p = 0 or Discrim(f,x) mod p = 0 do\n p := n extprime(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 pri me is better than what we had before:\n best_prime_so_far := p; \n number_factors := nops(v[2]);\n factors_so_far := [se q(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 0 {PARA 3 "" 0 "" {TEXT -1 9 "Examples." }}{EXCHG {PARA 0 "" 0 " " {TEXT -1 175 "In the following example we see that improvement 1 giv es a noticable speedup (about a factor 2 or 3 on my computer). So tryi ng a few primes instead of one prime is worthwhile." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "T := time(): my_factor(x^52-1): computati on_time = time()-T;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%1computation_ timeG$\"$!))!\"$" }}}{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^45+636717*x^42-69118164*x^39+626174940*x^3 6+9955088304*x^33+161139520824*x^30+82525678888*x^27+147157350765*x^24 -7164210121644*x^21-13937116483437*x^18-39511851689016*x^15-1761123309 3267*x^12-12897223914040*x^9-77448210738*x^6-47689721112*x^3-564926254 1:" }}}{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\"%(>#\"\"\"*&\"(L.Z \"F&)%\"xG\"#7F&F&*&\"&wh\"F&)F*\"\"$F&F&*&\"%toF&)F*\"\"'F&F&*&\"(,`' \\F&)F*\"\"*F&F&*&\"'uToF&)F*\"#:F&F&*&\"&3D#F&)F*\"#=F&!\"\"*&\"#@F&) F*\"#CF&F&*&\"$f(F&)F*FBF&F@*$)F*\"#FF&F&F&,6FHF&*&\"#LF&FCF&F@*&\"%D< F&FGF&F&*&\"&3k\"F&F>F&F&*&\"''[u\"F&F:F&F&*&\"'*)\\RF&F)F&F@*&\"':'G% F&F6F&F@*&\"(,*znF&F2F&F@*&\"(sWx#F&F.F&F@\"(`8d#F@F&" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%1computation_timeG$\"%+o!\"$" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 68 "I mprovement 2a, test one coefficient of g before computing all of g." } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 100 " The \"rootbound\" algorithm calculates a simple upper bound for the ab solute values of the roots of f." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 373 "In the following algorithm, when \+ we have a set S of factors in Z_p[x], before we compute all coefficien ts of the product g, we just compute one coefficient of g, and then co mpare that coefficient to a bound. If the coefficient is larger than t he bound, then \"test_some_coefficients\" returns \"true\" and then we can skip S without calculating the remaining coefficients of g." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 944 "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 := rootbound(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),ppow er);\n g := g/icontent(g);\n if divide(ff,g,'q') t hen\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 \"All factors should ha ve been removed by now\"\n fi;\n ff*result\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 345 "Here we'll implement the \"d-1 test\". That mean s that if g has degree d, then we calculate coeff(g,x,d-1). Now if g i s monic, then that coefficient is (-1) times the sum of the roots. So \+ it is bounded by: degree * rootbound. If g is not monic, but has lead ing coefficient l, then abs( coeff(g,x,d-1) ) must be no greater than \+ l*degree*rootbound." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 352 "tes t_some_coefficients := proc(S, x, ppower, l, tc, RootBound)\n local b d, l_times_coeff,i;\n # Bound for the second highest coefficient is: \n bd := RootBound * add(degree(i,x),i=S);\n # l times the actual se cond highest coefficient is:\n l_times_coeff := mods(l *\n add(co eff(i,x,degree(i,x)-1), i=S), ppower);\n evalb( l_times_coeff > l * b d)\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 23 "Try last example again." }}{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$\"%Sa!\"$" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 48 "This example runs faster now than it did before." }} }{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 43 "Improvement 2b, test two coefficients of g." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 386 "Now do the \"d-1 test\" as well as the \"constant term test\". If g is the product of the ele ments of S, then coeff(g,x,0) is easily calculated: just take the prod uct of the constant coefficients of the elements of S. Now if g is pri mitive (the gcd of the coeffs is 1) then coeff(g,x,0) must divide coef f(f,x,0) if g divides f. So abs(coeff(g,x,0)) can not be bigger than a bs(coeff(f,x,0))." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 566 "test_ some_coefficients := proc(S, x, ppower, l, tc, RootBound)\n local bd, l_times_coeff,i;\n # Bound for the second highest coefficient is:\n \+ bd := RootBound * add(degree(i,x),i=S);\n # l times the actual secon d highest coefficient is:\n l_times_coeff := mods(l *\n add(coeff (i,x,degree(i,x)-1), i=S), ppower);\n if l_times_coeff > l * bd then \n return true\n fi;\n i := mods(l * mul(coeff(i,x,0), i=S),ppow er):\n evalb(i <> 0 and irem(l*tc, i) <> 0)\n # So if lc*i is not di visible by i (when irem<>0)\n # then we return \"true\" which means: \+ skip this S.\nend:" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 23 "Try last \+ example 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$\"%!>$!\"$" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 101 "If we have a subset S whose product is g, and if th e degree of g is d, then 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 compare those numbers to some a-priori bounds, an d if they are bigger than the bound, then we can skip S." }}{PARA 0 " " 0 "" {TEXT -1 111 "In that case we save something, namely we save th e computation of all other coefficients of the product g of S." }}}}} {MARK "0 1 0" 1 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }