# MAS3008 -- SOLUTIONS TO MAPLE EXERCISES # # An executable MAPLE worksheet containing implementations of the # various algorithms can be downloaded from the Module Web page. # Alternatively, you might like to write your own procedures, as below > # EXERCISE 1 > # (a) > p:=nextprime(100000); p := 100003 > 2^(p-1) mod p; 1 > 2&^(p-1) mod p; 1 > # We can measure the CPU time taken (in seconds) using the time() # command: using the 2nd method of taking powers, the time is negligible > time(2^(p-1) mod p);time(2&^(p-1) mod p); .025 0. > # # Now try a much larger prime... > p:=nextprime(1000000); p := 1000003 # > 2^(p-1) mod p;2&^(p-1) mod p; time(2^(p-1) mod p);time(2&^(p-1) mod > p); 1 1 1.175 0. # (b) > powermod:=proc(a,k,n) local b,c,i;c:=a;b:=1;i:=k; > while i>0 do > if (i mod 2=1) then b:=b*c mod n; i:=i-1;fi; > c:=c*c mod n; i:=i/2; > od; > RETURN(b); > end: > > > powermod(29,405,1347); # 233 > p1:=nextprime(2^20);p2:=nextprime(2^30);p3:=nextprime(2^50); p1 := 1048583 p2 := 1073741827 p3 := 1125899906842679 > powermod(3,p1-1,p1);powermod(5,p2-1,p2);powermod(7,p3-1,p3); 1 1 1 > time(powermod(7,p3-1,p3)); .001 # (c) Now try same thing using a simple loop.... # > powermod2:=proc(a,k,n) local b,i; > b:=1; > for i from 1 to k do b:=b*a mod n; od; > RETURN(b); > end: > > > > > powermod2(2,256,257); 1 > p1:=nextprime(2^20); p1 := 1048583 > powermod2(2,p1-1,p1); 1 > time(powermod2(2,p1-1,p1)); 4.863 # This is vastly slower ! > > # Exercise 2. > > n:=2212374139; n := 2212374139 > powermod(2,n-1,n); 1965736235 # This shows n is not prime... > > isprime(n); false > ifactor(n); (9337) (236947) > n:=19707683773; n := 19707683773 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > powermod(11,n-1,n); 1 > powermod(13,n-1,n); 1 # These are consistent with n being prime: if n is not prime then it is # a pseudoprime to all the bases 2, 3, 5, 7, 11, 13. (This is unlikely # unless n was contrived to have this property.) > > isprime(n); true > # Now look for 10-digit numbers passing Fermat test to bases 2,3,5,7: # This can be done by trial and error (or alternatively, you could # programme a loop to search for suitable numbers). > > n:=1234567891; n := 1234567891 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > n1:=n; n1 := 1234567891 > n:=2395806251; n := 2395806251 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > n2:=n; n2 := 2395806251 > n:=2134642231; n := 2134642231 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > n3:=n; n3 := 2134642231 > n:=7232632703; n := 7232632703 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > n4:=n; n4 := 7232632703 > n:=3272052041; n := 3272052041 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > n5:=n; n5 := 3272052041 > n:=2425472839; n := 2425472839 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > n6:=n; n6 := 2425472839 > n:=2637564901; n := 2637564901 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > n7:=n; n7 := 2637564901 > n:=1362385309; n := 1362385309 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > n8:=n; n8 := 1362385309 > n:=3960384521; n := 3960384521 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > n9:=n; n9 := 3960384521 > n:=1524831907; n := 1524831907 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > n10:=n; n10 := 1524831907 > n1;n2;n3;n4;n5;n6;n7;n8;n9;n10; > 1234567891 2395806251 2134642231 7232632703 3272052041 2425472839 2637564901 1362385309 3960384521 1524831907 > isprime(n1); true > isprime(n2); true > isprime(n3); true > isprime(n4); true > isprime(n5); true > isprime(n6); true > isprime(n7); true > isprime(n8); true > isprime(n9); true > isprime(n10); true # We've now found ten 10-digit numbers passing the Fermat test to bases # 2, 3, 5, 7; and in fact all of them are primes. The Fermat test is # fairly good - most numbers that pass to base 2 are prime. But of # course, not all ... # # true > n:=1001152801; n := 1001152801 > powermod(2,n-1,n); 1 > powermod(3,n-1,n); 1 > powermod(5,n-1,n); 1 > powermod(7,n-1,n); 1 > ifactor(n); (11) (41) (61) (151) (241) # For each of the primes p dividing n, we have that p-1 divides n-1. # Hence n is a Carmichael number, and it will pass the Fermat primality # test to any base a with gcd(a,n)=1. It will fail if gcd(a,n)>1, e.g. # if a is any of the primes dividing n: > powermod(11,n-1,n);powermod(41,n-1,n); # # 910138911 781387553 > powermod(13,n-1,n);powermod(43,n-1,n); 1 1 > > # Exercise 3: # Miller-Rabin test for stong pseudoprimes. : # > spsp:=proc(n,a) local g,m,s,b,b1,i; > if a>=n or a<2 then RETURN(`Invalid input`) fi; > if n mod 2=0 then RETURN(`Fail: factor`,2) fi; > g:=gcd(a,n); if g>1 then RETURN(`Fail: factor`,g) fi; > m:=n-1; s:=0; > while m mod 2=0 do m:=m/2; s:=s+1 od; > b:=powermod(a,m,n); > if b=1 then RETURN(Pass) fi; > if b=n-1 then RETURN(Pass) fi; > for i from 1 to s-1 do > b1:=b*b mod n; > if b1=1 then RETURN(`Fail: factor`,gcd(b-1,n)) fi; > if b1=n-1 then RETURN(Pass) fi; > b:=b1; > od; > RETURN(Fail); > end: > # # Test detection of invaild input and special cases ... > spsp(561,1);spsp(561,561); > spsp(56,2); Invalid input Invalid input Fail: factor, 2 # Now try it "for real" .. > spsp(567,7);# Fail: factor, 7 > spsp(561,2); # Try it out on examples done in lectures ... Fail: factor, 33 > spsp(2047,2); Pass # > spsp(2047,3); Fail > n:=1373653; n := 1373653 > spsp(n,2);spsp(n,3);spsp(n,5); Pass Pass Fail # Now try numbers from Question 2 > n:=221237139; spsp(n,2); n := 221237139 Fail > n:=19707683773; spsp(n,2);spsp(n,3);spsp(n,5); n := 19707683773 Pass Pass Pass > # > n:=1001152801; n := 1001152801 > spsp(n,2); Fail: factor, 4154161 > ifactor(4154161); (11) (41) (61) (151) # So the Carmichael number in Question 2 is detected by the Miller-Rabin # test at the first attempt (base 2). Other bases work too.. > spsp(n,3); spsp(n,5); spsp(n,7); Fail: factor, 101321 Fail: factor, 4154161 Fail: factor, 4154161 > > # Question 4: Pollard p-1 method: # # (a) > pollard:=proc(n,a0,kmax) local a,g,k;a:=a0; > for k from 2 to kmax do > a:=powermod(a,k,n); > g:=gcd(a-1,n); > if g=n then RETURN(Fail1); fi; > if g>1 then RETURN(g,n,k) fi; > od; > RETURN(Fail2); > end: > > pollard(403,2,1000);pollard(1891,2,1000);pollard(5157437,2,1000); 13, 403, 4 Fail1 2269, 5157437, 9 > ifactor(1891); (31) (61) # So the 1st and 3rd values of n are easily factorised (using base 2). # For n=1891=31 x 61, both prime factors show up together at step 5, so # we can't factorise n by this method unless we are particularly lucky, # (as happens here for base 5): > pollard(1891,3,100); pollard(1891,5,100); Fail1 31, 1891, 3 # Now try it for the Fermat numbers. As these are closely related to # powers of 2, we wouldn't expect base 2 to work. > > F5:=2^32+1;F6:=2^64+1;F7:=2^128+1; F5 := 4294967297 F6 := 18446744073709551617 F7 := 340282366920938463463374607431768211457 > pollard(F5,2,100); Fail1 # However, other bases should be OK > pollard(F5,3,100); 641, 4294967297, 8 > pollard(F6,3,100); 274177, 18446744073709551617, 17 > pollard(F7,3,1000); Fail2 # So F5 and F6 are easily factorised by this method. F7 is much harder # to factorise, although it is composite: > isprime(F7); false > pollard(F7,5,1000);pollard(F7,11,1000); Fail2 Fail2 > # (b) First identify composite Mersenne numbers Mp for p<100: > for p from 2 to 99 do > if isprime(p) then if not isprime(2^p-1) then print(p);fi;fi;od; 11 23 29 37 41 43 47 53 59 67 71 73 79 83 97 # So there are 15 values of p to consider, viz 11, 23, 29, 37, 41, 43, # 47, 53, 59, 67, 71, 73, 79, 83, 97 > > M11:=2^11-1;M23:=2^23-1;M29:=2^29-1;M37:=2^37-1;M41:=2^41-1;M43:=2^43- > 1; M11 := 2047 M23 := 8388607 M29 := 536870911 M37 := 137438953471 M41 := 2199023255551 M43 := 8796093022207 > M47:=2^47-1;M53:=2^53-1;M59:=2^59-1;M67:=2^67-1;M71:=2^71-1;M73:=2^73- > 1; M47 := 140737488355327 M53 := 9007199254740991 M59 := 576460752303423487 M67 := 147573952589676412927 M71 := 2361183241434822606847 M73 := 9444732965739290427391 > M79:=2^79-1;M83:=2^83-1;M97:=2^97-1; M79 := 604462909807314587353087 M83 := 9671406556917033397649407 M97 := 158456325028528675187087900671 > pollard(M11,2,200); Fail1 # We cannot factor any of these using base 2 for Pollard's p-1 method, # since 2 must have order p for any prime factor q of Mp. # Try base 3 > pollard(M11,3,200);pollard(M23,3,200);pollard(M29,3,200); Fail1 47, 8388607, 23 Fail1 > pollard(M37,3,200);pollard(M41,3,200);pollard(M43,3,200); 223, 137438953471, 37 13367, 2199023255551, 163 431, 8796093022207, 43 > pollard(M47,3,200);pollard(M53,3,200);pollard(M59,3,200); 10610063, 140737488355327, 47 129728784761, 9007199254740991, 53 179951, 576460752303423487, 61 > pollard(M67,3,200);pollard(M71,3,200);pollard(M73,3,200); Fail2 Fail2 439, 9444732965739290427391, 73 > pollard(M79,3,200);pollard(M83,3,200);pollard(M97,3,200); 2687, 604462909807314587353087, 79 167, 9671406556917033397649407, 83 11447, 158456325028528675187087900671, 97 # Using kmax=200 we factorise all cases except p=11, 29, 67, 71. The # latter 2 cases can be handled by increasing kmax: > pollard(M67,3,5000); pollard(M71,3,5000); 193707721, 147573952589676412927, 2677 228479, 2361183241434822606847, 1609 # However p=11 and p-29 can't - notice the program halts with Fail1 in # these cases but Fail2 in the preceding two cases. # We can explain this by looking at the prime factorisation of M11 and # M29: > > ifactor(M11);ifactor(M29); (23) (89) (233) (1103) (2089) > ifactor(232);ifactor(1102);ifactor(2088); 3 (2) (29) (2) (19) (29) 3 2 (2) (3) (29) # For p=11, the prime factors q=23, 89 both have p-1 dividing k! for # the first time when k=22; and for p=29 all 3 factors q have q-1 # dividing k! for the first time when k=29. Thus these composite # numbers cannot be factorised by the p-1 method (unless there is a # "fluke" with the choice of base). # # # Now try again with a=5: > pollard(M11,5,200);pollard(M23,5,200);pollard(M29,5,200); Fail1 47, 8388607, 23 Fail1 > pollard(M37,5,200);pollard(M41,5,200);pollard(M43,5,200); 223, 137438953471, 37 13367, 2199023255551, 163 431, 8796093022207, 43 > pollard(M47,5,200);pollard(M53,5,200);pollard(M59,5,200); 10610063, 140737488355327, 47 129728784761, 9007199254740991, 53 179951, 576460752303423487, 61 > pollard(M67,5,200);pollard(M71,5,200);pollard(M73,5,200); Fail2 Fail2 439, 9444732965739290427391, 73 > pollard(M79,5,200);pollard(M83,5,200);pollard(M97,5,200); 2687, 604462909807314587353087, 79 167, 9671406556917033397649407, 83 11447, 158456325028528675187087900671, 97 # The results are the same, with the same 4 failures, and the same # outcome if we increase kmax: > pollard(M11,5,5000);pollard(M29,5,5000);pollard(M67,5,5000);pollard(M7 > 1,5,5000); Fail1 Fail1 193707721, 147573952589676412927, 2677 228479, 2361183241434822606847, 1609 # This is because the p-1 method is very insensitive to the choice of # base a: only in "fluke" situations will changing a affect the result # of the test. # # (c) Repunits: For k=1 we get 1, so start from k=2 > for k from 2 to 30 do > print("k"= k);N:=(10^k-1)/9; pollard(N,2,50); > od; > "k" = 2 N := 11 Fail1 "k" = 3 N := 111 3, 111, 2 "k" = 4 N := 1111 11, 1111, 5 "k" = 5 N := 11111 41, 11111, 5 "k" = 6 N := 111111 3, 111111, 2 "k" = 7 N := 1111111 239, 1111111, 17 "k" = 8 N := 11111111 11, 11111111, 5 "k" = 9 N := 111111111 3, 111111111, 2 "k" = 10 N := 1111111111 451, 1111111111, 5 "k" = 11 N := 11111111111 21649, 11111111111, 41 "k" = 12 N := 111111111111 3, 111111111111, 2 "k" = 13 N := 1111111111111 4187, 1111111111111, 13 "k" = 14 N := 11111111111111 11, 11111111111111, 5 "k" = 15 N := 111111111111111 3, 111111111111111, 2 "k" = 16 N := 1111111111111111 17, 1111111111111111, 4 "k" = 17 N := 11111111111111111 Fail2 "k" = 18 N := 111111111111111111 3, 111111111111111111, 2 "k" = 19 N := 1111111111111111111 Fail2 "k" = 20 N := 11111111111111111111 451, 11111111111111111111, 5 "k" = 21 N := 111111111111111111111 3, 111111111111111111111, 2 "k" = 22 N := 1111111111111111111111 11, 1111111111111111111111, 5 "k" = 23 N := 11111111111111111111111 Fail2 "k" = 24 N := 111111111111111111111111 3, 111111111111111111111111, 2 "k" = 25 N := 1111111111111111111111111 41, 1111111111111111111111111, 5 "k" = 26 N := 11111111111111111111111111 11, 11111111111111111111111111, 5 "k" = 27 N := 111111111111111111111111111 3, 111111111111111111111111111, 2 "k" = 28 N := 1111111111111111111111111111 11, 1111111111111111111111111111, 5 "k" = 29 N := 11111111111111111111111111111 3191, 11111111111111111111111111111, 11 "k" = 30 N := 111111111111111111111111111111 3, 111111111111111111111111111111, 2 # We fail to find a factor in the 4 cases k=2, 17, 19, 23. For k=2, we # have N=11 which is prime. # For k=17 we can factorise N by increasing kmax: > N:=(10^17-1)/9; pollard(N,2,2000); N := 11111111111111111 2071723, 11111111111111111, 1069 > ifactor(N); (5363222357) (2071723) > ifactor(5363222356);ifactor(2071722); 2 (2) (17) (389) (202753) (2) (3) (17) (19) (1069) # The prime factorisation of p-1 for the two prime factors p of N shows # that 1069 steps are needed. # # In the other two cases, the repunits are prime: # > isprime((10^19-1)/9);isprime((10^23-1)/9); true true > > # Question 5 -- Pollard's Rho Algorithm # # > rho:=proc(n,c,x0,imax) local x,y,i,g;x:=x0;y:=x0; > for i from 1 to imax do > x:=x^2+c mod n; > y:=y^2+c mod n; y:=y^2+c mod n; > g:=gcd(x-y,n); > if g=n then RETURN(Fail); fi; > if g>1 then RETURN(g,i); fi; > od; > RETURN(Fail); > end: > # (a) > n:=403;rho(n,1,1,100); n := 403 31, 3 > n:=1891;rho(n,1,1,100); n := 1891 31, 3 > n:=5157437;rho(n,1,1,100); n := 5157437 2273, 36 > rho(F5,1,1,100); 641, 22 > rho(F6,1,1,1000); 274177, 808 > rho(F7,1,1,10000); Fail > rho(F7,2,3,50000); Fail # We are able to factorise all the numbers in Q4(a) except F7 (which we # couldn't factorise by the p-1 method either) # # (b) Mersenne numbers (excluding the Mersenne primes) > rho(M11,1,1,100);rho(M23,1,1,100);rho(M29,1,1,100); 23, 6 47, 12 486737, 14 > rho(M37,1,1,100);rho(M41,1,1,100);rho(M43,1,1,100); 223, 10 13367, 42 431, 18 > rho(M47,1,1,100);rho(M53,1,1,100);rho(M59,1,1,500); 2351, 50 6361, 67 179951, 444 > rho(M67,1,1,10000);rho(M71,1,1,500);rho(M73,1,1,500); 193707721, 5528 228479, 386 439, 24 > rho(M79,1,1,500);rho(M83,1,1,500);rho(M97,1,1,500); 2687, 81 167, 10 11447, 100 # So we manage to factorise all these Mersenne numbers, although the # number of iterations required is unpredictable (notice we need 5528 # for M67). If we change either x0 or c, we get different results (the # number of iterations changes, and sometimes we get a different # factor). Usually the number of iterations is (very roughly) # proportional to the square root of the factor found. # > rho(M11,1,2,100);rho(M23,1,2,100);rho(M29,1,2,100); 23, 4 47, 8 233, 12 > rho(M37,1,2,100);rho(M41,1,2,100);rho(M43,1,2,100); 223, 10 13367, 42 431, 18 > rho(M47,1,2,100);rho(M53,1,2,100);rho(M59,1,2,500); 2351, 50 6361, 67 179951, 444 > rho(M67,1,2,10000);rho(M71,1,2,500);rho(M73,1,2,500); 193707721, 5528 228479, 386 439, 24 > rho(M79,1,2,500);rho(M83,1,2,500);rho(M97,1,2,500); 2687, 81 167, 5 11447, 100 > > > > > rho(M11,3,1,100);rho(M23,3,1,100);rho(M29,3,1,100); 23, 2 47, 7 1103, 8 > rho(M37,3,1,100);rho(M41,3,1,500);rho(M43,3,1,100); 223, 30 13367, 84 431, 22 > rho(M47,3,1,100);rho(M53,3,1,100);rho(M59,3,1,1000); 2351, 22 6361, 88 179951, 566 > rho(M67,3,1,10000);rho(M71,3,1,500);rho(M73,3,1,500); 193707721, 6919 228479, 96 439, 28 > rho(M79,3,1,500);rho(M83,3,1,500);rho(M97,3,1,500); 2687, 84 167, 10 11447, 140 > > > # (c) repunits: > > for k from 2 to 30 do > print("k"= k);N:=(10^k-1)/9; rho(N,1,1,50); > od; "k" = 2 N := 11 Fail "k" = 3 N := 111 3, 1 "k" = 4 N := 1111 11, 4 "k" = 5 N := 11111 41, 7 "k" = 6 N := 111111 3, 1 "k" = 7 N := 1111111 239, 6 "k" = 8 N := 11111111 11, 4 "k" = 9 N := 111111111 3, 1 "k" = 10 N := 1111111111 11, 4 "k" = 11 N := 11111111111 Fail "k" = 12 N := 111111111111 3, 1 "k" = 13 N := 1111111111111 79, 4 "k" = 14 N := 11111111111111 11, 4 "k" = 15 N := 111111111111111 3, 1 "k" = 16 N := 1111111111111111 11, 4 "k" = 17 N := 11111111111111111 Fail "k" = 18 N := 111111111111111111 3, 1 "k" = 19 N := 1111111111111111111 Fail "k" = 20 N := 11111111111111111111 11, 4 "k" = 21 N := 111111111111111111111 3, 1 "k" = 22 N := 1111111111111111111111 11, 4 "k" = 23 N := 11111111111111111111111 Fail "k" = 24 N := 111111111111111111111111 3, 1 "k" = 25 N := 1111111111111111111111111 41, 7 "k" = 26 N := 11111111111111111111111111 869, 4 "k" = 27 N := 111111111111111111111111111 3, 1 "k" = 28 N := 1111111111111111111111111111 11, 4 "k" = 29 N := 11111111111111111111111111111 16763, 14 "k" = 30 N := 111111111111111111111111111111 3, 1 # With x0=c=1 we fail to find a factor in 50 steps for the same cases as # in Q4(c), viz,k = 2, 17, 19, 23 and also for 11. For k=2, 19 and 23 we # know that N is prime. Try changing parameters in the other cases: # > > N:=(10^11-1)/9; N := 11111111111 > rho(N,1,1,200); 21649, 162 > rho(N,2,3,100); 21649, 29 # For k=11 we get a factorisation, either by increasing imax or changing # the pseudorandom sequence. For k=17, the prime factors of N are large, # so we are likely to need imax large to find a factorisation. The # number of iterations needed may still vary dramatically with the # particular sequence used... # > N:=(10^17-1)/9;ifactor(N); N := 11111111111111111 (5363222357) (2071723) > rho(N,1,1,10000); 2071723, 2862 > rho(N,2,3,500); 2071723, 327 >