# John Miller # hw24.txt - Homework assigned Mon, Apr 16, due Apr 22. Help := proc(): print(`PerrinPseudoPrimes(N), Perrin_power2(k), PseudoPrime(A,B,N)`): print(`SmallestPseudoPrime(A,B), SearchForLargePseudoPrimes(MaxA,MaxB,MinN)`): print(`SmallestPseudoPrimeModified(A,B), SearchForLargePseudoPrimesModified(MaxA,MaxB,MinN)`): print(`Mul(P,Q), Gp(S), CubeSubgroup()`): print(`PermToCyc(P), Poyla(G,c), CC(c)`): end: with(LinearAlgebra): ###################################################################### # Question #1 # Part(a) # PerrinPseudoPrimes(N): Find the Perrin pseudoprimes up to N. # PerrinPseudoPrimes(1000000); outputs {271441, 904631} PerrinPseudoPrimes := proc(N) local n, u0, u1, u2, u3, pseudoPrimes::set: pseudoPrimes := {}: u0 := 3: u1 := 0: u2 := 2: for n from 3 to N do: u3 := u0 + u1: if u3 mod n = 0 and not isprime(n) then pseudoPrimes := pseudoPrimes union {n}: fi: u0 := u1: u1:= u2: u2 := u3: od: pseudoPrimes: end: # end of PerrinPseudoPrimes(N) ###################################################################### # Question #1 # Part(b): Find P(2^20), where P is the Perrin sequence. # Answer is pretty big: # 4907486623903722805297725232245265827055249791214434138116795812716167879116566939532285784203887477[...127856 digits...]2821351663672505468646787355488245123653236152134178820592060343691155925212651423598055857578598490 # Perrin_power2(k): Inputs k and outputs the 2^k term of the Perrin # sequence. Perrin_power2 := proc(k::nonnegint) local i, M::Matrix: M := Matrix([[0,1,0],[0,0,1],[1,1,0]]): for i from 1 to k do: M := Multiply(M,M): od: Multiply(M,<3,0,2>)[1]: end: # end of Perrin_power2(k) ###################################################################### # Question 2 # We define the generalized Perrin sequence PAB(A,B,n) by: # PAB(A,B,0) = 3, PAB(A,B,1) = 0, P(A,B,2) = 2*A, and # PAB(A,B,n) = A*PAB(A,B,n-2) + B*PAB(A,B,n-3), for n >= 3. # PseudoPrime(A,B,N): Inputs positive integers A and B, and # a positive integer N, and finds all composite n such that # PAB(A,B,n)/n is an integer. PseudoPrime := proc(A::posint,B::posint,N::posint) local n, u0, u1, u2, u3, pseudoPrimes::set: pseudoPrimes := {}: u0 := 3: u1 := 0: u2 := 2*A: for n from 3 to N do: u3 := B*u0 + A*u1: if u3 mod n = 0 and not isprime(n) then pseudoPrimes := pseudoPrimes union {n}: fi: u0 := u1: u1:= u2: u2 := u3: od: pseudoPrimes: end: # end of PseudoPrime(A,B,N) # Some example output: # PseudoPrime(1,1,300000); outputs {271441} # PseudoPrime(2,1,300000); outputs: # {4, 8, 16, 32, 64, 128, 256, 368, 512, 705, 736, 1024, 1472, 2048, 2465, 2737, 2944, 3745, 4096, 4181, 5777, 5888, 6721, 8192, 8464, 10877, 11776, 13201, 15251, 16384, 16928, 23552, 24465, 29281, 32768, 33856, 34561, 35296, 35785, 47104, 49216, 50416, 51841, 54705, 64079, 64681, 65536, 67712, 67861, 68251, 70592, 75077, 80189, 90061, 94208, 96049, 97921, 98048, 98432, 100065, 100127, 100832, 105281, 113573, 118441, 131072, 135424, 141184, 146611, 161027, 162133, 163081, 179697, 186961, 188416, 194672, 194833, 196096, 196864, 197209, 201664, 202688, 209665, 219781, 228241, 229445, 231703, 252601, 254288, 254321, 257761, 262144, 268801, 270848, 272611, 282368, 283361} # PseudoPrime(1,2,300000); outputs: # {42, 1729, 6790, 15457, 19045, 26866, 32041, 214367} # PseudoPrime(1,3,300000); outputs: # {9, 27, 81, 135, 243, 729, 2187, 3159, 6561, 7155, 12177, 19683, 55971, 59049, 70443, 145881, 176661, 177147, 189999, 295245} # Note that these are all multiples of 3. # PseudoPrime(3,1,300000); outputs: # {114, 1474, 5833, 19345, 41154, 74670, 96139, 146081, 146611, 282133, 286903, 294409} ###################################################################### # Question 2 - EXTRA CREDIT! # SmallestPseudoPrime(A,B): Inputs positive integers A and B, and # a positive integer N, and finds the smallest composite n such that # PAB(A,B,n)/n is an integer. SmallestPseudoPrime := proc(A::posint,B::posint) local n, u0, u1, u2, u3: u0 := 3: u1 := 0: u2 := 2*A: for n from 3 do: u3 := B*u0 + A*u1: if u3 mod n = 0 and not isprime(n) then return(n): fi: u0 := u1: u1:= u2: u2 := u3: od: end: # end of SmallestPseudoPrime(A,B) ###################################################################### # Question 2 - EXTRA CREDIT! # SmallestPseudoPrimeModified(A,B): Inputs positive integers A and B, and # a positive integer N, and finds the smallest composite n such that # PAB(A,B,n)/n is an integer, and such that n is not divisible by A or B # (unless A or B is 1). SmallestPseudoPrimeModified := proc(A::posint,B::posint) local n, u0, u1, u2, u3: u0 := 3: u1 := 0: u2 := 2*A: for n from 3 do: u3 := B*u0 + A*u1: if u3 mod n = 0 and not isprime(n) and (A = 1 or n mod A <> 0) and (B = 1 or n mod B <> 0) then return(n): fi: u0 := u1: u1:= u2: u2 := u3: if n mod 100000 = 0 then print(`checking`,A,B,n): fi: od: end: # end of SmallestPseudoPrimeModified(A,B) ###################################################################### # Question 2 - EXTRA CREDIT! # SmallestPseudoPrimeModified_v2(A,B): Inputs positive integers A and B, and # a positive integer N, and finds the smallest composite n such that # PAB(A,B,n)/n is an integer, and such that n is not divisible by A or B # (unless A or B is 1). This version 2 uses binary powering of matrices # mod m in an attempt to avoid the huge memory usage of the original # version when searching for large pseudoprimes. It's also useful if you want # to start checking from n much bigger than 1. But even with the efficiency # of binary powering, it still runs pretty slowly. So the quickest approach is # to run the original version until the memory usage blows the machine up, # and then swith to this, starting from where the first one left off. SmallestPseudoPrimeModified_v2 := proc(A::posint,B::posint) local n, u0, u1, u2, u3, M::Matrix,Mpower::Matrix: u0 := 3: u1 := 0: u2 := 2*A: M := Matrix([[0,1,0],[0,0,1],[B,A,0]]): for n from 2000000 do: # was 3 if not isprime(n) and (A = 1 or n mod A <> 0) and (B = 1 or n mod B <> 0) then Mpower := LinearAlgebra[Modular][MatrixPower](n,M,n): u3 := LinearAlgebra[Modular][Multiply](n,Mpower,<3,0,2*A>)[1]: if u3 = 0 then return(n): fi: fi: if n mod 100000 = 0 then print(`checking`,A,B,n): fi: od: end: # end of SmallestPseudoPrimeModified(A,B) ###################################################################### # Question 2 - EXTRA CREDIT! # SearchForLargePseudoPrimes(MaxA,MaxB,MinN): For all A <= MaxA and # all B <= MaxB, outputs the smallest pseudoprimes if they are greater # than MinN. SearchForLargePseudoPrimes := proc(MaxA::posint,MaxB::posint,MinN::posint) local A,B,N,S::set: S := {}: for A from 1 to MaxA do: for B from 1 to MaxB do: N := SmallestPseudoPrime(A,B): if N >= MinN then S := S union {[A,B,N]}: print(A,B,N): fi: od: od: S: end: # end of SearchForLargePseudoPrimes(MaxA,MaxB,MinN) # Here is some output: # SearchForLargePseudoPrimes(100,100,1000) outputs: # {[1, 1, 271441], [1, 4, 3481], [1, 32, 1034], [1, 43, 1849], [1, 59, 2065], [1, 64, 2047], [1, 74, 1369], [1, 83, 2211], [9, 1, 1349], [27, 1, 7158], [41, 71, 1681], [47, 71, 1065], [61, 1, 1891], [61, 71, 1105], [97, 43, 1261], [97, 97, 1073]} # SearchForLargePseudoPrimes(1000,1000,20000) outputs: # {[1, 1, 271441], [1, 571, 62239], [1, 617, 37639], [1, 797, 84945], [1, 857, 35137], [1, 883, 49105], [1, 947, 22801], [1, 953, 27637], [283, 179, 32041], [283, 251, 43978], [467, 647, 62745], [467, 653, 50281], [509, 719, 25165], [521, 911, 27613], [587, 157, 20545], [587, 383, 20545], [587, 613, 20545], [643, 643, 168466], [691, 661, 44287], [691, 857, 55969], [719, 661, 25779], [719, 853, 22289], [719, 887, 23221], [729, 1, 107801], [773, 719, 27055], [787, 293, 31429], [859, 953, 22334], [929, 599, 28153], [929, 863, 30889]} # Two things to observe about this result: # All of the As and Bs are prime, except for 729=3^6, for which # A=729,B=1 yields a respectably large smallest pseudoprime of 107801. # The other large smallest pseudoprime is 168446, given by A=643,B=643. # SearchForLargePseudoPrimes(1000,10000,200000) outputs: # {[1,1,271441],[1,6829,293647],[1,9871,286259],[587,8093,299441],[677,8837,258529],[829,5153,376366],[839,4813,703921],[857,5087,340829]} # Note that a few of these are bigger than 271441. The largest # is 340829, given by A=839, B=4813. ###################################################################### # Question 2 - EXTRA CREDIT! # SearchForLargePseudoPrimesModified(MaxA,MaxB,MinN): For all A <= MaxA and # all B <= MaxB, outputs the smallest pseudoprimes if they are greater # than MinN, excluding pseudoprimes divisible by A or B (if A or B <> 1). SearchForLargePseudoPrimesModified := proc(MaxA::posint,MaxB::posint,MinN::posint) local A,B,N,S::set: S := {}: for A from 1 to MaxA do: for B from 1 to MaxB do: N := SmallestPseudoPrimeModified(A,B): if N >= MinN then S := S union {[A,B,N]}: print(A,B,N): fi: od: od: S: end: # end of SearchForLargePseudoPrimesModified(MaxA,MaxB,MinN) # SearchForLargePseudoPrimesModified(100,100,100000) outputs: # {[1, 1, 271441], [1, 3, 535931], [2, 79, 455533], [3, 61, 362881], [4, 2, 170039], [4, 8, 271441], [5, 11, 864434], [5, 41, 709171], [7, 13, 502681], [9, 81, 535931], [11, 4, 173049], [13, 23, 103361], [17, 59, 1605291], [19, 71, 111361], [25, 19, 226801], [29, 13, 234201], [31, 17, 1168883], [31, 59, 157807], [31, 97, 351649], [41, 8, 275887], [41, 43, 422801], [41, 64, 802105], [43, 41, 2976487], [49, 59, 136115], [59, 5, 185761], [59, 23, 140577], [79, 43, 409981], [79, 59, 288637], [97, 43, 162401]} # So using this search, we found some interesting "restricted" pseudoprimes: # 535931 (A= 1,B= 3) # 864434 (A= 5,B=11) # 1605291 (A=17,B=59) # 1167773 (A=31,B=17) # 802105 (A=41,B=64) (the only prominent composite among the As and Bs) # 2976487 (A=43,B=41) (2976497 = 863*3449) ###################################################################### # Question 3 # Mul(P,Q): Inputs two permutations (written as lists) # of the same lengs, and outputs their product. For example, # Mul([2,3,1,4],[2,1,4,3]) should output # [1,4,2,3] (i.e. 1 -> 2 -> 1, 2 -> 3 -> 4, 3 -> 1 -> 2, 4 -> 4 -> 3) Mul := proc(P::list,Q::list) local L::list, i: if nops(P) <> nops(Q) then return(FAIL): fi: [seq(Q[P[i]],i=1..nops(P))]: end: # end of Mul(P,Q) ###################################################################### # Question 4 # Gp(S): Inputs a set of permutaions of the same length, and outputs # the set of all permuations that belong to the group generated by # the members of S. For example, # Gp({[2,1]}); should output {[1,2],[2,1]}, and # Gp({[2,1,3,4],[1,2,4,3]}); should output # {[1,2,3,4],[2,1,3,4],[1,2,4,3],[2,1,4,3]}. Gp := proc(S::set) local i, oldGroup::set, newGroup::set, N::posint, s1::list, s2::list, sizes::set, IdPerm: sizes := {seq(nops(s),s in S)}: if nops(sizes) <> 1 then return(FAIL): else N := op(sizes): fi: IdPerm := [seq(i,i=1..N)]: oldGroup := S union {IdPerm}: newGroup := oldGroup: do: for s1 in oldGroup do: for s2 in oldGroup do: newGroup := newGroup union {Mul(s1,s2)}: od: od: if newGroup = oldGroup then return(newGroup): else oldGroup := newGroup: fi: od: end: # end of Gp(S) ###################################################################### # Question 5 # CubeSubgroup(): Finds the subgroup of S_6 generated by the # three rotations (about the x,y and z axes) of a cube. # Front:1; Back:6; Left:3; Right:4; Bottom:2; Top:5. # The output is: # {[1, 2, 3, 4, 5, 6], [1, 3, 5, 2, 4, 6], [1, 4, 2, 5, 3, 6], [1, 5, 4, 3, 2, 6], [2, 1, 4, 3, 6, 5], [2, 3, 1, 6, 4, 5], [2, 4, 6, 1, 3, 5], [2, 6, 3, 4, 1, 5], [3, 1, 2, 5, 6, 4], [3, 2, 6, 1, 5, 4], [3, 5, 1, 6, 2, 4], [3, 6, 5, 2, 1, 4], [4, 1, 5, 2, 6, 3], [4, 2, 1, 6, 5, 3], [4, 5, 6, 1, 2, 3], [4, 6, 2, 5, 1, 3], [5, 1, 3, 4, 6, 2], [5, 3, 6, 1, 4, 2], [5, 4, 1, 6, 3, 2], [5, 6, 4, 3, 1, 2], [6, 2, 4, 3, 5, 1], [6, 3, 2, 5, 4, 1], [6, 4, 5, 2, 3, 1], [6, 5, 3, 4, 2, 1]} CubeSubgroup := proc() local rotX::list, rotY::list, rotZ::list: rotX := [1,4,2,5,3,6]: rotY := [2,6,3,4,1,5]: rotZ := [4,2,1,6,5,3]: Gp({rotX,rotY,rotZ}): end: # end fo CubeSubgroup() ###################################################################### # Question 6 # PermToCyc(P): Inputs a permutation P and otuputs its cycle structure, # as a set of cycles, with the convention that the smallest entry # of the cycle is written first. For example, # PermToCyc([3,5,6,4,2,1]); should output {[1,3,6],[2,5],[4]}. PermToCyc := proc(P::list) local i, S::set, Pleft::set, L::list, firstUp, nextUp: S := {}: Pleft := convert(P,set): while Pleft <> {} do: firstUp := min(Pleft): L := [firstUp]: Pleft := Pleft minus {firstUp}: nextUp := P[firstUp]: while nextUp <> firstUp do: L := [op(L),nextUp]: Pleft := Pleft minus {nextUp}: nextUp := P[nextUp]: od: S := S union {L}: od: S: end: # end of PermToCyc(P) ###################################################################### # Question 7 # Polya(G,c): Inputs a subgroup G of the symmetric group S_n, and # outputs the number of ways to color a combinatorial structure on # n vertices with c colors, whose group of symmetry happens to be G. # For example, # Polya({[1,2],[2,1]},c); should output (c^2+c)/2. Polya := proc(G::set,c) local g::list: add(c^nops(PermToCyc(g)),g in G)/nops(G): end: # end of Poyla(G,c) ###################################################################### # Question 8 # CC(c): Outputs the # of ways of coloring a cube with c colors, # where two colorings are considered the same if one can be gotten # from the other by a sequence of rotations as above. # The output is: # (1/24)*c^6+(1/2)*c^3+(1/8)*c^4+(1/3)*c^2 # This is in Sloane as A047780 CC := proc(c) local G: G := CubeSubgroup(): Polya(G,c): end: # end of CC(c)