# Matthew Russell # Experimental Math # Homework 24 # I give permission to post this ##################################################################################### # Background: # The Perrin sequence P(n) is defined by # P(0)=3,P(1)=0,P(2)=2 # and for n ≥ 3, by the recurrence # P(n)=P(n-2)+P(n-3) # Perrin (and other people, e.g. E.B. Escott in 1909) proved that if n is prime, # then P(n)/n must be an integer. Perrin and Escott wondered whether the converse is true. # Escott believed that it was false, but did not find a counterexample. The first counterexample, # n=271441, was found in 1982 by Adams and Shanks, and again, in 2012, in a few minutes of programming # and a few seconds of computation, by Neil Lutz, and quite a few other people in today's class. # (a) Find all such n less than one million. # perrin(N): finds indices of all Perrin pseudoprimes up to N perrin:=proc(N) local L, j, new: L:=[3,0,2]: for j from 3 to N do new:=L[1]+L[2]: if evalb(0=new mod j) and not type(j,prime) then print(j,` gives Perrin pseudoprime`): fi: L:=[L[2],L[3],new]: od: return NULL: end: # The counterexamples are 271441=521^2 and 904631 (not a square). # (b) By using the Matrix formula given in the wiki article, find P(2^20). #GetPerrinN(N): finds Nth term of Perrin sequence, using matrix approach GetPerrinN:=proc(N) m:=Matrix([[0,1,0],[0,0,1],[1,1,0]]): m1:=LinearAlgebra[MatrixPower](m,N): v:=Vector([3,0,2]): v1:=LinearAlgebra[MatrixVectorMultiply](m1,v): return v1[1]: end: # GetPerrin(2^20); gives # 4.9074866239... x 10^127855 ##################################################################################### # 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 for n ≥ 3, by the recurrence # PAB(A,B,n)=A*PAB(A,B,n-2)+B*PAB(A,B,n-3) # Write a procedure, PseudoPrime(A,B,N) that inputs positive integers A and B # and a positive integer N # and finds all composite (i.e. non-prime) n such that PAB(A,B,n)/n is an integer. PseudoPrime:=proc(A,B,N) L:=[3,0,2*A]: Comp:={}: for j from 3 to N do new:=B*L[1]+A*L[2]: if evalb(0=new mod j) and not type(j,prime) then Comp:=Comp union {j}: fi: L:=[L[2],L[3],new]: od: return Comp: end: ##################################################################################### # Write a procedure Mul(P,Q) that inputs two permutations (written as lists) # of the same lengths and outputs their product. For example, # Mul([2,3,1,4],[2,1,4,3]); # should output # [1,4,2,3] Mul:=proc(P,Q) local i: return [seq(Q[P[i]],i=1..nops(P))]: end: ##################################################################################### # Write a procedure Gp(S) that inputs a set of permutations of the same length # (output FAIL if S is not a set or it is not a set of lists of the same length, etc.) # and outputs the set of all permutations that belong to the group generated by the members of S. # (Don't forget the identity permutation!). # 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) if not type(S,set) then return `FAIL`: fi: if nops(S)=0 then return `FAIL: S should be nonempty`: fi: length:=nops(S[1]): for i from 2 to nops(S) do if nops(S[i])<>length then return `FAIL`: fi: od: Done:={}: ToDo:={[seq(i,i=1..length)]}: while ToDo<>{} do T:=ToDo[1]: Done:=Done union {T}: ToDo:=ToDo minus {T}: new:=map2(Mul,T,S): ToDo:=ToDo union (new minus Done): od: return Done: end: ##################################################################################### # Find the subgroup of S_6 generated by the three rotations (about the x,y, and z axes) # of a cube. Using Nathaniel Shar's suggestion, let's stick to the usual convention: # Front: 1 ; Back: 6; Left:3; Right: 4; Bottom:2; Top=5. # The three rotations are: # S:={[4,2,1,6,5,3],[2,6,3,4,1,5],[1,4,2,5,3,6]}: # Then, # Gp(S); 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]} ##################################################################################### # Write a procedure PermToCyc(P) that inputs a permutation P and outputs 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) len:=nops(P): Left:={seq(i,i=1..len)}: cycs:={}: while Left<>{} do cy:=GetSingleCycle(P,Left[1]): Left:=Left minus (convert(cy,set)): cycs:=cycs union {cy}: od: return cycs: end: GetSingleCycle:=proc(P,i) L:=[i]: now:=P[i]: while now<>i do L:=[op(L),now]: now:=P[now]: od: return L: end: ##################################################################################### # Write a procedure Polya(G,c) that inputs a subgroup of the symmetric group S_n, G, # 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,c) return add(c^(nops(PermToCyc(g))),g in G)/nops(G): end: ##################################################################################### # Find a polynomial expression, CC(c), for the number 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. Is the sequence # [seq(CC(i),i=1..infinity)] # in Sloane? # CC(c) = Polya(Gp(S),c), where S:={[4,2,1,6,5,3],[2,6,3,4,1,5],[1,4,2,5,3,6]}. # = (1/24)*c^6+(1/2)*c^3+(1/8)*c^4+(1/3)*c^2 # seq(CC(i),i=1..10) = 1, 10, 57, 240, 800, 2226, 5390, 11712, 23355, 43450 # This is A047780. # I decided to consider the equivalent questions for: # # A regular tetrahedron: # S:={[1,3,4,2],[3,2,4,1],[2,4,3,1],[2,3,1,4]}: # Polya(Gp(S),c) = (1/12)*c^4+(11/12)*c^2, # which gives # 1, 5, 15, 36, 75, 141, 245, 400, 621, 925 # This is A006008. # # A regular octahedron: # S:={[1,4,5,6,7,2,3,8],[7,2,1,6,5,8,3,4],[5,4,3,8,7,6,1,2],[3,8,5,4,1,2,7,6],[2,7,8,3,4,1,6,5]}: # Polya(Gp(S),c); = (1/24)*c^8+(17/24)*c^4+(1/4)*c^2 # which gives # 1, 23, 333, 2916, 16725, 70911, 241913, 701968, 1798281, 4173775 # This is A000543.