#
# Main file for Gauss
# Share Library Version
#
# After reading in the code in this file, before using Gauss,
# the user must assign the global variable gauss the directory of where all the
# gauss codes are kept, e.g.  gauss := `/usr/local/maple/share/gauss`;
# Then the user should execute Gauss(); if s/he wants to use the standard
# abbreviations which are listed below in the code.  See ?Gauss for help
#
# Author: MBM 89-92
#

Gauss := proc(dir)

if not assigned(gauss) then 
    if nargs = 1 and type(dir,string) then gauss := dir
    elif assigned(sharename) then gauss := ``.sharename.`/gauss`;
    else
ERROR(`please assign "gauss" a string which is the path of the gauss directory`)
    fi;
fi;

macro(SM='SquareMatrix'):
macro(GF='GaloisField'):
macro(DUP='DenseUnivariatePolynomial'):
macro(MUP='MapleUnivariatePolynomial'):
macro(OUP='OrderedUnivariatePolynomial'):
macro(SAE='AlgebraicExtension'):
macro(DEV='DenseExponentVector'):
macro(PEV='PrimeExponentVector'):
macro(MEV='MapleExponentVector'):
macro(TEV='MacaulayExponentVector'):
macro(SDMP='SparseDistributedMultivariatePolynomial'):
macro(MMP='MapleMultivariatePolynomial'):
macro(TDMP='TableDistributedMultivariatePolynomial'):
macro(SMP='SparseMultivariatePolynomial'):
macro(UPS='UnivariatePowerSeries'):
macro(LUPS='LazyUnivariatePowerSeries'):
macro(QF='QuotientField'):
macro(RF='RationalFunction'):
macro(ENF='ExpandedNormalForm'):
macro(FNF='FactoredNormalForm'):
macro(GB='GrobnerBasis');

Z := Integer():
Q := Rational():

lprint(`----------------------- Gauss version 1.0 -----------------------`);
lprint(`Initially defined domains are Z and Q the integers and rationals.`);
lprint(`Abbreviations, e.g. DUP for DenseUnivariatePolynomial, also made.`);
interface(prompt=`Gauss >> `);

NULL

end:


#
# Basic routines for Guass
# This code defines a Domain or Category to be a table.
# Author: MBM 1989
#
newDomain := newCategory:
newCategory := proc(x)
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    if nargs = 1 then x else table([
    Domain=true,        # use for system checking
    Categories={},        # defines the algebraic views of a Domain
    Properties={},      # properties of the operations and views
    Signatures=table()  # types of the operations 
    ]) fi
end:

notImplemented := proc() 
option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
ERROR(`operation is not implemented`) end:
`type/Domain` := proc(d) 
option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
type(d,table) and d[Domain] end:
addCategory := proc(x,c) 
option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
x[Categories] := x[Categories] union {c} end:
addProperties := addProperty:
addProperty := proc(x) local p;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    if type(args[2],set) then p := args[2] else p := {args[2]} fi;
    x[Properties] := x[Properties] union p
end:

defOperations := defOperation:
defOperation := proc(operation,signature,D)
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    if type(operation,set) then map(defOperation,operation,signature,D)
    else # Note signature may be a set of signatures
        if not assigned(D[operation]) then D[operation] := notImplemented fi;
        D[Signatures][operation] := subs(D=`$`,signature)
    fi;
end:

hasCategory := proc(x,c) 
option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
member(c,x[Categories]) end:
hasCategories := proc(x,s) local c;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    for c in s do if not member(c,x[Categories]) then RETURN(false) fi od;
    true
end:

hasProperty := proc(x,a) 
option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
member(a,x[Properties]) end:
hasProperties := proc(x,a) local p;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    for p in a do if not member(p,x[Properties]) then RETURN(false) fi od;
    true
end:

hasOperation := proc(x,o) 
option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
assigned(x[o]) end:
hasOperations := proc(x,s) local o;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    for o in s do if not assigned(x[o]) then RETURN(false) fi od;
    true
end:

isDomain := proc(f,d) 
option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
evalb(f[DomainName] = d) end:
isImplemented := proc(d,f) 
option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
assigned(d[f]) and d[f] <> notImplemented end:

#
# Basic catefory definitions for Gauss
# Author: MBM 1989
#
Set := proc() local S;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    if nargs = 0 then S := newCategory() else S := args[1] fi;
    if hasCategory(S,Set) then RETURN(op(S)) fi;
    addCategory(S,Set);
    defOperations( {`=`,`<>`}, [S,S] &-> Boolean, S );
    defOperation( Random, [] &-> S, S );
    defOperation( Input, Expression &-> Union(S,FAIL), S );
    defOperation( Output, S &-> Expression, S );
    defOperation( Type, Expression &-> Boolean, S );
    S[`<>`] := subs(D = S,proc(x,y) not D[`=`](x,y) end);
    S[Output] := <x>; # use Maple by default
    op(S)
end:

Poset := PartiallyOrderedSet:
PartiallyOrderedSet := proc() local P;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    P := Set(args);
    addCategory(P,PartiallyOrderedSet);
    defOperations({`<`,`>`,`<=`,`>=`}, [P,P] &-> Union(Boolean,FAIL), P);
    addProperties( P, {
        Reflexive(`<=`),
        Transitive(`<=`),
        Antisymmetric(`<=`) } );
    op(P)
end:

OrderedSet := proc() local O;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    O := PartiallyOrderedSet(args);
    addCategory(O,OrderedSet);
    defOperations( {`<`,`>`,`<=`,`>=`}, [O,O] &-> Boolean, O );
    defOperations( {Max,Min}, [O,Nary(O)] &-> O, O );
    O[`>`] := subs('D' = O, proc(x,y) not (D[`=`](x,y) or D[`<`](x,y)) end);
    O[`>=`] := subs('D' = O, proc(x,y) not D[`<`](x,y) end);
    O[`<=`] := subs(`<` = `>`, O[`>=`]);
    O[Min] := subs('D' = O, proc(x) local k,m;
        m := x; # Min and Max are nary : $+ -> $
        for k from 2 to nargs do
            if D[`<`](args[k],m) then m := args[k] fi
        od;
        m
        end);
    O[Max] := subs(`<` = `>`, eval(O[Min]) );
    op(O)
end:

# Compute x1 o x2 o ... o xn
RepeatedSquaring := proc(Op,Inv,Identity,x,n) local e,r,y,z;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    y := Identity;
    if n = 1 then RETURN(x) elif n = 0 then RETURN(y) fi;
    if n < 0 then
        z := Inv(x);
        if z = FAIL then ERROR(`unable to compute inverse`) fi;
        e := -n
    else    z := x; e := n;
    fi;
    do  # binary exponentiation
        e := iquo(e,2,'r');
        if r = 1 then
        y := Op(z,y);
        if e = 0 then RETURN(y) fi
        fi;
        z := Op(z,z)
    od
end:

