Math Solutions 1995 October/Math_Solutions_CD-ROM_Walnut_Creek_October_1995.iso/pc/mac/discrete/lib/polynom.g

Text File | 1993-05-05 | 40.3 KB | 1,483 lines

#############################################################################
##
#A  polynom.g                   GAP library                      Frank Celler
##
#A  @(#)$Id: polynom.g,v 3.23 1993/02/11 16:47:56 fceller Rel $
##
#Y  Copyright 1990-1992,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
##
##  This file  contains those functions that  are  dispatchers for polynomial
##  functions  and  the  polynomial   arithmetic.    The  idea  to  represent
##  polynomials as records containing '.baseRing'  and '.coefficients' is due
##  to S.Linton.
##
#H  $Log: polynom.g,v $
#H  Revision 3.23  1993/02/11  16:47:56  fceller
#H  added <list> + <polynomial>
#H
#H  Revision 3.22  1993/02/08  11:56:05  fceller
#H  ground ring is now stored in polynomial record
#H
#H  Revision 3.21  1993/02/04  13:47:12  fceller
#H  added '/' and 'EuclideanQuotient'
#H
#H  Revision 3.20  1993/02/03  12:36:06  fceller
#H  fixed undeclared local variables
#H
#H  Revision 3.19  1993/01/25  13:18:26  fceller
#H  fixed '=' and '<' to handle iterated polynomials
#H
#H  Revision 3.18  1993/01/20  11:09:02  fceller
#H  changed 'DisplayPolynomial' again (in order to display
#H  polynomials over the cyclotmics)
#H
#H  Revision 3.17  1993/01/06  15:44:35  fceller
#H  changed 'LeadingCoefficient' to handle the trivial polynomial
#H
#H  Revision 3.16  1993/01/04  11:18:50  fceller
#H  added 'LeadingCoefficient'
#H
#H  Revision 3.15  1992/12/16  19:47:27  martin
#H  replaced quoted record names with escaped ones
#H
#H  Revision 3.14  1992/12/07  07:41:08  fceller
#H  fixed 'Quotient( 0, f )'
#H
#H  Revision 3.13  92/11/25  12:34:28  fceller
#H  added 'Degree'
#H
#H  Revision 3.12  1992/11/19  14:12:23  fceller
#H  changed 'DisplayPolynomial'
#H
#H  Revision 3.11  92/11/16  12:23:24  fceller
#H  added Laurent polynomials
#H
#H  Revision 3.10  1992/08/13  13:11:09  fceller
#H  added 'Indeterminate'
#H
#H  Revision 3.9  1992/07/24  07:08:07  fceller
#H  improved '*' to allow <nullpoly> * <int>
#H
#H  Revision 3.8  1992/07/23  09:04:18  fceller
#H  improved '*' to allow <int> * <polynomial>
#H
#H  Revision 3.7  1992/06/17  07:06:04  fceller
#H  moved '<somedomain>.operations.Polynomial' function to "<somedomain>.g"
#H
#H  Revision 3.6  1992/06/05  09:58:04  fceller
#H  improved 'Derivative'
#H
#H  Revision 3.5  1992/06/01  07:33:55  fceller
#H  added 'Value'
#H
#H  Revision 3.4  1992/05/25  09:10:05  fceller
#H  added 'EmbeddedPolynomial', fixed a bug
#H
#H  Revision 3.3  1992/05/07  10:29:16  fceller
#H  added 'CompanionMatrix', better pretty printing in 'DisplayPolynomial'
#H
#H  Revision 3.2  1992/04/21  10:37:47  fceller
#H  Initial GAP 3.2 release
##


#############################################################################
##
#F  InfoPoly1(...)  . . . . . . . . . . . infomation function for polynomials
#F  InfoPoly2(...)  . . . . . . . . . . . . .  debug function for polynomials
##
if not IsBound(InfoPoly1)    then InfoPoly1   := Ignore;  fi;
if not IsBound(InfoPoly2)    then InfoPoly2   := Ignore;  fi;


#############################################################################
##

#F  Polynomial( <R>, <coeffs>, <val> )  . . . . . . . . . polynomial over <R>
##
Polynomial := function( arg )
    local   R,  coeffs,  val;

    # check the argument list
    if 2 = Length(arg)  then
        R := arg[1];
        coeffs := arg[2];
        val := 0;
    elif 3 = Length(arg)  then
        R := arg[1];
        coeffs := arg[2];
        val := arg[3];
    else
        Error( "usage: Polynomial( <R>, <coeffs>, <val> )",
               "or Polynomial( <R>, <coeffs> )" );
    fi;

    # <R> must be a ring with one or a field
    if not IsField(R) and not IsRing(R)  then
        Error("<R> must be a field or ring");
    elif IsRing(R) and not IsBound(R.one)  then
        Error("<R> must be a ring with one");
    fi;

    # construct the polynomial
    return R.operations.Polynomial( R, coeffs, val );

end;


#############################################################################
##
#F  IsPolynomial( <obj> ) . . . . . . . . . . . . . . . is <obj> a polynomial
##
IsPolynomial := function( obj )
    return IsRec(obj) and IsBound(obj.isPolynomial) and obj.isPolynomial;
end;


