(**************************************************************************
*  The de Casteljau algorithm.
*                                                                         *
*     Version for polynomial curves  (11/25/2002                          *
*                                                                         *
*  Control polygon     cpoly = (b0,...,bm)
*  Affine frame on real line: (r, s)
*  pdecas[cpoly, r, s, t]: Computation of the current point F(t) 
*  To display pdecas, use showpic2D
*  poldecas[cpoly, r, s, l, m, n]: computes n + 1 points on the curve,
*  starting with t = l, and ending with t = m
*  To display poldecas, use showcur2D
*  subdiv[cpoly, r, s, n]: perform n steps  of subdivision (using t = 1/2)
*  To display subdiv, use showsub2D
***************************************************************************)

lerp[{p1__},{p2__},r_,s_,t_] := 
 Block[
 {p = {p1}, q = {p2}, res},
 res = (s - t)/(s - r)  p + (t - r)/(s - r) q;
 res
 ];

decas[{cpoly__}, r_, s_, t_] :=
 Block[
 {bb = {cpoly}, b = {},
 m, i, j, lseg = {}, res},
 (m = Length[bb] - 1;
 Do[
    Do[
       b = Append[b, lerp[bb[[i]], bb[[i+1]], r, s, t]];
       If[i > 1, lseg = Append[lseg, Line[{b[[i - 1]], b[[i]]}]]]
       , {i, 1, m - j + 1}
      ]; bb = b; b = {}, {j, 1, m}
   ];
 res := Append[lseg, bb[[1]]];
 res
 )
 ];

pdecas[{cpoly__}, r_, s_, t_] :=
 Block[
 {pt, ll, res},
 res = decas[{cpoly}, r, s, t];
 pt = Last[res]; ll = Drop[res, -1];
 res = Append[ll, {PointSize[0.01], Point[pt]}];
 res
 ];


contpoly[{cpoly__}] :=
 Block[
  {bb = {cpoly}, lst = {}, m, i},
   m = Length[bb] - 1;
   Do[lst = Append[lst, Line[{bb[[i - 1]], bb[[i]]}]], {i, 2, m + 1}
     ];
  lst
 ];

showpic2D[{poly__},r_,s_,lx_,mx_,ly_,my_,t_] :=   Show[
                       Graphics[{contpoly[{poly}], 
                       {Thickness[0.001], pdecas[{poly}, r, s, t]}}],
                       AspectRatio -> Automatic,
                       PlotRange -> {{lx, mx}, {ly, my}},
                       DisplayFunction -> $DisplayFunction];


poldecas[{cpoly__}, r_, s_, l_, m_, n_] :=
 Block[
 {bb = {cpoly}, lpt = {},
  i, lseg = {}, res, t},
 (t = l;
 Do[
    res = decas[bb, r, s, t]; lpt = Append[lpt, Last[res]];
    If[i > 1, lseg = Append[lseg, Line[{lpt[[i-1]], lpt[[i]]}]]]; 
    t = t + (m - l)/n, {i, 1, n + 1} 
   ];
 res = lseg;
 res
 )
 ];



subdecas[{cpoly__}, r_, s_, t_] :=
 Block[
 {bb = {cpoly}, b = {}, ud = {}, ld = {},
 m, i, j, lseg = {}, res},
 (m = Length[bb] - 1; ud = {bb[[1]]}; ld = {bb[[m + 1]]};
 Do[
    Do[
       b = Append[b, lerp[bb[[i]], bb[[i+1]], r, s, t]], {i, 1, m - j + 1}
      ];   
      ud = Append[ud, b[[1]]];
      ld = Prepend[ld, b[[m - j + 1]]];
      bb = b; b = {}, {j, 1, m} 
   ];
 res := Join[{ud},{ld}];
 res
 )
 ];


subdivstep[{poly__}, r_, s_] :=
 Block[
 {cpoly = {poly}, lpoly = {},
  t, l, i, res},
 (l = Length[cpoly]; t = (r + s)/2;
 Do[
    lpoly = Join[lpoly, subdecas[cpoly[[i]], r, s, t]] , {i, 1, l} 
   ];
 lpoly
 )
 ];

subdiv[{poly__}, r_, s_, n_] :=
 Block[
 {lpoly = {{poly}}, res = {}, sl,
  i, j, l1, l2},
 (
 Do[
    lpoly = subdivstep[lpoly, r, s], {i, 1, n} 
   ];
 l1 = Length[lpoly];
 Do[
    sl = lpoly[[i]]; l2 = Length[sl];            
    Do[If[j > 1, res = Append[res, Line[{sl[[j-1]], sl[[j]]}]]], {j, 1, l2}
      ], {i, 1, l1}
   ];
 res
 )
 ];


showcur2D[{poly__},r_,s_,lx_,mx_,ly_,my_,l_,m_,n_] :=   Show[
                       Graphics[{contpoly[{poly}], 
                       poldecas[{poly}, r, s, l, m, n]}],
                       AspectRatio -> Automatic,
                       PlotRange -> {{lx, mx}, {ly, my}},
                       DisplayFunction -> $DisplayFunction];

showsub2D[{poly__},r_,s_,lx_,mx_,ly_,my_,n_] :=   Show[
                       Graphics[{contpoly[{poly}], 
                       subdiv[{poly}, r, s, n]}],
                       AspectRatio -> Automatic,
                       PlotRange -> {{lx, mx}, {ly, my}},
                       DisplayFunction -> $DisplayFunction];