%!

% t y=[ y0 y1 ... yn ]
/bernstein { % t y 
  % constants y n t s=1-t
  % variables k C P
  dup length % t y n+1
  1 sub      % t y n
  3 -1 roll 1 % y n t 1
  1 index sub % y n t s
  % constants in place
  1           % y n t s k
  3 index 3 index mul % y n t s k C=nt
  5 index 0 get       % y n t s k C P=y0
  5 index {           % y n t s k C P
    % P -> s.P + C.y[k]
    % C -> C.t.(n-k)/(k+1) 
    % k -> k+1
    3 index mul       % y n t s k C P.s
    1 index           % y n t s k C P.s C
    7 index           % y n t s k C P.s C y 
    4 index get mul add  % y n t s k C P.s+C.y[k]=new P
    3 1 roll          % y n t s P* k C
    5 index           % y n t s P* k C n
    2 index sub mul   % y n t s P* k C.(n-k)
    1 index 1 add div % y n t s P* k C.(n-k)/(k+1)
    4 index mul       % y n t s P* k C*
    3 1 roll 1 add    % y n t s C* P* k*
    3 1 roll          % y n t s k* C* P*
  } repeat
  7 1 roll 6 { pop } repeat
} def