#############################################################################
##
#F  CompanionMatrix( <R>, <f> ) . . . . . . . . . . . . . .  companion matrix
##
CompanionMatrix := function( R, f )

    # check <R> and <f>
    if not IsPolynomial( f )  then
        Error( "<f> must be a polynomial" );
    elif f.valuation < 0  then
        Error( "<f> must not be a Laurent polynomial" );
    fi;

    # check that the coefficients of <f> are in the base ring of <R>
    if not ForAll( f.coefficients, x -> x in R.baseRing )  then
        Error( "cannot write <f> over <R>.baseRing" );
    fi;

    # return the companion matrix of <f>
    return R.operations.CompanionMatrix( R, f );

end;


#############################################################################
##
#F  EmbeddedPolynomial( <R>, <f> )  . . . . . . . embed polynomial <f> in <R>
##
EmbeddedPolynomial := function( R, f )

    # check <R> and <f>
    if not IsPolynomial( f )  then
        Error( "<f> must be a polynomial" );
    fi;

    # check that the coefficients of <f> are in the base ring of <R>
    if not ForAll( f.coefficients, x -> x in R.baseRing )  then
        Error( "cannot write <f> over <R>.baseRing" );
    fi;

    # convert
    return R.operations.EmbeddedPolynomial( R, f );

end;


#############################################################################
##
#F  RandomPol( <R>, <deg> ) . . . . . . . . . . . . . . . . random polynomial
##
##  'RandomPol' returns a  random polynomial of  degree less than or equal to
##  <deg>,  which must be a nonnegativ integer,  whose coefficients come from
##  a ring <R>.
##
RandomPol := function ( R, deg )
    local  rand,  i;

    rand := [];
    for i  in [ 1 .. deg+1 ]  do
        rand[i] := Random( R );
    od;
    return Polynomial( R, rand );

end;


#############################################################################
##
#F  Value( <f>, <x> ) . . . . . . . . . . . . . . . . . . . . return <f>(<x>)
##
Value := function( f, x )
    return f.operations.Value( f, x );
end;


#############################################################################
##
#F  Indeterminate( <R> )  . . . . . . . . . . . . . . indeterminate of a ring
##
Indeterminate := function( R )
    if not IsBound( R.indeterminate )  then
        R.indeterminate := R.operations.Indeterminate( R );
    fi;
    return R.indeterminate;
end;

X := Indeterminate;


#############################################################################
##

#F  LaurentPolynomialRing( <R> )  . . . .  full Laurent polynom ring over <R>
##
LaurentPolynomialRing := function( R )
    local   i,  B;

    # <R> must be a ring with one
    if not IsRing(R) and not IsField(R)  then
        Error( "<R> must be a ring or field" );
    elif IsRing(R) and not IsBound(R.one)  then
        Error( "<R> must be a ring with one" );
    fi;

    # check if the polynom ring is already known
    if not IsBound(R.laurentPolynomialRing)  then
        R.laurentPolynomialRing := R.operations.LaurentPolynomialRing(R);
        if IsBound(R.laurentPolynomialRing.generators)  then
            B := R.laurentPolynomialRing.generators;
            for i  in [ 1 .. Length(B) ]  do
                R.laurentPolynomialRing.(i) := B[i];
            od;
        fi;
    fi;
    return R.laurentPolynomialRing;

end;


#############################################################################
##
#F  IsLaurentPolynomialRing( <obj> )  . . . . . is <R> a Laurent polynom ring
##
IsLaurentPolynomialRing := function( obj )
    return     IsRec(obj)
           and IsBound(obj.isLaurentPolynomialRing)
           and obj.isLaurentPolynomialRing;
end;


#############################################################################
##
#V  LaurentPolynomialRingOps  . . . . . . . . . full Laurent polynom ring ops
##
LaurentPolynomialRingOps := Copy( RingOps );
LaurentPolynomialRingOps.name := "LaurentPolynomialRingOps";


#############################################################################
##
#F  LaurentPolynomialRingOps.\in . . . . . . . . . . . . . . membership test
##
LaurentPolynomialRingOps.\in := function( p, P )
    return     IsRec( p )
           and IsBound( p.isPolynomial )
           and p.isPolynomial
           and p.baseRing = P.baseRing;
end;


#############################################################################
##
#F  LaurentPolynomialRingOps.Factors( <R>, <f> )  . . . . . .  factors of <f>
##
LaurentPolynomialRingOps.Factors := function( R, f )
    local   g,  r;

    # factorize the standard associate of <f>
    g := StandardAssociate( R, f );
    r := Factors( PolynomialRing(R.baseRing), g );
    if 0 < Length(r)  then
    r[1] := Quotient(R,f,g) * r[1];
    fi;
    return r;

end;


#############################################################################
##
#F  LaurentPolynomialRingOps.Quotient( <R>, <l>, <r> )  . . . . . . <l> / <r>
##
LaurentPolynomialRingOps.Quotient := function ( R, f, g )
    local   m,  n,  i,  k,  c,  q,  R,  val;

    # get the base ring
    R := R.baseRing;

    # if <f> is zero return it
    if 0 = Length(f.coefficients)  then
    return f;
    fi;

    # get the value of the valuation of <f> / <g>
    val := f.valuation - g.valuation;

    # Try to divide <f> by <g>
    q := [];
    n := Length( g.coefficients );
    m := Length( f.coefficients ) - n;
    f := ShallowCopy( f.coefficients );
    if IsField(R)  then
        for i  in [0 .. m ]  do
            c := f[m-i+n] / g.coefficients[n];
            for k  in [ 1 .. n ]  do
                f[m-i+k] := f[m-i+k] - c * g.coefficients[k];
            od;
            q[m-i+1] := c;
        od;
    else
        for i  in [0 .. m ]  do
            c := Quotient( R, f[m-i+n], g.coefficients[n] );
            if c = false  then return false;  fi;
            for k  in [ 1 .. n ]  do
                f[m-i+k] := f[m-i+k] - c * g.coefficients[k];
            od;
            q[m-i+1] := c;
        od;
    fi;

    # Did the division work?
    for i  in [ 1 .. m+n ]  do
        if f[i] <> R.zero  then
            return false;
        fi;
    od;
    return Polynomial( R, q, val );