SemiGroup := proc() local G;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    G := Set(args);
    addCategory( G, SemiGroup );
    addProperty( G, Associative(`+`) );
    defOperation( `+`, [G,Nary(G)] &-> G, G );
    op(G)
end:

Monoid := proc() local M;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    M := SemiGroup(args);
    if hasCategory(M,Monoid) then RETURN(op(M)) fi;
    addCategory(M,Monoid);
    defOperation( 0, M, M );    # the additive identity
    addProperty( M, NormalForm );    # must be a unique constant
    op(M)
end:

AbelianMonoid := proc() local M;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    M := Monoid(args);
    addCategory(M,AbelianMonoid);
    addProperty(M,Commutative(`+`));
    op(M)
end:

OrderedAbelianMonoid := proc() local M;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    M := AbelianMonoid( OrderedSet(args) );
    addCategory(M,OrderedAbelianMonoid);
    op(M)
end:

Group := proc() local G;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    G := Monoid(args);
    addCategory(G,Group);
    defOperation( `-`, {G &-> G, [G,G] &-> G}, G );
    defOperation( `*`, [Integer,G] &-> G, G );
    G[`=`] := subs( 'D' = G, proc(a,b) evalb(D[`-`](a,b) = D[0]) end );
    G[`*`] := subs( 'D' = G, proc(n,x)
        RepeatedSquaring(D[`+`],D[`-`],D[0],x,n)
        end);
    op(G)
end:

AbelianGroup := proc() local G;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    G := AbelianMonoid(Group(args));
    addCategory(G,AbelianGroup);
    op(G)
end:

OrderedAbelianGroup := proc() local G;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    G := AbelianGroup( OrderedAbelianMonoid(args) );
    addCategory(G,OrderedAbelianGroup);
    op(G)
end:

Ring := proc() local R;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    R := AbelianGroup(args);
    defOperation( 1, R, R );
    defOperation( `*`, { [Integer,R] &-> R, [R,Nary(R)] &-> R }, R );
    defOperation( Inv, R &-> Union(R,FAIL), R );
    defOperation( `^`, [R,Integer] &-> R, R );
    addProperty( R, Associative(`*`) );
    addProperty( R, Distributes(`*`,`+`) );
    defOperation( Coerce, Integer &-> R, R );
    defOperation( Characteristic, Integer, R );
    addCategory(R,Ring);
    R[Coerce] := subs( 'D' = R, proc(n) D[`*`](n,D[1]) end );
    R[`-`] := subs( 'D' = R, proc(a,b)
        if nargs = 1 then D[`*`](-1,a) else D[`+`](a,D[`*`](-1,b)) fi
        end);
    R[`^`] := subs( 'D' = R, proc(x,n)
        RepeatedSquaring(D[`*`],D[Inv],D[1],x,n)
        end);
    op(R)
end:

CommutativeRing := proc() local R;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    R := Ring(args);
    addCategory(R,CommutativeRing);
    addProperty(R,Commutative(`*`));
    op(R);
end:

IntegralDomain := proc() local I;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    I := CommutativeRing(args);
    addCategory(I,IntegralDomain);
    addProperty(I,NoZeroDivisors);
    defOperations( {Unit,Normal}, I &-> I, I );
    defOperation( UnitNormal, [I] &-> [I,I,I], I );
    defOperation( Div, [I,I] &-> Union(I,FAIL), I );
    I[Normal] := subs('D' = I, proc(x)
        D[`*`](D[Inv](D[Unit](x)),x)
        end);
    I[UnitNormal] := subs('D' = I, proc(x) local u,a;
        u := D[Unit](x); a := D[Inv](u);
        [ u, D[`*`](a,x), a ]
        end);
    op(I)
end:


OrderedDomain := proc() local O;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    O := IntegralDomain( OrderedSet() );
    defOperation( Sign, O &-> UNION(1,-1,0), O );
    defOperation( Zero, O &-> Boolean, O );
    O[Zero] := subs( 'OO'=O, proc(x) evalb( OO[Sign](x) = 0 ) end );
    addCategory(O,OrderedDomain);
    op(O)
end:

GcdDomain := proc() local G;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    G := IntegralDomain(args);
    addCategory(G,GcdDomain);
    defOperations( {Gcd,Lcm}, Nary(G) &-> G, G );
    defOperation( RelativelyPrime, [G,G] &-> Boolean, G );
    G[RelativelyPrime] := subs( 'D' = G,
        proc(a,b) evalb( D[Gcd](a,b) = D[1] )
        end);
    G[Lcm] := subs( 'D' = G, proc(x) local l,y;
        if nargs = 0 then RETURN( D[0] ) else l := x fi;
        for y in {args[2..nargs]} do
            if y = D[0] then RETURN(y) fi;
            l := D[Div](D[`*`](y,l),D[Gcd](y,l))
        od;
        D[Normal](l)
        end);
    op(G)
end:

UniqueFactorizationDomain := proc() local U;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    U := GcdDomain(args);
    defOperation( Prime, U &-> Boolean, U );
    defOperations( {Factor,Sqrfree}, U &-> List(U,Nary(List(U,U))), U );
    addCategory(U,UniqueFactorizationDomain);
    op(U)
end:
    
EuclideanDomain := proc() local E;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    if nargs = 0 then E := newCategory() else E := args[1] fi;
    E := UniqueFactorizationDomain(E);
    addCategory(E,EuclideanDomain);
    defOperation( EuclideanNorm, E &-> Integer, E );
    defOperation( SmallerEuclideanNorm, [E,E] &-> Boolean, E );
    defOperations( {Rem,Quo}, {[E,E] &-> E,[E,E,Name] &-> E}, E );
    defOperation( Gcdex, {[E,E,Name,Name] &-> E,[E,E,Name] &-> E}, E );
    defOperation( Powmod, [E,Integer,E] &-> E, E );
    E[Quo] := subs('D' = E, proc(x,y,r) local t,q;
        t := D[Rem](x,y,q); if nargs = 3 then r := t fi; q
        end);
    E[Div] := subs('D' = E, proc(x,y) local q;
        if D[Rem](x,y,q) <> D[0] then FAIL else q fi
        end);
    E[SmallerEuclideanNorm] := subs('D' = E, proc(x,y)
        evalb( D[EuclideanNorm](x) < D[EuclideanNorm](y) )
        end);
    E[Powmod] := subs('D' = E, proc() PowerRemainder(D,args) end);
    E[Gcd] := subs('D' = E, proc() EuclideanAlgorithm(D,args) end);
    E[Gcdex] := subs('D' = E, proc() PrincipalIdeal(D,args) end);
    op(E)
end:

