end;
   [bits1, temp] = EncodeVector(x1, bits1, HL1, Maxx1, Meanx1);
   [bits2, temp] = EncodeVector(x2, bits2, HL2, Maxx2, Meanx2);
   bits=bits1+bits2+1;
else
   bits=bits+1;      % this is how many bits we are going to write
   if Debug
      disp(['EncodeVector: Level=',int2str(Level),'  ',int2str(L),...
            ' sybols stored in ',int2str(bits),' bits.']);
   end
   if Speed
      % advance Byte and BitPos without writing to y
      Byte=Byte+floor(bits/8);
      BitPos=BitPos-mod(bits,8);
      if (BitPos<=0); BitPos=BitPos+8; Byte=Byte+1; end;
   else
      % put the bits into y
      StartPos=Byte*8-BitPos;     % control variable
      PutBit(0);       % indicate that a sequence is coded
      PutVLIC(L);       
      PutHuffTab(HL);
      HK=huffcode(HL);
      for i=1:L;
         n=x(i)+1;    % symbol number (value 0 is first symbol, symbol 1)
         for k=1:HL(n)
            PutBit(HK(n,k));
         end
      end
      % check if one has used as many bits as calculated
      BitsUsed=Byte*8-BitPos-StartPos;
      if (BitsUsed~=bits)
         disp(['L=',int2str(L),'  max(x)=',int2str(max(x)),'  min(x)=',int2str(min(x))]);
         disp(['BitsUsed=',int2str(BitsUsed),'  bits=',int2str(bits)]);
         error(['Huff06-EncodeVector: Logical error, (BitsUsed~=bits).']); 
      end
   end
end
Level = Level + 1;
return    % end of EncodeVector

function x = DecodeVector
global y Byte BitPos
MaxL=50000;      % as in the EncodeVector function (line 216)
if GetBit
   x1=DecodeVector;
   x2=DecodeVector;
   L=length(x1)+length(x2);
   if (L>MaxL)
      x=[x1(:);x2(:)];
   else
      xm=median([x1;x2]);
      x=zeros(L,1);
      x(1)=x2(1);
      i1=0;i2=1;
      for i=2:L
         if (x(i-1) <= xm) 
            i1=i1+1; x(i)=x1(i1);
         else
            i2=i2+1; x(i)=x2(i2);
         end
      end
   end
else
   L=GetVLIC;
   if (L>1)
      x=zeros(L,1);
      HL=GetHuffTab;
      if length(HL)
         Htree=HuffTree(HL);
         root=1;pos=root;     
         l=0;  % number of symbols decoded so far
         while l<L
            if GetBit
               pos=Htree(pos,3);
            else
               pos=Htree(pos,2);
            end
            if Htree(pos,1)           % we have arrived at a leaf
               l=l+1;
               x(l)=Htree(pos,2)-1;   % value is one less than symbol number
               pos=root;              % start at root again
            end
         end
      else     % HL has length 0, that is empty Huffman table
         x=x+GetVLIC;
      end
   elseif L==0
      if GetBit
         % this is a Run + Value coded sequence
         x1=DecodeVector;
         x2=DecodeVector;
         % now build the actual sequence
         I=x1;      % runs  
         I=I+1;
         L=length(I);  % one more than the number of values in x
         for i=2:L;I(i)=I(i-1)+I(i); end;
         x=zeros(I(L)-1,1);
         x(I(1:(L-1)))=x2;  % values
      else
         x=[];    % this was really a length 0 sequence
      end
   elseif L==1
      x=GetVLIC;
   else
      error('DecodeVector: illegal length of sequence.');
   end
end
return    % end of DecodeVector

% Functions to write and read the Huffman Table Information
% The format is defined in HuffTabLen, we repeat it here
% Function assume that the table information is stored in the following format
%        previous symbol is set to the initial value 2, Prev=2
%        Then we have for each symbol a code word to tell its length 
%         '0'             - same length as previous symbol
%         '10'            - increase length by 1, and 17->1
%         '1100'          - decrease length by 1, and 0->16
%         '11010'         - increase length by 2, and 17->1, 18->2
%         '11011'         - One zero, unused symbol (twice for two zeros)
%         '111xxxx'       - set code length to CL=Prev+x (where 3 <= x <= 14)
%                           and if CL>16; CL=CL-16
%        we have 4 unused 7 bit code words, which we give the meaning
%         '1110000'+4bits - 3-18 zeros
%         '1110001'+8bits - 19-274 zeros, zeros do not change previous value
%         '1110010'+4bits - for CL=17,18,...,32, do not change previous value
%         '1111111'       - End Of Table