end;


#############################################################################
##
#F  LaurentPolynomialRingOps.Derivative( <R>, <f> ) . . . . . . .  derivative
##
Derivative := function( arg )
    local   R,  f;

    if Length(arg) = 2  then
        R := arg[1];
        f := arg[2];
    else
        f := arg[1];
        R := DefaultRing(f);
    fi;
    return R.operations.Derivative( R, f );

end;

LaurentPolynomialRingOps.Derivative := function( R, f )
    local  d,  i;

    d := [];
    for i  in [ 1 .. Length(f.coefficients) ]  do
        d[i] := (i+f.valuation-1) * f.coefficients[i];
    od;
    return Polynomial( R.baseRing, d, f.valuation-1 );

end;


#############################################################################
##
#F  LaurentPolynomialRingOps.EuclideanDegree( <R>, <f> )  . . . degree of <f>
##
LaurentPolynomialRingOps.EuclideanDegree := function( R, f )
    return Length(f.coefficients)-1;
end;


#############################################################################
##
#F  LaurentPolynomialRingOps.StandardAssociate( <R>, <f> ) standard associate
##
LaurentPolynomialRingOps.StandardAssociate := function( R, f )
    local   l,  a;

    # <f> should be nontrivial
    if 0 < Length(f.coefficients)  then

        # get standard associate of leading term
        l := f.coefficients[Length(f.coefficients)];
        a := Quotient( R.baseRing, StandardAssociate( R.baseRing, l ), l );
        f := Polynomial(R.baseRing,List(f.coefficients,x->a*x),0);
        f.valuation := 0;
    fi;
    return f;

end;


#############################################################################
##
#F  LaurentPolynomialRingOps.EmbeddedPolynomial( <R>, <f> )  embed polynomial
##
LaurentPolynomialRingOps.EmbeddedPolynomial := function( R, f )
    return Polynomial( R.baseRing, f.coefficients, f.valuation );
end;


#############################################################################
##
#F  LaurentPolynomialRingOps.Print( <P> ) . . . . . . . . . . pretty printing
##
LaurentPolynomialRingOps.Print := function( P )
    if IsBound( P.name )  then
        Print( P.name );
    else
        Print( "LaurentPolynomialRing( ", P.baseRing, " )" );
    fi;
end;


#############################################################################
##

#F  PolynomialRing( <R> ) . . . . . . . . . . . .  full polynom ring over <R>
##
PolynomialRing := function( R )
    local   i;

    # <R> must be a ring with one
    if not IsRing(R) and not IsField(R)  then
        Error( "<R> must be a ring or field" );
    elif IsRing(R) and not IsBound(R.one)  then
        Error( "<R> must be a ring with one" );
    fi;

    # check if the polynom ring is already known
    if not IsBound(R.polynomialRing)  then
        R.polynomialRing := R.operations.PolynomialRing(R);
        if IsBound(R.polynomialRing.generators)  then
            for i  in [ 1 .. Length(R.polynomialRing.generators) ]  do
                R.polynomialRing.(i) := R.polynomialRing.generators[i];
            od;
        fi;
    fi;
    return R.polynomialRing;

end;


#############################################################################
##
#F  IsPolynomialRing( <obj> ) . . . . . . . . . . . . . is <R> a polynom ring
##
IsPolynomialRing := function( obj )
    return     IsRec(obj)
           and IsBound(obj.isPolynomialRing)
           and obj.isPolynomialRing;
end;


#############################################################################
##
#V  PolynomialRingOps . . . . . . . . . . . . . . . . . full polynom ring ops
##
PolynomialRingOps := Copy( RingOps );
PolynomialRingOps.name := "PolynomialRingOps";


#############################################################################
##
#F  PolynomialRingOps.\in  . . . . . . membership test for full polynom ring
##
PolynomialRingOps.\in := function( p, P )
    return     IsRec( p )
           and IsBound( p.isPolynomial )
           and p.isPolynomial
           and 0 <= p.valuation
           and p.baseRing = P.baseRing;
end;