EuclideanAlgorithm := proc(E) local a,b,r,s;

    # Compute Gcd(x1,x2,...,xn)
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    s := {args[2..nargs]} minus {E[0]};
    if s = {} then RETURN(E[0]) fi;
    sort( [op(s)], E[SmallerEuclideanNorm] );

    a := E[Normal](s[1]);
    for b in subsop(1=NULL,s) while a <> E[1] do
        b := E[Normal](b);
        while b <> E[0] do
        r := E[Normal](E[Rem](a,b));
        a := b;
        b := r
        od;
    od;
    a
end:

PrincipalIdeal := proc(E,x,y,s,t) local c,d,c1,c2,d1,d2,g,r,r1,r2,q;
    # It is assumed that E is a Euclidean Domain
    # Given x and y solve g = s * a + t * b = Gcd(x,y) for s, t, and g
    # Algorithm 2.2 from Geddes textbook
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    c := x; c1 := E[1]; c2 := E[0];
    d := y; d1 := E[0]; d2 := E[1];
    while d <> E[0] do
        r := E[Rem](c,d,'q'); c := d; d := r;
        r1 := E[`-`](c1,E[`*`](q,d1)); c1 := d1; d1 := r1;
        if nargs > 4 then
        r2 := E[`-`](c2,E[`*`](q,d2)); c2 := d2; d2 := r2
        fi
    od;
    g := E[UnitNormal](c);
    s := E[`*`](g[3],c1);
    if nargs > 4 then t := E[`*`](g[3],c2) fi;
    g[2]
end:

PrincipalIdeal := proc(E,x,y,s,t) local ax,ay,c1,c2,c3,d1,d2,d3,r1,r2,r3,q;
    # Given x and y solve g = s * a + t * b = Gcd(x,y) for s, t, and g
    # Note: this is the normalized half extended Euclidean algorithm
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    ax := E[UnitNormal](x); c3 := ax[2]; ax := ax[3];
    ay := E[UnitNormal](y); d3 := ay[2]; ay := ay[3];
    c1 := E[1]; d1 := E[0]; c2 := E[0]; d2 := E[1];
    while d3 <> E[0] do
        r3 := E[UnitNormal]( E[Rem](c3,d3,'q') );
        r1 := E[`*`](r3[3],E[`-`](c1,E[`*`](q,d1)));
        if nargs > 4 then r2 := E[`*`](r3[3],E[`-`](c2,E[`*`](q,d2))) fi;
        c1 := d1; d1 := r1;
        c2 := d2; d2 := r2;
        c3 := d3; d3 := r3[2];
    od;
    s := E[`*`](ax,c1);
    if nargs > 4 then t := E[`*`](ay,c2) fi;
    c3
end:

PowerRemainder := proc(E,a,n,b) local d,e,y,z;
    # It is assumed that E is a Euclidean Domain
    # Compute Rem(a^n,b) using binary powering
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    if not type(n,integer) or n < 0 then
        ERROR(`2nd argument must be a non-negative integer`) fi;
    d := E[Normal](b);
    z := E[Rem](a,d);
    if n = 0 then
        if z = E[0] then ERROR(`0^0 is undefined`) fi;
        RETURN( E[1] )
    elif n = 1 then RETURN(z)
    elif n < 0 then
        if E[`=`]( E[Gcdex](z,d,'z'), E[1] ) then e := -n else
        ERROR(`unable to compute inverse`)
        fi;
    else e := n
    fi;
    y := E[1];
    do  # binary powering
        if irem(e,2,'e') = 1 then y := E[Rem]( E[`*`](z,y), d ) fi;
        if e = 0 then RETURN(y) fi;
        z := E[Rem]( E[`*`](z,z), d )
    od;
end:

Field := proc() local F,R;

    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    if nargs = 0 then F := newCategory() else F := args[1] fi;
    if hasCategory(F,Field) then RETURN(op(F)) fi;

    F := EuclideanDomain(F);
    addCategory( F, Field );
    defOperation( AbsoluteDegree, Integer, F );
    defOperation( Inv, F &-> F, F );
    defOperation( `/`, {[F,Integer] &-> F, [F,F] &-> F}, F );
    F[`/`] := subs('D' = F, proc(a,b)
        if type(b,integer) then D[`*`](a, D[Inv](D[Coerce](b)))
        else D[`*`](a, D[Inv](b)) fi
        end);
    F[Div] := subs('D' = F, proc(a,b) D[`/`](a,b) end);
    F[Norm] := proc() 
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    0 end;
    F[Rem] := subs('D' = F,proc(x,y,q)
        if nargs = 3 then q := D[`/`](x,y) fi; D[0]
        end);
    F[Quo] := subs('D' = F,proc(x,y,r)
        if nargs = 3 then r := D[0] fi; D[`/`](x,y)
        end);
    F[EuclideanNorm] := subs('D' = F, proc(x)
        if x = D[0] then D[0] else D[1] fi
        end);
    F[Unit] := subs('D' = F, proc(x) if x = D[0] then D[1] else x fi end);
    F[Normal] := subs('D' = F, proc(x) if x = D[0] then x else D[1] fi end);
    F[UnitNormal] := subs('D' = F, proc(x)
        if x = D[0] then [D[1],x,D[1]] else [x,D[1],D[Inv](x)] fi
        end);
    F[Factor] := <[x,[]]>;
    F[Sqrfree] := <[x,[]]>;
    F[Prime] := proc() 
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    false end;
    op(F)
end:

OrderedField := proc() local O;
        option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
        O := Field( OrderedDomain() );
        addCategory(O,OrderedField);
        op(O)
end:

