# 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)