#############################################################################
##
#F  PolynomialRingOps.Quotient( <R>, <f>, <g> ) . . . . . . . . . . <l> / <r>
##
PolynomialRingOps.Quotient := function ( R, f, g )
    local   m,  n,  i,  k,  c,  q,  R,  val;

    # <f> and <g> must have the same field
    if f.baseRing <> g.baseRing  then
        Error( "<f> and <g> must have the same base ring" );
    fi;
    R := f.baseRing;

    # if <f> is zero return it
    if 0 = Length(f.coefficients)  then
    return f;
    fi;

    # check the value of the valuation of <f> and <g>
    if f.valuation < g.valuation  then
        return false;
    fi;
    val := f.valuation - g.valuation;

    # Try to divide <f> by <g>
    q := [];
    n := Length( g.coefficients );
    m := Length( f.coefficients ) - n;
    f := ShallowCopy( f.coefficients );
    if IsField(R)  then
        for i  in [0 .. m ]  do
            c := f[m-i+n] / g.coefficients[n];
            for k  in [ 1 .. n ]  do
                f[m-i+k] := f[m-i+k] - c * g.coefficients[k];
            od;
            q[m-i+1] := c;
        od;
    else
        for i  in [0 .. m ]  do
            c := Quotient( R, f[m-i+n], g.coefficients[n] );
            if c = false  then return false;  fi;
            for k  in [ 1 .. n ]  do
                f[m-i+k] := f[m-i+k] - c * g.coefficients[k];
            od;
            q[m-i+1] := c;
        od;
    fi;

    # Did the division work?
    for i  in [ 1 .. m+n ]  do
        if f[i] <> R.zero  then
            return false;
        fi;
    od;
    return Polynomial( R, q, val );

end;


#############################################################################
##
#F  PolynomialRingOps.Derivative( <R>, <f> )  . .  derivative of a polynomial
##
PolynomialRingOps.Derivative := function( R, f )
    local  d,  i;

    if f.valuation+Length(f.coefficients) < 2  then
        return R.zero;
    else
        d := [];
        for i  in [ 1 .. Length(f.coefficients) ]  do
            d[i] := (i+f.valuation-1) * f.coefficients[i];
        od;
        return Polynomial( R.baseRing, d, f.valuation-1 );
    fi;

end;


#############################################################################
##
#F  PolynomialRingOps.EuclideanDegree( <R>, <f> ) . . . . . . . degree of <f>
##
PolynomialRingOps.EuclideanDegree := function( R, f )
    return ( Length(f.coefficients)-1 ) + f.valuation;
end;


#############################################################################
##
#F  PolynomialRingOps.CompanionMatrix( <R>, <f> ) . . companion matrix of <f>
##
##  Returns the companion matrix A of <f>.  If <f> has leading coefficient 1
##  the companion matrix has characteristic polynomial <f>.
##
PolynomialRingOps.CompanionMatrix := function( R, f )
    local   d,  C,  i;

    d := (Length(f.coefficients)-1) + f.valuation;
    C := NullMat( d, d, R.baseRing.zero );
    for i  in [ 1 .. d-1 ]  do
        C[i][i+1] := R.baseRing.one;
    od;
    for i  in [ f.valuation+1 .. d ]  do
        C[d][i] := -f.coefficients[i];
    od;
    return C;

end;


#############################################################################
##
#F  PolynomialRingOps.StandardAssociate( <R>, <f> ) standard associate of <f>
##
PolynomialRingOps.StandardAssociate := function( R, f )
    local   l,  a;

    # <f> should be nontrivial
    if 0 < Length(f.coefficients)  then

        # get standard associate of leading term
        l := f.coefficients[Length(f.coefficients)];
        a := Quotient( R.baseRing, StandardAssociate( R.baseRing, l ), l );
        f := Polynomial(R.baseRing,List(f.coefficients,x->a*x),f.valuation);
    fi;
    return f;

end;


#############################################################################
##
#F  PolynomialRingOps.EmbeddedPolynomial( <R>, <f> )  . . .  embed polynomial
##
PolynomialRingOps.EmbeddedPolynomial := function( R, f )
    return Polynomial( R.baseRing, f.coefficients, f.valuation );
end;


#############################################################################
##
#F  PolynomialRingOps.Print( <P> )  . . . . . . . . . . . . . pretty printing
##
PolynomialRingOps.Print := function( P )
    if IsBound( P.name )  then
        Print( P.name );
    else
        Print( "PolynomialRing( ", P.baseRing, " )" );
    fi;
end;


#############################################################################
##

#V  LaurentPolynomials  . . . . . . . . . . domain of all Laurent polynomials
##
LaurentPolynomials            := rec();

LaurentPolynomials.isDomain   := "true";

LaurentPolynomials.name       := "LaurentPolynomials";

LaurentPolynomials.isFinite   := false;
LaurentPolynomials.size       := "infinity";

LaurentPolynomials.operations := Copy( DomainOps );
LaurentPolynomialsOps         := LaurentPolynomials.operations;

LaurentPolynomialsOps.\in := function( p, Polynomials )
    return     IsRec( p )
           and IsBound( p.isPolynomial )
           and p.isPolynomial;
end;


#############################################################################
##
#F  LaurentPolynomialsOps.DefaultRing( <L> )  . . . . . . default ring of <L>
##
LaurentPolynomialsOps.DefaultRing := function( L )
    local   R,  l,  v;

    # get the ring
    R := L[1].baseRing;
    v := L[1].valuation;
    for l  in L  do
        if R <> l.baseRing  then
            Error( "sorry, all elements must have the same base ring" );
        fi;
        if l.valuation < v  then
            v := l.valuation;
        fi;
    od;

    # return the full (Laurent) polynom ring over <R>
    if v < 0  then
        return LaurentPolynomialRing(R);
    else
        return PolynomialRing(R);
    fi;

end;


#############################################################################
##

#V  Polynomials . . . . . . . . . . . . . . . . . . domain of all polynomials
##
Polynomials            := rec();

Polynomials.isDomain   := "true";

Polynomials.name       := "Polynomials";

Polynomials.isFinite   := false;
Polynomials.size       := "infinity";

Polynomials.operations := Copy( DomainOps );
PolynomialsOps         := Polynomials.operations;

