|| The definitive matrices program, including as a special treat.....
||  .....transitive closure on relations! [thrilling stuff...]

matrix == [[num]]
|| list of rows, each a list of num

element :: num -> num -> matrix -> num
|| element 2 3 [[147,432,865,124],[436,987,678,459]] = 678, for example
  element 1     j (first_row : other_rows) = first_row!(j-1)
  element (i+2) j (first_row : other_rows) = element (i+1) j other_rows

no_of_rows :: matrix -> num
  no_of_rows m = #m

no_of_columns :: matrix -> num
  no_of_columns (first_row : other_rows) = #first_row



add :: (num -> num -> num) -> matrix -> matrix -> matrix
|| add-function (eg, +) -> matrix -> matrix -> sum of these two matrices
  add f [] [] = []
  add f (first_row_of_A : other_rows_of_A) (first_row_of_B : other_rows_of_B)
    = add_rows f first_row_of_A first_row_of_B
      :
      add f other_rows_of_A other_rows_of_B
      where
      add_rows f [] [] = []
      add_rows f (e1_of_A : rest_of_A) (e1_of_B : rest_of_B)
        = f e1_of_A e1_of_B   :   add_rows f rest_of_A rest_of_B

multiply :: (num -> num -> num) -> (num -> num -> num) -> matrix -> matrix -> matrix
|| multiply-function -> add-function -> matrix -> matrix -> the product
  multiply f g m1 m2 
    = error "matrices are not the right shapes!", if no_of_columns m1 ~= no_of_rows m2
    = map do_row [1..no_of_rows m1], otherwise
      where
      do_row i = map do_element [1..no_of_columns m2]
                 where
                 do_element j = sum_using_g things_to_add_up
                                where
                                things_to_add_up = map contribution_to_sum [1..no_of_columns m1]
                                contribution_to_sum k = f (element i k m1) (element k j m2)
                                sum_using_g [e1] = e1
                                sum_using_g (e1:e2:rest) = g e1 (sum_using_g (e2:rest))





num_and :: num -> num -> num
|| for multiplying matrices representing relations
  num_and 0 0 = 0
  num_and 0 1 = 0
  num_and 1 0 = 0
  num_and 1 1 = 1

num_or :: num -> num -> num
|| for adding and multiplying matrices representing relations
  num_or 0 0 = 0
  num_or 0 1 = 1
  num_or 1 0 = 1
  num_or 1 1 = 1


compose_relations :: matrix -> matrix -> matrix
|| composes two relations by multiplying using num_and and num_or
  compose_relations = multiply num_and num_or


powers :: (matrix -> matrix -> matrix) -> matrix -> [matrix]
|| multiply function -> matrix -> [matrix, matrix^2, matrix^3, .....]
  powers f m = m : map (f m) (powers f m)


transitive_closure :: matrix -> matrix
|| computes the transitive closure of a matrix representing a relation
|| This is the one all the others have been building up to!!!
  transitive_closure m
    = sum_using_add (take (no_of_rows m) (powers compose_relations m))
      where
      sum_using_add [e1] = e1
      sum_using_add (e1:e2:rest) = add num_or e1 (sum_using_add (e2:rest))
