(* Helper function to determine if big is divided by factor k
 * times. *)
let divs big factor = (big mod factor) = 0

let is_odd n = (n mod 2) = 1

(* We are not interested in prime (and thus abelian) groups.  This
 * implements a naive check for primality. *)
let rec is_prime n = is_prime_h n 2
and is_prime_h n p =
  if n <= p then true
  else if (divs n p) then false
  else is_prime_h n (p + 1)

(* Compute the prime factorialization of n. *)
let rec prime_factors n = prime_factors_h n 2
and prime_factors_h n p =
  if n = 1 then []
  else if (is_prime p) && (divs n p) then 
    p::(prime_factors_h (n / p)) p
  else prime_factors_h n (p + 1)

let str_of_intl l =
  let (_, res) =
    List.fold_left
      (fun (pre, res) elt -> (", ", res ^ pre ^ (string_of_int elt)))
      ("", "")
      l in
      res

(* Generate a string representation of a prime factorialization. *)
let str_of_prime_fact factors =
  let str_of_pow base exp 
      = (string_of_int base) ^ "^" ^ (string_of_int exp) in
  let counts =
    List.fold_left
      (fun counts elt ->
         (match counts with
              (num, count)::rest -> 
                if num = elt then (num, count + 1)::rest
                else (elt, 1)::(num, count)::rest
            | _ -> [(elt, 1)]))
      []
      factors in
  let (_, res) =
    List.fold_right
      (fun (num, count) (pre, rep) ->
         (" * ", rep ^ pre ^ (str_of_pow num count)))
      counts
      ("", "") in
    res

(* Decide if n | p!.  A naive implementation would quickly lead to
 * 32-bit integer overflow, but this implementation is robust against
 * this issue. *)
let rec divs_fact p n =
  if n = 1 then true
  else
    let small_factors = prime_factors n in
    let rec fact_list n = if n = 1 then [] else n::(fact_list (n - 1)) in
    let big_factors = fact_list p in

    let (_, res) =
      List.fold_left
        (fun (big_facts, fits) lil_fact ->
           if not fits then (big_facts, false)
           else 
             let (new_acc, found) =
               List.fold_left
                 (fun (acc', found_div) big_fact -> 
                    if found_div then (big_fact::acc', true)
                    else if (divs big_fact lil_fact) then 
                      ((big_fact / lil_fact)::acc', true)
                    else (big_fact::acc', false))
                 ([], false)
                 big_facts in
               if not found then (new_acc, false)
               else (new_acc, true))
        (big_factors, true)
        small_factors in
      res

(* Decompose n into a product of the full power of p in it and the
 * remaining factor. *)
let rec pow_div n p =
  let rem = n mod p in
    if not (rem = 0) then (1, n)
    else
      let quot = n / p in
      let (pow, fac) = pow_div quot p in
        (pow * p, fac)

(* Sound checkers for if there does not exist a simple group of a
 * given size: *)

(* Decide if a number is the product of three distinct primes. *)
let rec is_prod_pqr n = is_prod_pqr_h n 2 false
and is_prod_pqr_h n p one_out =
  if n <= p then false
  else if (not (is_prime p)) or (not (divs n p)) then
    (is_prod_pqr_h n (p + 1) one_out)
  else 
    let fac = n / p in
      if one_out then (is_prime fac)
      else (is_prod_pqr_h fac (p + 1) true)

(* Decide if a number is the product of a power of a prime and a
 * prime. *)
let rec is_prod_paq n = is_prod_paq_h n 2
and is_prod_paq_h n p =
  if n <= p then false
  else if (not (is_prime p)) then (is_prod_paq_h n (p + 1))
  else
    let (pow, rem) = pow_div n p in
      if pow = 1 then is_prod_paq_h n (p + 1)
      else if pow = p then (rem = 1) or (is_pow_prime rem)
      else (rem = 1) or (is_prime rem)
and is_pow_prime n = is_pow_prime_h n 2
and is_pow_prime_h n p =
  if p <= n then is_prime p
  else if (not (is_prime p)) then is_pow_prime_h n (p + 1)
  else 
    let (pow, rem) = pow_div n p in
      if pow = 1 then is_pow_prime_h n (p + 1)
      else rem = 1

(* Compute the potential orders of normalizers. *)
let norm_orders n =
  (* Compute all valid values for |N(P)|. *)
  let pnorm_orders p =
    let (_, index) = pow_div n p in
    let rec index_exists i =
      if i >= index then []
      else if ((divs index i) && (i mod p = 1) && (divs_fact i n))
      then (index_exists (i + 1)) @ [n / i]
      else index_exists (i + 1) in
      index_exists 2 in

  (* Iterate over all primes that could divide n. *)
  let rec it_over_p p =
    if p >= n then []
    else if not ((is_prime p) && (divs n p)) then it_over_p (p + 1)
    else (p, pnorm_orders p)::(it_over_p (p + 1)) in
      it_over_p 2

(* Decide if the n! theorem invalidates a given order. *)
let invalid_by_n_fact n =
  List.exists
    (fun (p, orders) -> (List.length orders) = 0)
    (norm_orders n) 

(* Invalidate orders by a permutation group result that determines
 * that no simple group can have an order that is the product of 2 and
 * an odd number. *)
let is_prod_2odd n =
  if divs n 2 then is_odd (n / 2)
  else false

(* A complete checker to determine if there is a simple group of the
 * input size. *)
let may_exist_simple_grp_order n =
  let checkers = 
    [
     is_prime;
     is_prod_paq;
     is_prod_pqr; 
     is_prod_2odd;
     invalid_by_n_fact;
    ] in
    List.for_all (fun checker -> not (checker n)) checkers

(* Iterate over all numbers up to the upper bound. *)
let rec num_driver upper_bound n =
  if n >= upper_bound then []
  else
    let succs = num_driver upper_bound (n + 1) in
      if (may_exist_simple_grp_order n) then n::succs else succs

(* Top-level function. *)
let main () =
  let upper_bound = 
    if (Array.length Sys.argv) = 1 then 1000 
    else int_of_string Sys.argv.(1) in
  let candidates = num_driver upper_bound 2 in
  let res_str = 
    List.fold_left
      (fun acc n ->
         let num_str = 
           (string_of_int n) 
           ^ " = " 
           ^ (str_of_prime_fact (prime_factors n))
           ^ "\n" in
         let p_str =
           List.fold_left
             (fun res (p, orders) ->
                let order_str = str_of_intl orders in
                let norm_str = 
                  "\tLet P a Syl " 
                  ^ (string_of_int p) 
                  ^ "-subgp, |N_G(P)| in { " 
                  ^ order_str 
                  ^ " }\n" in
                  res ^ norm_str)
             ""
             (norm_orders n) in
           acc ^ num_str ^ p_str)
      ("There may exist simple groups of order less than "
       ^ (string_of_int upper_bound)
       ^ " of the following orders:\n")
      candidates in
    print_string res_str;;

(* Enter the program. *)
main () ;;