path: root/src/algebra/vector.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra vector.spad}
\author{The Axiom Team}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{category VECTCAT VectorCategory}
<<category VECTCAT VectorCategory>>=
)abbrev category VECTCAT VectorCategory
++ Author:
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors: DirectProductCategory, Vector
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ \spadtype{VectorCategory} represents the type of vector like objects,
++ i.e. finite sequences indexed by some finite segment of the
++ integers. The operations available on vectors depend on the structure
++ of the underlying components. Many operations from the component domain
++ are defined for vectors componentwise. It can by assumed that extraction or
++ updating components can be done in constant time.
 
VectorCategory(R:Type): Category == OneDimensionalArrayAggregate R with
    if R has AbelianSemiGroup then
      + : (%, %) -> %
        ++ x + y returns the component-wise sum of the vectors x and y.
        ++ Error: if x and y are not of the same length.
    if R has AbelianMonoid then
      zero: NonNegativeInteger -> %
        ++ zero(n) creates a zero vector of length n.
    if R has AbelianGroup then
      - : % -> %
        ++ -x negates all components of the vector x.
      - : (%, %) -> %
        ++ x - y returns the component-wise difference of the vectors x and y.
        ++ Error: if x and y are not of the same length.
      * : (Integer, %) -> %
        ++ n * y multiplies each component of the vector y by the integer n.
    if R has Monoid then
      * : (R, %) -> %
        ++ r * y multiplies the element r times each component of the vector y.
      * : (%, R) -> %
        ++ y * r multiplies each component of the vector y by the element r.
    if R has Ring then
      dot: (%, %) -> R
        ++ dot(x,y) computes the inner product of the two vectors x and y.
        ++ Error: if x and y are not of the same length.
      outerProduct: (%, %) -> Matrix R
        ++ outerProduct(u,v) constructs the matrix whose (i,j)'th element is
        ++ u(i)*v(j).
      cross: (%, %) -> %
        ++ vectorProduct(u,v) constructs the cross product of u and v.
        ++ Error: if u and v are not of length 3.
    if R has RadicalCategory and R has Ring then
      length: % -> R
        ++ length(v) computes the sqrt(dot(v,v)), i.e. the magnitude
      magnitude: % -> R
        ++ magnitude(v) computes the sqrt(dot(v,v)), i.e. the length
 add
    if R has AbelianSemiGroup then
      u + v ==
        (n := #u) ~= #v => error "Vectors must be of the same length"
        map(_+ , u, v)
 
    if R has AbelianMonoid then
      zero n == new(n, 0)
 
    if R has AbelianGroup then
      - u             == map(- #1, u)
      n:Integer * u:% == map(n * #1, u)
      u - v           == u + (-v)
 
    if R has Monoid then
      u:% * r:R       == map(#1 * r, u)
      r:R * u:%       == map(r * #1, u)
 
    if R has Ring then
      dot(u, v) ==
        #u ~= #v => error "Vectors must be of the same length"
        +/[qelt(u, i) * qelt(v, i) for i in minIndex u .. maxIndex u]
      outerProduct(u, v) ==
        matrix [[qelt(u, i) * qelt(v,j) for i in minIndex u .. maxIndex u] _
                for j in minIndex v .. maxIndex v]
      cross(u, v) ==
        #u ~= 3 or #v ~= 3 => error "Vectors must be of length 3"
        construct [qelt(u, 2)*qelt(v, 3) - qelt(u, 3)*qelt(v, 2) , _
                   qelt(u, 3)*qelt(v, 1) - qelt(u, 1)*qelt(v, 3) , _
                   qelt(u, 1)*qelt(v, 2) - qelt(u, 2)*qelt(v, 1) ]

    if R has RadicalCategory and R has Ring then
      length p ==
         sqrt(dot(p,p))
      magnitude p ==
         sqrt(dot(p,p))
 
@

\section{domain VECTOR Vector}
<<domain VECTOR Vector>>=
)abbrev domain VECTOR Vector
++ Author:
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors: DirectProduct
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This type represents vector like objects with varying lengths
++ and indexed by a finite segment of integers starting at 1.
 
Vector(R:Type): Exports == Implementation where
 Exports == VectorCategory R with 
   vector: List R -> %
     ++ vector(l) converts the list l to a vector.
 Implementation == OneDimensionalArray R add
     vector l == construct l
     -- We want maxIndex to be inlined.  Ideally, the definition should
     -- read
     --      maxIndex x == # rep x
     -- However, there is currently an infelicity in the compiler that
     -- prevents good uses of dependent domains.  So, we fall back to
     maxIndex x == sizeOfSimpleArray(x)$Foreign(Builtin)
     if R has ConvertibleTo InputForm then
       convert(x:%):InputForm ==
          convert [convert('vector)@InputForm,
                          convert(members x)@InputForm]

@

\section{package VECTOR2 VectorFunctions2}

<<package VECTOR2 VectorFunctions2>>=
)abbrev package VECTOR2 VectorFunctions2
++ Author:
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This package provides operations which all take as arguments
++ vectors of elements of some type \spad{A} and functions from \spad{A} to
++ another of type B. The operations all iterate over their vector argument
++ and either return a value of type B or a vector over B.
 
VectorFunctions2(A, B): Exports == Implementation where
  A, B: Type
 
  VA ==> Vector A
  VB ==> Vector B
  O2 ==> FiniteLinearAggregateFunctions2(A, VA, B, VB)
  UB ==> Union(B,"failed")
 
  Exports ==> with
    scan   : ((A, B) -> B, VA, B) -> VB
      ++ scan(func,vec,ident) creates a new vector whose elements are
      ++ the result of applying reduce to the binary function func,
      ++ increasing initial subsequences of the vector vec,
      ++ and the element ident.
    reduce : ((A, B) -> B, VA, B) -> B
      ++ reduce(func,vec,ident) combines the elements in vec using the
      ++ binary function func. Argument ident is returned if vec is empty.
    map    : (A -> B, VA) -> VB
      ++ map(f, v) applies the function f to every element of the vector v
      ++ producing a new vector containing the values.
    map : (A -> UB, VA) -> Union(VB,"failed")
      ++ map(f, v) applies the function f to every element of the vector v
      ++ producing a new vector containing the values or \spad{"failed"}.
 
  Implementation ==> add
    scan(f, v, b)   == scan(f, v, b)$O2
    reduce(f, v, b) == reduce(f, v, b)$O2
    map(f:(A->B), v:VA):VB == map(f, v)$O2

    map(f:(A -> UB), a:VA):Union(VB,"failed") ==
     res : List B  := []
     for u in entries(a) repeat
       r := f u
       r = "failed" => return "failed"
       res := [r::B,:res]
     vector reverse! res

@
\section{category DIRPCAT DirectProductCategory}
<<category DIRPCAT DirectProductCategory>>=
)abbrev category DIRPCAT DirectProductCategory
 
--% DirectProductCategory
 
++ Author:
++ Date Created:
++ Date Last Updated: June 17, 2010
++ Basic Functions:
++ Related Constructors: DirectProduct
++ Also See: VectorCategory
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This category represents a finite cartesian product of a given type.
++ Many categorical properties are preserved under this construction.
 
DirectProductCategory(dim:NonNegativeInteger, R:Type): Category ==
  Join(IndexedAggregate(Integer,R),FiniteAggregate R,CoercibleTo Vector R) with
         directProduct: Vector R -> %
           ++ directProduct(v) converts the vector v to become
           ++ a direct product. Error: if the length of v is
           ++ different from dim.
         if R has SetCategory then FullyRetractableTo R
         if R has Ring then
           BiModule(R, R)
           DifferentialExtension R
           FullyLinearlyExplicitRingOver R
           unitVector: PositiveInteger -> %
             ++ unitVector(n) produces a vector with 1 in position n and
             ++ zero elsewhere.
           dot: (%, %) -> R
             ++ dot(x,y) computes the inner product of the vectors x and y.
         if R has AbelianSemiGroup then AbelianSemiGroup
         if R has CancellationAbelianMonoid then CancellationAbelianMonoid
         if R has AbelianMonoid then AbelianMonoid
         if R has AbelianGroup then AbelianGroup
         if R has Monoid then LinearSet R
         if R has Finite then Finite
         if R has CommutativeRing then Module R
         if R has OrderedSet then OrderedSet
         if R has OrderedAbelianMonoidSup then OrderedAbelianMonoidSup
         if R has Field then VectorSpace R
 add
      if R has Ring then
        equation2R: Vector % -> Matrix R
 
        coerce(n:Integer):%          == n::R::%
        characteristic == characteristic$R
        differentiate(z:%, d:R -> R) == map(d, z)
 
        equation2R v ==
          ans:Matrix(R) := new(dim, #v, 0)
          for i in minRowIndex ans .. maxRowIndex ans repeat
            for j in minColIndex ans .. maxColIndex ans repeat
              qsetelt!(ans, i, j, qelt(qelt(v, j), i))
          ans
 
        reducedSystem(m:Matrix %):Matrix(R) ==
          empty? m => new(0, 0, 0)
          reduce(vertConcat, [equation2R row(m, i)
                 for i in minRowIndex m .. maxRowIndex m])$List(Matrix R)
 
        reducedSystem(m:Matrix %, v:Vector %):
          Record(mat:Matrix R, vec:Vector R) ==
            vh:Vector(R) :=
              empty? v => empty()
              rh := reducedSystem(v::Matrix %)@Matrix(R)
              column(rh, minColIndex rh)
            [reducedSystem(m)@Matrix(R), vh]
 
      if R has Finite then size == size()$R ** dim
 
      if R has Field then
        x / b       == x * inv b
        dimension() == dim::CardinalNumber
 
@
\section{domain DIRPROD DirectProduct}
<<domain DIRPROD DirectProduct>>=
)abbrev domain DIRPROD DirectProduct
++ Author:
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors: Vector
++ Also See: OrderedDirectProduct
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This type represents the finite direct or cartesian product of an
++ underlying component type. This contrasts with simple vectors in that
++ the members can be viewed as having constant length. Thus many
++ categorical properties can by lifted from the underlying component type.
++ Component extraction operations are provided but no updating operations.
++ Thus new direct product elements can either be created by converting
++ vector elements using the \spadfun{directProduct} function
++ or by taking appropriate linear combinations of basis vectors provided
++ by the \spad{unitVector} operation.
 
DirectProduct(dim:NonNegativeInteger, R:Type):
  DirectProductCategory(dim, R) == Vector R add
      coerce(z:%):Vector(R)        == copy rep z
      coerce(r:R):%                == per new(dim, r)$Vector(R)
 
      members x == VEC2LIST(x)$Lisp
 
      directProduct z ==
        #z = dim => per copy z
        error "Not of the correct length"
 
 
      if R has SetCategory then
        same?: % -> Boolean
        same? z == every?(#1 = z(minIndex z), z)
 
        x = y == and/[qelt(rep x,i) = qelt(rep y,i) for i in 1..dim]
 
        retract(z:%):R ==
          same? z => z(minIndex z)
          error "Not retractable"
 
        retractIfCan(z:%):Union(R, "failed") ==
          same? z => z(minIndex z)
          "failed"
 
 
      if R has AbelianSemiGroup then
        u:% + v:% == per map(_+ , rep u, rep v)
 
      if R has AbelianMonoid then
        0 == per zero(dim)$Vector(R)
 
      if R has Monoid then
        1 == per new(dim, 1)$Vector(R)
        u:% * r:R       == per map(#1 * r, rep u)
        r:R * u:%       == per map(r * #1, rep u)
        x:% * y:%       == per [x.i * y.i for i in 1..dim]$Vector(R)
 
 
      if R has CancellationAbelianMonoid then
        subtractIfCan(u:%, v:%) ==
          w := new(dim,0)$Vector(R)
          for i in 1..dim repeat
            (c := subtractIfCan(qelt(rep u, i), qelt(rep v,i))) case nothing =>
                    return nothing
            qsetelt!(w, i, c::R)
          just per w
 
      if R has Ring then
 
        u:% * v:%                    == per map(_* ,rep u,rep v)
 
        recip z ==
          w := new(dim,0)$Vector(R)
          for i in minIndex w .. maxIndex w repeat
            (u := recip qelt(z, i)) case "failed" => return "failed"
            qsetelt!(w, i, u::R)
          per w
 
        unitVector i ==
          v:= new(dim,0)$Vector(R)
          v.i := 1
          per v
 
      if R has OrderedSet then
        x < y ==
          for i in 1..dim repeat
             qelt(x,i) < qelt(y,i) => return true
             qelt(x,i) > qelt(y,i) => return false
          false

      if R has OrderedAbelianMonoidSup then
        sup(x,y) == per map(sup, rep x, rep y)
 
@
\section{package DIRPROD2 DirectProductFunctions2}
<<package DIRPROD2 DirectProductFunctions2>>=
)abbrev package DIRPROD2 DirectProductFunctions2
++ Author:
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This package provides operations which all take as arguments
++ direct products of elements of some type \spad{A} and functions from \spad{A} to another
++ type B. The operations all iterate over their vector argument
++ and either return a value of type B or a direct product over B.
 
DirectProductFunctions2(dim, A, B): Exports == Implementation where
  dim : NonNegativeInteger
  A, B: Type
 
  DA ==> DirectProduct(dim, A)
  DB ==> DirectProduct(dim, B)
  VA ==> Vector A
  VB ==> Vector B
  O2 ==> FiniteLinearAggregateFunctions2(A, VA, B, VB)
 
  Exports ==> with
    scan   : ((A, B) -> B, DA, B) -> DB
      ++ scan(func,vec,ident) creates a new vector whose elements are
      ++ the result of applying reduce to the binary function func,
      ++ increasing initial subsequences of the vector vec,
      ++ and the element ident.
    reduce : ((A, B) -> B, DA, B) -> B
      ++ reduce(func,vec,ident) combines the elements in vec using the
      ++ binary function func. Argument ident is returned if the vector is empty.
    map    : (A -> B, DA) -> DB
      ++ map(f, v) applies the function f to every element of the vector v
      ++ producing a new vector containing the values.
 
  Implementation ==> add
    import FiniteLinearAggregateFunctions2(A, VA, B, VB)
 
    map(f, v)       == directProduct map(f, v::VA)
    scan(f, v, b)   == directProduct scan(f, v::VA, b)
    reduce(f, v, b) == reduce(f, v::VA, b)

@

\section{Vector space of finite dimension given by a basis}


<<domain LINBASIS LinearBasis>>=
++ Author: Gabriel Dos Reis
++ Date Created: July 2, 2010
++ Date Last Modified: July 2, 2010
++ Descrption:
++   Representation of a vector space basis, given by symbols.
)abbrev domain LINBASIS LinearBasis
LinearBasis(vars: List Symbol): Public == Private where
  Public == Join(OrderedFinite,CoercibleFrom OrderedVariableList vars) with
    dual: DualBasis vars -> %
      ++ \spad{dual f} constructs the dual vector of a linear form
      ++ which is part of a basis.
  Private == OrderedVariableList vars add
    coerce(s: Rep): % == per s
    dual f == index(lookup f)@%

@

<<domain DBASIS DualBasis>>=
++ Author: Gabriel Dos Reis
++ Date Created: July 2, 2010
++ Date Last Modified: July 2, 2010
++ Descrption:
++   Representation of a dual vector space basis, given by symbols.
)abbrev domain DBASIS DualBasis
DualBasis(vars: List Symbol): Public == Private where
  Public == Join(OrderedFinite) with
    dual: LinearBasis vars -> %
      ++ \spad{dual x} constructs the dual vector of a linear element
      ++ which is part of a basis.
  Private == OrderedVariableList vars add
    dual v == index(lookup v)@%
    coerce(x: %): OutputForm ==
      superscript(convert(rep x)@Symbol,['_*::OutputForm])::OutputForm

@

<<domain LINELT LinearElement>>=
)abbrev domain LINELT LinearElement
++ Author: Gabriel Dos Reis
++ Date Created: June 30, 2010
++ Date Last Modified: July 2, 2010
++ Description:
++   A simple data structure for elements that form a
++   vector space of finite dimension over a given field,
++   with a given symbolic basis.
LinearElement(K,B): Public == Private where
  K: Field
  B: List Symbol
  macro Basis == LinearBasis B
  Public == Join(VectorSpace K,CoercibleFrom Basis,_
                 IndexedDirectProductCategory(K,Basis)) with
    linearElement: List K -> %
      ++ \spad{linearElement [x1,..,xn]} returns a linear element
      ++  with coordinates \spad{[x1,..,xn]} with respect to
      ++ the basis elements \spad{B}.
    coordinates: % -> Vector K
      ++ \spad{coordinates x} returns the coordinates of the linear
      ++ element with respect to the basis \spad{B}.
  Private == FreeModule(K,Basis) add
    coerce(b: Basis): % ==
      per monomial(1$K,b)

    dimension() ==
      size()$Basis::CardinalNumber

    linearElement cs ==
      reduce(_+@(%,%)->%,
         [per monomial(c,index(i)$Basis) for c in cs for i in 1..],0)

    coordinates x ==
      n := #B
      v: Vector K := new(n,0$K)
      ts := terms rep x
      for i in 1..n repeat
        t := first ts
        lookup index t = i => 
          v.i := coefficient t
          ts := rest ts
      v
@


<<domain LINFORM LinearForm>>=
)abbrev domain LINFORM LinearForm
++ Author: Gabriel Dos Reis
++ Date Created: July 2, 2010
++ Date Last Modified: July 2, 2010
++ Description:
++   A simple data structure for linear forms on a vector space of
++   finite dimension over a given field, with a given symbolic basis.
LinearForm(K,B): Public == Private where
  K: Field
  B: List Symbol
  macro Basis == DualBasis B
  Public == Join(VectorSpace K,Eltable(LinearElement(K,B),K)) with
    linearForm: List K -> %
      ++ \spad{linearForm [x1,..,xn]} constructs
      ++ a linear form with coordinates \spad{[x1,..,xn]} with
      ++ respect to the basis elements \spad{DualBasis B}.
    coordinates: % -> Vector K
      ++ \spad{coordinates x} returns the coordinates of the linear
      ++ form with respect to the basis \spad{DualBasis B}.
  Private == FreeModule(K,Basis) add
    dimension() ==
      size()$Basis::CardinalNumber

    linearForm cs ==
      reduce(_+@(%,%)->%,
         [per monomial(c,index(i)$Basis) for c in cs for i in 1..],0)

    coordinates f ==
      n := #B
      v: Vector K := new(n,0$K)
      ts := terms rep f
      for i in 1..n repeat
        t := first ts
        lookup index t = i => 
          v.i := coefficient t
          ts := rest ts
      v

    elt(f,x) ==
      dot(coordinates f, coordinates x)$Vector(K)

@


\section{License}

<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--Copyright (C) 2007-2013, Gabriel Dos Reis.
--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>>
 
<<category VECTCAT VectorCategory>>
<<domain VECTOR Vector>>
<<package VECTOR2 VectorFunctions2>>
<<category DIRPCAT DirectProductCategory>>
<<domain DIRPROD DirectProduct>>
<<package DIRPROD2 DirectProductFunctions2>>

<<domain LINBASIS LinearBasis>>
<<domain DBASIS DualBasis>>

<<domain LINELT LinearElement>>
<<domain LINFORM LinearForm>>
 
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}