--* From BMT%WATSON.vnet.ibm.com@yktvmh.watson.ibm.com  Fri Oct 22 00:36:48 1993
--* Received: from yktvmh.watson.ibm.com by radical.watson.ibm.com (AIX 3.2/UCB 5.64/900524)
--*           id AA12200; Fri, 22 Oct 1993 00:36:48 -0400
--* Received: from watson.vnet.ibm.com by yktvmh.watson.ibm.com (IBM VM SMTP V2R3)
--*    with BSMTP id 4018; Fri, 22 Oct 93 00:43:17 EDT
--* Received: from YKTVMH by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id <A.BMT.NOTE.VAGENT2.1977.Oct.22.00:43:16.-0400>
--*           for asbugs@watson; Fri, 22 Oct 93 00:43:17 -0400
--* Received: from YKTVMH by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id 1974; Fri, 22 Oct 1993 00:43:15 EDT
--* Received: from cyst.watson.ibm.com by yktvmh.watson.ibm.com (IBM VM SMTP V2R3)
--*    with TCP; Fri, 22 Oct 93 00:43:14 EDT
--* Received: from spadserv.watson.ibm.com by cyst.watson.ibm.com (AIX 3.2/UCB 5.64/900528)
--*   id AA29794; Fri, 22 Oct 1993 00:42:46 -0400
--* Received: by spadserv.watson.ibm.com (AIX 3.2/UCB 5.64/900524)
--*           id AA26390; Fri, 22 Oct 1993 00:44:49 -0400
--* Date: Fri, 22 Oct 1993 00:44:49 -0400
--* From: bmt@spadserv.watson.ibm.com
--* X-External-Networks: yes
--* Message-Id: <9310220444.AA26390@spadserv.watson.ibm.com>
--* To: asbugs@watson.ibm.com
--* Subject: import in default clause for category causes syntax errors [algcat.as][32.1 (current)]

--@ Fixed  by:  PI   Fri Mar 11 15:44:14 EST 1994 
--@ Tested by:  none 
--@ Summary:    see file 
		of the "then"/"else" must be wrapped with braces. It's a syntax
		constraint.


#include "aslib.as"
FINRALG ==> FiniteRankAlgebra

FINRALG_- ==> FiniteRankAlgebra_&

FiniteRankAlgebra(R: CommutativeRing,UP: UnivariatePolynomialCategory R)== with
        {
          Algebra R;
          rank: () -> PositiveInteger;
                  ++ `rank()' returns the rank of the algebra.

          regularRepresentation: (%,Vector %) -> Matrix R;
                  ++ `regularRepresentation(a,basis)' returns the matrix of the
                  ++ linear map defined by left multiplication by `a' with
                  ++ respect to the `basis' `basis'.

          trace: % -> R;
                  ++ `trace(a)' returns the trace of the regular representation
                  ++ of `a' with respect to any basis.

          norm: % -> R;
                  ++ `norm(a)' returns the determinant of the regular
                  ++ representation of `a' with respect to any basis.

          coordinates: (%,Vector %) -> Vector R;
                  ++ `coordinates(a,basis)' returns the coordinates of `a' with
                  ++ respect to the `basis' `basis'.

          coordinates: (Vector %,Vector %) -> Matrix R;
                  ++ `coordinates([v1,...,vm], basis)' returns the coordinates
                  ++ of the `vi'`'s' with to the basis `basis'. The coordinates
                  ++ of `vi' are contained in the `i'th row of the matrix
                  ++ returned by this function.

          represents: (Vector R,Vector %) -> %;
                  ++ `represents([a1,..,an],[v1,..,vn])' returns `a1*v1 + ... +
                  ++ an*vn'.

          discriminant: Vector % -> R;
                  ++ `discriminant([v1,..,vn])' returns
                  ++ `determinant(traceMatrix([v1,..,vn]))'.

          traceMatrix: Vector % -> Matrix R;
                  ++ `traceMatrix([v1,..,vn])' is the `n'-by-`n' matrix (
                  ++ `Tr'(`vi' * `vj') )

          characteristicPolynomial: % -> UP;
                  ++ `characteristicPolynomial(a)' returns the characteristic
                  ++ polynomial of the regular representation of `a' with
                  ++ respect to any basis.

          if R has Field then
            minimalPolynomial: % -> UP;
                    ++ `minimalPolynomial(a)' returns the minimal polynomial of
                    ++ `a'.

          if R has CharacteristicZero then
            CharacteristicZero;
          if R has CharacteristicNonZero then
            CharacteristicNonZero;
          default {
                  import UP;
                  import R;
                  import SingleInteger;
                  import Segment Integer;
                  import Vector %;
                  import List List R;
                  discriminant (v:Vector %):R == {
                          import Matrix R;
                          determinant traceMatrix v;
                  }

                  coordinates(v:Vector %,b:Vector %):Matrix R == {
                          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:Vector R,b:Vector %):% == {
                          local i:SingleInteger;
                          m := minIndex v - 1;
                          reduce(+,
                            [v (i::Integer + m) * b (i::Integer + m)
                              for (free i) in 1..retract rank]);
                  }

                  traceMatrix (v:Vector %):Matrix R == {
                          import List R;
                          matrix (
                            [[trace(v(i) * v(j))
                              for j in minIndex(v)..maxIndex(v)]$List(R) for i in minIndex(v)..maxIndex(v)]$List(List(R)));
                  }

                  regularRepresentation(x:%,b:Vector %):Matrix R == {
                          import Integer;
                          import Vector R;
                          local i:SingleInteger;
                          m := minIndex b - 1;
                          matrix (
                            [parts(coordinates(x * b(i::Integer + m),b))
                              for (free i) in 1..retract(rank)]$List(List(R)));
                  }

          }
}

