jdh30

let branches m =
  { { Vec2.rotate(π/2 - asin(4/5));
      Vec2.scale(4/5, 4/5);
      Vec2.translate(0, -1) };
    { Vec2.translate(1, -1);
      Vec2.rotate(-π/2 + asin(3/5));
      Vec2.scale(3/5, 3/5);
      Vec2.translate(1/3, 0) } }
  @ Array.map [m2 → Array.concat{m; m2}]

jdh30

let toHtml f =
  let x = Html "𝑥" in
  let pow = (* xⁿ *)
    [ 0 → {Html "1"}
    | 1 → {x}
    | n → {x; Html.tag "sup" {"style", "font-size:70%"} (Html.ofNumber n)} ] in
  let append(xs, x) = Array.Unsafe.append xs x in
  let appends(xs, xss) = Array.fold append xs xss in
  let minus = Html.of " - " in
  let ellipses = Html " + …" in

jdh30

let safe (x1, y1) (x2, y2) =
  x1 <> x2 && y1 <> y2 && x2 - x1 <> y2 - y1 && x1 - y2 <> x2 - y1

let rec next n s =
  Stack.pop s @
  [ None -> None
  | Some((qs, ps), s) ->
      Stack.pop ps @
      [ None -> if Stack.length qs = n then Some(qs, s) else next n s
      | Some(q, ps) ->

jdh30

let pearsonsCorrelationCoefficient xys =
  let n = Array.length xys in
  let ∑xy, ∑x, ∑y, ∑x2, ∑y2 =
    Array.fold [(∑xy, ∑x, ∑y, ∑x2, ∑y2), (x, y) →
      ∑xy+x*y, ∑x+x, ∑y+y, ∑x2+x*x, ∑y2+y*y]
      (0, 0, 0, 0, 0) xys in
  (n*∑xy-∑x*∑y)/(√(n*∑x2-∑x*∑x)*√(n*∑y2-∑y*∑y))

jdh30

let primes =
  let a = {2} in
  let grow() =
    let p0 = Array.get a (Array.length a - 1) + 1 in
    let b = Array.init p0 [_ → True] in
    let () =
      Array.fold [(), di →
        let rec loop i =
          if i >= Array.length b then () else
            let () = Array.Unsafe.set b i False in

jdh30

let primeFactors x =
  let rec loop fs c p =
    if p < c * c then Array.Unsafe.append fs p else
      if mod(p, c) = 0 then
        loop (Array.Unsafe.append fs c) c (p / c)
      else loop fs (c + 1) p in
  loop {} 2 x
  

jdh30

let rec fixedPoint f x =
  let fx = f x in
  if fx = x then x else fixedPoint f fx

let ∂ fxs f xs (i, xi) =
  let () = xi + δ @ Array.Unsafe.set xs i in
  let f2 = (f xs - fxs) / δ in
  let () = xi @ Array.Unsafe.set xs i in
  f2

jdh30

let integrate (x0, x1) f =
  let integrateAux f (x0, x2) =
    let δ = ε in
    let simpsonsRule dx fx0 fx1 fx2 = dx * (fx0 + 4*fx1 + fx2) / 6 in
    let rec aux x0 fx0 x2 fx2 x4 fx4 =
      let dx = x2 - x0 in
      let x1 = (x0 + x2) / 2 in
      let x3 = (x2 + x4) / 2 in
      if x0 = x1 || x2 = x3 then 0 else
        let fx1 = f x1 in

jdh30

let lagrange xys =
  let n = # xys in
  let xys = Array.sortBy [x, _ → x] xys in
  let lagrange j x =
    for 0 1 (n-1) 1 [y, i →
      if i = j then y else
        let xi, _ = Array.get xys i in
        let xj, _ = Array.get xys j in
        y * (x - xi) / (xj - xi)] in
  [ x →

jdh30

let δ = ∛ ε

let d f x = (f(x+δ) - f(x-δ))/(2*δ)

let rec fixedPoint f x =
  let fx = f x in
  if fx = x then x else fixedPoint f fx

let newtonRaphson f df x = fixedPoint [x → x - f x / df x] x