function y=xor(a,b,c,d)
  [lhs,rhs] = argn(0);
  if rhs<=3 then d=0; end
  if rhs<=2 then c=0; end
  if rhs<=1 then b=0; end
    
  y=modulo(a+b+c+d,2)
endfunction

function [bitseries,nextstnum] = mlbs(log2N,bitno,startnum)
  
  // get default arguments
  [lhs,rhs] = argn(0);
  if rhs<=2 then startnum=2^log2N-1, warning('No start number specified!'); end
  if rhs<=1 then bitno=log2N, warning('No series length specified!'); end
  // get taps for Maximum Length LFSR
  select log2N
    case 1 then 
      error('log2N=1 is not supported yet.');
    case 2 then 
      tapa=2,  tapb=1,  tapc=log2N+1, tapd=log2N+1;
    case 3 then 
      tapa=3,  tapb=2,  tapc=log2N+1, tapd=log2N+1;
    case 4 then 
      tapa=4,  tapb=3,  tapc=log2N+1, tapd=log2N+1;
    case 5 then 
      tapa=5,  tapb=3,  tapc=log2N+1, tapd=log2N+1;
    case 6 then 
      tapa=6,  tapb=5,  tapc=log2N+1, tapd=log2N+1;
    case 7 then 
      tapa=7,  tapb=6,  tapc=log2N+1, tapd=log2N+1;
    case 8 then 
      tapa=8,  tapb=6,  tapc=5, tapd=4;
    case 9 then 
      tapa=9,  tapb=5,  tapc=log2N+1, tapd=log2N+1;
    case 10 then 
      tapa=10, tapb=7,  tapc=log2N+1, tapd=log2N+1;
    case 11 then 
      tapa=11, tapb=9,  tapc=log2N+1, tapd=log2N+1;
    case 12 then 
      tapa=12, tapb=6,  tapc=4, tapd=1;
    case 13 then 
      tapa=13, tapb=4,  tapc=3, tapd=1;
    case 14 then 
      tapa=14, tapb=5,  tapc=3, tapd=1;
    case 15 then 
      tapa=15, tapb=14, tapc=log2N+1, tapd=log2N+1;
    case 16 then 
      tapa=16, tapb=15, tapc=13, tapd=4;
    case 17 then 
      tapa=17, tapb=14, tapc=log2N+1, tapd=log2N+1;
    case 18 then 
      tapa=18, tapb=11, tapc=log2N+1, tapd=log2N+1;
    case 19 then 
      tapa=19, tapb=6,  tapc=2,  tapd=1;
    case 20 then 
      tapa=20, tapb=17, tapc=log2N+1, tapd=log2N+1;
    case 21 then 
      tapa=21, tapb=19, tapc=log2N+1, tapd=log2N+1;
    case 22 then 
      tapa=22, tapb=21, tapc=log2N+1, tapd=log2N+1;
    case 23 then 
      tapa=23, tapb=18, tapc=log2N+1, tapd=log2N+1;
    case 24 then 
      tapa=24, tapb=23, tapc=22, tapd=17;
    case 25 then 
      tapa=25, tapb=22, tapc=log2N+1, tapd=log2N+1;
    case 26 then 
      tapa=26, tapb=6,  tapc=2, tapd=1;
    case 27 then 
      tapa=27, tapb=5,  tapc=2, tapd=1;
    case 28 then 
      tapa=28, tapb=25, tapc=log2N+1, tapd=log2N+1;
    case 29 then 
      tapa=29, tapb=27, tapc=log2N+1, tapd=log2N+1;
    case 30 then 
      tapa=30, tapb=6,  tapc=4, tapd=1;
    case 31 then 
      tapa=31, tapb=28, tapc=log2N+1, tapd=log2N+1;
    case 32 then 
      tapa=32, tapb=22, tapc=2, tapd=1;
    else
      error('log2N>32 is not supported yet.');
  end
  // parse start number
  if startnum<0 then 
    error('startnum out of range.'); 
  elseif startnum==0 then 
    warning('startnum is 0.');
  elseif startnum>=2^log2N then 
    error('startnum out of range.'); 
  else
    startnum_col=zeros(log2N,1);
    startnum_str=dec2bin(startnum,log2N);
    for i=1:log2N
      //startnum_col(i,1)=sscanf(part(startnum_str,i),'%d')
      startnum_col(i,1)=sscanf(part(startnum_str,(log2N+1-i)),'%d')
    end
  end
  // initialize data structure   
  newbits    = zeros(log2N,1);
  oldbits    = zeros(log2N+1,1);
  mlbseries  = zeros(ceil(bitno/log2N)*log2N,1);
  nextstbits = zeros(log2N,1);
  nextstnum  = 0;
  // commence generation
  newbits = startnum_col;
  for i=1:ceil(bitno/log2N)
    for j=1:log2N,
      oldbits(1:log2N,1) = newbits;
      newbit  = xor(oldbits(tapa,1),oldbits(tapb,1),oldbits(tapc,1),oldbits(tapd,1));
      newbits(1,1)= newbit;
      for k=1:log2N-1,
        newbits(k+1,1) = oldbits(k,1);
      end
      mlbseries((i-1)*log2N+j,1)= newbit;
    end
  end
  // format results
  bitseries  = mlbseries(1:bitno,1);
  if bitno<=log2N then
    nextstbits = mlbseries(1:log2N,1);
  else
    nextstbits = mlbseries(bitno-log2N+1:bitno,1);
  end
  for i=1:log2N
    if nextstbits(i,1)==1 then
      nextstnum=nextstnum+2^(i-1);
    else
      nextstnum=nextstnum;
    end 
  end
  
endfunction

//// mlbs cycle test procedure, 2 to 17 verified.
//m=2;
//n=17;
//c=10;
//testCycle = zeros(n-m+1,1);
//for i=m:n
//  i
//  log2N=i;
//  bitno=c*(2^log2N-1);
//  [bitseries,nextstnum]=mlbs(log2N,bitno,2^log2N-1);
//  for j=1:bitno-log2N+1
//    if bitseries(j:j+log2N-1,1)==bitseries(1:log2N,1) then
//      if modulo(j-1,2^log2N-1)==0 then
//        testCycle(i-m+1,1)=testCycle(i-m+1,1)+1;
//      else
//        j
//      end
//    end
//  end
//end
//
//testCycle

    

