factorization.ml

let factors max_n =
  let open Bigarray in
  let size = (max_n - 3) / 2 in
  let f = Array1.create int c_layout (size + 1) in
  for n = 0 to size do Array1.set f n 0 done;
  for n = 0 to size do
    if Array1.get f n = 0 then (
      let nn = 2 * n + 3 in
      Array1.set f n nn;
      let m_times_n = ref (2 * n * n + 6 * n + 3) in
      while !m_times_n <= size do
        Array1.set f !m_times_n nn;
        m_times_n := !m_times_n + nn;
      done)
  done;
  f

let dec f n =
  let rec log_p acc n p =
    if n mod p <> 0 then acc, n
    else log_p (acc + 1) (n / p) p
  in
  let rec aux acc n =
    if n = 1 then acc
    else
      let p = Bigarray.Array1.get f ((n-3)/2) in
      let q, n = log_p 0 n p in
      aux ((p, q) :: acc) n
  in
  let q2, n = log_p 0 n 2 in
  List.sort compare (aux (if q2 = 0 then [] else [2, q2]) n)

let print_dec d =
  let print_pair (p, q) = Printf.printf "(%d, %d) " p q in
  List.iter print_pair d;
  print_newline ()

let () =
  let factors = factors 100_000_000 in
  print_dec (dec factors 12_345);
  print_dec (dec factors 123_456);
  print_dec (dec factors 1_234_567);
  print_dec (dec factors 12_345_678);
  print_dec (dec factors 256);
  print_dec (dec factors 10000);
  print_dec (dec factors 99999)