using Oscar

spor_names = all_character_table_names(is_sporadic_simple => true,
  is_duplicate_table => false; ordered_by = order);
println(spor_names)

normindices = Dict("Co1" => 1, "J4" => 1, "F3+" => 1, "B" => 1, "M" => 1);
for name in spor_names
  n = character_table("$(name)N2")
  if n == nothing
    println(name)
  else
    two_part = 1
    ord = order(n)
    while iseven(ord)
      ord = div(ord, 2)
      two_part = 2 * two_part
    end
    normindices[name] = div(order(n), two_part)
  end
end

maxbound = [];
for name in filter(x -> x != "M" && x != "B", spor_names)
  t = character_table(name)
  mx = [character_table(m) for m in maxes(t)]
  odd = filter(s -> isodd(div(order(t), order(s))), mx)
  primperm = [trivial_character(s)^t for s in odd]
  wsum = normindices[name] * sum(pi -> [ZZ(x.data) // degree(pi) for x in values(pi)], primperm)
  push!(maxbound, [name, maximum(wsum[2:length(wsum)])])
end
sort!(maxbound, by = x -> -x[2]);
println(maxbound[1])

filename = string(GAP.Globals.Filename(GAP.Globals.DirectoriesPackageLibrary(GAP.Obj("ctbllib"),GAP.Obj("data")), GAP.Obj("prim_perm_M.json")));
str = read(filename, String);
monster_prim_data = Oscar.JSON.parse(str; inttype = BigInt)[2];
length(monster_prim_data)

t = character_table("M");
monstermaxindices = [];
monstermaxtables = [];
for entry in monster_prim_data
  if length(entry) == 1
    s = character_table(entry[1])
    push!(monstermaxtables, s)
    push!(monstermaxindices, div(order(t), order(s)))
  else
    push!(monstermaxtables, nothing)
    push!(monstermaxindices, entry[2][1])
  end
end;
odd_prim = Vector{ZZRingElem}[];
for i in 1:length(monster_prim_data)
  if isodd(monstermaxindices[i])
    if monstermaxtables[i] != nothing
      push!(odd_prim, [ZZ(x.data) for x in trivial_character(monstermaxtables[i])^t])
    else
      push!(odd_prim, Vector{ZZRingElem}(monster_prim_data[i][2]))
    end
  end
end;
length(odd_prim)

wsum = normindices["M"] * sum(pi -> [x // pi[1] for x in values(pi)], odd_prim);
mx = maximum(wsum[2:length(wsum)])

function upper_bound( tblS, tblM, p )
  facts = collect(factor(order(tblS)))
  pos = findfirst(x -> x[1] == p, facts)

  ppart = p^facts[pos][2]
  ppartposS = filter(i -> orders_class_representatives(tblS)[i] == ppart, 1:nrows(tblS))
  if 0 < length(ppartposS)
    # P is cyclic.
    if tblM == nothing
      return (orders_centralizers(tblS)[ppartposS[1]] * euler_phi(ppart)
              // length(ppartposS)) // ppart
    else
      ppartposM = filter(i -> orders_class_representatives(tblM)[i] == ppart, 1:nrows(tblM))
      return (orders_centralizers(tblS)[ppartposS[1]] * euler_phi(ppart)
              // length(ppartposS)) //
             (orders_centralizers(tblM)[ppartposM[1]] * euler_phi(ppart)
              // length(ppartposM))
    end
  end

  # Compute |N_S(P)|.
  n = character_table("$(identifier(tblS))N$p")
  if n != nothing
    N_S = order(n)
  elseif p == 2
    N_S = ppart * normindices[identifier(tblS)]
  elseif identifier(tblS) == "M" && p == 3
    # The Sylow 3-normalizer is contained in 3^(3+2+6+6):(L3(3)xSD16)
    N_S = ppart * 2^6
  elseif identifier(tblS) == "F3+" && p == 3
    N_S = ppart * 8
  else
    error("cannot compute |N_S(P)|")
  end

  if tblM == nothing
    return N_S // ppart
  elseif sum(class_lengths(tblM)[class_positions_of_pcore(tblM, p)]) == ppart
    # P is normal in M.
    return N_S // order(tblM)
  end
    
  # Inspect known character tables of subgroups of M.
  f = N_S // ppart
  for subname in Oscar.names_of_fusion_sources(tblM)
    u = character_table(subname)
    if class_positions_of_kernel(known_class_fusion(u, tblM)[2]) == [1] &&
       sum([order(u) // x for x in orders_centralizers(u)[class_positions_of_pcore(u, p)]]) == ppart
      f = min(f, N_S // order(u));
    end
  end

  return f;
end;

badcases_thm = [];
badcases_strong = [];
for name in filter(x -> x != "M" && x != "B", spor_names)
  t = character_table(name)
  orders = orders_class_representatives(t)
  n = nrows(t)
  mx = [character_table(m) for m in maxes(t)]
  for p in prime_divisors(order(t))
    primperm = []
    good = filter(s -> mod(order(t) // order(s), p) != 0, mx)
    for s in good
      if known_class_fusion(s, t)[1]
        push!(primperm, [ZZ(x.data) for x in trivial_character(s)^t])  
      else
        ind = Set([induce(trivial_character(s), t, map)
                   for map in possible_class_fusions(s, t)])
        length(ind) == 1 || error("permutation character not uniquely determined")
        push!(primperm, [ZZ(x.data) for x in pop!(ind)])    
      end
    end
    weighted = []
    for i in 1:length(good)
      bd = upper_bound(t, good[i], p)
      pi = primperm[i]
      push!(weighted, [bd * pi[j] // pi[1] for j in 1:length(pi)])
    end
    sumx = sum(weighted)
    badpos = filter(i -> sumx[i] >= 1, 2:length(sumx))
    if length(badpos) != 0
      # We have to deal with this case for the stronger variant of Theorem 1.7.
      push!(badcases_strong, (name, p, copy(badpos)))
      for i in copy(badpos)
        facts = collect(factor(orders[i]))
        if length(facts) == 1
          q = facts[1][1]
          # If we know some `q`-element for which the bound is good enough
          # then we need not deal with this class for Theorem 1.7.
          if any(j -> mod(orders[j], q) == 0 &&
                      length(collect(factor(orders[j]))) == 1 &&
                      !(j in badpos), 2:n)
            setdiff!(badpos, i)
          end
        end
      end
      length(badpos) != 0 && push!(badcases_thm, (name, p, badpos))
    end
  end
end;

badcases_thm

badcases_strong

function prove_generation(name, prg, p)
  S = atlas_group(name)
  P = sylow_subgroup(S, p)[1]
  reps = evaluate(prg, gens(S))
  good = elem_type(S)[]
  for x in filter(!isone, reps)
    U = P
    g = x
    while order(U) < order(S)
      g = x^rand(S)
      U = sub(S, vcat(gens(P), [g]))[1]
    end
    push!(good, g)
  end
  return (P, good)
end;

for (name, p, badpos) in badcases_strong
  prg = atlas_program(name, :classes)
  if prg == nothing
    println("no classes script for ", (name, p, badpos))
  else
    prove_generation(name, prg, p)
  end
end

S = atlas_group("He");
x = rand(S);
while order(x) != 10 x = rand(S); end
x = x^5;
P5 = sylow_subgroup(S, 5)[1];
g = x^rand(S);
U = sub(S, vcat(gens(P5), [g]))[1]
while order(U) != order(S)
  g = x^rand(S)
  U = sub(S, vcat(gens(P5), [g]))[1]
end

S = atlas_group("Fi22");
x = rand(S);
while order(x) != 30 x = rand(S); end
x = x^15;
P5 = sylow_subgroup(S, 5)[1];
g = x^rand(S);
U = sub(S, vcat(gens(P5), [g]))[1]
while order(U) != order(S)
  g = x^rand(S)
  U = sub(S, vcat(gens(P5), [g]))[1]
end
P7 = sylow_subgroup(S, 7)[1];
g = x^rand(S);
U = sub(S, vcat(gens(P7), [g]))[1]
while order(U) != order(S)
  g = x^rand(S);
  U = sub(S, vcat(gens(P7), [g]))[1]
end

S = atlas_group("Fi23");
x = rand(S);
while order(x) != 42 x = rand(S); end
x = x^21;
P5 = sylow_subgroup(S, 5)[1];
g = x^rand(S);
U = sub(S, vcat(gens(P5), [g]))[1]
while order(U) != order(S)
  g = x^rand(S)
  U = sub(S, vcat(gens(P5), [g]))[1]
end
P7 = sylow_subgroup(S, 7)[1];
g = x^rand(S);
U = sub(S, vcat(gens(P7), [g]))[1]
while order(U) != order(S)
  g = x^rand(S);
  U = sub(S, vcat(gens(P7), [g]))[1]
end

x60 = rand(S);
while order(x60) != 60 x60 = rand(S); end
x = x60^30;
g = x^rand(S);
U = sub(S, vcat(gens(P5), [g]))[1]
while order(U) != order(S)
  g = x^rand(S)
  U = sub(S, vcat(gens(P5), [g]))[1]
end

x = rand(S);
while order(x) != 60 x = rand(S); end
x = x60^20;
g = x^rand(S);
U = sub(S, vcat(gens(P5), [g]))[1]
while order(U) != order(S)
  g = x^rand(S)
  U = sub(S, vcat(gens(P5), [g]))[1]
end

t = character_table("B");
ords = orders_class_representatives(t);
pos7 = filter(i -> ords[i] == 7, 1:length(ords))
pos47 = filter(i -> ords[i] == 47, 1:length(ords))
class_multiplication_coefficient(t, 2, pos7[1], pos47[1]) != 0
filter(x -> mod(order(character_table(x)), 47) == 0, maxes(t))

t = character_table("M");
orders = orders_class_representatives(t);
for p in setdiff(prime_divisors(order(t)), [19])
  goodpos = filter(i -> mod(monstermaxindices[i], p) != 0, 1:length(monster_prim_data))
  sumx = zeros(ZZRingElem, ncols(t))
  for i in goodpos
    if length(monster_prim_data[i]) == 2
      # We know the permutation character but not the subgroup table.
      bd = upper_bound(t, nothing, p)
      sumx = sumx + [bd * ZZ(x) // monstermaxindices[i] for x in monster_prim_data[i][2]]
    else
      s = monstermaxtables[i]
      if known_class_fusion(s, t)[1]
        # We can compute the permutation character.
        bd = upper_bound(t, s, p)
        sumx = sumx + [bd * ZZ(x.data) // monstermaxindices[i] for x in trivial_character(s)^t]
      else
        # We get only candidates for the permutation character.
        cand = Set([[ZZ(x.data) for x in induce(trivial_character(s), t, map)]
                    for map in possible_class_fusions(s, t)])
        # For each column, take the maximum of the possible values.        
        mat = matrix(ZZ, collect(cand))
        mx = [maximum(mat[1:nrows(mat), i]) for i in 1:ncols(mat)]
        bd = upper_bound(t, s, p)
        sumx = sumx + [bd * x // monstermaxindices[i] for x in mx]
      end
    end
  end
  badpos = filter(i -> sumx[i] >= 1, 2:length(sumx));
  if length(badpos) > 0
    error("p = $p: check open cases in $badpos");
  end
end;

pos19 = filter(i -> orders[i] == 19, 1:length(orders))

pos71 = filter(i -> orders[i] == 71, 1:length(ords))

class_multiplication_coefficient(t, pos19[1], pos19[1], pos71[1]) != 0

any(x -> mod(order(t) // x, 19 * 71) == 0, monstermaxindices)

any(x -> mod(order(character_table(x)), 71) == 0, ["L2(13)", "U3(4)", "U3(8)", "Sz(8)"])