#!dotnet fsi
// See www.temmel.me/tuin/nearly_equilateral.html

let pi : float = 3.141592653

let max_alpha : float = pi / 12.0

let max_alpha_sine : float = sin max_alpha

let equilateral_angle : float = pi / 3.0

let ceil (x : float) : int =
  x |> ceil |> int

let floor (x : float) : int =
  x |> floor |> int

type point = int * int

let origin : point = (0, 0)

let difference ((a_x, a_y): point) ((b_x, b_y): point) : point =
  (a_x - b_x, a_y - b_y)

let same ((a_x, a_y): point) ((b_x, b_y): point) : bool =
  (a_x = b_x) && (a_y = b_y)

let inner_product ((a_x, a_y): point) ((b_x, b_y): point) : int =
  a_x * b_x + a_y * b_y

let norm_square (A : point) : int =
  inner_product A A

let norm (a: point) : float =
  norm_square a |> float |> sqrt

// See https://en.wikipedia.org/wiki/Inner_product#/media/File:Inner-product-angle.png
let cos_angle (A: point) (B: point) : float =
  let acos angle = System.Math.Acos angle
  let numerator = inner_product A B |> float
  let denumerator = (norm A) * (norm B)
  acos (numerator / denumerator)

let alpha ((a_x, a_y): point) =
  cos_angle (a_x, a_y) (1, 0)

type c_approximation = double * double

let b_to_c_approx (B : point) : c_approximation =
  let alpha = alpha B
  let beta = equilateral_angle - alpha
  let l = norm B
  let B_x = fst B |> float
  let B_y = snd B |> float
  (B_x - l * cos beta, B_y + l * sin beta)

let c_approx_to_cs ((l_x, l_y): c_approximation) : list<point> =
  let d_x = floor l_x
  let u_x = ceil l_x
  let d_y = floor l_y
  let u_y = ceil l_y
  [(d_x, d_y); (d_x, u_y); (u_x, d_y); (u_x, u_y)]

let b_to_cs (B : point) : list<point> =
  let cs = B |> b_to_c_approx |> c_approx_to_cs
  let filter_c C =
    not (same origin C || same B C)
  List.filter filter_c cs

type solution = point * point

let b_to_solutions (B : point) : list<solution> =
  B |> b_to_cs |> List.map (fun C -> (B, C))

let min_max_of_three (a: float, b: float, c: float) : float * float =
  let values = [a; b; c]
  let m = List.min values
  let M = List.max values
  (m, M)

let min_max_diff (m : float, M : float) =
  M - m

let min_max_ratio (m : float, M : float) =
  m / M - 1.0 |> abs

let min_max_diff_of_three (triple : float * float * float) : float =
  min_max_of_three triple |> min_max_diff

let min_max_ratio_of_three (triple : float * float * float) : float =
  min_max_of_three triple |> min_max_ratio

let sides ((B, C) : solution) : float * float * float =
  let a = difference B C |> norm
  let b = C |> norm
  let c = B |> norm
  (a, b, c)

let angles ((B, C) : solution) : float * float * float =
  let alpha = cos_angle B C
  let beta = cos_angle (difference origin B) (difference C B)
  let gamma = cos_angle (difference origin C) (difference B C)
  (alpha, beta, gamma)

let radii (a, b, c) : float * float =
  let s = (a + b + c) / 2.0
  let prod = a * b * c
  let area = sqrt(s * (s - a) + (s - b) * (s - c))
  let r = area / s
  let R = prod / 4.0 / area
  (r, R)

let radius_circumscribed_equilateral : float -> float =
  fun a -> a / (sqrt 3.0)

let error_absolute_sides (s : solution) : float =
  sides s |> min_max_diff_of_three

let error_relative_sides (s : solution) : float =
  sides s |> min_max_ratio_of_three

let error_absolute_angles (s : solution) : float =
  angles s |> min_max_diff_of_three

let error_relative_angles (s : solution) : float =
  angles s |> min_max_ratio_of_three

let error_absolute_radii (s : solution) : float =
  s |> sides |> radii |> min_max_diff