PolynomialsOps.\in := function( p, Polynomials )
    return     IsRec( p )
           and IsBound( p.isPolynomial )
           and p.isPolynomial
           and 0 <= p.valuation;
end;


#############################################################################
##
#F  PolynomialsOps.DefaultRing( <L> ) . . . . . . . . . . default ring of <L>
##
PolynomialsOps.DefaultRing := LaurentPolynomialsOps.DefaultRing;


#############################################################################
##

#V  PolynomialOps   . . . . . . . . . . . . . . . operations for a polynomial
##
PolynomialOps := rec();
PolynomialOps.name := "PolynomialOps";


#############################################################################
##
#F  PolynomialOps.Print( <f> )  . . . . . .  print a polynomial in a nice way
##
PolynomialOps.Print := function( pol )
    local   R;

    R := pol.baseRing;
    if IsBound(R.indeterminate) and IsBound(R.indeterminate.name)  then
        DisplayPolynomial( pol, [R.indeterminate.name], false );
    else
    DisplayPolynomial( pol, ["X(",R,")"], false );
    fi;
end;


#############################################################################
##
#F  PolynomialOps.Value( <f>, <x> ) . . . . . . . . . . evaluate a polynomial
##
PolynomialOps.Value := function( f, x )
    local  val, i, id;

    id  := x ^ 0;
    val := 0 * id;
    i := Length(f.coefficients);
    while 0 < i  do
        val := val*x + id*f.coefficients[i];
        i := i-1;
    od;
    if 0 <> f.valuation  then
        val := val * x^f.valuation;
    fi;
    return val;

end;


#############################################################################
##
#F  PolynomialOps.Depth( <f> )  . . . . . . . . . ((R[x_1])[x_2])..)[x_depth]
##
PolynomialOps.Depth := function( f )
    local    R,  d;

    if not IsBound(f.depth)  then
        R := f.baseRing;
        d := 1;
        while IsPolynomialRing(R) or IsLaurentPolynomialRing(R)  do
            R := R.baseRing;
            d := d + 1;
        od;
        f.depth := d;
    fi;
    return f.depth;

end;


#############################################################################
##
#F  PolynomialOps.GroundRing( <f> ) . . . . . . . . . . .  ground ring of <f>
##
PolynomialOps.GroundRing := function( f )
    local    R;

    if not IsBound(f.groundRing)  then
        R := f.baseRing;
        while IsPolynomialRing(R) or IsLaurentPolynomialRing(R)  do
            R := R.baseRing;
        od;
        f.groundRing := R;
    fi;
    return f.groundRing;
end;


#############################################################################
##
#F  PolynomialOps.Degree( <f> ) . . . . . . . . . . .  degree of a polynomial
##
Degree := function( obj )
    if not IsBound(obj.degree)  then
    obj.degree := obj.operations.Degree(obj);
    fi;
    return obj.degree;
end;

PolynomialOps.Degree := function( f )
    return ( Length(f.coefficients)-1 ) + f.valuation;
end;


#############################################################################
##
#F  PolynomialOps.LeadingCoefficient( <f> ) . . .  leading coefficient of <f>
##
LeadingCoefficient := function( f )
    if not IsBound(f.leadingCoefficient)  then
    f.leadingCoefficient := f.operations.LeadingCoefficient(f);
    fi;
    return f.leadingCoefficient;
end;

PolynomialOps.LeadingCoefficient := function( f )
    if 0 < Length(f.coefficients)  then
        return f.coefficients[Length(f.coefficients)];
    else
        return f.baseRing.zero;
    fi;
end;


#############################################################################
##
#F  PolynomialOps.\= . . . . . . . . . . . .  equaltity test for polynomials
##
PolynomialOps.\= := function( l, r )
    local   R,  sum,  i;

    # handle the case <something> = <polynomial>
    if not IsPolynomial( l )  then
        return false;

    # handle the case <polynomial> = <something>
    elif not IsPolynomial( r )  then
        return false;

    # handle special case of iterated polynomials
    elif l.operations.Depth(l) <> r.operations.Depth(r)  then
        return false;

    # return 'false' if <l> and <r> have different rings
    elif l.operations.GroundRing(l) <> r.operations.GroundRing(r)  then
        return false;

    # compare valuation
    elif l.valuation <> r.valuation  then
        return false;

    # compare coefficients
    else
        return l.coefficients = r.coefficients;
    fi;

end;


#############################################################################
##
#F  PolynomialOps.\< . . . . . . . . . . . . . . . comparison of polynomials
##
PolynomialOps.\< := function( l, r )

    if not IsPolynomial( l )  then
        return l < r.coefficients;
    elif not IsPolynomial( r )  then
        return l.coefficients < r;
    elif l.operations.Depth(l) <> r.operations.Depth(r)  then
        return l.operations.Depth(l) < r.operations.Depth(r);
    elif l.valuation <> r.valuation  then
        return l.valuation < r.valuation;
    elif Length( l.coefficients ) = Length( r.coefficients )  then
        return l.coefficients < r.coefficients;
    else
        return Length( l.coefficients ) < Length( r.coefficients );
    fi;

end;