function PutHuffTab(HL)
global y Byte BitPos

HL=HL(:);
% if (max(HL) > 32) 
%    disp(['PutHuffTab: To large value in HL, max(HL)=',int2str(max(HL))]); 
% end
% if (min(HL) < 0)
%    disp(['PutHuffTab: To small value in HL, min(HL)=',int2str(min(HL))]); 
% end
Prev=2;
ZeroCount=0;
L=length(HL);

for l=1:L
   if HL(l)==0
      ZeroCount=ZeroCount+1;
   else
      while (ZeroCount > 0)
         if ZeroCount<3
            for i=1:ZeroCount
               PutBit(1);PutBit(1);PutBit(0);PutBit(1);PutBit(1);
            end
            ZeroCount=0; 
         elseif ZeroCount<19 
            PutBit(1);PutBit(1);PutBit(1);PutBit(0);PutBit(0);PutBit(0);PutBit(0);
            for (i=4:-1:1); PutBit(bitget(ZeroCount-3,i)); end;
            ZeroCount=0; 
         elseif ZeroCount<275
            PutBit(1);PutBit(1);PutBit(1);PutBit(0);PutBit(0);PutBit(0);PutBit(1);
            for (i=8:-1:1); PutBit(bitget(ZeroCount-19,i)); end;
            ZeroCount=0; 
         else 
            PutBit(1);PutBit(1);PutBit(1);PutBit(0);PutBit(0);PutBit(0);PutBit(1);
            for (i=8:-1:1); PutBit(1); end;
            ZeroCount=ZeroCount-274; 
         end
      end
      if HL(l)>16
         PutBit(1);PutBit(1);PutBit(1);PutBit(0);PutBit(0);PutBit(1);PutBit(0);
         for (i=4:-1:1); PutBit(bitget(HL(l)-17,i)); end;
      else
         Inc=HL(l)-Prev;
         if Inc<0; Inc=Inc+16; end;
         if (Inc==0) 
            PutBit(0);
         elseif (Inc==1) 
            PutBit(1);PutBit(0);
         elseif (Inc==2) 
            PutBit(1);PutBit(1);PutBit(0);PutBit(1);PutBit(0);
         elseif (Inc==15) 
            PutBit(1);PutBit(1);PutBit(0);PutBit(0);
         else 
            PutBit(1);PutBit(1);PutBit(1);
            for (i=4:-1:1); PutBit(bitget(Inc,i)); end;
         end
         Prev=HL(l);
      end
   end
end
for (i=7:-1:1); PutBit(1); end;       % the EOT codeword

return;  % end of PutHuffTab

function HL=GetHuffTab
global y Byte BitPos

Debug=0;
Prev=2;
ZeroCount=0;
HL=zeros(10000,1);
HLi=0;
EndOfTable=0;

while ~EndOfTable
   if GetBit
      if GetBit
         if GetBit
            Inc=0;
            for (i=1:4); Inc=Inc*2+GetBit; end;
            if Inc==0
               ZeroCount=0;
               for (i=1:4); ZeroCount=ZeroCount*2+GetBit; end;
               HLi=HLi+ZeroCount+3;
            elseif Inc==1
               ZeroCount=0;
               for (i=1:8); ZeroCount=ZeroCount*2+GetBit; end;
               HLi=HLi+ZeroCount+19;
            elseif Inc==2           % HL(l) is large, >16
               HLi=HLi+1;
               HL(HLi)=0;
               for (i=1:4); HL(HLi)=HL(HLi)*2+GetBit; end;
               HL(HLi)=HL(HLi)+17;
            elseif Inc==15
               EndOfTable=1;
            else
               Prev=Prev+Inc;
               if Prev>16; Prev=Prev-16; end;
               HLi=HLi+1;HL(HLi)=Prev;
            end
         else
            if GetBit
               if GetBit
                  HLi=HLi+1;
               else
                  Prev=Prev+2;
                  if Prev>16; Prev=Prev-16; end;
                  HLi=HLi+1;HL(HLi)=Prev;
               end
            else
               Prev=Prev-1;
               if Prev<1; Prev=16; end;
               HLi=HLi+1;HL(HLi)=Prev;
            end
         end
      else
         Prev=Prev+1;
         if Prev>16; Prev=1; end;
         HLi=HLi+1;HL(HLi)=Prev;
      end
   else
      HLi=HLi+1;HL(HLi)=Prev;
   end