let error_relative_radii (s : solution) : float =
  s |> sides |> radii |> min_max_ratio

let error_absolute_circumcircle (s : solution) : float =
  let (m, M) = s |> sides |> min_max_of_three
  let r = radius_circumscribed_equilateral m
  let R = radius_circumscribed_equilateral M
  R - r

type solution_with_errors =
  {
  Solution : solution
  Absolute_Angles : float
  Absolute_Radii : float
  Absolute_Sides : float
  Relative_Angles : float
  Relative_Radii : float
  Relative_Sides : float
  }

// works only for n in {0, .., 99}
let point_fixed_with ((x,y) : point) : string =
  if x < 0 || 100 < x || y < 0 || 100 < y
  then failwith "Arguments not in [0, 99]."
  else sprintf "(%2i, %2i)" x y

let solution_fixed_with ((B, C) : solution) : string * string =
  (point_fixed_with B, point_fixed_with C)

type swe_fixed_with =
  string * string * string * string * string * string * string * string

let swe_fixed_with (swe : solution_with_errors) : swe_fixed_with =
  let (b, c) = solution_fixed_with swe.Solution
  let error_fw e = sprintf "%.4f" e
  (b, c
  , error_fw swe.Absolute_Angles
  , error_fw swe.Absolute_Radii
  , error_fw swe.Absolute_Sides
  , error_fw swe.Relative_Angles
  , error_fw swe.Relative_Radii
  , error_fw swe.Relative_Sides
  )

let join_swe_fixed_with (sep : string) ((b, c, a_a, a_r, a_s, r_a, r_r, r_s) : swe_fixed_with) : string =
  String.concat sep [b; c; a_a; a_r; a_s; r_a; r_r; r_s]

let swe_fw_bars (swe : solution_with_errors) : string =
  swe |> swe_fixed_with |> join_swe_fixed_with " | "

let swes_fw_bars (swes : list<solution_with_errors>) : string =
  let header : string =
    ("B       ", "C       ", "AbsAng", "AbsRad", "AbsSid", "RelAng", "RelRad", "RelSid")
    |> join_swe_fixed_with " | "
  let lines : list<string>= List.map swe_fw_bars swes
  String.concat "\n" (header::lines)

let swes_cli (swes : list<solution_with_errors>) : string =
  let header =
    ("B       ", "C       ", "AbsAng", "AbsRad", "AbsSid", "RelAng", "RelRad", "RelSid")
    |> join_swe_fixed_with " | "
  let lines = List.map swe_fw_bars swes
  String.concat "\n" (header::lines)

let swes_html (swes : list<solution_with_errors>) : string =
  let tablebars strings =
    "| " + (String.concat " | " strings) + " |"
  let swe_html swe =
    let point (x, y) = sprintf "$(%2i,%2i)$" x y
    let error e = sprintf "$%.4f$  " e
    [ point (fst swe.Solution)
    ; point (snd swe.Solution)
    ; error swe.Absolute_Angles
    ; error swe.Absolute_Radii
    ; error swe.Absolute_Sides
    ; error swe.Relative_Angles
    ; error swe.Relative_Radii
    ; error swe.Relative_Sides
    ]
  let header : string =
    ["$B$      "; "$C$      "; "$A_W(B,C)$"; "$A_C(B,C)$"; "$A_S(B,C)$"; "$R_A(B,C)$"; "$R_C(B,C)$"; "$R_S(B,C)$"]
    |> tablebars
  let lines =
    List.map (fun swe -> swe |> swe_html |> tablebars) swes
  String.concat "\n" (header::lines)

let with_errors (s : solution) =
  { Solution = s
    Absolute_Angles = error_absolute_angles s
    Absolute_Radii = error_absolute_radii s
    Absolute_Sides = error_absolute_sides s
    Relative_Angles = error_relative_angles s
    Relative_Radii = error_relative_radii s
    Relative_Sides = error_relative_sides s
  }

let bs_for_n (n : int) : list<point> =
  let y_max : int = floor ((float n) * max_alpha_sine)
  [ 0..y_max ]
  |> List.map (fun y -> (n, y))