#############################################################################
##
#F  PolynomialOps.\+ . . . . . . . . . . . . . . . .  sum of two polynomials
##
PolynomialOps.\+ := function( l, r )
    local   R,  sum,  val,  vdf,  i;

    # handle the case that one argument is a list
    if IsList(l)  then
        return List( l, x -> x+r );
    elif IsList(r)  then
        return List( r, x -> l+x );
    fi;

    # handle the case <scalar> + <polynomial>
    if not IsPolynomial(l)  then
        if IsInt(l)  then
            l := l * r.baseRing.one;
        elif not l in r.baseRing  then
            Error( "<l> must lie in the base ring of <r>" );
        fi;
        R := r.baseRing;
        l := R.operations.Polynomial( R, [l], 0 );
    fi;

    # handle the case <polynomial> + <scalar>
    if not IsPolynomial(r)  then
        if IsInt(r)  then
            r := r * l.baseRing.one;
        elif not r in l.baseRing  then
            Error( "<r> must lie in the base ring of <l>" );
        fi;
        R := l.baseRing;
    r := R.operations.Polynomial( R, [r], 0 );
    fi;

    # handle special case of iterated polynomials
    if l.operations.Depth(l) < r.operations.Depth(r)  then
        R := DefaultRing(r);
        l := l * R.one;
    elif r.operations.Depth(r) < l.operations.Depth(l)  then
        R := DefaultRing(l);
        r := R.one * r;
    fi;

    # give up if we have different rings
    if l.operations.GroundRing(l) <> r.operations.GroundRing(r)  then
        Error( "polynomials must have the same ring" );

    # if <l> is the null polynomial return <r>
    elif Length(l.coefficients) = 0  then
    return r;

    # if <r> is the null polynomial return <l>
    elif Length(r.coefficients) = 0  then
        return l;

    # sum of two polynomials
    else

        # get a common ring
        R := l.baseRing;

        # add both coefficients lists
        if l.valuation < r.valuation  then
        val := l.valuation;
            sum := Copy( l.coefficients );
            vdf := r.valuation - l.valuation;
            for i  in [ Length(sum)+1 .. Length(r.coefficients)+vdf ] do
            sum[i] := R.zero;
            od;
            for i  in [ 1 .. Length(r.coefficients) ]  do
            sum[i+vdf] := sum[i+vdf] + r.coefficients[i];
            od;
        else
        val := r.valuation;
            sum := Copy( r.coefficients );
            vdf := l.valuation - r.valuation;
            for i  in [ Length(sum)+1 .. Length(l.coefficients)+vdf ] do
            sum[i] := R.zero;
            od;
            for i  in [ 1 .. Length(l.coefficients) ]  do
            sum[i+vdf] := l.coefficients[i] + sum[i+vdf];
            od;
        fi;
    fi;

    # return the sum
    return R.operations.Polynomial( R, sum, val );

end;


#############################################################################
##
#F  PolynomialOps.\- . . . . . . . . . . . . . difference of two polynomials
##
PolynomialOps.\- := function( l, r )
    local   R,  dif,  val,  vdf,  i;

    # handle the case that one argument is a list
    if IsList(l)  then
        return List( l, x -> x-r );
    elif IsList(r)  then
        return List( r, x -> l-x );
    fi;

    # handle the case <scalar> - <polynomial>
    if not IsPolynomial(l)  then
        if IsInt(l)  then
            l := l * r.baseRing.one;
        elif not l in r.baseRing  then
            Error( "<l> must lie in the base ring of <r>" );
        fi;
        R := r.baseRing;
        l := R.operations.Polynomial( R, [l], 0 );
    fi;

    # handle the case <polynomial> - <scalar>
    if not IsPolynomial(r)  then
        if IsInt(r)  then
            r := r * l.baseRing.one;
        elif not r in l.baseRing  then
            Error( "<r> must lie in the base ring of <l>" );
        fi;
        R := l.baseRing;
    r := R.operations.Polynomial( R, [r], 0 );
    fi;

    # handle special case of iterated polynomials
    if l.operations.Depth(l) < r.operations.Depth(r)  then
        R := DefaultRing(r);
        l := l * R.one;
    elif r.operations.Depth(r) < l.operations.Depth(l)  then
        R := DefaultRing(l);
        r := R.one * r;
    fi;

    # give up if we have different rings
    if l.operations.GroundRing(l) <> r.operations.GroundRing(r)  then
        Error( "polynomials must have the same ring" );

    # if <l> is the null polynomial return -<r>
    elif Length(l.coefficients) = 0  then
    return -r;

    # if <r> is the null polynomial return <l>
    elif Length(r.coefficients) = 0  then
        return l;

    # difference of two polynomials
    else

        # get a common ring
        R := l.baseRing;

        # add both coefficients lists
        if l.valuation < r.valuation  then
        val := l.valuation;
            dif := Copy( l.coefficients );
            vdf := r.valuation - l.valuation;
            for i  in [ Length(dif)+1 .. Length(r.coefficients)+vdf ] do
            dif[i] := R.zero;
            od;
            for i  in [ 1 .. Length(r.coefficients) ]  do
            dif[i+vdf] := dif[i+vdf] - r.coefficients[i];
            od;
        else
        val := r.valuation;
            dif := Copy( r.coefficients ) * (R.zero-R.one);
            vdf := l.valuation - r.valuation;
            for i  in [ Length(dif)+1 .. Length(l.coefficients)+vdf ] do
            dif[i] := R.zero;
            od;
            for i  in [ 1 .. Length(l.coefficients) ]  do
            dif[i+vdf] := l.coefficients[i] + dif[i+vdf];
            od;
        fi;
    fi;

    # return the dif
    return R.operations.Polynomial( R, dif, val );

