\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra pseudolin.spad}
\author{Bruno Zuercher}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package PSEUDLIN PseudoLinearNormalForm}
<<package PSEUDLIN PseudoLinearNormalForm>>=
)abbrev package PSEUDLIN PseudoLinearNormalForm
++ Normal forms of pseudo-linear operators
++ Author: Bruno Zuercher
++ Date Created: November 1993
++ Date Last Updated: 12 April 1994
++ Description:
++   PseudoLinearNormalForm provides a function for computing a block-companion
++   form for pseudo-linear operators.
PseudoLinearNormalForm(K:Field): Exports == Implementation where
  ER  ==> Record(C: Matrix K, g: Vector K)
  REC ==> Record(R: Matrix K, A: Matrix K, Ainv: Matrix K)

  Exports ==> with
    normalForm: (Matrix K, Automorphism K, K -> K) -> REC
      ++ normalForm(M, sig, der) returns \spad{[R, A, A^{-1}]} such that
      ++ the pseudo-linear operator whose matrix in the basis \spad{y} is
      ++ \spad{M} had matrix \spad{R} in the basis \spad{z = A y}.
      ++ \spad{der} is a \spad{sig}-derivation.
    changeBase: (Matrix K, Matrix K, Automorphism K, K -> K) -> Matrix K
      ++ changeBase(M, A, sig, der): computes the new matrix of a pseudo-linear
      ++ transform given by the matrix M under the change of base A
    companionBlocks: (Matrix K, Vector K) -> List ER
      ++ companionBlocks(m, v) returns \spad{[[C_1, g_1],...,[C_k, g_k]]}
      ++ such that each \spad{C_i} is a companion block and
      ++ \spad{m = diagonal(C_1,...,C_k)}.

  Implementation ==> add
    normalForm0: (Matrix K, Automorphism K, Automorphism K, K -> K) -> REC
    mulMatrix: (Integer, Integer, K) -> Matrix K
      -- mulMatrix(N, i, a): under a change of base with the resulting matrix of
      -- size N*N the following operations are performed:
      -- D1: column i will be multiplied by sig(a)
      -- D2: row i will be multiplied by 1/a
      -- D3: addition of der(a)/a to the element at position (i,i)
    addMatrix: (Integer, Integer, Integer, K) -> Matrix K
      -- addMatrix(N, i, k, a): under a change of base with the resulting matrix
      -- of size N*N the following operations are performed:
      -- C1: addition of column i multiplied by sig(a) to column k
      -- C2: addition of row k multiplied by -a to row i
      -- C3: addition of -a*der(a) to the element at position (i,k)
    permutationMatrix: (Integer, Integer, Integer) -> Matrix K
      -- permutationMatrix(N, i, k): under a change of base with the resulting
      -- permutation matrix of size N*N the following operations are performed:
      -- P1: columns i and k will be exchanged
      -- P2: rows i and k will be exchanged
    inv: Matrix K -> Matrix K
      -- inv(M): computes the inverse of a invertable matrix M.
      -- avoids possible type conflicts

    inv m                      == inverse(m) :: Matrix K
    changeBase(M, A, sig, der) == inv(A) * (M * map(sig #1, A) + map(der, A))
    normalForm(M, sig, der)    == normalForm0(M, sig, inv sig, der)

    companionBlocks(R, w) ==
      -- decomposes the rational matrix R into single companion blocks
      -- and the inhomogenity w as well
      i:Integer := 1
      n := nrows R
      l:List(ER) := empty()
      while i <= n repeat
        j := i
        while j+1 <= n and R(j,j+1) = 1 repeat j := j+1
        --split block now
        v:Vector K := new((j-i+1)::NonNegativeInteger, 0)
        for k in i..j repeat v(k-i+1) := w k
        l := concat([subMatrix(R,i,j,i,j), v], l)
        i := j+1
      l

    normalForm0(M, sig, siginv, der) ==
      -- the changes of base will be incremented in B and Binv,
      -- where B**(-1)=Binv; E defines an elementary matrix
      B, Binv, E    : Matrix K
      recOfMatrices : REC
      N := nrows M
      B := diagonalMatrix [1 for k in 1..N]
      Binv := copy B
      -- avoid unnecessary recursion
      if diagonal?(M) then return [M, B, Binv]
      i : Integer := 1
      while i < N repeat
        j := i + 1
        while j <= N and M(i, j) = 0 repeat  j := j + 1
        if j <= N then
          -- expand companionblock by lemma 5
          if j ~= i+1 then
            -- perform first a permutation
            E := permutationMatrix(N, i+1, j)
            M := changeBase(M, E, sig, der)
            B := B*E
            Binv := E*Binv
          -- now is M(i, i+1) ~= 0
          E := mulMatrix(N, i+1, siginv inv M(i,i+1))
          M := changeBase(M, E, sig, der)
          B := B*E
          Binv := inv(E)*Binv
          for j: local in 1..N repeat
            if j ~= i+1 then
              E := addMatrix(N, i+1, j, siginv(-M(i,j)))
              M := changeBase(M, E, sig, der)
              B := B*E
              Binv := inv(E)*Binv
          i := i + 1
        else
          -- apply lemma 6
          for j: local in i..2 by -1 repeat
            for k in (i+1)..N repeat
              E := addMatrix(N, k, j-1, M(k,j))
              M := changeBase(M, E, sig, der)
              B := B*E
              Binv := inv(E)*Binv
          j := i + 1
          while j <= N and M(j,1) = 0 repeat  j := j + 1
          if j <= N then
            -- expand companionblock by lemma 8
            E := permutationMatrix(N, 1, j)
            M := changeBase(M, E, sig, der)
            B := B*E
            Binv := E*Binv
            -- start again to establish rational form
            i := 1
          else
            -- split a direct factor
            recOfMatrices :=
              normalForm(subMatrix(M, i+1, N, i+1, N), sig, der)
            setsubMatrix!(M, i+1, i+1, recOfMatrices.R)
            E := diagonalMatrix [1 for k in 1..N]
            setsubMatrix!(E, i+1, i+1, recOfMatrices.A)
            B := B*E
            setsubMatrix!(E, i+1, i+1, recOfMatrices.Ainv)
            Binv := E*Binv
            -- M in blockdiagonalform, stop program
            i := N
      [M, B, Binv]

    mulMatrix(N, i, a) ==
      M : Matrix K := diagonalMatrix [1 for j in 1..N]
      M(i, i) := a
      M

    addMatrix(N, i, k, a) ==
      A : Matrix K := diagonalMatrix [1 for j in 1..N]
      A(i, k) := a
      A

    permutationMatrix(N, i, k) ==
      P : Matrix K := diagonalMatrix [1 for j in 1..N]
      P(i, i) := P(k, k) := 0
      P(i, k) := P(k, i) := 1
      P

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<package PSEUDLIN PseudoLinearNormalForm>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}