Generator := proc(n) option remember,
    `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
 rand(n) end:
Finite := proc() local F,p,G;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    F := NULL;
    if nargs = 0 then
    elif type([args],[integer]) then p := args[1];
    elif type([args],[Domain]) then F := args[1];
    elif type([args],[integer,Domain]) then p := args[1]; F := args[2]
    elif type([args],[Domain,integer]) then p := args[2]; F := args[3]
    else ERROR(`bad arguments`)
    fi;
    F := Set(F);
    defOperation( Size, Integer, F );
    defOperation( Index, Integer &-> F, F );
    defOperation( Universe, [] &-> Nary(F), F );
    defOperation( Lookup, F &-> Integer, F );
    if assigned(p) then F[Size] := p fi;
    F[Index] := subs('D' = F, proc(k) op(k,[D[Universe]()]) end);
    F[Lookup] := subs('D' = F, proc(x) local p;
        if member(x,[D[Universe]()],p) then p
            else ERROR(`Finite: bad argument to Lookup`) fi
        end);
    F[Random] := subs('D' = F,proc() D[Index](Generator(D[Size])()) end);
    F[Universe] := subs('D' = F, proc() local k; option remember;
        seq(D[Index](k), k=0..D[Size]-1)
        end);
    addCategory(F,Finite);
    op(F)
end:

TranscendentalFunctions := proc() local T;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    T := Field(args);
    defOperations( {Exp,Ln,Log,Sin,Cos,Tan,Sinh,Cosh,Tanh},
        T &-> Union(T,FAIL), T );
    T[Characteristic] := 0;
    T[Log] := subs( 'G'=T, proc(x) G[Ln](x) end );
    T[Tan] := subs( 'G'=T, proc(x) G[Sin](x)/G[Cos](x) end );
    T[Sinh] := subs( 'G'=T, proc(x) local t;
        t := G[Exp](x);
        G[`/`](G[`-`](t,G[Inv](t)),2)
        end);
    T[Cosh] := subs( 'G'=T, proc(x) local t;
        t := G[Exp](x);
        G[`/`](G[`+`](t,G[Inv](t)),2)
        end);
    T[Tanh] := subs( 'G'=T, proc(x) local t;
                t := G[`^`](G[Exp](x),2);
                G[`/`](G[`-`](t,G[1]),G[`+`](t,G[1]))
            end);
        op(T)
end:

DifferentialField := proc() local D;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    D := Field(args);
    addCategory(D, DifferentialField);
    defOperation( Diff, D &-> D, D );

    if D[Characteristic]=0 then
       defOperation( Integrate, D &-> Union(D, FAIL), D );
    fi;
    op(D)
end:

FiniteField := proc(p,k) local F;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    F := Field(Finite(p^k,args[3..nargs]));
    F[Characteristic] := p;
    F[AbsoluteDegree] := k;
    addCategory(F,FiniteField);
    op(F)
end:

#
# Initial definitions for Gauss
# Author: MBM 89 - 92
#

# These definitions should go internal in Maple V release III
`+` := proc(a,b) convert([args],`+`) end:
`*` := proc(a,b) convert([args],`*`) end:
`-` := proc(x,y) if nargs = 1 then -x else x-y fi end:
`^` := proc(x,y) x^y end:
`/` := proc(x,y) x/y end:
`=` := proc(x,y) evalb(x=y) end:
`<` := proc(x,y) evalb(x<y) end:
`>` := proc(x,y) evalb(x>y) end:

# Polynomial definitions
UnivariatePolynomial := 'readlib('UnivariatePolynomial',
    ``.gauss.`/UP/UP.m`)':
DenseUnivariatePolynomial := 'readlib('DenseUnivariatePolynomial',
    ``.gauss.`/UP/DUP.m`)':
MapleUnivariatePolynomial := 'readlib('MapleUnivariatePolynomial',
    ``.gauss.`/UP/MUP.m`)':
OrderedUnivariatePolynomial := 'readlib('OrderedUnivariatePolynomial',
    ``.gauss.`/UP/OUP.m`)':

UnivariatePowerSeries := 'readlib('UnivariatePowerSeries',
    ``.gauss.`/UPS/UPS.m`)':
LazyUnivariatePowerSeries := 'readlib('LazyUnivariatePowerSeries',
    ``.gauss.`/UPS/LUPS.m`)':
ODESolve := 'readlib('ODESolve',
    ``.gauss.`/UPS/PSODES.m`)':


# Rational function and Quotient field definitions
QuotientField := 'readlib('QuotientField',``.gauss.`/QF/ENF.m`)':
ExpandedNormalForm := 'readlib('QuotientField',``.gauss.`/QF/ENF.m`)':
FactoredNormalForm := 'readlib('FNQuotientField',``.gauss.`/QF/FNF.m`)':
RationalFunction := 'readlib('RationalFunction',``.gauss.`/QF/RF.m`)':

# Matrix definitions
Matrix := 'readlib('Matrix', ``.gauss.`/MX/Matrix.m`)':
SquareMatrix := 'readlib('SquareMatrix', ``.gauss.`/SM/SM.m`)':

# Basic number domains
Integer := 'readlib('Integer',``.gauss.`/Z.m`)':
Rational := 'readlib('Rational',``.gauss.`/Q.m`)':
Zmod := 'readlib('Zmod', ``.gauss.`/FF/Zmod.m`)':
Gaussian := 'readlib('Gaussian',``.gauss.`/G.m`)':
GaloisField := 'readlib('GaloisField',``.gauss.`/FF/GF.m`)':

AlgebraicExtension := 'readlib('AlgebraicExtension', ``.gauss.`/SAE/AE.m`)':

Maple := 'readlib('Maple',``.gauss.`/Maple.m`)':

# Multivariate polynomial and related definitions
DenseExponentVector := 'readlib('DenseExponentVector',
    ``.gauss.`/EV/Dense.m`)':
PrimeExponentVector := 'readlib('PrimeExponentVector',
    ``.gauss.`/EV/Prime.m`)':
MacaulayExponentVector := 'readlib('MacaulayExponentVector',
    ``.gauss.`/EV/Macaulay.m`)':
MapleExponentVector := 'readlib('MapleExponentVector',
    ``.gauss.`/EV/Maple.m`)':
ExponentVector := 'readlib('ExponentVector',
    ``.gauss.`/EV/EV.m`)':
DistributedMultivariatePolynomial :=
    'readlib('DistributedMultivariatePolynomial', ``.gauss.`/DMP/DMP.m`)':
SparseDistributedMultivariatePolynomial :=
    'readlib('SparseDistributedMultivariatePolynomial',
    ``.gauss.`/DMP/SDMP.m`)':
TableDistributedMultivariatePolynomial :=
    'readlib('TableDistributedMultivariatePolynomial',
    ``.gauss.`/DMP/TDMP.m`)':
MapleMultivariatePolynomial := 
    'readlib('MapleMultivariatePolynomial',
    ``.gauss.`/DMP/MMP.m`)':
SparseMultivariatePolynomial := proc(R,X)
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    SparseDistributedMultivariatePolynomial(R,DenseExponentVector(X)):
end:

# Groebner Bases computation
GrobnerBasis := 'readlib('GB', ``.gauss.`/GB/GB.m`)':


#
#--> show(D); show(D,operations); show(D,categories); show(D,properties);
# Routine to print out operations in a domain D
# Author MBM: 1989
#

macro(format=`gauss/show/format`):
macro(format1=`gauss/show/format1`):
macro(formatSignature=`gauss/show/formatSignature`):
macro(formatSequence=`gauss/show/formatSequence`):
show := proc(D) local c,f,n,o,x;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    if not type(D,Domain) then ERROR(`domain expected`) fi;
    if type(D,name) then n := D
    elif assigned(D[DomainName]) then n := D[DomainName]
    else n := `$`
    fi;
    if nargs = 1 then
	c := D[Categories];
	o := map(op,[indices(D[Signatures])]);
	o := subs( {`0`=0,`1`=1}, sort( subs({0=`0`,1=`1`},o), lexorder ) );
	lprint();
	lprint(`Categories:`,D[Categories]);
	if D[Properties] <> [] then lprint(`Properties:`,D[Properties]) fi;
	lprint(`Operations:`,o);
	lprint()
    elif args[2] = 'categories' then D[Categories]
    elif args[2] = 'properties' then D[Properties]
    elif args[2] = 'operations' then
	f := map(op,[indices(D[Signatures])]);
	f := subs( {`0`=0,`1`=1}, sort( subs({0=`0`,1=`1`},f), lexorder ) );
	f := map(format,f,D,n);
	c := cat(`     Signatures for constructor `,n);
	n := `     note: operations prefixed by  --  are not available`;
	lprint(); lprint(c); lprint(n); lprint();
	for x in f do lprint(x) od;
	lprint();
    else ERROR(`dont know how to show`,args[2])
    fi
end:

format := proc(operation,D,n) local signature;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    signature := eval(D[Signatures][operation],1);
    if n <> `$` then
        signature := subs( `$` = n, Union(n,n) = n, signature );
    fi;
    format1(signature,D,operation)
end:

format1 := proc(signature,D,operation) local prefix;
    option `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    if type(signature,function) or type(signature,string) then
        if eval(D[operation],1) = notImplemented
            then prefix := '`  --  `'
            else prefix := '`      `'
        fi;
        cat( prefix, operation, '` : `', formatSignature(signature) )
    elif type(signature,set) then
        op(map(format1,signature,D,operation))
    else ERROR(`bad signature`)
    fi
