#!/usr/local/bin/maple #-*- maplev -*- # Nathaniel Shar # HW 24 # Experimental Mathematics # It is okay to link to this assignment on the course webpage. Help := proc(): print(`PerrinPseudoprimes(n), PerrinMatrix(n), Perrin(n), PABPseudoprimes(A,B,n), Mul(P, Q), Gp(S), PermToCyc(P), Polya(G, c)`): end: ############# # Problem 1 # ############# PerrinPseudoprimes := proc(n) local a,b,c, i, t, rv: a := 3: b := 0: c := 2: i := 2: rv := {}: while i <= n do: i := i+1: t := a+b: a := b: b := c: c := t: if c mod i = 0 and not isprime(i) then: rv := rv union {i}: fi: od: return rv: end: # The only two less than 10^6 are 271441 and 904631. PerrinMatrix := proc(n) option remember: if n = 0 then: return Matrix(3, shape=identity): elif n mod 2 = 0 then: return PerrinMatrix(n/2)^2: else: return Matrix([[0,1,0], [0,0,1], [1,1,0]]) . PerrinMatrix((n-1)/2)^2: fi: end: Perrin := proc(n): return (PerrinMatrix(n).Matrix(3,1,[3,0,2]))[1,1]: end: # Perrin(2^20) is approximately 0.491*10^128056. ############# # Problem 2 # ############# PABPseudoprimes := proc(A,B,n) local a,b,c, i, t, rv: a := 3: b := 0: c := 2*A: i := 2: rv := {}: while i <= n do: i := i+1: t := B*a+A*b: a := b: b := c: c := t: if c mod i = 0 and not isprime(i) then: rv := rv union {i}: fi: od: return rv: end: ############# # Problem 3 # ############# Mul := proc(P, Q): return map(i -> Q[i], P): end: ############# # Problem 4 # ############# Gp := proc(S) local todo, G, product, rho, pi: G := S: todo := S: while todo <> {} do: rho := todo[1]: todo := todo minus {rho}: for pi in G do: product := Mul(pi, rho): if not member(product, G) then: todo := todo union {product}: G := G union {product}: fi: od: od: return G: end: ############# # Problem 5 # ############# # [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]} ############# # Problem 6 # ############# PermToCyc := proc(P) local rv, c, i, j, seen: seen := {}: rv := {}: for i from 1 to nops(P) do: if not member(i, seen) then: seen := seen union {i}: c := [i]: j := P[i]: while j <> i do: c := [op(c),j]: seen := seen union {j}: j := P[j]: od: rv := rv union {c}: fi: od: return rv: end: ############# # Problem 7 # ############# Polya := proc(G, c): return add(c^nops(PermToCyc(pi)), pi=G)/nops(G): end: ############# # Problem 8 # ############# # CC(c) = (1/24)*(c^6 + 3c^4 + 12c^3 + 8c^2). # Of course the sequence CC(1), CC(2), ... is in OEIS (A047780).