function  H = graystep(d, B)
%  H = graystep(d, B)  takes bit-row  B  (every element is 0 or 1)
%   through one step forward,  if  d > 0 ,  or otherwise backward,
%   through the  Gray Cyclic Binary Codes.  Only one element of  H
%   differs from the corresponding element of  BĘ;  and repeated 
%   invocations  B = graystep(d, B)  step through all  2^length(B)
%   different  Gray  codes represented as bit-rows.  Each step 
%   takes time proportional to  length(B)  but far less than,  say,
%       H = int2btr( graynext(d, btr2int(B), length(B)) ) .
%  Actually  H == graystep(d, (B(:)' ~= 0)) ,  as if array  B  had
%   been converted to a bit-row whose every nonzero element is  1 .
%  See also  grays,  btr2gray,  gray2btr,  graynext,  grayndcs.

%  Graystep  was adapted from  J. Boothroyd's Algorithm #246  on
%  p. 701 of Comm. ACM vol. 7 (1964)  by  W. Kahan,  14 Feb. 2009  

H = (B(:)' ~= 0) ;  n = length(H) ;  d = (d > 0) ;
B = cumsum(H) ;  d = abs( (floor(B(n)*0.5)*2 - B(n)) + d )  ;
j = n+1 - sum(B ~= 0) ;  %  H(j) is first nonzero element if any.
if  d ,  H(1) = ~H(1) ;
  else
    if  (j < n) ,  H(j+1) = ~H(j+1) ;  else  H(n) = ~H(n) ;  end  
  end