end;


#############################################################################
##
#F  PolynomialOps.\* . . . . . . . . . . . . . .  product of two polynomials
##
PolynomialOps.\* := function( l, r )
    local   R,  prd,  z,  m,  n,  i,  j,  val;

    # handle the case that one argument is a list
    if IsList(l)  then
        return List(l, x -> x*r);
    elif IsList(r)  then
        return List(r, x -> l*x);

    # handle the case <scalar> * <polynomial>
    elif not IsPolynomial(l)  then
        if IsInt(l) and IsBound(r.baseRing.one)  then
            l := l * r.baseRing.one;
        elif not l in r.baseRing  then
        Error("<l> must lie in the base ring of <r>");
    fi;
    R := r.baseRing;
    if l = r.baseRing.zero or r.coefficients = []  then
        prd := [];
            val := 0;
        else
        prd := l * r.coefficients;
            val := r.valuation;
    fi;

    # handle the case <polynomial> * <scalar>
    elif not IsPolynomial(r)  then
        if IsInt(r) and IsBound(l.baseRing.one)  then
            r := l.baseRing.one * r;
    elif not r in l.baseRing  then
        Error("<r> must lie in the base ring of <l>");
    fi;
    R := l.baseRing;
    if r = l.baseRing.zero or l.coefficients = []  then
        prd := [];
            val := 0;
        else
        prd := l.coefficients * r;
            val := l.valuation;
    fi;

    # the ground fields must be equal
    elif l.operations.GroundRing(l) <> r.operations.GroundRing(r)  then
        Error( "polynomials must have the same ring" );

    # handle special case of iterated polynomials
    elif l.operations.Depth(l) <> r.operations.Depth(r)  then
        if l.operations.Depth(l) < r.operations.Depth(r)  then
            prd := List( r.coefficients, x -> l*x );
            val := r.valuation;
            R   := r.baseRing;
        else
            prd := List( l.coefficients, x -> x*r );
            val := l.valuation;
            R   := l.baseRing;
        fi;

    # multiply two polynomials
    else

        # get a common ring
        R := l.baseRing;

        # fold the product
        m := Length(l.coefficients);
        n := Length(r.coefficients);
        prd := [];
        for i  in [ 1 .. m+n-1 ]  do
            z := R.zero;
            for j  in [ Maximum(i+1-n,1) .. Minimum(m,i) ]  do
                z := z + l.coefficients[j]*r.coefficients[i+1-j];
            od;
            prd[i] := z;
        od;
        val := l.valuation + r.valuation;
    fi;

    # return the product
    return R.operations.Polynomial( R, prd, val );

end;


#############################################################################
##
#F  PolynomialOps.\/ . . . . . . . . . . . . . . quotient of two polynomials
##
PolynomialOps.\/ := function( l, r )
    local   R,  quo,  z,  m,  n,  i,  j,  val;

    # handle the case that <l> argument is a list
    if IsList(l)  then
        return List(l, x -> x/r);

    # handle the case <scalar> / <polynomial>
    elif not IsPolynomial(l)  then
        if IsInt(l) and IsBound(r.baseRing.one)  then
            l := l * r.baseRing.one;
        elif not l in r.baseRing  then
        Error("<l> must lie in the base ring of <r>");
    fi;
        r := r^-1;
    R := r.baseRing;
    if l = r.baseRing.zero or r.coefficients = []  then
        quo := [];
            val := 0;
        else
        quo := l * r.coefficients;
            val := r.valuation;
    fi;

    # handle the case <polynomial> / <scalar>
    elif not IsPolynomial(r)  then
        if IsInt(r) and IsBound(l.baseRing.one)  then
            r := l.baseRing.one * r;
    elif not r in l.baseRing  then
        Error("<r> must lie in the base ring of <l>");
    fi;
    R := l.baseRing;
    if r = l.baseRing.zero or l.coefficients = []  then
            Error( "<r> must be non-zero" );
        else
        quo := l.coefficients * (1/r);
            val := l.valuation;
    fi;

    # the ground fields must be equal
    elif l.operations.GroundRing(l) <> r.operations.GroundRing(r)  then
        Error( "polynomials must have the same ring" );

    # handle special case of iterated polynomials
    elif l.operations.Depth(l) <> r.operations.Depth(r)  then
        if l.operations.Depth(l) < r.operations.Depth(r)  then
            quo := List( r.coefficients, x -> l/x );
            val := r.valuation;
            R   := r.baseRing;
        else
            quo := List( l.coefficients, x -> x/r );
            val := l.valuation;
            R   := l.baseRing;
        fi;

    # compute the quotient and return in this case
    else
        quo := Quotient( l, r );
        if quo = false  then
            Error( "cannot divide <l> by <r>" );
        fi;
        return quo;
    fi;

    # return the product
    return R.operations.Polynomial( R, quo, val );

end;