FRAMALG ==> FramedAlgebra

FRAMALG_- ==> FramedAlgebra_&

FramedAlgebra(R: CommutativeRing,UP: UnivariatePolynomialCategory R)== with {
        FiniteRankAlgebra(R,UP);
        basis: () -> Vector %;
                ++ `basis()' returns the fixed `R'-module basis.

        coordinates: % -> Vector R;
                ++ `coordinates(a)' returns the coordinates of `a' with respect
                ++ to the fixed `R'-module basis.

        coordinates: Vector % -> Matrix R;
                ++ `coordinates([v1,...,vm])' returns the coordinates of the
                ++ `vi'`'s' with to the fixed basis. The coordinates of `vi' are
                ++ contained in the `i'th row of the matrix returned by this
                ++ function.

        represents: Vector R -> %;
                ++ `represents([a1,..,an])' returns `a1*v1 + ... + an*vn', where
                ++ `v1', ..., `vn' are the elements of the fixed basis.

        convert: % -> Vector R;
                ++ `convert(a)' returns the coordinates of `a' with respect to
                ++ the fixed `R'-module basis.

        convert: Vector R -> %;
                ++ `convert([a1,..,an])' returns `a1*v1 + ... + an*vn', where
                ++ `v1', ..., `vn' are the elements of the fixed basis.

        traceMatrix: () -> Matrix R;
                ++ `traceMatrix()' is the `n'-by-`n' matrix ( `Tr(\spad{vi' *
                ++ vj)} ), where `v1', ..., `vn' are the elements of the fixed
                ++ basis.

        discriminant: () -> R;
                ++ `discriminant()' = determinant(traceMatrix()).

        regularRepresentation: % -> Matrix R;
                ++ `regularRepresentation(a)' returns the matrix of the linear
                ++ map defined by left multiplication by `a' with respect to the
                ++ fixed basis.

        default {
                import UP;
                import R;
                convert (x:%):Vector R == coordinates x;
                convert (v:Vector R):% == represents v;
                traceMatrix:Matrix R == traceMatrix basis;
                discriminant:R == discriminant basis;
                regularRepresentation (x:%):Matrix R ==
                  regularRepresentation(x,basis);

                coordinates (x:Vector %):Matrix R == coordinates(x,basis);
                represents (x:Vector R):% == represents(x,basis);
                coordinates (v:Vector %):Matrix R == {
                        import Vector %;
                        import Segment Integer;
                        m := new(#v,rank::NonNegativeInteger,0)$Matrix(R);
                        for i in minIndex v..maxIndex v for j in minRowIndex m..
                           repeat setRow!(m,j,coordinates qelt(v,i));
                        m;
                }

                regularRepresentation (x:%):Matrix R == {
                        import Vector %;
                        import Segment Integer;
                        m :=
                          new(n := rank::NonNegativeInteger,n::
                            NonNegativeInteger,0)$Matrix(R);
                        b := basis;
                        for i in minIndex b..maxIndex b for j in minRowIndex m..
                           repeat setRow!(m,j,coordinates (x * qelt(b,i)));
                        m;
                }

                characteristicPolynomial (x:%):UP == {
                        import
                          MatrixCategoryFunctions2(R,Vector R,Vector R,Matrix R,
                            UP,Vector UP,Vector UP,Matrix UP);
                        import NonNegativeInteger;
                        mat00 := regularRepresentation x;
                        mat0 :=
                          map(((%1:R):UP) +-> %1::UP,mat00)$_
                              MatrixCategoryFunctions2(R,Vector(R),Vector(R),
                                Matrix(R),UP,Vector(UP),Vector(UP),Matrix(UP));
                        mat1:Matrix UP :=
                          scalarMatrix(rank::NonNegativeInteger,monomial(1,1)$UP
                            );
                        determinant (mat1 - mat0);
                }

        }
}

MONOGEN ==> MonogenicAlgebra

MONOGEN_- ==> MonogenicAlgebra_&

MonogenicAlgebra(R: CommutativeRing,UP: UnivariatePolynomialCategory R)== with {
        FramedAlgebra(R,UP);
        CommutativeRing;
        ConvertibleTo UP;
        FullyRetractableTo R;
        FullyLinearlyExplicitRingOver R;
        generator: () -> %;
                ++ `generator()' returns the generator for this domain.

        definingPolynomial: () -> UP;
                ++ `definingPolynomial()' returns the minimal polynomial which
                ++ `generator()' satisfies.

        reduce: UP -> %;
                ++ `reduce(up)' converts the univariate polynomial `up' to an
                ++ algebra element, reducing by the `definingPolynomial()' if
                ++ necessary.

        convert: UP -> %;
                ++ `convert(up)' converts the univariate polynomial `up' to an
                ++ algebra element, reducing by the `definingPolynomial()' if
                ++ necessary.

        lift: % -> UP;
                ++ `lift(z)' returns a minimal degree univariate polynomial up
                ++ such that `z=reduce up'.

        if R has Finite then
          Finite;
        if R has Field then
          Field;
          DifferentialExtension R;
          reduce: Fraction UP -> Union(%,Enumeration failed);
                  ++ `reduce(frac)' converts the fraction `frac' to an algebra
                  ++ element.

          derivationCoordinates: (Vector %,R -> R) -> Matrix R;
                  ++ `derivationCoordinates(b, ')' returns `M' such that `b' = M
                  ++ b'.

        if R has FiniteFieldCategory then
          FiniteFieldCategory;
        default {
                import UP;
                import R;
                convert (x:%):UP == lift x;
                convert (p:UP):% == reduce p;
                generator:% == {
                        import NonNegativeInteger;
                        reduce monomial(1,1)$UP;
                }

                norm (x:%):R == resultant(definingPolynomial,lift x);
                retract (x:%):R == retract lift x;
                retractIfCan (x:%):Union(R,Enumeration failed) ==
                  retractIfCan lift x;

                basis:Vector % == {
                        import Integer;
                        import SingleInteger;
                        import Vector %;
                        local i:SingleInteger;
                        [reduce monomial(1,i::NonNegativeInteger)$UP
                          for (free i) in 0..
                            retract ((rank - 1)::NonNegativeInteger)];
                }

                characteristicPolynomial (x:%):UP == {
                        import
                          CharacteristicPolynomialInMonogenicalAlgebra(R,UP,%);
                        characteristicPolynomial(x)$_
                            CharacteristicPolynomialInMonogenicalAlgebra(R,UP,%);
                }

                if R has Finite then
                  size:NonNegativeInteger == size$R ^ rank;
                  random:% == {
                          import SingleInteger;
                          import Vector R;
                          local i:SingleInteger;
                          represents (
                            [random$R for (free i) in 1..retract(rank)]$Vector(R));
                  }

                if R has Field then
                  reduce (x:Fraction UP):Union(%,Enumeration failed) ==
                    exquo(reduce numer x,reduce denom x);

                  differentiate(x:%,d:R -> R):% == {
                          p := definingPolynomial;
                          yprime := - reduce map(d,p) / reduce differentiate p;
                          reduce map(d,lift x) + yprime *
                            reduce differentiate lift x;
                  }

                  derivationCoordinates(b:Vector %,d:R -> R):Matrix R == {
                          import Vector %;
                          coordinates(
                            map((local (local %1:%):%) +-> differentiate(%1,d),b
                              ),b);
                  }

                  recip (x:%):Union(%,Enumeration failed) == {
                          import Record(coef1: UP,coef2: UP);
                          (bc := extendedEuclidean(lift x,definingPolynomial,1))
                             case failed => failed::Union(%,Enumeration failed);
                          reduce (bc::Record(coef1: UP,coef2: UP)).coef1::
                            Union(%,Enumeration failed);
                  }

        }
}

CPIMA ==> CharacteristicPolynomialInMonogenicalAlgebra

Pol ==> SparseUnivariatePolynomial

CharacteristicPolynomialInMonogenicalAlgebra(R: CommutativeRing,PolR:
  UnivariatePolynomialCategory R,E: MonogenicAlgebra(R,PolR)): with {
        characteristicPolynomial: E -> PolR;
}
== add {
        import E;
        import PolR;
        import R;
        import
          UnivariatePolynomialCategoryFunctions2(R,PolR,PolR,
            SparseUnivariatePolynomial PolR);
        import NonNegativeInteger;
        import SparseUnivariatePolynomial PolR;
        import
          UnivariatePolynomialCategoryFunctions2(R,PolR,PolR,
            SparseUnivariatePolynomial PolR);
        XtoY (Q:PolR):SparseUnivariatePolynomial PolR ==
          map(((%1:R):PolR) +-> monomial(%1,0),Q);

        P := XtoY definingPolynomial$E;
        X := monomial(monomial(1,1)$PolR,0);
        characteristicPolynomial (x:E):PolR == {
                Qx:PolR := lift x;
                return resultant(P,X - XtoY Qx);
        }

}