/* IO(1) */

void FSM_inc( vector _x, int _n)
{
/* increment x, base n */

  int i, s, c = 1;	/* index, sum, carry */

  for( i = cols(_x); i > 0; --i) {

    s = c + _x(i);

    if( s > _n) {
      s = 1;
      c = 1;
    } else
      c = 0;

    _x(i) = s;
  }
}

matrix FSM( vector _s, int _L)
{
/* produce state matrix for L-state FSM with symbols s */

  matrix states;
  vector c;
  int nsyms, nstates;
  int i, j;

  nstates = nsyms = cols(_s);	/* nstates = nsyms ** L */
  for( i = 1; i < _L; ++i)
    nstates *= nsyms;
  
  states = zeros( nstates, _L);
  c = ones( 1, _L);

  for( i = 1; i <= nstates; ++i) {
    states(i) = _s(c);
    FSM_inc( c, nsyms);
  }

  return states;
}

matrix FSM_trans( matrix _states, vector _f) 
{
/* produce state transition matrix */

  matrix trans;
  int r = rows(_states), c = cols(_states);
  int i, j;

  trans = HUGE_VAL * ones( r, r);

  for( i = 1; i <= r; ++i)
    for( j = 1; j <= r; ++j)
    if( c == 1 || all( _states(i,1:c-1) == _states(j,2:c))) {
      trans(i,j) = [_states(j,:) _states(i,c)] * _f;
    }

  return trans;
}

vector FSM_out( vector _x, matrix _states, int _zstate, vector _f)
{
/* produce FSM output for input vector x */

  vector y;
  vector xx;
  int m = rows(_f);
  int n = cols(_x);
  int i;

  y = _f(1) * _x;

  for( i = 1; i < m; ++i) {
    xx = [_states(_zstate,1:i), _x(1:n-i)];
    y += _f(i+1) * xx;
  }

  return y;
}
