Initial population.

1 files changed, 208 insertions, 0 deletions

diff --git a/src/algebra/pseudolin.spad.pamphlet b/src/algebra/pseudolin.spad.pamphlet
new file mode 100644
index 00000000..3eb384b4
--- /dev/null
+++ b/src/algebra/pseudolin.spad.pamphlet
@@ -0,0 +1,208 @@
+\documentclass{article}
+\usepackage{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 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 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}