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 * 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 FUNC y = Sqrt(x) IS 
  y ~ 1 AND y * y = x 
END;

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 := 
          Sqrt(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) / 1000) REL 
    (u, v) 
END;

PRIVATE CONST Steps = 10;

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