end
if HLi>0
   HL=HL(1:HLi);
else
   HL=[];
end

if Debug
   % check if this is a valid Huffman table
   temp=sum(2.^(-nonzeros(HL)));
   if temp ~=1
      error(['GetHuffTab: HL table is no good, temp=',num2str(temp)]);
   end
end

return;  % end of GetHuffTab

% Functions to write and read a Variable Length Integer Code word
% This is a way of coding non-negative integers that uses fewer 
% bits for small integers than for large ones. The scheme is:
%   '00'   +  4 bit  - integers from 0 to 15
%   '01'   +  8 bit  - integers from 16 to 271
%   '10'   + 12 bit  - integers from 272 to 4367
%   '110'  + 16 bit  - integers from 4368 to 69903
%   '1110' + 20 bit  - integers from 69940 to 1118479
%   '1111' + 24 bit  - integers from 1118480 to 17895695
%   not supported  - integers >= 17895696 (=2^4+2^8+2^12+2^16+2^20+2^24)
function PutVLIC(N)
global y Byte BitPos
if (N<0)
   error('Huff06-PutVLIC: Number is negative.'); 
elseif (N<16)
   PutBit(0);PutBit(0);
   for (i=4:-1:1); PutBit(bitget(N,i)); end;
elseif (N<272)
   PutBit(0);PutBit(1);
   N=N-16;
   for (i=8:-1:1); PutBit(bitget(N,i)); end;
elseif (N<4368)
   PutBit(1);PutBit(0);
   N=N-272;
   for (i=12:-1:1); PutBit(bitget(N,i)); end;
elseif (N<69940)
   PutBit(1);PutBit(1);PutBit(0);
   N=N-4368;
   for (i=16:-1:1); PutBit(bitget(N,i)); end;
elseif (N<1118480)
   PutBit(1);PutBit(1);PutBit(1);PutBit(0);
   N=N-69940;
   for (i=20:-1:1); PutBit(bitget(N,i)); end;
elseif (N<17895696)
   PutBit(1);PutBit(1);PutBit(1);PutBit(1);
   N=N-1118480;
   for (i=24:-1:1); PutBit(bitget(N,i)); end;
else
   error('Huff06-PutVLIC: Number is too large.'); 
end
return

function N=GetVLIC
global y Byte BitPos
N=0;
if GetBit
   if GetBit
      if GetBit
         if GetBit
            for (i=1:24); N=N*2+GetBit; end;
            N=N+1118480;
         else
            for (i=1:20); N=N*2+GetBit; end;
            N=N+69940;
         end
      else
         for (i=1:16); N=N*2+GetBit; end;
         N=N+4368;
      end
   else
      for (i=1:12); N=N*2+GetBit; end;
      N=N+272;
   end
else
   if GetBit
      for (i=1:8); N=N*2+GetBit; end;
      N=N+16;
   else
      for (i=1:4); N=N*2+GetBit; end;
   end
end
return

% Functions to write and read a Bit
function PutBit(Bit)
global y Byte BitPos
BitPos=BitPos-1;
if (~BitPos); Byte=Byte+1; BitPos=8; end; 
y(Byte) = bitset(y(Byte),BitPos,Bit);
return
   
function Bit=GetBit
global y Byte BitPos
BitPos=BitPos-1;
if (~BitPos); Byte=Byte+1; BitPos=8; end; 
Bit=bitget(y(Byte),BitPos);
return;
   
% this function is a variant of the standard hist function
function Hi=IntHist(W,i1,i2);
W=W(:);
%if (rem(i1,1) | rem(i2,1));   error('Non integers'); end;
L=length(W);
Hi=zeros(i2-i1+1,1);
if (i2-i1)>50
   for l=1:L
      i=W(l)-i1+1;
      Hi(i)=Hi(i)+1;
   end
else
   for i=i1:i2
      I=find(W==i);
      Hi(i-i1+1)=length(I);
   end
end
return;