let bs_up_to_n (n: int) : list<point> =
  seq {1..n}
  |> Seq.map bs_for_n
  |> Seq.toList
  |> List.concat

let solutions_up_to (n: int) : list<solution> =
  let bs = bs_up_to_n n
  bs
  |> List.map b_to_solutions
  |> List.concat
  |> Set.ofList
  |> Set.toList

let take n l =
  l |> Seq.truncate n |> Seq.toList

let pr_nl o =
  printfn "%A" o

let pr_swe_fw_nl swes =
  swes |> swes_fw_bars |> printfn "%s"

printfn "---------------------------------"
printfn "Errors of (2, 0), (1, 2) solution."
with_errors ((2, 0), (1, 2)) |> swe_fw_bars |> pr_nl

printfn "---------------------------------"
printfn "Errors of (4, 1), (1, 4) solution."
with_errors ((4, 1), (1, 4)) |> swe_fw_bars |> pr_nl

printfn "---------------------------------"
printfn "Bs up to n = 10."
bs_up_to_n 10 |> pr_nl

printfn "---------------------------------"
printfn "C approximations up to n = 10."
bs_up_to_n 10 |> List.map (fun B -> (B, b_to_c_approx B)) |> pr_nl

printfn "---------------------------------"
printfn "Solutions up to n = 10."
solutions_up_to 10 |> pr_nl

let sol_10 = solutions_up_to 10
let rated_10 = sol_10 |> List.map with_errors
let abs_angles_10 = rated_10 |> List.sortBy (fun swe -> swe.Absolute_Angles)
let abs_sides_10 = rated_10 |> List.sortBy (fun swe -> swe.Absolute_Sides)

printfn "---------------------------------"
printfn "Angles up to n = 10."
abs_angles_10 |> pr_swe_fw_nl

printfn "---------------------------------"
printfn "Sides up to n = 10."
abs_sides_10 |> pr_swe_fw_nl

printfn "---------------------------------"
printfn "Top 7 of angles up to n = 10."
abs_angles_10 |> take 7 |> pr_swe_fw_nl

printfn "---------------------------------"
printfn "Top 7 of sides up to n = 10."
abs_sides_10 |> take 7 |> pr_swe_fw_nl

printfn "---------------------------------"
let sol_50 : list<solution> = solutions_up_to 50
let rated_50 = sol_50 |> List.map with_errors
let abs_angles_50 = rated_50 |> List.sortBy (fun swe -> swe.Absolute_Angles)
let abs_radii_50 = rated_50 |> List.sortBy (fun swe -> swe.Absolute_Radii)
let abs_sides_50 = rated_50 |> List.sortBy (fun swe -> swe.Absolute_Sides)
let rel_angle_50 = rated_50 |> List.sortBy (fun swe -> swe.Relative_Angles)
let rel_radii_50 = rated_50 |> List.sortBy (fun swe -> swe.Relative_Radii)
let rel_sides_50 = rated_50 |> List.sortBy (fun swe -> swe.Relative_Sides)

printfn "Top 3 of absolute angles up to n = 50."
abs_angles_50 |> take 3 |> pr_swe_fw_nl

printfn "Top 3 of absolute radii up to n = 50."
abs_radii_50 |> take 3 |> pr_swe_fw_nl

printfn "Top 3 of absolute sides up to n = 50."
abs_sides_50 |> take 3 |> pr_swe_fw_nl

printfn "Top 3 of relative angles up to n = 50."
rel_angle_50 |> take 3 |> pr_swe_fw_nl

printfn "Top 3 of relative radii up to n = 50."
rel_radii_50 |> take 3 |> pr_swe_fw_nl

printfn "Top 3 of relative sides up to n = 50."
rel_sides_50 |> take 3 |> pr_swe_fw_nl

rel_sides_50 |> take 3 |> swes_html |> printf "%s\n"

sol_50 |> List.map (fun s -> (s, error_absolute_circumcircle s)) |> List.sortBy snd |> take 5 |> List.map pr_nl

let sol_41_sym = ((4, 1), (1, 4))
(sol_41_sym, error_absolute_circumcircle sol_41_sym) |> pr_nl