# Compute linear recursive sequences using galois stepping (I think --
# or maybe I have it backwards)
# FIXME: Check names for this stuff!

# Galois Matrix
# Given a linear combining rule (a_1*x_+...+a_n*x_n=x_(n+1)),
# gives the galois stepping matrix
SetHelp ("GaloisMatrix", "combinatorics", "Galois matrix given a linear combining rule (a_1*x_+...+a_n*x_n=x_(n+1))")
function GaloisMatrix(combining_rule) = (
 [[0;I(columns(combining_rule)-1)],combining_rule.']
)
protect ("GaloisMatrix")

# Linear recursive sequence
SetHelp ("LinearRecursiveSequence", "combinatorics", "Compute linear recursive sequence using galois stepping")
function LinearRecursiveSequence(seed_values,combining_rule,n) =
(
 k=columns(seed_values);
 if (k > n >= 0) then                      # If asks for one of the seed values, return it
   seed_values@(n+1)
 else                                      # otherwise...
 (
  G=galois_matrix(combining_rule);         # form the galois matrix
  if (n >= k) then
   (seed_values*G^(n-k+1))@(k)       # ...and step it enough times
  else if (IsInvertible(G)) then
   (seed_values*G^n)@(1)             # ...or step it backwards
  else null                         # (if sequence is reversible)
 )
)
protect ("LinearRecursiveSequence")