girth.rb

require "matrix"
require "prime"
require "set"

Infinity = Float::INFINITY

CHARS = [+1, -1, +2, -2]
A = Matrix[[1, 2], [0, 1]]
B = Matrix[[1, 0], [2, 1]]

Ainv = Matrix[[1, -2], [0, 1]]
Binv = Matrix[[1, 0], [-2, 1]]

MAP_SHOW = {+1 => "A", -1 => "a", +2 => "B", -2 => "b"}
MAP_MAT = {+1 => A, -1 => Ainv, +2 => B, -2 => Binv}

def reverse(chars, mat)
  a, b, c, d = mat.to_a.flatten
  [chars.reverse, Matrix[[d, b], [c, a]]]
end

def swap_ab(chars, mat)
  a, b, c, d = mat.to_a.flatten
  [chars.map{|c| {+1 => +2, +2 => +1, -1 => -2, -2 => -1}[c] }, Matrix[[d, c], [b, a]]]
end

def swap_inv(chars, mat)
  a, b, c, d = mat.to_a.flatten
  [chars.map{|c| -c }, Matrix[[a, -b], [-c, d]]]
end

# 同値であって文字cが+1になるようなものを返す
def replace(chars, mat, c)
  case c
  when +1
    [chars, mat]
  when -1
    swap_inv(chars, mat)
  when +2
    swap_ab(chars, mat)
  when -2
    swap_inv(*swap_ab(chars, mat))
  end
end

def generate(n)
  i = +1
	generate0(n-1, i) {|x| yield [i] + x }
end

def generate0(n, last)
  if n == 0
    yield []
  else
    CHARS.each do |i|
      if i != -last
        generate0(n-1, i) {|x| yield [i] + x }
      end
    end
  end
end

def show(chars)
  chars.map{|x| MAP_SHOW[x] }.join("")
end

def subst(chars)
  chars.inject(Matrix.I(2)) {|r,c| r * MAP_MAT[c] }
end

def find_divisor(mat)
  gcd(mat[0,0] - 1, mat[0,1], mat[1,0], mat[1,1] - 1)
end

def gcd(*xs)
  xs.inject(0) {|g, x| gcd0(g, x) }
end

def gcd0(a, b)
  if b == 0
    a.abs
  else
    gcd0(b, a % b)
  end
end

def main()
  girth = {}
  memo = [[[], Matrix.I(2)]]
  (1..14).each do |n|
    seen = Set.new
    puts "length = #{n}"
    new_memo = []
    memo.each do |chars0, mat0|
      [+1, +2, -2].each do |c|
        chars, mat = replace(chars0, mat0, c)
        chars = [+1] + chars
        mat = A * mat
        next if seen.include?(chars)
        charsd, matd = replace(*reverse(chars, mat), chars.last)
        seen << charsd
        new_memo << [chars, mat]
        new_memo << [charsd, matd] if chars != charsd
        
        d = find_divisor(mat)
        factors = Prime.prime_division(d)
        ps = []
        factors.each do |p,e|
          if n <= (girth[p] || Infinity)
            ps << p
            girth[p] = n
          end 
        end
        if ps != []
          puts "#{show(chars)}: #{mat} = E (mod #{ps.join(",")})"
        end
      end
    end
    memo = new_memo
  end
  girth.keys.sort.each do |p|
    puts "G(#{p}) = #{girth[p]}"
  end
end

main() if $0 == __FILE__