// Initialization variables:
a = p+1;
b = a+1;
nb = 1;
c1 = 0;
c1 (nt+1) = 1;
mbar = p + 2 + nt + m;
ki = 1;
si = 1;
while b < mbar do
  cnb+1 = 0; // Initialize the next extraction operator row.
  // Count multiplicity of the knot at location b.
  add = 0;
  if si <= m && spans(si) = ki do
    mult = 0;
    add = 1;
    // Add the new knot to the knot vector.
    Ubar(b+1:mbar+p-m+si) = Ubar(b:mbar+p-m+si-1);
    Ubar(b) = U(si);
    si = si + 1;
  else
    ki = ki+1;
    i = b;
    while b < m && Ubar(b+1) == Ubar(b) do b = b+1;
    mult = b-i+1;
  end
  if mult < p do
    numer = Ubar(b)-Ubar(a);
    for j = p,p-1,. . .,mult+1 do
      alphas(j-mult) = numer / (Ubar(a+j+add)-Ubar(a));
    end
    r = p-mult;
    // Update the matrix coefficients for r new knots
    for j=1,2,. . .,r do
      save = r-j+1;
      s = mult+j;
      for k=p+1,p,. . .,s+1 do
      alpha = alphas(k-s);
      cnb (k) = alpha*cnb (k) + (1.0-alpha)*cnb (k-1);
    end
    if b < m do
      // Update overlapping coefficients of the next operator row.
      cnb+1 (save) = cnb (p+1);
    end
  end
  nb = nb + 1; // Finished with the current operator.
  if b < m do
    // Update indices for the next operator.
    a = b;
    b = b+1;
  end
end