open-axiom-source 1.4.1+svn~2626-2ubuntu2 / usr / lib / open-axiom / src / algebra / naalgc.spad

--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.

)abbrev category MONAD Monad
++ Authors: J. Grabmeier, R. Wisbauer
++ Date Created: 01 March 1991
++ Date Last Updated: 11 June 1991
++ Basic Operations: *, **
++ Related Constructors: SemiGroup, Monoid, MonadWithUnit
++ Also See:
++ AMS Classifications:
++ Keywords: Monad,  binary operation
++ Reference:
++  N. Jacobson: Structure and Representations of Jordan Algebras
++  AMS, Providence, 1968
++ Description:
++  Monad is the class of all multiplicative monads, i.e. sets
++  with a binary operation.
Monad(): Category == SetCategory with
    --operations
      *: (%,%) -> %
        ++ a*b is the product of \spad{a} and b in a set with
        ++ a binary operation.
      rightPower: (%,PositiveInteger) -> %
        ++ rightPower(a,n) returns the \spad{n}-th right power of \spad{a},
        ++ i.e. \spad{rightPower(a,n) := rightPower(a,n-1) * a} and
        ++ \spad{rightPower(a,1) := a}.
      leftPower: (%,PositiveInteger) -> %
        ++ leftPower(a,n) returns the \spad{n}-th left power of \spad{a},
        ++ i.e. \spad{leftPower(a,n) := a * leftPower(a,n-1)} and
        ++ \spad{leftPower(a,1) := a}.
      **: (%,PositiveInteger) -> %
        ++ a**n returns the \spad{n}-th power of \spad{a},
        ++ defined by repeated squaring.
    add
      import RepeatedSquaring(%)
      x:% ** n:PositiveInteger == expt(x,n)
      rightPower(a,n) ==
        one? n => a
        res := a
        for i in 1..(n-1) repeat res := res * a
        res
      leftPower(a,n) ==
        one? n => a
        res := a
        for i in 1..(n-1) repeat res := a * res
        res

)abbrev category MONADWU MonadWithUnit
++ Authors: J. Grabmeier, R. Wisbauer
++ Date Created: 01 March 1991
++ Date Last Updated: 11 June 1991
++ Basic Operations: *, **, 1
++ Related Constructors: SemiGroup, Monoid, Monad
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Keywords: Monad with unit, binary operation
++ Reference:
++  N. Jacobson: Structure and Representations of Jordan Algebras
++  AMS, Providence, 1968
++ Description:
++  MonadWithUnit is the class of multiplicative monads with unit,
++  i.e. sets with a binary operation and a unit element.
++ Axioms
++    leftIdentity("*":(%,%)->%,1)   \tab{30} 1*x=x
++    rightIdentity("*":(%,%)->%,1)  \tab{30} x*1=x
++ Common Additional Axioms
++    unitsKnown---if "recip" says "failed", that PROVES input wasn't a unit
MonadWithUnit(): Category == Monad with
    --constants
      1: constant ->  %
        ++ 1 returns the unit element, denoted by 1.
    --operations
      one?: % -> Boolean
        ++ one?(a) tests whether \spad{a} is the unit 1.
      rightPower: (%,NonNegativeInteger) -> %
        ++ rightPower(a,n) returns the \spad{n}-th right power of \spad{a},
        ++ i.e. \spad{rightPower(a,n) := rightPower(a,n-1) * a} and
        ++ \spad{rightPower(a,0) := 1}.
      leftPower: (%,NonNegativeInteger) -> %
        ++ leftPower(a,n) returns the \spad{n}-th left power of \spad{a},
        ++ i.e. \spad{leftPower(a,n) := a * leftPower(a,n-1)} and
        ++ \spad{leftPower(a,0) := 1}.
      "**": (%,NonNegativeInteger) -> %
        ++ \spad{a**n} returns the \spad{n}-th power of \spad{a},
        ++ defined by repeated squaring.
      recip: % -> Union(%,"failed")
        ++ recip(a) returns an element, which is both a left and a right
        ++ inverse of \spad{a},
        ++ or \spad{"failed"} if such an element doesn't exist or cannot
        ++ be determined (see unitsKnown).
      leftRecip: % -> Union(%,"failed")
        ++ leftRecip(a) returns an element, which is a left inverse of \spad{a},
        ++ or \spad{"failed"} if such an element doesn't exist or cannot
        ++ be determined (see unitsKnown).
      rightRecip: % -> Union(%,"failed")
        ++ rightRecip(a) returns an element, which is a right inverse of
        ++ \spad{a}, or \spad{"failed"} if such an element doesn't exist
        ++ or cannot be determined (see unitsKnown).
    add
      import RepeatedSquaring(%)
      one? x == x = 1
      x:% ** n:NonNegativeInteger ==
         zero? n => 1
         expt(x,n pretend PositiveInteger)
      rightPower(a: %,n: NonNegativeInteger) ==
        zero? n => 1
        res: % := 1
        for i in 1..n repeat res := res * a
        res
      leftPower(a: %,n: NonNegativeInteger) ==
        zero? n => 1
        res: % := 1
        for i in 1..n repeat res := a * res
        res

)abbrev category NARNG NonAssociativeRng
++ Author: J. Grabmeier, R. Wisbauer
++ Date Created: 01 March 1991
++ Date Last Updated: 03 July 1991
++ Basic Operations: +, *, -, **
++ Related Constructors: Rng, Ring, NonAssociativeRing
++ Also See:
++ AMS Classifications:
++ Keywords: not associative ring
++ Reference:
++  R.D. Schafer: An Introduction to Nonassociative Algebras
++  Academic Press, New York, 1966
++ Description:
++  NonAssociativeRng is a basic ring-type structure, not necessarily
++  commutative or associative, and not necessarily with unit.
++  Axioms
++    x*(y+z) = x*y + x*z
++    (x+y)*z = x*z + y*z
++  Common Additional Axioms
++    noZeroDivisors  ab = 0 => a=0 or b=0
NonAssociativeRng(): Category == Join(AbelianGroup,Monad)  with
    associator: (%,%,%) -> %
      ++ associator(a,b,c) returns \spad{(a*b)*c-a*(b*c)}.
    commutator: (%,%) -> %
      ++ commutator(a,b) returns \spad{a*b-b*a}.
    antiCommutator: (%,%) -> %
      ++ antiCommutator(a,b) returns \spad{a*b+b*a}.
  add
    associator(x,y,z) == (x*y)*z - x*(y*z)
    commutator(x,y) == x*y - y*x
    antiCommutator(x,y) == x*y + y*x

)abbrev category NASRING NonAssociativeRing
++ Author: J. Grabmeier, R. Wisbauer
++ Date Created: 01 March 1991
++ Date Last Updated: 11 June 1991
++ Basic Operations: +, *, -, **
++ Related Constructors: NonAssociativeRng, Rng, Ring
++ Also See:
++ AMS Classifications:
++ Keywords: non-associative ring with unit
++ Reference:
++  R.D. Schafer: An Introduction to Nonassociative Algebras
++  Academic Press, New York, 1966
++ Description:
++  A NonAssociativeRing is a non associative rng which has a unit,
++  the multiplication is not necessarily commutative or associative.
NonAssociativeRing(): Category == Join(NonAssociativeRng,MonadWithUnit) with
    --operations
      characteristic: NonNegativeInteger
        ++ characteristic() returns the characteristic of the ring.
        --we can not make this a constant, since some domains are mutable
      coerce: Integer -> %
        ++ coerce(n) coerces the integer n to an element of the ring.
   add
      n:Integer
      coerce(n) == n * 1$%

)abbrev category NAALG NonAssociativeAlgebra
++ Author: J. Grabmeier, R. Wisbauer
++ Date Created: 01 March 1991
++ Date Last Updated: 11 June 1991
++ Basic Operations: +, -, *, **
++ Related Constructors: Algebra
++ Also See:
++ AMS Classifications:
++ Keywords: nonassociative algebra
++ Reference:
++  R.D. Schafer: An Introduction to Nonassociative Algebras
++  Academic Press, New York, 1966
++ Description:
++   NonAssociativeAlgebra is the category of non associative algebras
++   (modules which are themselves non associative rngs).
++   Axioms
++      r*(a*b) = (r*a)*b = a*(r*b)
NonAssociativeAlgebra(R:CommutativeRing): Category == _
  Join(NonAssociativeRng, Module R) with
    --operations
    plenaryPower : (%,PositiveInteger) -> %
      ++ plenaryPower(a,n) is recursively defined to be
      ++ \spad{plenaryPower(a,n-1)*plenaryPower(a,n-1)} for \spad{n>1}
      ++ and \spad{a} for \spad{n=1}.
  add
    plenaryPower(a,n) ==
      one? n => a
      n1 : PositiveInteger := (n-1)::NonNegativeInteger::PositiveInteger
      plenaryPower(a,n1) * plenaryPower(a,n1)

)abbrev category FINAALG FiniteRankNonAssociativeAlgebra
++ Author: J. Grabmeier, R. Wisbauer
++ Date Created: 01 March 1991
++ Date Last Updated: 12 June 1991
++ Basic Operations: +,-,*,**, someBasis
++ Related Constructors: FramedNonAssociativeAlgebra, FramedAlgebra,
++   FiniteRankAssociativeAlgebra
++ Also See:
++ AMS Classifications:
++ Keywords: nonassociative algebra, basis
++ References:
++   R.D. Schafer: An Introduction to Nonassociative Algebras
++   Academic Press, New York, 1966
++
++   R. Wisbauer: Bimodule Structure of Algebra
++   Lecture Notes Univ. Duesseldorf 1991
++ Description:
++   A FiniteRankNonAssociativeAlgebra is a non associative algebra over
++   a commutative ring R which is a free \spad{R}-module of finite rank.
FiniteRankNonAssociativeAlgebra(R:CommutativeRing):
 Category == NonAssociativeAlgebra R with
    someBasis : () -> Vector %
      ++ someBasis() returns some \spad{R}-module basis.
    rank : () -> PositiveInteger
      ++ rank() returns the rank of the algebra as \spad{R}-module.
    conditionsForIdempotents: Vector % -> List Polynomial R
      ++ conditionsForIdempotents([v1,...,vn]) determines a complete list
      ++ of polynomial equations for the coefficients of idempotents
      ++ with respect to the \spad{R}-module basis \spad{v1},...,\spad{vn}.
    structuralConstants: Vector % -> Vector Matrix R
      ++ structuralConstants([v1,v2,...,vm]) calculates the structural
      ++ constants \spad{[(gammaijk) for k in 1..m]} defined by
      ++ \spad{vi * vj = gammaij1 * v1 + ... + gammaijm * vm},
      ++ where \spad{[v1,...,vm]} is an \spad{R}-module basis
      ++ of a subalgebra.
    leftRegularRepresentation: (% , Vector %) -> Matrix R
      ++ leftRegularRepresentation(a,[v1,...,vn]) returns the matrix of
      ++ the linear map defined by left multiplication by \spad{a}
      ++ with respect to the \spad{R}-module basis \spad{[v1,...,vn]}.
    rightRegularRepresentation: (% , Vector %) -> Matrix R
      ++ rightRegularRepresentation(a,[v1,...,vn]) returns the matrix of
      ++ the linear map defined by right multiplication by \spad{a}
      ++ with respect to the \spad{R}-module basis \spad{[v1,...,vn]}.
    leftTrace: %  -> R
      ++ leftTrace(a) returns the trace of the left regular representation
      ++ of \spad{a}.
    rightTrace: %  -> R
      ++ rightTrace(a) returns the trace of the right regular representation
      ++ of \spad{a}.
    leftNorm: %  -> R
      ++ leftNorm(a) returns the determinant of the left regular representation
      ++ of \spad{a}.
    rightNorm: %  -> R
      ++ rightNorm(a) returns the determinant of the right regular
      ++ representation of \spad{a}.
    coordinates: (%, Vector %) -> Vector R
      ++ coordinates(a,[v1,...,vn]) returns the coordinates of \spad{a}
      ++ with respect to the \spad{R}-module basis \spad{v1},...,\spad{vn}.
    coordinates: (Vector %, Vector %) -> Matrix R
      ++ coordinates([a1,...,am],[v1,...,vn]) returns a matrix whose
      ++ i-th row is formed by the coordinates of \spad{ai}
      ++ with respect to the \spad{R}-module basis \spad{v1},...,\spad{vn}.
    represents: (Vector R, Vector %) -> %
      ++ represents([a1,...,am],[v1,...,vm]) returns the linear
      ++ combination \spad{a1*vm + ... + an*vm}.
    leftDiscriminant: Vector % -> R
      ++ leftDiscriminant([v1,...,vn]) returns  the determinant of the
      ++ \spad{n}-by-\spad{n} matrix whose element at the \spad{i}-th row
      ++ and \spad{j}-th column is given by the left trace of the product
      ++ \spad{vi*vj}.
      ++ Note: the same as \spad{determinant(leftTraceMatrix([v1,...,vn]))}.
    rightDiscriminant: Vector % -> R
      ++ rightDiscriminant([v1,...,vn]) returns  the determinant of the
      ++ \spad{n}-by-\spad{n} matrix whose element at the \spad{i}-th row
      ++ and \spad{j}-th column is given by the right trace of the product
      ++ \spad{vi*vj}.
      ++ Note: the same as \spad{determinant(rightTraceMatrix([v1,...,vn]))}.
    leftTraceMatrix: Vector % -> Matrix R
      ++ leftTraceMatrix([v1,...,vn]) is the \spad{n}-by-\spad{n} matrix
      ++ whose element at the \spad{i}-th row and \spad{j}-th column is given
      ++ by the left trace of the product \spad{vi*vj}.
    rightTraceMatrix: Vector % -> Matrix R
      ++ rightTraceMatrix([v1,...,vn]) is the \spad{n}-by-\spad{n} matrix
      ++ whose element at the \spad{i}-th row and \spad{j}-th column is given
      ++ by the right trace of the product \spad{vi*vj}.
    leftCharacteristicPolynomial: % -> SparseUnivariatePolynomial R
      ++ leftCharacteristicPolynomial(a) returns the characteristic
      ++ polynomial of the left regular representation of \spad{a}
      ++ with respect to any basis.
    rightCharacteristicPolynomial: % -> SparseUnivariatePolynomial R
      ++ rightCharacteristicPolynomial(a) returns the characteristic
      ++ polynomial of the right regular representation of \spad{a}
      ++ with respect to any basis.

    --we not necessarily have a unit
    --if R has CharacteristicZero then CharacteristicZero
    --if R has CharacteristicNonZero then CharacteristicNonZero

    commutative?:()-> Boolean
      ++ commutative?() tests if multiplication in the algebra
      ++ is commutative.
    antiCommutative?:()-> Boolean
      ++ antiCommutative?() tests if \spad{a*a = 0}
      ++ for all \spad{a} in the algebra.
      ++ Note: this implies \spad{a*b + b*a = 0} for all \spad{a} and \spad{b}.
    associative?:()-> Boolean
      ++ associative?() tests if multiplication in algebra
      ++ is associative.
    antiAssociative?:()-> Boolean
      ++ antiAssociative?() tests if multiplication in algebra
      ++ is anti-associative, i.e. \spad{(a*b)*c + a*(b*c) = 0}
      ++ for all \spad{a},b,c in the algebra.
    leftAlternative?: ()-> Boolean
      ++ leftAlternative?() tests if \spad{2*associator(a,a,b) = 0}
      ++ for all \spad{a}, b in the algebra.
      ++ Note: we only can test this; in general we don't know
      ++ whether \spad{2*a=0} implies \spad{a=0}.
    rightAlternative?: ()-> Boolean
      ++ rightAlternative?() tests if \spad{2*associator(a,b,b) = 0}
      ++ for all \spad{a}, b in the algebra.
      ++ Note: we only can test this; in general we don't know
      ++ whether \spad{2*a=0} implies \spad{a=0}.
    flexible?: ()->  Boolean
      ++ flexible?() tests if \spad{2*associator(a,b,a) = 0}
      ++ for all \spad{a}, b in the algebra.
      ++ Note: we only can test this; in general we don't know
      ++ whether \spad{2*a=0} implies \spad{a=0}.
    alternative?: ()-> Boolean
      ++ alternative?() tests if
      ++ \spad{2*associator(a,a,b) = 0 = 2*associator(a,b,b)}
      ++ for all \spad{a}, b in the algebra.
      ++ Note: we only can test this; in general we don't know
      ++ whether \spad{2*a=0} implies \spad{a=0}.
    powerAssociative?:()-> Boolean
      ++ powerAssociative?() tests if all subalgebras
      ++ generated by a single element are associative.
    jacobiIdentity?:() -> Boolean
      ++ jacobiIdentity?() tests if \spad{(a*b)*c + (b*c)*a + (c*a)*b = 0}
      ++ for all \spad{a},b,c in the algebra. For example, this holds
      ++ for crossed products of 3-dimensional vectors.
    lieAdmissible?: () -> Boolean
      ++ lieAdmissible?() tests if the algebra defined by the commutators
      ++ is a Lie algebra, i.e. satisfies the Jacobi identity.
      ++ The property of anticommutativity follows from definition.
    jordanAdmissible?: () -> Boolean
      ++ jordanAdmissible?() tests if 2 is invertible in the
      ++ coefficient domain and the multiplication defined by
      ++ \spad{(1/2)(a*b+b*a)} determines a
      ++ Jordan algebra, i.e. satisfies the Jordan identity.
      ++ The property of \spadatt{commutative("*")}
      ++ follows from by definition.
    noncommutativeJordanAlgebra?: () -> Boolean
      ++ noncommutativeJordanAlgebra?() tests if the algebra
      ++ is flexible and Jordan admissible.
    jordanAlgebra?:() -> Boolean
      ++ jordanAlgebra?() tests if the algebra is commutative,
      ++ characteristic is not 2, and \spad{(a*b)*a**2 - a*(b*a**2) = 0}
      ++ for all \spad{a},b,c in the algebra (Jordan identity).
      ++ Example:
      ++ for every associative algebra \spad{(A,+,@)} we can construct a
      ++ Jordan algebra \spad{(A,+,*)}, where \spad{a*b := (a@b+b@a)/2}.
    lieAlgebra?:() -> Boolean
      ++ lieAlgebra?() tests if the algebra is anticommutative
      ++ and \spad{(a*b)*c + (b*c)*a + (c*a)*b = 0}
      ++ for all \spad{a},b,c in the algebra (Jacobi identity).
      ++ Example:
      ++ for every associative algebra \spad{(A,+,@)} we can construct a
      ++ Lie algebra \spad{(A,+,*)}, where \spad{a*b := a@b-b@a}.

    if R has IntegralDomain then
      -- we not neccessarily have a unit, hence we don't inherit
      -- the next 3 functions anc hence copy them from MonadWithUnit:
      recip: % -> Union(%,"failed")
        ++ recip(a) returns an element, which is both a left and a right
        ++ inverse of \spad{a},
        ++ or \spad{"failed"} if there is no unit element, if such an
        ++ element doesn't exist or cannot be determined (see unitsKnown).
      leftRecip: % -> Union(%,"failed")
        ++ leftRecip(a) returns an element, which is a left inverse of \spad{a},
        ++ or \spad{"failed"} if there is no unit element, if such an
        ++ element doesn't exist or cannot be determined (see unitsKnown).
      rightRecip: % -> Union(%,"failed")
        ++ rightRecip(a) returns an element, which is a right inverse of
        ++ \spad{a},
        ++ or \spad{"failed"} if there is no unit element, if such an
        ++ element doesn't exist or cannot be determined (see unitsKnown).
      associatorDependence:() -> List Vector R
        ++ associatorDependence() looks for the associator identities, i.e.
        ++ finds a basis of the solutions of the linear combinations of the
        ++ six permutations of \spad{associator(a,b,c)} which yield 0,
        ++ for all \spad{a},b,c in the algebra.
        ++ The order of the permutations is \spad{123 231 312 132 321 213}.
      leftMinimalPolynomial : % -> SparseUnivariatePolynomial R
        ++ leftMinimalPolynomial(a) returns the polynomial determined by the
        ++ smallest non-trivial linear combination of left powers of \spad{a}.
        ++ Note: the polynomial never has a constant term as in general
        ++ the algebra has no unit.
      rightMinimalPolynomial : % -> SparseUnivariatePolynomial R
        ++ rightMinimalPolynomial(a) returns the polynomial determined by the
        ++ smallest non-trivial linear
        ++ combination of right powers of \spad{a}.
        ++ Note: the polynomial never has a constant term as in general
        ++ the algebra has no unit.
      leftUnits:() -> Union(Record(particular: %, basis: List %), "failed")
        ++ leftUnits() returns the affine space of all left units of the
        ++ algebra, or \spad{"failed"} if there is none.
      rightUnits:() -> Union(Record(particular: %, basis: List %), "failed")
        ++ rightUnits() returns the affine space of all right units of the
        ++ algebra, or \spad{"failed"} if there is none.
      leftUnit:() -> Union(%, "failed")
        ++ leftUnit() returns a left unit of the algebra
        ++ (not necessarily unique), or \spad{"failed"} if there is none.
      rightUnit:() -> Union(%, "failed")
        ++ rightUnit() returns a right unit of the algebra
        ++ (not necessarily unique), or \spad{"failed"} if there is none.
      unit:() -> Union(%, "failed")
        ++ unit() returns a unit of the algebra (necessarily unique),
        ++ or \spad{"failed"} if there is none.
      -- we not necessarily have a unit, hence we can't say anything
      -- about characteristic
      -- if R has CharacteristicZero then CharacteristicZero
      -- if R has CharacteristicNonZero then CharacteristicNonZero
      unitsKnown
        ++ unitsKnown means that \spadfun{recip} truly yields reciprocal
        ++ or \spad{"failed"} if not a unit,
        ++ similarly for \spadfun{leftRecip} and
        ++ \spadfun{rightRecip}. The reason is that we use left, respectively
        ++ right, minimal polynomials to decide this question.

  add
    --n := rank()
    --b := someBasis()
    --gamma : Vector Matrix R := structuralConstants b
    -- here is a problem: there seems to be a problem having local
    -- variables in the capsule of a category, furthermore
    -- see the commented code of conditionsForIdempotents, where
    -- we call structuralConstants, which also doesn't work
    -- at runtime, i.e. is not properly inherited, hence for
    -- the moment we put the code for
    -- conditionsForIdempotents, structuralConstants, unit, leftUnit,
    -- rightUnit into the domain constructor ALGSC
    V  ==> Vector
    M  ==> Matrix
    REC  ==> Record(particular: Union(V R,"failed"),basis: List V R)
    LSMP ==> LinearSystemMatrixPackage(R,V R,V R, M R)


    SUP ==>  SparseUnivariatePolynomial
    NNI ==>  NonNegativeInteger
    -- next 2 functions: use a general characteristicPolynomial
    leftCharacteristicPolynomial a ==
       n := rank()$%
       ma : Matrix R := leftRegularRepresentation(a,someBasis()$%)
       mb : Matrix SUP R := zero(n,n)
       for i in 1..n repeat
         for j in 1..n repeat
           mb(i,j):=
             i=j => monomial(ma(i,j),0)$SUP(R) - monomial(1,1)$SUP(R)
             monomial(ma(i,j),1)$SUP(R)
       determinant mb

    rightCharacteristicPolynomial a ==
       n := rank()$%
       ma : Matrix R := rightRegularRepresentation(a,someBasis()$%)
       mb : Matrix SUP R := zero(n,n)
       for i in 1..n repeat
         for j in 1..n repeat
           mb(i,j):=
             i=j => monomial(ma(i,j),0)$SUP(R) - monomial(1,1)$SUP(R)
             monomial(ma(i,j),1)$SUP(R)
       determinant mb



    leftTrace a ==
      t : R := 0
      ma : Matrix R := leftRegularRepresentation(a,someBasis()$%)
      for i in 1..rank()$% repeat
        t := t + elt(ma,i,i)
      t

    rightTrace a ==
      t : R := 0
      ma : Matrix R := rightRegularRepresentation(a,someBasis()$%)
      for i in 1..rank()$% repeat
        t := t + elt(ma,i,i)
      t

    leftNorm a == determinant leftRegularRepresentation(a,someBasis()$%)

    rightNorm a == determinant rightRegularRepresentation(a,someBasis()$%)


    antiAssociative?() ==
      b := someBasis()
      n := rank()
      for i in 1..n repeat
        for j in 1..n repeat
          for k in 1..n repeat
            not zero? ( (b.i*b.j)*b.k + b.i*(b.j*b.k) )  =>
              messagePrint("algebra is not anti-associative")$OutputForm
              return false
      messagePrint("algebra is anti-associative")$OutputForm
      true


    jordanAdmissible?() ==
      b := someBasis()
      n := rank()
      recip(2 * 1$R) case "failed" =>
        messagePrint("this algebra is not Jordan admissible, as 2 is not invertible in the ground ring")$OutputForm
        false
      for i in 1..n repeat
       for j in 1..n repeat
        for k in 1..n repeat
         for l in 1..n repeat
           not zero? ( _
             antiCommutator(antiCommutator(b.i,b.j),antiCommutator(b.l,b.k)) + _
             antiCommutator(antiCommutator(b.l,b.j),antiCommutator(b.k,b.i)) + _
             antiCommutator(antiCommutator(b.k,b.j),antiCommutator(b.i,b.l))   _
                      ) =>
               messagePrint("this algebra is not Jordan admissible")$OutputForm
               return false
      messagePrint("this algebra is Jordan admissible")$OutputForm
      true

    lieAdmissible?() ==
      n := rank()
      b := someBasis()
      for i in 1..n repeat
       for j in 1..n repeat
        for k in 1..n repeat
          not zero? (commutator(commutator(b.i,b.j),b.k) _
                  + commutator(commutator(b.j,b.k),b.i) _
                  + commutator(commutator(b.k,b.i),b.j))   =>
            messagePrint("this algebra is not Lie admissible")$OutputForm
            return false
      messagePrint("this algebra is Lie admissible")$OutputForm
      true

    -- conditionsForIdempotents b  ==
    --   n := rank()
    --   gamma : Vector Matrix R := structuralConstants b
    --   listOfNumbers : List String :=  [string(q)$String for q in 1..n]
    --   symbolsForCoef : Vector Symbol :=
    --     [concat("%", concat("x", i))::Symbol  for i in listOfNumbers]
    --   conditions : List Polynomial R := []
    --  for k in 1..n repeat
    --    xk := symbolsForCoef.k
    --    p : Polynomial R :=  monomial( - 1$Polynomial(R), [xk], [1] )
    --    for i in 1..n repeat
    --      for j in 1..n repeat
    --        xi := symbolsForCoef.i
    --        xj := symbolsForCoef.j
    --        p := p + monomial(_
    --          elt((gamma.k),i,j) :: Polynomial(R), [xi,xj], [1,1])
    --    conditions := cons(p,conditions)
    --  conditions

    structuralConstants b ==
      --n := rank()
      -- be careful with the possibility that b is not a basis
      m : NonNegativeInteger := (maxIndex b) :: NonNegativeInteger
      sC : Vector Matrix R := [new(m,m,0$R) for k in 1..m]
      for i in 1..m repeat
        for j in 1..m repeat
          covec : Vector R := coordinates(b.i * b.j, b)
          for k in 1..m repeat
             setelt( sC.k, i, j, covec.k )
      sC

    if R has IntegralDomain then

      leftRecip x ==
        zero? x => "failed"
        lu := leftUnit()
        lu case "failed" => "failed"
        b := someBasis()
        xx : % := (lu :: %)
        k  : PositiveInteger := 1
        cond : Matrix R := coordinates(xx,b) :: Matrix(R)
        listOfPowers : List % := [xx]
        while rank(cond) = k repeat
          k := k+1
          xx := xx*x
          listOfPowers := cons(xx,listOfPowers)
          cond := horizConcat(cond, coordinates(xx,b) :: Matrix(R) )
        vectorOfCoef : Vector R := (nullSpace(cond)$Matrix(R)).first
        invC := recip vectorOfCoef.1
        invC case "failed" => "failed"
        invCR : R :=  - (invC :: R)
        reduce(_+,[(invCR*vectorOfCoef.i)*power for i in _
         2..maxIndex vectorOfCoef for power in reverse listOfPowers])


      rightRecip x ==
        zero? x => "failed"
        ru := rightUnit()
        ru case "failed" => "failed"
        b := someBasis()
        xx : % := (ru :: %)
        k  : PositiveInteger := 1
        cond : Matrix R := coordinates(xx,b) :: Matrix(R)
        listOfPowers : List % := [xx]
        while rank(cond) = k repeat
          k := k+1
          xx := x*xx
          listOfPowers := cons(xx,listOfPowers)
          cond := horizConcat(cond, coordinates(xx,b) :: Matrix(R) )
        vectorOfCoef : Vector R := (nullSpace(cond)$Matrix(R)).first
        invC := recip vectorOfCoef.1
        invC case "failed" => "failed"
        invCR : R :=  - (invC :: R)
        reduce(_+,[(invCR*vectorOfCoef.i)*power for i in _
         2..maxIndex vectorOfCoef for power in reverse listOfPowers])


      recip x ==
        lrx := leftRecip x
        lrx case "failed" => "failed"
        rrx := rightRecip x
        rrx case "failed" => "failed"
        (lrx :: %) ~= (rrx :: %)  => "failed"
        lrx :: %


      leftMinimalPolynomial x ==
        zero? x =>  monomial(1$R,1)$(SparseUnivariatePolynomial R)
        b := someBasis()
        xx : % := x
        k  : PositiveInteger := 1
        cond : Matrix R := coordinates(xx,b) :: Matrix(R)
        while rank(cond) = k repeat
          k := k+1
          xx := x*xx
          cond := horizConcat(cond, coordinates(xx,b) :: Matrix(R) )
        vectorOfCoef : Vector R := (nullSpace(cond)$Matrix(R)).first
        res : SparseUnivariatePolynomial R := 0
        for i in 1..k repeat
          res := res+monomial(vectorOfCoef.i,i)$(SparseUnivariatePolynomial R)
        res

      rightMinimalPolynomial x ==
        zero? x =>  monomial(1$R,1)$(SparseUnivariatePolynomial R)
        b := someBasis()
        xx : % := x
        k  : PositiveInteger := 1
        cond : Matrix R := coordinates(xx,b) :: Matrix(R)
        while rank(cond) = k repeat
          k := k+1
          xx := xx*x
          cond := horizConcat(cond, coordinates(xx,b) :: Matrix(R) )
        vectorOfCoef : Vector R := (nullSpace(cond)$Matrix(R)).first
        res : SparseUnivariatePolynomial R := 0
        for i in 1..k repeat
          res := res+monomial(vectorOfCoef.i,i)$(SparseUnivariatePolynomial R)
        res



      associatorDependence() ==
        n := rank()
        b := someBasis()
        cond : Matrix(R) := new(n**4,6,0$R)$Matrix(R)
        z : Integer := 0
        for i in 1..n repeat
         for j in 1..n repeat
          for k in 1..n repeat
           a123 : Vector R := coordinates(associator(b.i,b.j,b.k),b)
           a231 : Vector R := coordinates(associator(b.j,b.k,b.i),b)
           a312 : Vector R := coordinates(associator(b.k,b.i,b.j),b)
           a132 : Vector R := coordinates(associator(b.i,b.k,b.j),b)
           a321 : Vector R := coordinates(associator(b.k,b.j,b.i),b)
           a213 : Vector R := coordinates(associator(b.j,b.i,b.k),b)
           for r in 1..n repeat
            z:= z+1
            setelt(cond,z,1,elt(a123,r))
            setelt(cond,z,2,elt(a231,r))
            setelt(cond,z,3,elt(a312,r))
            setelt(cond,z,4,elt(a132,r))
            setelt(cond,z,5,elt(a321,r))
            setelt(cond,z,6,elt(a213,r))
        nullSpace(cond)

    jacobiIdentity?()  ==
      n := rank()
      b := someBasis()
      for i in 1..n repeat
       for j in 1..n repeat
        for k in 1..n repeat
          not zero? ((b.i*b.j)*b.k + (b.j*b.k)*b.i + (b.k*b.i)*b.j) =>
            messagePrint("Jacobi identity does not hold")$OutputForm
            return false
      messagePrint("Jacobi identity holds")$OutputForm
      true

    lieAlgebra?()  ==
      not antiCommutative?() =>
        messagePrint("this is not a Lie algebra")$OutputForm
        false
      not jacobiIdentity?() =>
        messagePrint("this is not a Lie algebra")$OutputForm
        false
      messagePrint("this is a Lie algebra")$OutputForm
      true




    jordanAlgebra?()  ==
      b := someBasis()
      n := rank()
      recip(2 * 1$R) case "failed" =>
        messagePrint("this is not a Jordan algebra, as 2 is not invertible in the ground ring")$OutputForm
        false
      not commutative?() =>
        messagePrint("this is not a Jordan algebra")$OutputForm
        false
      for i in 1..n repeat
       for j in 1..n repeat
        for k in 1..n repeat
         for l in 1..n repeat
           not zero? (associator(b.i,b.j,b.l*b.k)+_
               associator(b.l,b.j,b.k*b.i)+associator(b.k,b.j,b.i*b.l)) =>
             messagePrint("not a Jordan algebra")$OutputForm
             return false
      messagePrint("this is a Jordan algebra")$OutputForm
      true

    noncommutativeJordanAlgebra?() ==
      b := someBasis()
      n := rank()
      recip(2 * 1$R) case "failed" =>
        messagePrint("this is not a noncommutative Jordan algebra, as 2 is not invertible in the ground ring")$OutputForm
        false
      not flexible?()$% =>
        messagePrint("this is not a noncommutative Jordan algebra, as it is not flexible")$OutputForm
        false
      not jordanAdmissible?()$% =>
        messagePrint("this is not a noncommutative Jordan algebra, as it is not Jordan admissible")$OutputForm
        false
      messagePrint("this is a noncommutative Jordan algebra")$OutputForm
      true

    antiCommutative?() ==
      b := someBasis()
      n := rank()
      for i in 1..n repeat
        for j in i..n repeat
          not zero? (i=j => b.i*b.i; b.i*b.j + b.j*b.i) =>
            messagePrint("algebra is not anti-commutative")$OutputForm
            return false
      messagePrint("algebra is anti-commutative")$OutputForm
      true

    commutative?() ==
      b := someBasis()
      n := rank()
      for i in 1..n repeat
       for j in i+1..n repeat
         not zero? commutator(b.i,b.j) =>
           messagePrint("algebra is not commutative")$OutputForm
           return false
      messagePrint("algebra is commutative")$OutputForm
      true


    associative?() ==
      b := someBasis()
      n := rank()
      for i in 1..n repeat
       for j in 1..n repeat
        for k in 1..n repeat
         not zero? associator(b.i,b.j,b.k) =>
           messagePrint("algebra is not associative")$OutputForm
           return false
      messagePrint("algebra is associative")$OutputForm
      true

    leftAlternative?() ==
      b := someBasis()
      n := rank()
      for i in 1..n repeat
       for j in 1..n repeat
        for k in 1..n repeat
         not zero? (associator(b.i,b.j,b.k) + associator(b.j,b.i,b.k)) =>
           messagePrint("algebra is not left alternative")$OutputForm
           return false
      messagePrint("algebra satisfies 2*associator(a,a,b) = 0")$OutputForm
      true

    rightAlternative?() ==
      b := someBasis()
      n := rank()
      for i in 1..n repeat
       for j in 1..n repeat
        for k in 1..n repeat
         not zero? (associator(b.i,b.j,b.k) + associator(b.i,b.k,b.j)) =>
           messagePrint("algebra is not right alternative")$OutputForm
           return false
      messagePrint("algebra satisfies 2*associator(a,b,b) = 0")$OutputForm
      true

    flexible?() ==
      b := someBasis()
      n := rank()
      for i in 1..n repeat
       for j in 1..n repeat
        for k in 1..n repeat
         not zero? (associator(b.i,b.j,b.k) + associator(b.k,b.j,b.i)) =>
           messagePrint("algebra is not flexible")$OutputForm
           return false
      messagePrint("algebra satisfies 2*associator(a,b,a) = 0")$OutputForm
      true

    alternative?() ==
      b := someBasis()
      n := rank()
      for i in 1..n repeat
       for j in 1..n repeat
        for k in 1..n repeat
         not zero? (associator(b.i,b.j,b.k) + associator(b.j,b.i,b.k)) =>
           messagePrint("algebra is not alternative")$OutputForm
           return false
         not zero? (associator(b.i,b.j,b.k) + associator(b.i,b.k,b.j)) =>
           messagePrint("algebra is not alternative")$OutputForm
           return false
      messagePrint("algebra satisfies 2*associator(a,b,b) = 0 =  2*associator(a,a,b) = 0")$OutputForm
      true

    leftDiscriminant v == determinant leftTraceMatrix v
    rightDiscriminant v == determinant rightTraceMatrix v

    coordinates(v:Vector %, b:Vector %) ==
      m := new(#v, #b, 0)$Matrix(R)
      for i in minIndex v .. maxIndex v for j in minRowIndex m .. repeat
        setRow!(m, j, coordinates(qelt(v, i), b))
      m

    represents(v, b) ==
      m := minIndex v - 1
      reduce(_+,[v(i+m) * b(i+m) for i in 1..maxIndex b])

    leftTraceMatrix v ==
      matrix [[leftTrace(v.i*v.j) for j in minIndex v..maxIndex v]$List(R)
               for i in minIndex v .. maxIndex v]$List(List R)

    rightTraceMatrix v ==
      matrix [[rightTrace(v.i*v.j) for j in minIndex v..maxIndex v]$List(R)
               for i in minIndex v .. maxIndex v]$List(List R)

    leftRegularRepresentation(x, b) ==
      m := minIndex b - 1
      matrix
       [parts coordinates(x*b(i+m),b) for i in 1..rank()]$List(List R)

    rightRegularRepresentation(x, b) ==
      m := minIndex b - 1
      matrix
       [parts coordinates(b(i+m)*x,b) for i in 1..rank()]$List(List R)

)abbrev category FRNAALG FramedNonAssociativeAlgebra
++ Author: J. Grabmeier, R. Wisbauer
++ Date Created: 01 March 1991
++ Date Last Updated: 11 June 1991
++ Basic Operations: +,-,*,**,basis
++ Related Constructors: FiniteRankNonAssociativeAlgebra, FramedAlgebra,
++   FiniteRankAssociativeAlgebra
++ Also See:
++ AMS Classifications:
++ Keywords: nonassociative algebra, basis
++ Reference:
++  R.D. Schafer: An Introduction to Nonassociative Algebras
++  Academic Press, New York, 1966
++ Description:
++   FramedNonAssociativeAlgebra(R) is a
++   \spadtype{FiniteRankNonAssociativeAlgebra} (i.e. a non associative
++   algebra over R which is a free \spad{R}-module of finite rank)
++   over a commutative ring R together with a fixed \spad{R}-module basis.
FramedNonAssociativeAlgebra(R:CommutativeRing): Category == _
  Join(FiniteRankNonAssociativeAlgebra(R),Eltable(Integer,R)) with
  --operations
    basis: () -> Vector %
      ++ basis() returns the fixed \spad{R}-module basis.
    coordinates: % -> Vector R
      ++ coordinates(a) returns the coordinates of \spad{a}
      ++ with respect to the
      ++ fixed \spad{R}-module basis.
    coordinates: Vector % -> Matrix R
      ++ coordinates([a1,...,am]) returns a matrix whose i-th row
      ++ is formed by the coordinates of \spad{ai} with respect to the
      ++ fixed \spad{R}-module basis.
    structuralConstants:() -> Vector Matrix R
      ++ structuralConstants() calculates the structural constants
      ++ \spad{[(gammaijk) for k in 1..rank()]} defined by
      ++ \spad{vi * vj = gammaij1 * v1 + ... + gammaijn * vn},
      ++ where \spad{v1},...,\spad{vn} is the fixed \spad{R}-module basis.
    conditionsForIdempotents: () -> List Polynomial R
      ++ conditionsForIdempotents() determines a complete list
      ++ of polynomial equations for the coefficients of idempotents
      ++ with respect to the fixed \spad{R}-module basis.
    represents: Vector R -> %
      ++ represents([a1,...,an]) returns \spad{a1*v1 + ... + an*vn},
      ++ where \spad{v1}, ..., \spad{vn} are the elements of the
      ++ fixed \spad{R}-module basis.
    convert: % -> Vector R
      ++ convert(a) returns the coordinates of \spad{a} with respect to the
      ++ fixed \spad{R}-module basis.
    convert: Vector R -> %
      ++ convert([a1,...,an]) returns \spad{a1*v1 + ... + an*vn},
      ++ where \spad{v1}, ..., \spad{vn} are the elements of the
      ++ fixed \spad{R}-module basis.
    leftDiscriminant : () -> R
      ++ leftDiscriminant() returns the
      ++ determinant of the \spad{n}-by-\spad{n}
      ++ matrix whose element at the \spad{i}-th row and \spad{j}-th column is
      ++ given by the left trace of the product \spad{vi*vj}, where
      ++ \spad{v1},...,\spad{vn} are the
      ++ elements of the fixed \spad{R}-module basis.
      ++ Note: the same as \spad{determinant(leftTraceMatrix())}.
    rightDiscriminant : () -> R
      ++ rightDiscriminant() returns the determinant of the \spad{n}-by-\spad{n}
      ++ matrix whose element at the \spad{i}-th row and \spad{j}-th column is
      ++ given by the right trace of the product \spad{vi*vj}, where
      ++ \spad{v1},...,\spad{vn} are the elements of
      ++ the fixed \spad{R}-module basis.
      ++ Note: the same as \spad{determinant(rightTraceMatrix())}.
    leftTraceMatrix : () -> Matrix R
      ++ leftTraceMatrix() is the \spad{n}-by-\spad{n}
      ++ matrix whose element at the \spad{i}-th row and \spad{j}-th column is
      ++ given by left trace of the product \spad{vi*vj},
      ++ where \spad{v1},...,\spad{vn} are the
      ++ elements of the fixed \spad{R}-module
      ++ basis.
    rightTraceMatrix : () -> Matrix R
      ++ rightTraceMatrix() is the \spad{n}-by-\spad{n}
      ++ matrix whose element at the \spad{i}-th row and \spad{j}-th column is
      ++ given by the right trace of the product \spad{vi*vj}, where
      ++ \spad{v1},...,\spad{vn} are the elements
      ++ of the fixed \spad{R}-module basis.
    leftRegularRepresentation : % -> Matrix R
      ++ leftRegularRepresentation(a) returns the matrix of the linear
      ++ map defined by left multiplication by \spad{a} with respect
      ++ to the fixed \spad{R}-module basis.
    rightRegularRepresentation : % -> Matrix R
      ++ rightRegularRepresentation(a) returns the matrix of the linear
      ++ map defined by right multiplication by \spad{a} with respect
      ++ to the fixed \spad{R}-module basis.
    if R has Field then
      leftRankPolynomial : () -> SparseUnivariatePolynomial Polynomial R
        ++ leftRankPolynomial() calculates the left minimal polynomial
        ++ of the generic element in the algebra,
        ++ defined by the same structural
        ++ constants over the polynomial ring in symbolic coefficients with
        ++ respect to the fixed basis.
      rightRankPolynomial : () -> SparseUnivariatePolynomial Polynomial R
        ++ rightRankPolynomial() calculates the right minimal polynomial
        ++ of the generic element in the algebra,
        ++ defined by the same structural
        ++ constants over the polynomial ring in symbolic coefficients with
        ++ respect to the fixed basis.
    apply: (Matrix R, %) -> %
      ++ apply(m,a) defines a left operation of n by n matrices
      ++ where n is the rank of the algebra in terms of matrix-vector
      ++ multiplication, this is a substitute for a left module structure.
      ++ Error: if shape of matrix doesn't fit.
    --attributes
    --attributes
      --separable <=> discriminant() ~= 0
  add

    V  ==> Vector
    M  ==> Matrix
    P  ==> Polynomial
    F  ==> Fraction
    REC  ==> Record(particular: Union(V R,"failed"),basis: List V R)
    LSMP ==> LinearSystemMatrixPackage(R,V R,V R, M R)
    CVMP ==> CoerceVectorMatrixPackage(R)

    --GA ==> GenericNonAssociativeAlgebra(R,rank()$%,_
    -- [random()$Character :: String :: Symbol for i in 1..rank()$%], _
    -- structuralConstants()$%)
    --y : GA := generic()
    if R has Field then
      leftRankPolynomial() ==
        n := rank()
        b := basis()
        gamma : Vector Matrix R := structuralConstants b
        listOfNumbers : List String :=  [string(q)$String for q in 1..n]
        symbolsForCoef : Vector Symbol :=
          [concat("%", concat("x", i))::Symbol  for i in listOfNumbers]
        xx : M P R := new(1,n,0)
        mo : P R
        x : M P R := new(1,n,0)
        for i in 1..n repeat
          mo := monomial(1, [symbolsForCoef.i], [1])$(P R)
          qsetelt!(x,1,i,mo)
        y : M P R := copy x
        k  : PositiveInteger := 1
        cond : M P R := copy x
        -- multiplication in the generic algebra means using
        -- the structural matrices as bilinear forms.
        -- left multiplication by x, we prepare for that:
        genGamma : V M P R :=  coerceP$CVMP gamma
        x := reduce(horizConcat,[x*genGamma(i) for i in 1..#genGamma])
        while rank(cond) = k repeat
          k := k+1
          for i in 1..n repeat
            setelt(xx,[1],[i],x*transpose y)
          y := copy xx
          cond := horizConcat(cond, xx)
        vectorOfCoef : Vector P R := (nullSpace(cond)$Matrix(P R)).first
        res : SparseUnivariatePolynomial P R := 0
        for i in 1..k repeat
         res := res+monomial(vectorOfCoef.i,i)$(SparseUnivariatePolynomial  P R)
        res

      rightRankPolynomial() ==
        n := rank()
        b := basis()
        gamma : Vector Matrix R := structuralConstants b
        listOfNumbers : List String :=  [string(q)$String for q in 1..n]
        symbolsForCoef : Vector Symbol :=
          [concat("%", concat("x", i))::Symbol  for i in listOfNumbers]
        xx : M P R := new(1,n,0)
        mo : P R
        x : M P R := new(1,n,0)
        for i in 1..n repeat
          mo := monomial(1, [symbolsForCoef.i], [1])$(P R)
          qsetelt!(x,1,i,mo)
        y : M P R := copy x
        k  : PositiveInteger := 1
        cond : M P R := copy x
        -- multiplication in the generic algebra means using
        -- the structural matrices as bilinear forms.
        -- left multiplication by x, we prepare for that:
        genGamma : V M P R :=  coerceP$CVMP gamma
        x := reduce(horizConcat,[genGamma(i)*transpose x for i in 1..#genGamma])
        while rank(cond) = k repeat
          k := k+1
          for i in 1..n repeat
            setelt(xx,[1],[i],y * transpose x)
          y := copy xx
          cond := horizConcat(cond, xx)
        vectorOfCoef : Vector P R := (nullSpace(cond)$Matrix(P R)).first
        res : SparseUnivariatePolynomial P R := 0
        for i in 1..k repeat
         res := res+monomial(vectorOfCoef.i,i)$(SparseUnivariatePolynomial  P R)
        res

      leftUnitsInternal : () -> REC
      leftUnitsInternal() ==
        n := rank()
        b := basis()
        gamma : Vector Matrix R := structuralConstants b
        cond : Matrix(R) := new(n**2,n,0$R)$Matrix(R)
        rhs : Vector(R) := new(n**2,0$R)$Vector(R)
        z : Integer := 0
        addOn : R := 0
        for k in 1..n repeat
         for i in 1..n repeat
           z := z+1   -- index for the rows
           addOn :=
             k=i => 1
             0
           setelt(rhs,z,addOn)$Vector(R)
           for j in 1..n repeat  -- index for the columns
             setelt(cond,z,j,elt(gamma.k,j,i))$Matrix(R)
        solve(cond,rhs)$LSMP


      leftUnit() ==
        res : REC := leftUnitsInternal()
        res.particular case "failed" =>
          messagePrint("this algebra has no left unit")$OutputForm
          "failed"
        represents (res.particular :: V R)

      leftUnits() ==
        res : REC := leftUnitsInternal()
        res.particular case "failed" =>
          messagePrint("this algebra has no left unit")$OutputForm
          "failed"
        [represents(res.particular :: V R)$%, _
          map(represents, res.basis)$ListFunctions2(Vector R, %) ]

      rightUnitsInternal : () -> REC
      rightUnitsInternal() ==
        n := rank()
        b := basis()
        gamma : Vector Matrix R := structuralConstants b
        condo : Matrix(R) := new(n**2,n,0$R)$Matrix(R)
        rhs : Vector(R) := new(n**2,0$R)$Vector(R)
        z : Integer := 0
        addOn : R := 0
        for k in 1..n repeat
         for i in 1..n repeat
           z := z+1   -- index for the rows
           addOn :=
             k=i => 1
             0
           setelt(rhs,z,addOn)$Vector(R)
           for j in 1..n repeat  -- index for the columns
             setelt(condo,z,j,elt(gamma.k,i,j))$Matrix(R)
        solve(condo,rhs)$LSMP

      rightUnit() ==
        res : REC := rightUnitsInternal()
        res.particular case "failed" =>
          messagePrint("this algebra has no right unit")$OutputForm
          "failed"
        represents (res.particular :: V R)

      rightUnits() ==
        res : REC := rightUnitsInternal()
        res.particular case "failed" =>
          messagePrint("this algebra has no right unit")$OutputForm
          "failed"
        [represents(res.particular :: V R)$%, _
          map(represents, res.basis)$ListFunctions2(Vector R, %) ]

      unit() ==
        n := rank()
        b := basis()
        gamma : Vector Matrix R := structuralConstants b
        cond : Matrix(R) := new(2*n**2,n,0$R)$Matrix(R)
        rhs : Vector(R) := new(2*n**2,0$R)$Vector(R)
        z : Integer := 0
        u : Integer := n*n
        addOn : R := 0
        for k in 1..n repeat
         for i in 1..n repeat
           z := z+1   -- index for the rows
           addOn :=
             k=i => 1
             0
           setelt(rhs,z,addOn)$Vector(R)
           setelt(rhs,u,addOn)$Vector(R)
           for j in 1..n repeat  -- index for the columns
             setelt(cond,z,j,elt(gamma.k,j,i))$Matrix(R)
             setelt(cond,u,j,elt(gamma.k,i,j))$Matrix(R)
        res : REC := solve(cond,rhs)$LSMP
        res.particular case "failed" =>
          messagePrint("this algebra has no unit")$OutputForm
          "failed"
        represents (res.particular :: V R)
    apply(m:Matrix(R),a:%) ==
      v : Vector R := coordinates(a)
      v := m *$Matrix(R) v
      convert v


    structuralConstants()   == structuralConstants basis()
    conditionsForIdempotents() == conditionsForIdempotents basis()
    convert(x:%):Vector(R)  == coordinates(x, basis())
    convert(v:Vector R):%   == represents(v, basis())
    leftTraceMatrix()       == leftTraceMatrix basis()
    rightTraceMatrix()      == rightTraceMatrix basis()
    leftDiscriminant()      == leftDiscriminant basis()
    rightDiscriminant()     == rightDiscriminant basis()
    leftRegularRepresentation x == leftRegularRepresentation(x, basis())
    rightRegularRepresentation x == rightRegularRepresentation(x, basis())
    coordinates(x: %)       == coordinates(x, basis())
    represents(v:Vector R):%== represents(v, basis())

    coordinates(v:Vector %) ==
      m := new(#v, rank(), 0)$Matrix(R)
      for i in minIndex v .. maxIndex v for j in minRowIndex m .. repeat
        setRow!(m, j, coordinates qelt(v, i))
      m