(** Generating primes *)

open Builtin
open Basic_arithmetics

(* Initializing list of integers for eratosthenes's sieve. Naive
   version.
*)

(** List composed of 2 and then odd integers starting at 3.
    @param n number of elements in the list of integers.
 *)
let init_eratosthenes n =
  let rec init_eratosthenes_rec i =
  match i with
      i when i > n -> []
    | i when modulo i 2 = 1 -> i::init_eratosthenes_rec (i+2)
    | i -> init_eratosthenes_rec (i+1)
  in if n > 1 then
      2::init_eratosthenes_rec 3
    else [];;

(* Eratosthenes sieve. *)

(** Eratosthene sieve.
    @param n limit of list of primes, starting at 2.
*)
let rec remove_multiples_from_list e l =
  match l with
      [] -> []
    | e1::l1 -> begin
      match modulo e1 e with
          0 -> remove_multiples_from_list e l1
        | _ -> e1::remove_multiples_from_list e l1
    end;;

let eratosthenes n =
  let list = init_eratosthenes n in
  let rec eratosthenes_rec l =
    match l with
        [] -> []
      | e::l1 -> begin
        let l1 = remove_multiples_from_list e l1 in
        e::eratosthenes_rec l1
      end
  in eratosthenes_rec list

(* Write and read into file functions for lists. *)

(** Write a list into a file. Element seperator is newline.
    @param file path to write to.
 *)
let write_list li file =
  let oc = open_out file in
  let rec write l =
    match l with
        [] -> close_out oc
      | e::l1 -> Printf.fprintf oc "%d\n" e; write l1
  in write li;;

(** Write a list of prime numbers up to limit into a txt file.
    @param n limit of prime numbers up to which to build up a list of primes.
    @param file path to write to.
*)
let write_list_primes n file =
  let list = eratosthenes n in
  write_list list file;;

(** Read file safely ; catch End_of_file exception.
    @param in_c input channel.
 *)
let input_line_opt in_c =
  try Some (input_line in_c)
  with End_of_file -> None

(** Create a list out of reading a line per line channel.
    @param in_c input channel.
 *)
let create_list in_c =
  let rec read acc =
    match input_line_opt in_c with
        None -> acc
      | Some l -> l::read acc
  in read [];;

(** Load list of primes into OCaml environment.
    @param file path to load from.
 *)
let read_list_primes file =
  let list = create_list (open_in file) in
  let rec put list =
    match list with
        [] -> []
      | e::l -> (int_of_string e)::put l
  in put list;;

(* Auxiliary functions to extract big prime numbers for testing
   purposes.
 *)

(** Get biggest prime.
    @param l list of prime numbers.
 *)
let rec last_element l = match l with
  | [] -> failwith "Builtin.generate_primes.last_element: Your list \
                    is empty. "
  | e::[] -> e
  | h::t -> last_element t

(** Get two biggest primes.
    @param l list of prime numbers.
 *)
let rec last_two l = match l with
  | [] | [_] -> failwith "Builtin.generate_primes.last_two: List has \
                          to have at least two prime numbers."
  | e::g::[] -> (e, g)
  | h::t -> last_two t
;;

(* Generating couples of prime numbers for specific or fun
   purposes.
 *)

(** Finding couples of primes where second entry is twice the first
    plus 1.
    @param limit positive integer bounding searched for primes.
    @param isprime function testing for (pseudo)primality.
*)


let double_primes limit isprime =
  let list = eratosthenes limit in
  let rec double_primes_rec l =
    match l with
        [] -> []
      | e::l1 -> begin
        let double = e * 2 + 1 in
        if isprime double then
          (e, double)::double_primes_rec l1
        else
          double_primes_rec l1
      end
  in double_primes_rec list;;

(** Finding twin primes.
    @param limit positive integer bounding searched for primes.
    @param isprime function testing for (pseudo)primality.
 *)
let twin_primes limit isprime =
  let list = eratosthenes limit in
  let rec double_primes_rec l =
    match l with
        [] -> []
      | e::l1 -> begin
        let double = e + 2 in
        if isprime double then
          (e, double)::double_primes_rec l1
        else
          double_primes_rec l1
      end
  in double_primes_rec list;;