#############################################################################
##
#F  PolynomialOps.\^ . . . . . . . . . . . . . . . . . power of a polynomial
##
PolynomialOps.\^ := function( l, r )
    local   R,  pow;

    # <l> must be a polynomial and <r> a non-negative integer
    if not IsPolynomial(l)  then
        Error("<l> must be a polynomial");
    fi;
    if not IsInt(r)  then
        Error("<r> must be an integer");
    fi;

    # invert <l> if necessary and possible
    if r < 0  then
        R := LaurentPolynomialRing( l.baseRing );
        l := R.operations.Quotient( R, R.one, l );
        if l = false  then
            Error( "cannot invert <l> in the laurent polynomial ring" );
        fi;
        r := -r;
    fi;

    # if <r> is zero, return x^0, if <r> is one return <l>
    R := l.baseRing;
    if r = 0  then
        return Polynomial( R, [ R.one ] );
    elif r = 1  then
        return l;
    fi;

    # if <l> is of degree less than 2, return
    if Length(l.coefficients) = 0  then
        return l;
    elif Length(l.coefficients) = 1  then
        return Polynomial( R, [ l.coefficients[1]^r ], l.valuation*r );
    fi;

    # use repeated squaring
    pow := Polynomial( R, [ R.one ] );
    while 0 < r  do
        if r mod 2 = 1  then
        pow := pow * l;
            r   := r - 1;
        fi;
        if 1 < r  then
        l := l * l;
            r := r / 2;
        fi;
    od;

    # return the power
    return pow;

end;


#############################################################################
##

#F  DisplayPolynomial( <p>, <str> ) . . . . . . . . . .  display a polynomial
##
DisplayPolynomial := function( arg )

    local   i,            # loop variable
        p,              # the polynomial
            wantNl,         # try if we want a newline at the end
        str,            # string list for indeterminate
        nam,            # loop variable for <str>
            cof,            # coefficients of <p>
            len,            # length of coefficients of <p>
            isMinus,        # if 'true' print "x^2 - x"
            isPrimeField;   # if 'true' print "Z(3)*(2*x^2+1)"

    # get arguments, <wantNl> is 'false' by default
    p := arg[1];
    str := arg[2];
    if IsString(str)  then str := [ str ];  fi;
    wantNl := 3 = Length(arg) and arg[3];

    # check if the base ring of <p> is a finite prime field
    isPrimeField := IsFiniteField(p.baseRing) and (p.baseRing.degree = 1);

    # if the polynomial is of length one, ignore <isPrimeField>
    if  Length(p.coefficients) < 2  then isPrimeField := false;  fi;

    # check if can print a minus sign
    isMinus := IsInt(p.baseRing.zero);

    # catch some trivial cases
    if p.valuation = 0 and Length(p.coefficients) < 2  then
    if Length(p.coefficients) = 0  then
        Print("0*");
        for nam  in str  do
        Print( nam );
        od;
        Print("^0");
    elif p.coefficients[1] = p.baseRing.one  then
        for nam  in str  do
        Print( nam );
        od;
        Print("^0");
    elif isMinus and p.coefficients[1] = -p.baseRing.one  then
        Print("-");
        for nam  in str  do
        Print( nam );
        od;
        Print("^0");
    else
        Print( p.coefficients[1], "*" );
        for  nam  in str  do
        Print( nam );
        od;
        Print("^0");
    fi;
    return;
    fi;

    # if <isPrimeField> is 'true' print the one in Zp
    if isPrimeField  then  Print( p.baseRing.one, "*(" );  fi;

    # print the polynomial
    len := Length(p.coefficients);
    for i  in Reversed([ 1 .. len ])  do
        cof := p.coefficients[i];
        if cof <> p.baseRing.zero  then
            if i < len and isMinus and IsInt(cof) and cof < 0  then
                Print(" - ");
                cof := -cof;
            elif i = len and isMinus and cof = -p.baseRing.one  then
                Print("-");
                cof := -cof;
            elif i < len  then
                Print(" + ");
            fi;
            if 1 = i + p.valuation  then
                if isPrimeField  then
                    Print(Int(cof));
                elif IsFFE(cof) or IsInt(cof) then
                    Print(cof);
        else
            Print( "(", cof, ")" );
                fi;
            elif 2 = i + p.valuation  then
                if cof <> p.baseRing.one  then
                       if isPrimeField  then
                        Print( Int(cof), "*" );
            for nam in str  do
                Print( nam );
            od;
            elif IsFFE(cof) or IsInt(cof)  then
                        Print( cof, "*" );
            for nam in str  do
                Print( nam );
            od;
                    else
                        Print( "(", cof, ")*" );
            for nam in str  do
                Print( nam );
            od;
                    fi;
                else
            for nam in str  do
            Print( nam );
            od;
                fi;
            else
                if cof <> p.baseRing.one  then
                    if isPrimeField  then
                        Print( Int(cof), "*" );
            for nam in str  do
                Print( nam );
            od;
            elif IsFFE(cof) or IsInt(cof)  then
                        Print( cof, "*" );
            for nam in str  do
                Print( nam );
            od;
                    else
                        Print( "(", cof, ")*" );
            for nam in str  do
                Print( nam );
            od;
                    fi;
                else
            for nam in str  do
            Print( nam );
            od;
                fi;
            if i+p.valuation-1 < 0  then
            Print( "^(", i+p.valuation-1, ")" );
        else
                    Print( "^", i+p.valuation-1 );
        fi;
            fi;
        fi;
    od;
    if isPrimeField  then Print(")");  fi;
    if wantNl  then Print("\n");  fi;

end;


#############################################################################
##

#E  Emacs . . . . . . . . . . . . . . . . . . . . . . . local emacs variables
##
##  Local Variables:
##  mode:               outline
##  outline-regexp:     "#F\\|#V\\|#E"
##  eval:               (local-set-key "\t" 'tab-to-tab-stop)
##  eval:               (hide-body)
##  End:
##