PRIVATE FUNC t = EstimateT0(p, abcd, bcd) IS 
  (E x, y :: 
    p = (x, y) REL (abcd, bcd) AND t = 0.5 + x / 6) 
END;

PRIVATE FUNC p = FHat(t, f0, f1, f2) IS 
  p = 
    R2.Plus(f0, 
            R2.Plus(R2.Times(t, f1), 
                    R2.Times(t * t, f2))) 
END;

PRIVATE CONST Steps = 10;

PRIVATE PROC DrawFHat(f0, f1, f2) IS 
  VAR t = 0, cnt = Steps, delta = 1 / Steps IN 
    PS.MoveTo(FHat(t, f0, f1, f2)); 
    DO 
      cnt > 0 -> 
        t := t + delta; 
        PS.LineTo(FHat(t, f0, f1, f2)); 
        cnt := cnt - 1 
    OD 
  END 
END;

PROC DrawApprox(a, b, c, d) IS 
  IF 
    VAR 
      ab = (0.5, 0) REL (a, b), 
      bc = (0.5, 0) REL (b, c), 
      cd = (0.5, 0) REL (c, d), 
      abc = (0.5, 0) REL (ab, bc), 
      bcd = (0.5, 0) REL (bc, cd), 
      abcd = (0.5, 0) REL (abc, bcd), 
      alpha = R2.Minus(bcd, abcd), 
      beta = R2.Minus(cd, bcd), 
      delta = R2.Minus(beta, alpha), 
      gamma = R2.Minus(delta, alpha), 
      a0 = R2.Plus(abcd, R2.Times(3, gamma)), 
      delta2 = R2.Times(2, delta), 
      a1 = R2.Times(6, R2.Minus(alpha, delta2)), 
      a2 = R2.Times(12, delta) 
    IN 
      DrawFHat(a0, a1, a2) 
    END 
  FI 
END;

PRIVATE FUNC c = C0(a0, a1, p) IS 
  c = 2 * a1 * (a0 - p) 
END;

PRIVATE FUNC c = C1(a0, a1, a2, p) IS 
  c = 2 * (a1 + 2 * a2 * (a0 - p)) 
END;

PRIVATE FUNC c = C2(a1, a2) IS 
  c = 6 * a1 * a2 
END;

PRIVATE FUNC c = C3(a2) IS c = 4 * a2 * a2 END;

PRIVATE FUNC t = EstimateT(p, a, b, c, d) IS 
  (E ab = (0.5, 0) REL (a, b), 
     bc = (0.5, 0) REL (b, c), 
     cd = (0.5, 0) REL (c, d), 
     abc = (0.5, 0) REL (ab, bc), 
     bcd = (0.5, 0) REL (bc, cd), 
     abcd = (0.5, 0) REL (abc, bcd), 
     alpha = R2.Minus(bcd, abcd), 
     beta = R2.Minus(cd, bcd), 
     delta = R2.Minus(beta, alpha), 
     gamma = R2.Minus(delta, alpha), 
     a0 = R2.Plus(abcd, R2.Times(3, gamma)), 
     delta2 = R2.Times(2, delta), 
     a1 = R2.Times(6, R2.Minus(alpha, delta2)), 
     a2 = R2.Times(12, delta), a0x, a0y, a1x, 
     a1y, a2x, a2y, px, py, 
     c0 = C0(a0x, a1x, px) + C0(a0y, a1y, py), 
     c1 
   = C1(a0x, a1x, a2x, px) + C1(a0y, a1y, a2y, py), 
     c2 = C2(a1x, a2x) + C2(a1y, a2y), 
     c3 = C3(a2x) + C3(a2y), t2 = t * t :: 
    a0 = (a0x, a0y) AND 
    a1 = (a1x, a1y) AND 
    a2 = (a2x, a2y) AND 
    p = (px, py) AND 
    t ~ EstimateT0(p, abcd, bcd) AND 
    0 = c3 * t2 * t + c2 * t2 + c1 * t + c0) 
END;

PRED On(p, a, b, c, d) IS 
  (E t ~ EstimateT(p, a, b, c, d), 
     ab ~ (t, 0) REL (a, b), 
     bc ~ (t, 0) REL (b, c), 
     cd ~ (t, 0) REL (c, d), 
     abc ~ (t, 0) REL (ab, bc), 
     bcd ~ (t, 0) REL (bc, cd) :: 
    ab = (t, 0) REL (a, b) AND 
    bc = (t, 0) REL (b, c) AND 
    cd = (t, 0) REL (c, d) AND 
    abc = (t, 0) REL (ab, bc) AND 
    bcd = (t, 0) REL (bc, cd) AND 
    p = (t, 0) REL (abc, bcd)) 
END;

