function [mk,blaises,mp1,pmp1] = mpnext(blaises,m,pm)
%MPNEXT next section of MM
%
%        [MK,BLAISES,MP1,PMP1] = MPNEXT(BLAISES,M,PM)
%
% input:  BLAISES := ( (k+h-1 choose h-1) : h=1:d+1 ), hence k+1 = BLAISES(2);
%               M := the last column of the matrix containing in its rows the
%                     elements of A_k(d) := { a in Z_+^d : |a| = k } in 
%                     lexicographic order;
%              PM := the result of  [ignored, PM] = SORT(M) ;
% output:  MK  for which  MK(i,n(a)) = n(a+e_i), i=1:d , with  n(a)  the place
%             at which  a  occurs in the lexicographic ordering of  A_{|a|}(d);
%          as well as  BLAISES , M , PM  for  k+1  rather than  k .
%       
% This is of use in working with multivariate polynomials in power form.
%
%   See also: MPMAK, MPAPI

% cb: 9aug97
%   Copyright (c) 1997 by C. de Boor

if nargin==1  % initialize the process
   d = blaises(1);
   blaises = ones(1,d+1);
   mk = zeros(d,0);
   mp1 = [0]; pmp1 = [1];
else          % update it
   d = length(blaises)-1;
   % Ready to generate, in MP1, the last column of the matrix which has
   % in its rows the elements of A_{k+1}(d), in lexicographic order. This
   % column starts off with  k+1 , followed by increasingly longer initial
   % segments of the last column of the same matrix for  A_k(d) , presumably
   % available in  M .
   mp1 = [blaises(2)]; % recall that  blaises(2) equals  k+1 
   for jj=2:d
      mp1 =[mp1 m(1:blaises(jj))];
   end
   [ignored,pmp1] = sort(mp1);  %  gets the permutation that orders  MP1 
   nakd = blaises(d);    % = # A_k(d) 
   blaises = cumsum(blaises); % raises BLAISES to the next level
   rangek = (blaises(d)-nakd) + [1:nakd]; % the tail-end of interest
   % Ready to generate  MK  for which MK(i,n(a)) = n(a+e_i), i=1:d, a in A_k(d).
   % This is done as follows. It is easily seen that the map  a |--> a+e_1
   % carries  A_k(d)  without any order change to the tail-end of the ordered
   %  A_{k+1(d) , hence  MK(1,:) = RANGEK := (BLAISES(d)-NAKD) + [1:NAKD] .
   % For  a+e_j , let  Ak  be the matrix whose rows contain the elements of
   %  A_k(d)  in lexicographic order. Then, in particular,  Ak(:,d) = M
   % and  Ak(PM,[d 1:d-1]) = Ak  (since  PM  is the permutation which puts  
   %  M  into increasing order without unnecessary interchanges). Hence  
   %        Ak(PM^{-1},[2:d 1]) = Ak  and  Akp1(PMP1^{-1},[2:d,1]) = Akp1 .
   % Therefore,
   %          MK(2,PM^{-1}) = PMP1^{-1}(RANGEK) ,
   % and, in general
   %    MK(i,PM^{-i+1}) = PMP1^{-i+1}(RANGEK), i=1:d
   % and, of course,  PM^{-d} = id, PMP1^{-d} = id .
   % This is realized below by starting with  NEXT = 1:BLAISES(d), and
   % repeatedly permuting the columns of  MK  by  PM^{-1}  and the columns of
   % NEXT, in concert, by PMP1^{-1}, and then setting the appropriate row of
   % the permuted  MK  equal to the tail-end of the analogously permuted  NEXT .
   
   mm = zeros(d,nakd); next = 1:blaises(d);
   for i=1:d
      mk(i,:) = next(rangek);
      mk(:,pm) = mk; next(pmp1) = next;
   end
end