Infinite power series

InfinitePowerSeries.tq

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
  let rec loop needsSign n xs =
    let plus = if needsSign then {Html.of " + "} else {} in
    let recur xss = loop True (n+1) (Array.fold appends xs xss) in
    let recurp p q t = recur{plus; {Q(p, q) @ Html.ofRational; Html " "}; t} in
    let recurn p q t = recur{{minus; Q(-p, q) @ Html.ofRational; Html " "}; t} in
    if n > 8 then append(xs, ellipses) else
      (f n, n) @
      [ Q(0, 1), n → loop needsSign (n+1) xs (* 0 + … *)
      | Q(p, q), 0 → if p >= 0 then recurp p q {} else recurn p q {}
      | Q(1, 1), 1 → recur{plus; {x}}
      | Q(p, q), 1 → if p >= 0 then recurp p q {x} else recurn p q {x}
      | Q(1, 1), n → recur {plus; pow n}
      | Q(p, q), n →
          if p >= 0 then recurp p q (pow n) else recurn p q (pow n) ] in
  loop False 0 {} @ Html.concat

let print str series =
  {String.concat{str; " = "} @ Html.of; toHtml series}
  @ Html.concat
  @ Html.div {}
  @ yield

let constant k = [0 → k | _ → Q.z 0]
let one = constant (Q.z 1)
let x = [1 → Q.z 1 | _ → Q.z 0]
let rest f n = f(n+1)
let xTimes f = [0 → Q.z 0 | n → f(n-1)]
let add f g n = Q.add (f n) (g n)
let sub f g n = Q.sub (f n) (g n)
let cons(x, xs) n = if n=0 then x else xs(n-1)
let differentiate f n = Q.mul (Q.z(n+1)) (f(n+1))
let integrate f n = if n=0 then Q.z 0 else Q.div (f(n-1)) (Q.z n)
let scale x xs n = Q.mul x (xs n)

let rec mul f g n =
  Q.add
    (Q.mul (f 0) (g n))
    (xTimes(mul (rest f) g) n)

let rec div f g = memoize [divfg →
  cons(Q.div (f 0) (g 0),
       scale
        (Q.div (Q.z 1) (g 0))
           (sub (rest f) (mul divfg (rest g))))]

let logOnePlusX = integrate(div one (add x one))

let o f g =
  let rec o f g n = cons(f 0, mul (rest g) (o (rest f) g)) n in
  o f g

let inverse f =
  if f 0 <> Q.zero then
    panic "Can only invert series with a zero constant"
  else
    memoize [inversef → cons(Q.z 0, div one (o (rest f) inversef))]

let exp = add (inverse logOnePlusX) one
let atan = integrate(div one (add one (xTimes x)))
let tan = inverse atan
let sin = memoize [sin → integrate (sub one (integrate sin))]
let cos = differentiate sin