PROC ShowT(p, a, b, c, d) IS 
  VAR t = EstimateT(p, a, b, c, d) IN 
    PS.Type(p, Unparse.Value(t)) 
  END | PS.Type(p, "None") 
END;

PRIVATE PROC res := DPrime(t, p, a, b, c, d) IS 
  IF 
    VAR 
      ab = (0.5, 0) REL (a, b), 
      bc = (0.5, 0) REL (b, c), 
      cd = (0.5, 0) REL (c, d), 
      abc = (0.5, 0) REL (ab, bc), 
      bcd = (0.5, 0) REL (bc, cd), 
      abcd = (0.5, 0) REL (abc, bcd), 
      alpha = R2.Minus(bcd, abcd), 
      beta = R2.Minus(cd, bcd), 
      delta = R2.Minus(beta, alpha), 
      gamma = R2.Minus(delta, alpha), 
      a0 = R2.Plus(abcd, R2.Times(3, gamma)), 
      delta2 = R2.Times(2, delta), 
      a1 = R2.Times(6, R2.Minus(alpha, delta2)), 
      a2 = R2.Times(12, delta), a0x, a0y, a1x, 
      a1y, a2x, a2y, px, py, 
      c0 = C0(a0x, a1x, px) + C0(a0y, a1y, py), 
      c1 
    = C1(a0x, a1x, a2x, px) + 
      C1(a0y, a1y, a2y, py), 
      c2 = C2(a1x, a2x) + C2(a1y, a2y), 
      c3 = C3(a2x) + C3(a2y), t2 = t * t 
    IN 
      a0 = (a0x, a0y) AND 
      a1 = (a1x, a1y) AND 
      a2 = (a2x, a2y) AND 
      p = (px, py) -> 
        res := c3 * t2 * t + c2 * t2 + c1 * t + c0 
    END 
  FI 
END;

PRIVATE FUNC b = B0(a0, p) IS 
  (E diff = (a0 - p) :: b = diff * diff) 
END;

PRIVATE FUNC b = B1(a0, a1, p) IS 
  b = 2 * a1 * (a0 - p) 
END;

PRIVATE FUNC b = B2(a0, a1, a2, p) IS 
  b = a1 + 2 * a2 * (a0 - p) 
END;

PRIVATE FUNC b = B3(a1, a2) IS 
  b = 2 * a1 * a2 
END;

PRIVATE FUNC b = B4(a2) IS b = a2 * a2 END;

PRIVATE PROC res := Dist(t, p, a, b, c, d) IS 
  IF 
    VAR 
      ab = (0.5, 0) REL (a, b), 
      bc = (0.5, 0) REL (b, c), 
      cd = (0.5, 0) REL (c, d), 
      abc = (0.5, 0) REL (ab, bc), 
      bcd = (0.5, 0) REL (bc, cd), 
      abcd = (0.5, 0) REL (abc, bcd), 
      alpha = R2.Minus(bcd, abcd), 
      beta = R2.Minus(cd, bcd), 
      delta = R2.Minus(beta, alpha), 
      gamma = R2.Minus(delta, alpha), 
      a0 = R2.Plus(abcd, R2.Times(3, gamma)), 
      delta2 = R2.Times(2, delta), 
      a1 = R2.Times(6, R2.Minus(alpha, delta2)), 
      a2 = R2.Times(12, delta), a0x, a0y, a1x, 
      a1y, a2x, a2y, px, py, 
      b0 = B0(a0x, px) + B0(a0y, py), 
      b1 = B1(a0x, a1x, px) + B1(a0y, a1y, py), 
      b2 
    = B2(a0x, a1x, a2x, px) + 
      B2(a0y, a1y, a2y, py), 
      b3 = B3(a1x, a2x) + B3(a1y, a2y), 
      b4 = B4(a2x) + B4(a2y), t2 = t * t 
    IN 
      a0 = (a0x, a0y) AND 
      a1 = (a1x, a1y) AND 
      a2 = (a2x, a2y) AND 
      p = (px, py) -> 
        res := 
          b4 * t2 * t2 + b3 * t2 * t + b2 * t2 + 
          b1 * t + b0 
    END 
  FI 
END;

PRIVATE PROC 
  res := PlotP(t, p, a, b, c, d, u, v) IS 
  res := 
    (t, Dist(t, p, a, b, c, d) / 50) REL (u, v) 
END;

PROC Plot(p, a, b, c, d, u, v) IS 
  VAR t = 0 IN 
    PS.MoveTo(PlotP(t, p, a, b, c, d, u, v)); 
    DO 
      t <= 1 -> 
        PS.LineTo(PlotP(t, p, a, b, c, d, u, v)); 
        t := t + 0.1 
    OD 
  END 
END;