end:

formatSignature := proc(x) local first, rest; option remember,system,
    `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    if type(x,function) then
        if op(0,x) = `&->` then
            cat(    formatSignature(op(1,x)), '` -> `',
                formatSignature(op(2,x)) )
        elif op(0,x) = List then cat('`[`',formatSequence(x),'`]`')
        elif op(0,x) = Nary then
            if nops(x) = 1
            then cat(formatSignature(op(1,x)),'`*`')
            else cat(formatSignature([op(x)]),'`*`')
            fi
        else cat(op(0,x),formatSignature([op(x)]))
        fi
    elif type(x,list) then cat('`(`',formatSequence(x),'`)`')
    elif type(x,integer) then convert(x,name)
    elif type(x,name) then x
    elif type(x,range) then
        first := formatSignature(op(1,x));
        rest := formatSignature(op(2,x));
        cat(first,` .. `,rest)
    else ERROR(`bad signature`)
    fi
end:

formatSequence := proc(x) local first, rest, comma; options remember,system,
    `Copyright 1992 Wissenschaftliches Rechnen, ETH Zurich`;
    if nops(x) = 0 then RETURN() fi;
    comma := '`,`';
    first := formatSignature(op(1,x));
    rest := map(formatSignature,subsop(1=NULL,x));
    first, op(map(proc(x,c) c,x end,rest,comma))
end:

`help/text/Gauss` := TEXT(
`HELP FOR: Gauss version 1.0`,
`   `,
`This is the share library version of Gauss.  You have loaded it into Maple.`,
`To use the share library version of Gauss, you must now set the gauss variable`,
`to the directory where the subroutines for Gauss library are kept so that Maple`,
`can find all the routines.  You should find a directory called gauss in`,
`the share library.  You need to do something like`,
`   `,
`> gauss := ``/usr/monagan/maple/share/gauss``;`,
`   `,
`Then you should execute Gauss();  This assigns the variables Z and Q to the`,
`domain of integers and rationals respectively, by executing`,
`   `,
`> Z := Integer():`,
`> Q := Rational():`,
`   `,
`and defines a number of other standard abbreviations for the domains known`,
`to Gauss -- see ?Gauss,domains.  Having defined the integers Z, you can`,
`compute the remainder of two integers m,n by doing`,
`   `,
`> Z[Rem](m,n);`,
`   `,
`Gauss is intended to be used as a tool for developing code for complicated`,
`algorithms.  The computational facilities offered by Gauss differ from Maple`,
`in that the user can parameterize a domain, e.g. a polynomial or matrix ring,`,
`by the coefficient ring.  This means that the code is more general as it`,
`will work in principle for any coefficient ring.  It is expected that some`,
`Maple library codes will be replaced by Gauss codes over the next versions`,
`of Maple because of this greater generality that Gauss offers.`,
`   `,
`When computing in Gauss, one first constructs a so called "domain of`,
`computation".  For example, if we want to compute with univariate polynomials`,
`with rational coefficients, one first creates the Q[x] domain as follows`,
`   `,
`> P := DenseUnivariatePolynomial(Q,x):`,
`   `,
`The value assigned to P is a Maple table of the operations that are`,
`available for computing with polynomials in Q[x].  For example, once could`,
`then compute the remainder of two polynomials a and b by doing`,
`   `,
`> P[Rem](a,b);`,
`   `,
`We suggest the reader read the help file ?Gauss,example to see a Gauss`,
`sample session before proceeding.`,
`   `,
`The principle advantage Gauss offers over Maple is in allowing one to`,
`parameterize domains, e.g. polynomial, series, matrix rings by a coefficient`,
`ring R.  The idea is, if you know how to compute in a coefficient ring R,`,
`and you write code for polynomials, then that code should work for any ring R.`,
`You should not have to rewrite the polynomial code for a new R.  Consider`,
`another example.  Suppose you have implement Gaussian elimination to solve`,
`a linear system of equations over a field F.  What Gauss allows you to do`,
`is to parameterize your subroutine by F so that your routine works for all`,
`fields, not just the field of rational numbers that you may have had in mind`,
`when you first wanted to write the routine.`,
`   `,
`This is accomplished by grouping together all operations in a ring or a field`,
`into a data structure called a domain.  In Gauss, a domain is just a Maple`,
`table of operations.  See ?Gauss,domain for some more details on domains`,
`and a list of known domains.`,
`   `,
`Coding in Gauss is quite straightforward though a little cumbersome due to`,
`having to package call every operation as you have seen in the examples.`,
`See ?Gauss,coding for an example of how to code a simple function which`,
`computes with values of a domain.`,
`Coding domains is somewhat more difficult and also cumbersome because nested`,
`lexical scopes and closures are not supported in Maple and must be simulated.`,
`We suggest that the reader study the code of existing domains.`,
`   `,
`Some of the main ideas behind Gauss come from the AXIOM system, which was`,
`formerly called Scratchpad II.  AXIOM also has the notion of parameterized`,
`types which it calls domains.  The author of Gauss would like to acknowledge`,
`the use of the primary idea behind AXIOM, namely that of of passing as a`,
`parameter a collection of functions as a single unit which the authors`,
`of AXIOM have termed a "domain".`,
`   `,
`For further information, please contact Michael Monagan: monagan@inf.ethz.ch`,
`Department of Computer Science, ETH Zentrum, CH 8092 Zurich, Switzerland.`,
`   `,
`SEE ALSO: Gauss[domain], Gauss[example], Gauss[coding]`
):
`help/text/gauss` := ":
`help/Gauss/text/coding` := TEXT(
`HELP FOR: Gauss[coding] - writing functions in Gauss`,
`   `,
`The basic idea for writing code in Gauss for computing with elements`,
`of a domain(s) is to pass the domain(s) as an argument(s) to the procedure.`,
`Essentially passing a collection of routines for manipulating elements`,
`of the domain.  E.g., let us write a routine to evaluate a univariate`,
`polynomial a(x) at x=b.  Our routine would look like this`,
`   `,
`    Evaluate := proc(P,a,b) local R,k,d,r;`,
`        if not hasCategory(P,UnivariatePolynomial) then ERROR(``...``) fi;`,
`        R := P[CoefficientRing];`,
`   `,
`We pass the domain P as the first argument and check that it is a univariate`,
`polynomial domain then since we need to do coefficient operations, we get`,
`the coefficient ring and call it R.`,
`   `,
`Next, we check the argument types of a and x as follows`,
`   `,
`        if not P[Type](a) then ERROR(``2nd argument must be of type P``) fi;`,
`        if not R[Type](b) then ERROR(``3rd argument must be of type R``) fi;`,
`   `,
`Now we can do the polynomial evaluation using Horners rule in the normal way.`,
`We need to use the Degree and Coeff functions from the univariate polynomial`,
`domain P, and the arithmetic operations ``+`` and ``*`` from the coefficient`,
`domain R.   `,
`   `,
`        d := P[Degree](a);`,
`        r := P[Coeff](a,d);`,
`        for k from d-1 by -1 to 0 do r := R[``+``](R[``*``](r,b),P[Coeff](a,k)) od;`,
`        r   `,
`    end:   `,
`   `,
`Note: the overhead of the table referencing in this example is quite small.`,
`I.e. the time to access the procedures P[Degree] and R[``+``] etc.  One might`,
`think to optimize this by factoring the table referencing operations out of`,
`the inner loop as follows.  Define three local variables cof, add, mul, and`,
`   `,
`	cof := eval(P[Coeff]);`,
`	mul := eval(R[``*``]);`,
`	add := eval(R[``+``]);`,
`   `,
`Note: P[Coeff], R[``*``] and R[``+``] are assigned to Maple procedures so eval`,
`must be used to evaluate to the procedure.  Then recoding the main loop as`,
`   `,
`        for k from d-1 by -1 to 0 do r := add(mul(r,b),cof(a,k)) od;`,
`   `,
`But this saves very little time.  The main overhead in Gauss comes not from`,
`this table subscripting but from the fact that almost every function in Gauss`,
`does a Maple procedure call.  I.e. the operations P[Degree], P[Coeff], R[``+``],`,
`R[``*``] are in general Maple procedure calls, which execute slower than`,
`builtin Maple functions.  E.g. P[Degree] and P[Coeff] will execute slower`,
`on the polynomial data structure than the builtin Maple functions degree`,
`and coeff do on the builtin Maple sum-of-products data structure.`
):
`help/gauss/text/coding` := ":
`help/gauss/text/domain` := TEXT(
`HELP FOR: Gauss[domain] - domains (parameterized types)`,
`   `,
`Domains in Gauss are Maple functions which return Maple tables of`,
`operations for manipulating objects in the domain.  E.g. Integer() returns`,
`a table of operations for computing with integers including`,
```+`` addition, ``-`` subtraction, ``*`` multiplication etc.`,
`   `,
`Domains may be parameterized by other domains and values, e.g. the`,
`domain DenseUnivariatePolynomial(R,x) takes a coefficient ring R and a`,
`variable x as a parameter.  The coefficient ring must be a Gauss domain`,
`which belongs to the category Ring, i.e. supports all the operations`,
`of a ring.  The variable x must be a Maple name.`,
`   `,
`All domains support belong to the category Set which supports the operations`,
`   `,
`1: =, <> -- boolean equality of domains elements`,
`2: Input -- for converting Maple expressions into the domain data representation`,
`3: Output -- for converting from the domain representation to an output form`,
`4: Random -- for generating a pseudo-random value from the domain`,
`5: Type -- for testing if a value is a valid domain element`,
`   `,
`The command show(D,operations) can be used to print out all the operations`,
`that are defined for a domain.  Operations marked by -- are not implemented.`,
`A list of the domains constructors in Gauss is`,
`   `,
`Z	Integer()						`,
`Q	Rational()					`,
`G	Gaussian(R:Ring)`,
`Zmod	Zmod(n:posint)`,
`GF	GaloisField(p:primeint,k:posint)`,
`   `,
`DUP	DenseUnivariatePolynomial(R:Ring,x:name)`,
`OUP	OrderedUnivariatePolynomial(P:UnivariatePolynomial(R),`,
`				    f:(R,R) -> Boolean)`,
`   `,
`DEV	DenseExponentVector(X:list(name))`,
`PEV	PrimeExponentVector(X:list(name))`,
`MEV	MapleExponentVector(X:list(name))`,
`TEV	MacaulayExponentVector(X:list(name))`,
`TDMP	TableDistributedMultivariatePolynomial(R:Ring,E:ExponentVector)`,
`SDMP	SparseDistributedMultivariatePolynomial(R:Ring,E:ExponentVector)`,
`   `,
`QF	QuotientField(D:GcdDomain)`,
`ENF	ExpandedNormalForm(D:GcdDomain)`,
`FNF	FactoredNormalForm(D:GcdDomain)`,
`RF	RationalFunction(D:GcdDomain,X:list(name))`,
`LUPS	LazyUnivariatePowerSeries(R:Ring,x:name)`,
`	Matrix(R:Ring)`,
`SM	SquareMatrix(n:posint,R:Ring)`,
`SAE	AlgebraicExtension(D:UnvivaraitePolynomial,m:D)	`,
`   `,
`In addition, there are some special domains that use the Maple representation`,
`for polynomials to try to get back some efficiency for integer and rational`,
`coefficients.`,
`   `,
`MUP	MapleUnivariatePolynomial(R:{Z,Q,Zmod}, x:name)`,
`MMP	MapleMultivariatePolynomial(R:{Z,Q,Zmod}, X:list(name))`,
`   `,
`And there are some packages of routines for computing with these`,
`objects, but they are very much in an experimental status.`
):
`help/Gauss/text/domain` := ":
`help/gauss/text/example` := TEXT(
`HELP FOR: Gauss[example]`,
`   `,
`This is Maple sample session for showing how to use Gauss.`,
`Note that all Gauss functions begin with an upper case letter.`,
`   `,
`    |\\^/|     Maple V Release 2`,
`._|\\|   |/|_. Copyright (c) 1981-1991 by the University of Waterloo.`,
` \\  MAPLE  /  All rights reserved. MAPLE is a registered trademark of`,
` <____ ____>  Waterloo Maple Software.`,
`      |       Type ? for help.`,
`> with(Gauss);`,
`----------------------- Gauss version 1.0 ----------------------`,
`Initially defined domains are Z and Q the integers and rationals`,
`   `,
`                                     [init]`,
`   `,
`# The domains Z (the integers) and Q (the rationals) have been defined.`,
`# Lets do some operations`,
`> Z[Gcd](8,12);`,
`                                       4`,
`   `,
`> Q[``+``](1/2,1/3,1/4);`,
`                                      13`,
`                                     ----`,
`                                      12`,
`   `,
`# What is the object Z and Q?  Z and Q are Maple tables`,
`> type(Z,table);`,
`                                     true`,
`   `,
`# The contains operations (Maple procedures) for computing in Z and Q`,
`# What operations are available?`,
`> show(Z,operations);`,
`   `,
`     Signatures for constructor Z`,
`     note: operations prefixed by  --  are not available`,
`   `,
`      * : (Z,Z*) -> Z`,
`      * : (Integer,Z) -> Z`,
`      + : (Z,Z*) -> Z`,
`      - : Z -> Z`,
`      - : (Z,Z) -> Z`,
`      0 : Z`,
`      1 : Z`,
`      < : (Z,Z) -> Boolean`,
`      <= : (Z,Z) -> Boolean`,
`      <> : (Z,Z) -> Boolean`,
`      = : (Z,Z) -> Boolean`,
`      > : (Z,Z) -> Boolean`,
`      >= : (Z,Z) -> Boolean`,
`      Abs : Z -> Z`,
`      Characteristic : Integer`,
`      Coerce : Integer -> Z`,
`      Div : (Z,Z) -> Union(Z,FAIL)`,
`      EuclideanNorm : Z -> Integer`,
`      Factor : Z -> [Z,[Z,Z]*]`,
`      Gcd : Z* -> Z`,
`      Gcdex : (Z,Z,Name) -> Z`,
`      Gcdex : (Z,Z,Name,Name) -> Z`,
`      Input : Expression -> Union(Z,FAIL)`,
`      Inv : Z -> Union(Z,FAIL)`,
`      Lcm : Z* -> Z`,
`      Max : (Z,Z*) -> Z`,
`      Min : (Z,Z*) -> Z`,
`      Modp : (Z,Z) -> Z`,
`      Mods : (Z,Z) -> Z`,
`      ModularHomomorphism : () -> (Z -> Z,Z)`,
`      Normal : Z -> Z`,
`      Output : Z -> Expression`,
`      Powmod : (Z,Integer,Z) -> Z`,
`      Prime : Z -> Boolean`,
`      Quo : (Z,Z) -> Z`,
`      Quo : (Z,Z,Name) -> Z`,
`      Random : () -> Z`,
`      RelativelyPrime : (Z,Z) -> Boolean`,
`      Rem : (Z,Z) -> Z`,
`      Rem : (Z,Z,Name) -> Z`,
`      Sign : Z -> UNION(1,-1,0)`,
`      SmallerEuclideanNorm : (Z,Z) -> Boolean`,
`      Sqrfree : Z -> [Z,[Z,Z]*]`,
`      Type : Expression -> Boolean`,
`      Unit : Z -> Z`,
`      UnitNormal : (Z) -> (Z,Z,Z)`,
`      Zero : Z -> Boolean`,
`      ^ : (Z,Integer) -> Z`,
`   `,
`# Lets do some operations in Q[x], i.e. univariate polynomials over Q`,
`# First we have to create the domain Q[x], lets call it D`,
`> D := DenseUnivariatePolynomial(Q,x):`,
`# The name DenseUnivaritePolynomial indicates that the data structure being`,
`# used is a dense one.  Lets input a polyomial`,
`> m := D[Input](x^4-10*x^2+1);`,
`   `,
`                            m := [1, 0, -10, 0, 1]`,
`   `,
`# Notice that the polynomial is displayed as a Maple list of coefficients`,
`> D[Output](m);`,
`                                 4       2`,
`                                x  - 10 x  + 1`,
`   `,
`# What can we do with the polynomial m?`,
`> D[Degree](m); # compute it's degree`,
`                                       4`,
`   `,
`> D[``^``](m,2); # square it`,
`   `,
`                       [1, 0, -20, 0, 102, 0, -20, 0, 1]`,
`   `,
`# What other operations are there for univariate polynomials?`,
`# A list of the operations is given by  show(D,operations);`,
`> a := D[Random]();`,
`                                    35    50    55`,
`                           a := [- ----, ----, ----]`,
`                                    97    79    37`,
`   `,
`# Every Gauss domain has a Random function which returns a randomly generated`,
`# element from a domain, which is useful for writing documentation examples!`,
`> D[Output](a);`,
`                             55   2    50       35`,
`                            ---- x  + ---- x - ----`,
`                             37        79       97`,
`   `,
`# Gauss can also compute with matrices and other objects`,
`# Lets compute the inverse of a 3 by 3 Hilbert matrix using Gauss`,
`# First, we must create the matrix domain`,
`> M := SquareMatrix(3,Q):`,
`> A := M[Input]([[1,1/2,1/3],[1/2,1/3,1/4],[1/3,1/4,1/5]]);`,
`   `,
`            A := [[1, 1/2, 1/3], [1/2, 1/3, 1/4], [1/3, 1/4, 1/5]]`,
`   `,
`> M[Det](A);`,
`                                    1/2160`,
`   `,
`> M[Inv](A);`,
`               [[9, -36, 30], [-36, 192, -180], [30, -180, 180]]`,
`   `,
`# Lets now compute with matrices of polynomials in Q[x].`,
`# Lets use 2 by 2 matrices to keep the examples small.`,
`> M := SquareMatrix(2,D):`,
`> A := M[Input]([[x,x^2],[1-x^2,1+x^2]]);`,
`   `,
`              A := [[[0, 1], [0, 0, 1]], [[1, 0, -1], [1, 0, 1]]]`,
`   `,
`> D[Output]( M[Det](A) );  # Compute and output det(A)`,
`   `,
`                                4    3    2`,
`                               x  + x  - x  + x`,
`   `,
`> M[Output]( M[Inv](A) );  # Compute and output A^(-1)`,
`Error, (in notImplemented) operation is not implemented`,
`# That's right, you can't compute the inverse of a matrix of polynomials,`,
`# because the result is in general a rational function.`,
`# Now let us demonstrate the power of Gauss by doing the same Matrix`,
`# problems with a Matrix of different entries, this time algebraic numbers.`,
`# For example, suppose we want to compute with the roots of the polyomial`,
`# m = x^4-10*x^2+1.  In Maple one uses the RootOf function as follows`,
`> alias(alpha=RootOf(x^4-10*x^2+1=0,x)):`,
`> simplify(1/alpha);`,
`                                     3`,
`                              - alpha  + 10 alpha`,
`   `,
`# In Gauss, we will create a simple algebraic extension using the polynomial`,
`# m as the minimal polynomial thus`,
`> F := SAE(D,m):`,
`> a := F[Input](x);`,
`                                  a := [0, 1]`,
`   `,
`> ai := F[Inv](a);`,
`                             ai := [0, 10, 0, -1]`,
`   `,
`> F[Output](ai);`,
`                                     3`,
`                                  - x  + 10 x`,
`   `,
`# Gauss generalizes very naturally to compute over F`,
`> M := SquareMatrix(2,F):`,
`> A := M[Input]([[1,x],[x^3,x^2]]);`,
`   `,
`                A := [[[1], [0, 1]], [[0, 0, 0, 1], [0, 0, 1]]]`,
`   `,
`> F[Output]( M[Det](A) );`,
`                                       2`,
`                                  - 9 x  + 1`,
`   `,
`> M[Output]( M[Inv](A) );`,
`   `,
`             2             3                  3                 2`,
`    [[- 1/8 x  + 9/8, 9/8 x  - 89/8 x], [1/8 x  - 9/8 x, - 9/8 x  + 89/8]]`,
`   `,
`# For our last examples, some calculations with univariate power series`,
`# First lets create a univariate power series domain in x over Q`,
`> PS := LazyUnivariatePowerSeries(Q,x):`,
`> show(PS,operations);`,
`   `,
`     Signatures for constructor PS`,
`     note: operations prefixed by  --  are not available`,
`   `,
`      * : (PS,PS*) -> PS`,
`      * : (Integer,PS) -> PS`,
`      + : (PS,PS*) -> PS`,
`      - : (PS,PS) -> PS`,
`      - : PS -> PS`,
`      / : (PS,Integer) -> PS`,
`      / : (PS,PS) -> PS`,
`      0 : PS`,
`      1 : PS`,
`      <> : (PS,PS) -> Boolean`,
`      = : (PS,PS) -> Boolean`,
`  --  AbsoluteDegree : Integer`,
`      Characteristic : Integer`,
`      Coeff : (PS,Integer) -> Q`,
`      CoefficientRing : Ring`,
`      Coerce : Integer -> PS`,
`      Constant : Q -> PS`,
`      Cos : PS -> Union(PS,FAIL)`,
`      Cosh : PS -> Union(PS,FAIL)`,
`      Diff : PS -> PS`,
`      Div : (PS,PS) -> Union(PS,FAIL)`,
`      EuclideanNorm : PS -> Integer`,
`      Exp : PS -> Union(PS,FAIL)`,
`      Factor : PS -> [PS,[PS,PS]*]`,
`      Gcd : PS* -> PS`,
`      Gcdex : (PS,PS,Name) -> PS`,
`      Gcdex : (PS,PS,Name,Name) -> PS`,
`      Input : Expression -> Union(PS,FAIL)`,
`      Integrate : PS -> Union(PS,FAIL)`,
`      Inv : PS -> PS`,
`      Lcm : PS* -> PS`,
`      Ln : PS -> Union(PS,FAIL)`,
`      Log : PS -> Union(PS,FAIL)`,
`      Lorder : PS -> Integer`,
`      Monomial : () -> PS`,
`      Monomial : Integer -> PS`,
`      Normal : PS -> PS`,
`      Output : PS -> Expression`,
`      Powmod : (PS,Integer,PS) -> PS`,
`      Prime : PS -> Boolean`,
`      Quo : (PS,PS,Name) -> PS`,
`      Quo : (PS,PS) -> PS`,
`      R* : (Q,PS) -> PS`,
`      R/ : (PS,Q) -> PS`,
`      Random : () -> PS`,
`      RelativelyPrime : (PS,PS) -> Boolean`,
`      Rem : (PS,PS) -> PS`,
`      Rem : (PS,PS,Name) -> PS`,
`      Series : [Q,Integer] -> PS`,
`      Shift : (PS,Integer) -> PS`,
`      Sin : PS -> Union(PS,FAIL)`,
`      Sinh : PS -> Union(PS,FAIL)`,
`      SmallerEuclideanNorm : (PS,PS) -> Boolean`,
`      Sqrfree : PS -> [PS,[PS,PS]*]`,
`      Tan : PS -> Union(PS,FAIL)`,
`      Tanh : PS -> Union(PS,FAIL)`,
`      Type : PS -> Boolean`,
`      Unit : PS -> PS`,
`      UnitNormal : (PS) -> (PS,PS,PS)`,
`      Variable : Name`,
`      ^ : (PS,Integer) -> PS`,
`      ^ : (PS,Rational) -> PS`,
`      ^ : (PS,PS) -> PS`,
`      order : PS -> Integer`,
`   `,
`# Lets compute the series for exp(x)`,
`> a := PS[Input](x);`,
`                                    a := x`,
`   `,
`> e := PS[Exp](a);`,
`   `,
`                             2        3         4          5      6`,
`           e := 1 + x + 1/2 x  + 1/6 x  + 1/24 x  + 1/120 x  + O(x )`,
`   `,
`> PS[``^``](e,3); # exp(x)^3`,
`   `,
`                            2        3         4    81   5      6`,
`             1 + 3 x + 9/2 x  + 9/2 x  + 27/8 x  + ---- x  + O(x )`,
`                                                    40`,
`   `,
`# Lazy univariate power series are "lazy" which means coefficients are`,
`# computed on demand, i.e. we can compute a series to higher order without`,
`# haveing to recompute any previously computed coefficients.`,
`> PS[Output](e,10);`,
`   `,
`                   2        3         4          5          6           7`,
`      1 + x + 1/2 x  + 1/6 x  + 1/24 x  + 1/120 x  + 1/720 x  + 1/5040 x`,
`   `,
`                      8             9              10      11`,
`           + 1/40320 x  + 1/362880 x  + 1/3628800 x   + O(x  )`,
`   `,
`# Here are the Legendre polynomials computed from their generating function:`,
`#            1/sqrt(1-2*x*t+t^2) = sum( L[n](x)*t^k, k=0.. )`,
`# First create the domain: LUPS == LazyUnivariatePowerSeries`,
`> P := LUPS(LUPS(Q,x),t):`,
`> p := P[Input](1-2*x*t+t^2);`,
`                                                2`,
`                              p := 1 - 2 x t + t`,
`   `,
`> P[``^``](p, -1/2);`,
`   `,
`                         2   2                   3   3`,
` 1 + x t + (- 1/2 + 3/2 x ) t  + (- 3/2 x + 5/2 x ) t`,
`   `,
`                     2         4   4                   3         5   5      6`,
`      + (3/8 - 15/4 x  + 35/8 x ) t  + (15/8 x - 35/4 x  + 63/8 x ) t  + O(t )`,
`   `,
`# Compare the output from Gauss with that from Maple`,
`> seq(orthopoly[P](i,x), i=0..5);`,
`   `,
`                          2                 3              2         4`,
`       1, x, - 1/2 + 3/2 x , - 3/2 x + 5/2 x , 3/8 - 15/4 x  + 35/8 x ,`,
`   `,
`                          3         5`,
`           15/8 x - 35/4 x  + 63/8 x`,
`   `,
`> quit   `
):
`help/Gauss/text/example` := ":

save `Gauss.m`;
quit
