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

Text File | 1993-05-05 | 46.3 KB | 1,582 lines

#############################################################################
##
#A  polyrat.g                   GAP library                      Frank Celler
##
#A  @(#)$Id: polyrat.g,v 3.19 1993/02/19 15:50:51 fceller Exp $
##
#Y  Copyright 1990-1992,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
##
##  This file contains functions for polynomials over finite fields.
##
#H  $Log: polyrat.g,v $
#H  Revision 3.19  1993/02/19  15:50:51  fceller
#H  fixed a bug in '-'
#H
#H  Revision 3.18  1993/02/12  12:00:08  martin
#H  multiplication must check for zero before using 'FastPolynomial'
#H
#H  Revision 3.17  1993/02/11  16:47:56  fceller
#H  added <list> + <polynomial>
#H
#H  Revision 3.16  1993/02/09  14:25:55  martin
#H  made undefined globals local
#H
#H  Revision 3.15  1993/02/08  11:57:07  fceller
#H  some speed ups in polynomial arithmetic
#H
#H  Revision 3.14  1993/01/22  12:03:48  fceller
#H  added missing 'isPolynomialRing' in 'RationalsPolynomials'
#H
#H  Revision 3.13  1993/01/08  09:13:08  fceller
#H  added Hensel 'Factors' for polynomials over the rationals
#H
#H  Revision 3.12  1993/01/04  11:19:14  fceller
#H  added 'Gcd'
#H
#H  Revision 3.11  1992/12/16  19:47:27  martin
#H  replaced quoted record names with escaped ones
#H
#H  Revision 3.10  1992/12/07  07:40:04  fceller
#H  added <f> mod <n>
#H
#H  Revision 3.9  1992/11/25  12:34:54  fceller
#H  added 'Degree'
#H
#H  Revision 3.8  1992/11/16  17:48:31  fceller
#H  added 'String' support
#H
#H  Revision 3.7  92/11/16  12:23:40  fceller
#H  added Laurent polynomials
#H
#H  Revision 3.6  92/10/28  09:15:08  fceller
#H  fixed 'mod' for zero polynomial
#H
#H  Revision 3.5  1992/07/24  07:08:07  fceller
#H  improved '*' to allow <nullpoly> * <int>
#H
#H  Revision 3.4  1992/07/23  09:20:48  fceller
#H  improved '*' to allow <list> * <polynomial>
#H
#H  Revision 3.3  1992/06/17  07:06:04  fceller
#H  moved '<somedomain>.operations.Polynomial' function to "<somedomain>.g"
#H
#H  Revision 3.2  1992/06/01  07:32:24  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;
if not IsBound(InfoPoly3)    then InfoPoly3   := Ignore;  fi;


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

#V  RationalsPolynomialOps  . . . . . . . . . . polynomial over the rationals
##
RationalsPolynomialOps := Copy( PolynomialOps );
RationalsPolynomialOps.name := "RationalsPolynomialOps";


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

    # 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 (IsRec(l) and IsBound(l.isPolynomial) and l.isPolynomial)  then

        # <r> must have the rationals as base ring
        if Rationals <> r.baseRing  then
            Error( "<r> must have the rationals as base ring" );
        fi;

        # <l> must lie in the base ring of <r>
        if not IsRat(l)  then
            Error( "<l> must lie in the base ring of <r>" );
        fi;

    # if <l> is trivial return <r>
    if l = 0  then
        return r;
    fi;

    # otherwise convert <l> into a polynomial
        l := RationalsOps.FastPolynomial( Rationals, [l], 0 );
    fi;

    # handle the case <polynomial> + <scalar>
    if not (IsRec(r) and IsBound(r.isPolynomial) and r.isPolynomial)  then

        # <l> must have the rationals as base ring
        if Rationals <> l.baseRing  then
            Error( "<l> must have the rationals as base ring" );
        fi;

        # <r> must lie in the base ring of <l>
        if not IsRat(r)  then
            Error( "<r> must lie in the base ring of <l>" );
        fi;

    # if <r> is trivial return <l>
    if r = 0  then
        return l;
    fi;

    # otherwise convert <r> into a polynomial
    r := RationalsOps.FastPolynomial( Rationals, [r], 0 );
    fi;

    # depth greater than one are handle by our superclass
    if 1 <> l.operations.Depth(l) or 1 <> r.operations.Depth(r)  then
    return PolynomialOps.\+( l, r );

    # give up if we have rings other then the rationals
    elif Rationals <> l.baseRing  then
        Error( "<l> must have the rationals as base ring" );
    elif Rationals <> r.baseRing  then
        Error( "<r> must have the rationals as base 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 and the valuation minimum;
        vdf := r.valuation - l.valuation;

        # if <r>.valuation is the minimum shift <l>
        if r.valuation < l.valuation  then
            val := r.valuation;
            sum := ShiftedCoeffs( l.coefficients, -vdf );

            # the rationals are commutative
            AddCoeffs( sum, r.coefficients );

        # if <l>.valuation is the minimum shift <r>
        elif l.valuation < r.valuation  then
            val := l.valuation;
            sum := ShiftedCoeffs( r.coefficients, vdf );

            # the rationals are commutative
            AddCoeffs( sum, l.coefficients );

        # otherwise they are equal
        else
            sum := SumCoeffs( l.coefficients, r.coefficients );
            val := l.valuation;
        fi;

        # return the sum
    if 0 < Length(sum) and 0 <> sum[1]  then
            return RationalsOps.FastPolynomial( Rationals, sum, val );
    else
            return RationalsOps.Polynomial( Rationals, sum, val );
    fi;

    fi;

end;


#############################################################################
##
#F  RationalsPolynomialOps.\- . . . . . . . . . . . . diff of two polynomials
##
RationalsPolynomialOps.\- := function( l, r )
    local   dif,  val,  vdf;

    # 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 (IsRec(l) and IsBound(l.isPolynomial) and l.isPolynomial)  then

        # <r> must have the rationals as base ring
        if Rationals <> r.baseRing  then
            Error( "<r> must have the rationals as base ring" );
        fi;

        # <l> must lie in the base ring of <r>
        if not IsRat(l)  then
            Error( "<l> must lie in the base ring of <r>" );
        fi;

    # if <l> is trivial return -<r>
    if l = 0  then
            return RationalsOps.FastPolynomial(
               Rationals,
               (-1) * r.coefficients,
               r.valuation );
    fi;

    # otherwise convert <l> into a polynomial
        l := RationalsOps.FastPolynomial( Rationals, [l], 0 );
    fi;

    # handle the case <polynomial> - <scalar>
    if not (IsRec(r) and IsBound(r.isPolynomial) and r.isPolynomial)  then

        # <l> must have the rationals as base ring
        if Rationals <> l.baseRing  then
            Error( "<l> must have the rationals as base ring" );
        fi;

        # <r> must lie in the base ring of <l>
        if not IsRat(r)  then
            Error( "<r> must lie in the base ring of <l>" );
        fi;

    # if <r> is trivial return <l>
    if r = 0  then
        return l;
    fi;

    # otherwise convert <r> into a polynomial
    r := RationalsOps.FastPolynomial( Rationals, [r], 0 );
    fi;

    # depth greater than one are handle by our superclass
    if 1 <> l.operations.Depth(l) or 1 <> r.operations.Depth(r)  then
    return PolynomialOps.\-( l, r );

    # give up if we have rings other then the rationals
    elif Rationals <> l.baseRing  then
        Error( "<l> must have the rationals as base ring" );
    elif Rationals <> r.baseRing  then
        Error( "<r> must have the rationals as base 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 and the valuation minimum;
        vdf := r.valuation - l.valuation;

        # if <r>.valuation is the minimum shift <l>
        if r.valuation < l.valuation  then
            val := r.valuation;
            dif := ShiftedCoeffs( l.coefficients, -vdf );

            # the rationals are commutative
            AddCoeffs( dif, r.coefficients, -1 );

        # if <l>.valuation is the minimum shift <r>
        elif l.valuation < r.valuation  then
            val := l.valuation;
            dif := (-1)*ShiftedCoeffs( r.coefficients, vdf );

            # the rationals are commutative
            AddCoeffs( dif, l.coefficients );

        # otherwise they are equal
        else
            dif := Copy(l.coefficients);
            AddCoeffs( dif, r.coefficients, -1 );
            val := l.valuation;
        fi;

        # return the difference
    if 0 < Length(dif) and 0 <> dif[1]  then
            return RationalsOps.FastPolynomial( Rationals, dif, val );
    else
            return RationalsOps.Polynomial( Rationals, dif, val );
    fi;

    fi;

end;


#############################################################################
##
#F  RationalsPolynomialOps.\*  . . . . . . . . .  product of two polynomials
##
RationalsPolynomialOps.\* := function( l, r )
    local   R,  prd,  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 (IsRec(l) and IsBound(l.isPolynomial) and l.isPolynomial)  then

        # <r> must have the rationals as base ring
        if Rationals <> r.baseRing  then
            Error( "<r> must have the rationals as base ring" );
        fi;

        # <l> must lie in the base ring of <r>
        if not IsRat(l)  then
        Error( "<l> must lie in the base ring of <r>" );
    fi;
    if l = 0 or r.coefficients = []  then
        prd := [];
        val := 0;
        else
        prd := l * r.coefficients;
        val := r.valuation;
    fi;

        # compute the product
        prd            := RationalsOps.FastPolynomial( Rationals, prd, val );
        prd.depth      := 1;
        prd.groundRing := Rationals;

        # and return
        return prd;

    # handle the case <polynomial> * <scalar>
    elif not (IsRec(r) and IsBound(r.isPolynomial) and r.isPolynomial)  then

        # <l> must have the rationals as base ring
        if Rationals <> l.baseRing  then
            Error( "<l> must have the rationals as base ring" );
        fi;

        # <r> must lie in the base ring of <r>
    if not IsRat(r)  then
        Error( "<r> must lie in the base ring of <l>" );
    fi;
    if r = l.baseRing.zero or l.coefficients = []  then
        prd := [];
        val := 0;
        else
        prd := l.coefficients * r;
        val := l.valuation;
    fi;

        # compute the product
        prd            := RationalsOps.FastPolynomial( Rationals, prd, val );
        prd.depth      := 1;
        prd.groundRing := Rationals;

        # and return
        return prd;

    # our superclass will handle different depth
    elif 1 <> l.operations.Depth(l) or 1 <> r.operations.Depth(r) then
        return PolynomialOps.\*( l, r );

    # give up if we have rings other then the rationals
    elif Rationals <> l.baseRing  then
        Error( "<l> must have the rationals as base ring" );
    elif Rationals <> r.baseRing  then
        Error( "<r> must have the rationals as base ring" );

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

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

    # multiply two polynomials
    else

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

        # use 'ProductCoeffs' in order to fold product
        prd := ProductCoeffs( l.coefficients, r.coefficients );
    val := l.valuation + r.valuation;

        # compute the product
        prd            := R.operations.FastPolynomial( R, prd, val );
        prd.depth      := 1;
        prd.groundRing := Rationals;

        # and return
        return prd;

    fi;

end;


#############################################################################
##
#F  RationalsPolynomialOps.\mod  . . . . . . .  remainder of two polynomials
##
RationalsPolynomialOps.\mod := function( l, r )
    local   R,  rem,  val,  vdf;

    # <l> must be a polynomial
    if not IsPolynomial(l)  then
        Error( "<l> must be a polynomial" );
    fi;

    # if <r> is a integer reduce the coefficients of <l>
    if IsInt(r)  then
    rem := List( l.coefficients, x -> x mod r );
    return Polynomial( l.baseRing, rem, l.valuation );
    fi;


    # otherwise <r> must be a non-zero polynomial
    if not IsPolynomial(r)  then
        Error( "<r> must be a polynomial" );
    fi;
    if Length(r.coefficients) = 0  then
        Error( "<r> must be non zero" );
    fi;

    # our superclass will handle different depth
    if l.operations.Depth(l) <> r.operations.Depth(r)  then
        return PolynomialOps.\mod( l, r );

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

    # multiply two polynomials
    else

        # if one is a Laurent polynomial use 'Mod'
        if l.valuation < 0 or r.valuation < 0  then
            return Mod( DefaultRing(l,r), l, r );
        fi;

        # get a common ring and the value difference
        R   := l.baseRing;
        vdf := r.valuation - l.valuation;

        # if <r>.valuation is the minimum shift <l>
        if r.valuation < l.valuation  then
            val := r.valuation;
            rem := ShiftedCoeffs( l.coefficients, -vdf );
            ReduceCoeffs( rem, r.coefficients );

        # if <l>.valuation is the minimum shift <r>
        elif l.valuation < r.valuation  then
            r   := ShiftedCoeffs( r.coefficients, vdf );
            rem := RemainderCoeffs( l.coefficients, r );
            val := l.valuation;

        # otherwise they are equal
        else
            rem := RemainderCoeffs( l.coefficients, r.coefficients );
            val := l.valuation;
        fi;
    fi;

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

end;


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

    # <l> must be a polynomial over the rationals and <r> an integer
    if not IsPolynomial(l) or l.baseRing <> Rationals  then
        Error( "<l> must be a polynomial over the rationals" );
    fi;
    if not IsInt(r)  then
        Error( "<r> must be an integer" );
    fi;

    # invert <l> if necessary
    if r < 0  then
        R := LaurentPolynomialRing( l.baseRing );
        l := R.operations.Quotient( R, R.one, l );
        r := -r;
    fi;
    R := l.baseRing;

    # if <r> is zero, return x^0
    if r = 0  then
        return RationalsOps.FastPolynomial( Rationals, [R.one], 0 );

    # if <r> is one return <l>
    elif r = 1  then
        return l;
    fi;

    # if <l> is trivial return
    if Length(l.coefficients) = 0  then
        return l;

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

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

    # return the power
    return RationalsOps.FastPolynomial( Rationals, pow, val );

end;


#############################################################################
##
#F  RationalsPolynomialOps.String( <f> )  . . . . . construct a pretty string
##
RationalsPolynomialOps.String := function( f )
    local   x,  i,  d,  v,  s,  l;

    # find a name for the indeterminate
    x := Indeterminate(f.baseRing);
    if IsBound(x.name)  then x := x.name;  else x := "x";  fi;

    # run through the coefficients of <f>
    v := f.valuation-1;
    l := Length(f.coefficients);
    for i  in Reversed([ 1 .. l ])  do
        d := f.coefficients[i];
        if 0 <> d  then
            if i = l and d = 1 and i+v <> 0  then
                s := "";
            elif i = l and d = 1  then
                s := "1";
            elif i = l and d = -1 and i+v <> 0  then
                s := "-";
            elif i = l and d = -1  then
                s := "-1";
            elif i = l  then
                s := String(d);
            elif d = 1 and i+v <> 0  then
                s := ConcatenationString( s, "+" );
            elif d = 1  then
                s := ConcatenationString( s, "+1" );
            elif d = -1 and i+v <> 0  then
                s := ConcatenationString( s, "-" );
            elif d = -1  then
                s := ConcatenationString( s, "-1" );
            elif d < 0  then
                s := ConcatenationString( s, String(d) );
            elif 0 < d  then
                s := ConcatenationString( s, "+", String(d) );
            else
                Error( "internal error in 'RationalsPolynomialOps.String'" );
            fi;
            if i+v < 0 or 1 < i+v then
                s := ConcatenationString( s, x, "^", String(i+v) );
            elif i+v = 1  then
                s := ConcatenationString( s, x );
            fi;
    fi;
    od;

    # catch a special case
    if l = 0  then s := "0";  fi;
    return s;

end;



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

#V  RationalsPolynomials  . . . . .  domain of polynomials over the rationals
##
RationalsPolynomials             := Copy( Polynomials );
RationalsPolynomials.name        := "RationalsPolynomials";
RationalsPolynomials.operations  := Copy( FieldPolynomialRingOps );
RationalsPolynomialsOps          := RationalsPolynomials.operations;

# show that this a polynomial ring
RationalsPolynomials.isPolynomialRing := true;

# rationals polynomials form a ring
RationalsPolynomials.isDomain := true;
RationalsPolynomials.isRing   := true;

# set known properties
RationalsPolynomials.isFinite := false;
RationalsPolynomials.size     := "infinity";

# add properties of polynom ring over a field
RationalsPolynomials.isCommutativeRing         := true;
RationalsPolynomials.isIntegralRing            := true;
RationalsPolynomials.isUniqueFactorizationRing := true;
RationalsPolynomials.isEuclideanRing           := true;

# set one, zero and base ring
RationalsPolynomials.one  := Polynomial( Rationals, [ 1 ] );
RationalsPolynomials.zero := Polynomial( Rationals, [] );
RationalsPolynomials.baseRing := Rationals;


#############################################################################
##
#F  RationalsPolynomialsOps.\in  . . . . . . . . . . . . . . membership test
##
RationalsPolynomialsOps.\in := function( p, RationalsPolynomials )
    return     IsRec( p )
           and IsBound( p.isPolynomial )
           and p.isPolynomial
           and IsField( p.baseRing )
       and 0 <= p.valuation
           and p.baseRing = Rationals;
end;


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


#############################################################################
##
#F  RationalsPolynomialsOps.PowerMod( <R>, <g>, <e>, <m> )  . . . . . . power
##
RationalsPolynomialsOps.PowerMod := function( R, g, e, m )

    # if <m> is of degree one return the zero polynomial
    if Degree(m) = 1  then
    return R.zero;

    # if <e> is zero return one
    elif e = 0  then
    return R.one;
    fi;

    # reduce polynomial
    g := R.operations.Mod( R, g, m );

    # and invert if necessary
    if e < 0  then
    g := R.operations.QuotientMod( R, R.one, g, m );
        if g = false  then
        Error( "<g> must be invertable module <m>" );
        fi;
        e := -e;
    fi;

    # use 'PowerModCoeffs' to power polynomial
    g := ShiftedCoeffs( g.coefficients, g.valuation );
    m := ShiftedCoeffs( m.coefficients, m.valuation );
    return Polynomial( R.baseRing, PowerModCoeffs( g, e, m ) );

end;


#############################################################################
##
#F  RationalsPolynomialsOps.IntegerPolynomial( <R>, <f> ) . . . . convert <f>
##
RationalsPolynomialsOps.IntegerPolynomial := function( R, f )
    local   lcm,  c;

    # compute lcm of denominator
    lcm := 1;
    for c  in f.coefficients  do
    lcm := LcmInt( lcm, Denominator(c) );
    od;

    # remove all denominators
    f := f * lcm;

    # remove gcd of coefficients
    return f * (1/Gcd(f.coefficients));

end;


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

#F  LandauMignotteBoundGcd( <f>, <g> )    . . . .  bound for gcd of <f> and <g>
##
LandauMignotteBoundGcd := function( f, g )
    local  A,  B,  sum,  srt;

    # compute norm of <f>
    sum := Sum( f.coefficients, x -> x^2 );
    srt := RootInt( sum, 2 );
    if srt*srt < sum  then srt := srt+1;  fi;
    A := srt / AbsInt(LeadingCoefficient(f));

    # compte norm of <g>
    sum := Sum( g.coefficients, x -> x^2 );
    srt := RootInt( sum, 2 );
    if srt*srt < sum  then srt := srt+1;  fi;
    B := srt / AbsInt(LeadingCoefficient(g));

    # return the bound
    return   2^Minimum( Degree(f), Degree(g) )
           * GcdInt( LeadingCoefficient(f), LeadingCoefficient(g) )
           * Minimum( A, B );

end;


#############################################################################
##
#F  RationalsPolynomialsOps.GcdModPrime( <f>, <g>, <p>, <a> ) . . gcd mod <p>
##
RationalsPolynomialsOps.GcdModPrime := function( f, g, p, a )
    local   F,  c,  tab,  cof,  log,  x;

    # compute in the finite field F_<p>
    f := f mod p;
    g := g mod p;
    F := GF(p);
    c := Gcd( Polynomial( F, f.coefficients*F.one, f.valuation ),
              Polynomial( F, g.coefficients*F.one, g.valuation ) );

    # convert back to integers
    tab := IntegerTable(F);
    cof := [];
    for x  in [ 1 .. Length(c.coefficients) ]  do
        log := LogVecFFE(c.coefficients,x);
        if log = false  then
            Add( cof, 0 );
        else
            Add( cof, (tab[log+1]*a) mod p );
        fi;
    od;
    return Polynomial( f.baseRing, cof, c.valuation );

end;


#############################################################################
##
#F  RationalsPolynomialsOps.Gcd( <R>, <f>, <g> )  . . . . gcd of <f> and <g>
##
RationalsPolynomialsOps.Gcd := function( R, f, g )

    # check trivial cases
    if -1 = Degree(f)  then
        return g;
    elif -1 = Degree(g)  then
        return f;
    elif 0 = Degree(f) or 0 = Degree(g)  then
        return R.one;
    fi;

    # convert polynomials into integer polynomials
    f := R.operations.IntegerPolynomial( R, f );
    g := R.operations.IntegerPolynomial( R, g );
    InfoPoly2( "#I  <f> = ", f, "\n" );
    InfoPoly2( "#I  <g> = ", g, "\n" );

    # return the standard associate
    return StandardAssociate( R, R.operations.IGcd( R, f, g ) );

end;

RationalsPolynomialsOps.IGcd := function( R, f, g )
    local   a,  t;

    # compute the Landau Mignotte bound for the gcd
    t := rec( prime := 1 );
    t.bound := 2 * Int(LandauMignotteBoundGcd( f, g )+1);
    InfoPoly2( "#I  Landau-Mignotte bound = ", t.bound/2, "\n" );

    # avoid gcd of leading coefficients
    a := GcdInt( LeadingCoefficient(f), LeadingCoefficient(g) );
    repeat

        # start with first prime avoiding gcd of leading coefficients
        repeat t.prime := NextPrimeInt(t.prime);  until a mod t.prime <> 0;

        # compute modular gcd with leading coefficient <a>
        t.gcd := RationalsPolynomialsOps.GcdModPrime( f, g, t.prime, a );
        InfoPoly2( "#I  gcd mod ", t.prime, " = ", t.gcd, "\n" );

        # loop until we have success
        repeat
            if 0 = Degree(t.gcd)  then
                InfoPoly2( "#I  <f> and <g> are relative prime\n" );
                return R.one;
            fi;
        until RationalsPolynomialsOps.Gcd1( R, t, a, f, g );
    until t.correct;

    # return the gcd
    return t.gcd;

end;

RationalsPolynomialsOps.Gcd1 := function( R, t, a, f, g )
    local   G,  P;

    # <P> will hold the product of primes use so far
    t.modulo := t.prime;

    # <G> will hold the approximation of the gcd
    G := t.gcd;

    # use next prime until we reach the Landau-Mignotte bound
    while t.modulo < t.bound  do
        repeat t.prime := NextPrimeInt(t.prime);  until a mod t.prime <> 0;

        # compute modular gcd
        t.gcd := RationalsPolynomialsOps.GcdModPrime( f, g, t.prime, a );
        InfoPoly2( "#I  gcd mod ", t.prime, " = ", t.gcd, "\n" );

        # if the degree of <C> is smaller we started with wrong <p>
        if Degree(t.gcd) < Degree(G)  then
            InfoPoly2( "#I  found lower degree, restarting\n" );
            return false;
        fi;

        # if the degrees of <C> and <G> are equal use chinese remainder
        if Degree(t.gcd) = Degree(G)  then
            P := G;
            G := RationalsPolynomialsOps.GcdCRT(G,t.modulo,t.gcd,t.prime);
            t.modulo := t.modulo * t.prime;
            InfoPoly2( "#I  gcd mod ", t.modulo, " = ", G, "\n" );
            if G = P  then
                t.correct :=     Quotient(R,f,G)<>false
                             and Quotient(R,g,G)<>false;
                if t.correct  then
                    InfoPoly2( "#I  found correct gcd\n" );
                    t.gcd := G;
                    return true;
                fi;
            fi;
        fi;
    od;

    # check if <G> is correct but return 'true' in any case
    t.correct := Quotient(R,f,G)<>false and Quotient(R,g,G)<>false;
    t.gcd := G;
    return true;

end;

RationalsPolynomialsOps.GcdCRT := function( f, p, g, q )
    local   min,  cf,  lf,  cg,  lg,  i,  P,  m,  r;

    # remove valuation
    min := Minimum( f.valuation, g.valuation );
    if f.valuation <> min  then
        cf := ShiftedCoeffs( f.coefficients, f.valuation - min );
    else
        cf := ShallowCopy(f.coefficients);
    fi;
    lf := Length(cf);
    if g.valuation <> min  then
        cg := ShiftedCoeffs( g.coefficients, g.valuation - min );
    else
        cg := ShallowCopy(g.coefficients);
    fi;
    lg := Length(cg);

    # use chinese remainder
    r := [ p, q ];
    P := p * q;
    m := P/2;
    for i  in [ 1 .. Minimum(lf,lg) ]  do
        cf[i] := ChineseRem( r, [ cf[i], cg[i] ] );
        if m < cf[i]  then cf[i] := cf[i] - P;  fi;
    od;
    if lf < lg  then
        for i  in [ lf+1 .. lg ]  do
            cf[i] := ChineseRem( r, [ 0, cg[i] ] );
            if m < cf[i]  then cf[i] := cf[i] - P;  fi;
        od;
    elif lg < lf  then
        for i  in [ lg+1 .. lf ]  do
            cf[i] := ChineseRem( r, [ cf[i], 0 ] );
            if m < cf[i]  then cf[i] := cf[i] - P;  fi;
        od;
    fi;

    # return the polynomial
    return Polynomial( f.baseRing, cf, min );

end;


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

#F  RationalsPolynomialsOps.QuotientModPrime(<R>,<f>,<g>,<p>) . . .  quotient
##
RationalsPolynomialsOps.QuotientModPrime := function( R, f, g, p )
    local   m,  n,  i,  k,  c,  q,  R,  val;

    # get base ring
    R := R.baseRing;

    # reduce <f> and <g> mod <p>
    f := f mod p;
    g := g mod p;

    # 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>, compute mod <p>
    q := [];
    n := Length( g.coefficients );
    m := Length( f.coefficients ) - n;
    f := ShallowCopy( f.coefficients );
    for i  in [0 .. m ]  do
        c := f[m-i+n] / g.coefficients[n] mod p;
        for k  in [ 1 .. n ]  do
            f[m-i+k] := ( f[m-i+k] - c * g.coefficients[k] ) mod p;
        od;
        q[m-i+1] := c;
    od;

    # 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  LandauMignotteBound( <f> )  . . . . . . . . . . bound for factors of <f>
##
LandauMignotteBound := function( f )
    local   sum,  srt;

    sum := Sum( f.coefficients, x -> x^2 );
    srt := RootInt( sum, 2 );
    if srt*srt < sum  then srt := srt + 1;  fi;
    return 2^Degree(f) * srt;
end;


#############################################################################
##
#F  RationalsPolynomialsOps.SquareHensel( <R>, <f>, <lp> )
##
RationalsPolynomialsOps.SquareHensel := function( R, f, lp )

    local   p,              # prime of <lp>
            q,              # current modulus
            q1,             # last modulus
            q2,             # <q> / 2
            lc,             # leading coefficient of <f>
            k,              # Landau-Mignotte bound
            prd,            # product of <lp> or <l>
            pd2,            # product of <l>
            rp,             # representation of gcd(<lp>)
            tab,            # integer table of <fp>
            rep,            # lifted representation of gcd(<lp>)
            fac,            # true factor found so far
            fcn,            # index of true factor in <l>
            dis,            # distance of <f> and <l>
            cor,            # correction
            rcr,            # inverse corrections
            quo,            # quotient
            sum,            # temp
            aa, bb,         # left and right subproducts
            l,              # factors mod <q>
            lq1,            # factors mod <q1>
            g, cof, log,    # temps
            i, j,  x;       # loop

    # get the leading coefficient of <f>
    lc := LeadingCoefficient(f);

    # compute the Landau-Mignotte bound
    k := 2 * AbsInt(lc) * LandauMignotteBound(f);
    InfoPoly2( "#I  Landau-Mignotte bound = ", k/2/AbsInt(lc), "\n" );

    # get the prime <p>
    p := lp[1].baseRing.char;

    # compute the representation of 1 mod <p>
    prd := Product(lp);
    rp  := GcdRepresentation( List( lp, x -> Quotient( prd, x ) ) );

    # convert representation <rp> back to the integers
    tab := IntegerTable( lp[1].baseRing );
    rep := [];
    for g  in rp  do
        cof := [];
        for x  in [ 1 .. Length(g.coefficients) ]  do
            log := LogVecFFE(g.coefficients,x);
            if log = false  then
                Add( cof, 0 );
            else
                Add( cof, tab[log+1] );
            fi;
        od;
        Add( rep, Polynomial( Rationals, cof, g.valuation ) );
    od;

    # convert factors <lp> back to the integers
    l := [];
    for g  in lp  do
        cof := [];
        for x  in [ 1 .. Length(g.coefficients) ]  do
            log := LogVecFFE(g.coefficients,x);
            if log = false  then
                Add( cof, 0 );
            else
                Add( cof, tab[log+1] );
            fi;
        od;
        Add( l, Polynomial( Rationals, cof, g.valuation ) );
    od;

    # and check it
    InfoPoly3( "#C  difference mod ", p, " = ",
               ((1/lc mod p)*f-Product(l)) mod p, " (should be 0)\n" );

    # start Hensel until <q> is greater than k
    q   := p^2;
    q1  := p;
    fac := [];
    while q1 < k  do
        InfoPoly2( "#I  computing mod ", q, "\n" );

        # compute discrepancies for factors
        dis := ( ( ( 1/lc mod q ) * f - Product(l) ) mod q ) * (1/q1);
        InfoPoly2( "#I  discrepancy = ", dis, "\n" );

        # compute corrections for factors
        lq1 := ShallowCopy(l);
        cor := [];
        for i  in [ 1 .. Length(l) ]  do
            cor[i] := rep[i] * dis mod l[i] mod q1;
            InfoPoly2( "#I  ", i, ".th correction = ", cor[i], "\n" );
            if 0 < Length(cor[i].coefficients)  then
                l[i] := l[i] + q1 * cor[i];
            fi;
            InfoPoly2( "#I  ", i, ".th factor = ", l[i], "\n" );
        od;

        # if this is not the last step update <rep>
        if q < k  then

            # compute the products Product(<lq1>)/<lq1>[i]
            aa := [ R.one ];
            for i  in [ 1 .. Length(l)-1 ]  do
                aa[i+1] := aa[i] * lq1[i] mod q;
            od;
            bb := [ R.one ];
            for i  in [ 1 .. Length(l)-1 ]  do
                bb[i+1] := bb[i] * lq1[Length(l)+1-i] mod q;
            od;
            bb := Reversed(bb);

            # compute discrepancies for inverses
            rcr := [];
            for i  in [ 1 .. Length(l) ]  do
                prd := aa[i] * bb[i] mod q;
                dis := ( ( 1 - rep[i] * prd ) mod l[i] mod q ) * (1/q1);
                InfoPoly2( "#I  ", i, ".th rep discrepancy = ", dis, "\n" );
                sum := 0 * dis;
                prd := prd mod q1;
                for j  in [ 1 .. Length(l) ]  do
                    if j <> i  then
                      quo := R.operations.QuotientModPrime(R,prd,lq1[j],q1);
                       sum := sum + quo * cor[j];
                    fi;
                od;
                rcr[i] := (rep[i]*(dis-rep[i]*sum)) mod lq1[i] mod q1;
                InfoPoly2("#I  ", i, ".th rep correction = ", rcr[i], "\n");
            od;

            # correct the inverses
            for i  in [ 1 .. Length(l) ]  do
                rep[i] := ( rep[i] + q1 * rcr[i] ) mod q;
                InfoPoly2( "#I  ", i, ".th representation = ", rep[i], "\n" );
            od;
        fi;

        # check for true factors
        fcn := [];
        q2  := q/2;
        for i  in [ 1 .. Length(l) ]  do
            if 0 = Length(cor[i].coefficients)  then
                g := [];
                for j  in [ 1 .. Length(l[i].coefficients) ]  do
                    g[j] := ( l[i].coefficients[j] * lc ) mod q;
                    if q2 < g[j]  then
                        g[j] := g[j] - q;
                    fi;
                od;
                g := Polynomial( Rationals, g * (1/Gcd(g)), l[i].valuation );
                if f.coefficients[1] mod g.coefficients[1] = 0
                   and LeadingCoefficient(f)
                       mod LeadingCoefficient(g) = 0
                then
                    quo := Quotient( R, f, g );
                    if quo <> false  then
                        InfoPoly2( "#I  found true factor = ", l[i], "\n" );
                        Add( fcn, i );
                        f := quo;
                    fi;
                fi;
            fi;
        od;

        # if we have found a true factor update everything
        if 0 < Length(fcn)  then

            # append factors to <fac>
            Append( fac, Sublist( l, fcn ) );

            # compute new Landau-Mignotte bound
            lc := LeadingCoefficient(f);
            k  := 2 * AbsInt(lc) * LandauMignotteBound(f);
            InfoPoly2( "#I  new L-M bound = ", k/2/AbsInt(lc), "\n" );

            # remove true factors from <l> and corresponding <rep>
            prd := Product( Sublist( l, fcn ) ) mod q;
            fcn := Filtered( [ 1 .. Length(l) ], x -> not x in fcn );
            l   := Sublist( l, fcn );
            rep := List( Sublist( rep, fcn ), x -> prd * x );

            # reduce <rep>[i] mod <l>[i]
            for i  in [ 1 .. Length(l) ]  do
                rep[i] := rep[i] mod l[i] mod q;
            od;
        fi;

        # and check it
        if q < k and InfoPoly3 <> Ignore  then
            prd := Product(l);
            for i  in [ 1 .. Length(l) ]  do
                InfoPoly3( "#C  rep[", i, "] * P[", i, "] = ",
                           (rep[i]*Quotient(R,prd,l[i])) mod l[i] mod q,
                           " (should be 1)\n" );
            od;
        fi;
        InfoPoly3( "#C  difference mod ", q, " = ",
                   ((1/lc mod q)*f-Product(l)) mod q,
                   " (should be 0)\n" );

        # square modulus
        q1 := q;
        q  := q^2;
    od;

    # convert the coefficients c_i such that |c_i| <= k/2
    rep := [];
    q2  := q1 / 2;
    for g  in l  do
        cof := ShallowCopy(g.coefficients);
        for x  in [ 1 .. Length(cof) ]  do
            cof[x] := cof[x]*lc mod q1;
            if q2 < cof[x]  then
                cof[x] := cof[x] - q1;
            fi;
        od;

        # make polynomial primitive
        g := Polynomial( Rationals, cof * (1/Gcd(cof)), g.valuation );

        # check if <cof> is a true factor
        if f.coefficients[1] mod g.coefficients[1] = 0
           and LeadingCoefficient(f) mod LeadingCoefficient(g) = 0
        then
            quo := Quotient( R, f, g );
            if quo <> false  then
                InfoPoly2( "#I  found a true factor = ", g, "\n" );
                Add( fac, g );
                f  := quo;
                lc := LeadingCoefficient(f);
            else
                Add( rep, g );
            fi;
        else
            Add( rep, g );
        fi;
    od;

    # return the result
    return rec( polynomial     := f,
                factors        := fac,
                remaining      := rep,
                landauMignotte := LandauMignotteBound(f) * AbsInt(lc),
                primePower     := q1,
                primePower2    := q1/2 );

end;


#############################################################################
##
#F  RationalsPolynomialsOps.FactorsModPrime( <R>, <f> )      find suitable prime
##
##  <f> must be squarefree.  We test 3 "small" and 2 "big" primes.
##
RationalsPolynomialsOps.FactorsModPrime := function( R, f )

    local   i,          # loop
            lc,         # leading coefficient of <f>
            p,          # current prime
            PR,         # polynomial ring over F_<p>
            fp,         # <f> in <R>
            lp,         # factors of <fp>
            min,        # minimal number of factors so far
            P,          # best prime so far
            LP,         # factorization of <f> mod <P>
            deg,        # possible degrees of factors
            t,          # return record
            tmp;

    # set minimal number of factors to the degree of <f>
    min := Degree(f)+1;
    lc  := LeadingCoefficient(f);

    # find a suitable prime
    t := rec();
    p := 1;
    for i  in [ 1 .. 5 ]  do

        # reset <p> to big prime after first 3 test
        if i = 4  then p := Maximum( p, 1000 );  fi;

        # find a prime not dividing <lc> and <f>_<p> squarefree
        repeat
            repeat
                p  := NextPrimeInt(p);
            until lc mod p <> 0;
            PR  := PolynomialRing(GF(p));
            tmp := 1/lc mod p;
            fp  := List( f.coefficients, x->(tmp*x mod p)*PR.baseRing.one );
            fp  := Polynomial( PR.baseRing, fp, f.valuation );
        until 0 = Degree(Gcd(fp,Derivative(fp)));

        # factorise <f> modulo <p>
        lp := PR.operations.Factors( PR, fp );

        # if <fp> is irreducible so is <f>
        if 1 = Length(lp)  then
            InfoPoly2( "#I  <f> mod ", p, " is irreducible\n" );
            t.isIrreducible := true;
            return t;
        else
            InfoPoly2( "#I  found ", Length(lp), " factors mod ", p, "\n" );
        fi;

        # choose a minimal prime with minimal number of factors
        if Length(lp) <= min  then
            min := Length(lp);
            P   := p;
            LP  := lp;
        fi;

        # compute the possible degrees
        tmp := Set( List( Combinations( List( lp, Degree ) ), g -> Sum(g) ) );
        if 1 = i  then
            deg := tmp;
        else
            deg := Intersection( deg, tmp );
        fi;

        # if there is only one possible degree != 0 then <f> is irreducible
        if 2 = Length(deg)  then
            InfoPoly2( "#I  <f> must be irreducible, only one degree left\n" );
            t.isIrreducible := true;
            return t;
        fi;

    od;

    # return the chosen prime
    InfoPoly2( "#I  choosing prime ", P, " with ", Length(LP), " factors\n" );
    t.isIrreducible := false;
    t.prime         := P;
    t.factors       := LP;
    t.degrees       := deg;
    return t;

end;


#############################################################################
##
#F  RationalsPolynomialsOps.FactorsSquarefree( <R>, <f> ) . .  factors of <f>
##
##  <f> must be square free and must have a constant term.
##
RationalsPolynomialsOps.FactorsSquarefree := function( R, f )
    local   t,  h,  fac,  rem,  c,  cbs,  i,  j,  sl,  prd,  cof,  q;

    # find a suitable prime, if <f> is irreducible return
    t := R.operations.FactorsModPrime( R, f );
    if t.isIrreducible  then return [f];  fi;
    InfoPoly2( "#I  using prime ", t.prime, " for factorization\n" );

    # start Hensel
    h := R.operations.SquareHensel( R, f, t.factors );

    # check for any direct factors
    f   := h.polynomial;
    fac := h.factors;
    rem := h.remaining;
    InfoPoly2( "#I  found ", Length(fac), " true factors, ",
               Length(rem), " remaining\n" );

    # try all combinations
    c := 2;
    while 2*c <= Length(rem)  do
        cbs := Combinations( [ 1 .. Length(rem) ], c );
        j   := 0;
        repeat
            j  := j + 1;
            sl := Sublist( rem, cbs[j] );

            # the degree of the combination must be possible
            if Sum(List(sl,Degree)) in t.degrees  then

                # compute the product and reduce
                prd := Product( sl );
                cof := [];
                for i  in [ 1 .. Length(prd.coefficients) ]  do
                    cof[i] := prd.coefficients[i] mod h.primePower;
                    if h.primePower2 < cof[i]  then
                        cof[i] := cof[i] - h.primePower;
                    fi;
                od;

                # make the product primitive
                cof := cof * (1/Gcd(cof));
                prd := Polynomial( Rationals, cof, prd.valuation );

                # make sure that the quotient has a chance
                if f.coefficients[1] mod prd.coefficients[1] <> 0
                   or LeadingCoefficient(f)
                      mod LeadingCoefficient(prd) <> 0
                then
                    InfoPoly2( "#I  ignoring combination ", cbs[j], "\n" );
                    q := false;
                else
                    InfoPoly2( "#I  testing combination ", cbs[j], "\n" );
                    q := Quotient( f, prd );
                fi;
            else
                InfoPoly2( "#I  skipping combination ", cbs[j], "\n" );
                q := false;
            fi;
        until q <> false or j = Length(cbs);

        # if the quotient is ok, remove the factors, set <f> and <fac>
        if q <> false  then
            f := q;
            InfoPoly2( "#I  found true factor = ", prd, "\n" );
            Add( fac, prd );
            rem := List( Filtered( [ 1 .. Length(rem) ],
                             x -> not x in cbs[j] ),
                         y -> rem[y] );
        else
            c := c + 1;
        fi;
    od;

    # the product of the remaing one is a true factor
    if 0 < Length(rem)  then
        prd := Product(rem);
        cof := [];
        for i  in [ 1 .. Length(prd.coefficients) ]  do
            cof[i] := prd.coefficients[i] mod h.primePower;
            if h.primePower2 < cof[i]  then
                cof[i] := cof[i] - h.primePower;
            fi;
        od;
        cof := cof * (1/Gcd(cof));
        Add( fac, Polynomial( Rationals, cof, prd.valuation ) );
    fi;

    # return the factors
    return fac;

end;


#############################################################################
##
#F  RationalsPolynomialsOps.IFactors( <R>, <f> )  . . . . . .  factors of <f>
##
RationalsPolynomialsOps.IFactors := function( R, f )
    local   l,  v,  g,  q,  s,  r,  x;

    # if <f> is trivial return
    if 0 = Length( f.coefficients )  then
        InfoPoly2( "#I  <f> is trivial\n" );
        return [ f ];
    fi;

    # remove a valuation
    v := f.valuation;
    f := Polynomial( R.baseRing, f.coefficients );
    x := Indeterminate(R.baseRing);

    # if <f> is constant return
    if 0 = Degree(f)  then
        InfoPoly2( "#I  <f> is a power of x\n" );
        s    := List( [ 1 .. v ], f -> x );
        s[1] := s[1] * f.coefficients[1];
        return s;
    fi;

    # if <f> linear return
    if 1 = Degree(f)  then
        InfoPoly2( "#I  <f> is a linear\n" );
        s := List( [ 1 .. v ], f -> x );
        Add( s, f );
        return s;
    fi;

    # make <f> integral, primitive and square free
    g := R.operations.IGcd( R, f, Derivative(f) );
    q := R.operations.IntegerPolynomial( R, Quotient( R, f, g ) );
    q := q * SignInt(LeadingCoefficient(q));
    InfoPoly2( "#I  factorizing ", q, "\n" );

    # and factorize <q>
    if Degree(q) < 2  then
        InfoPoly2( "#I  <f> is a linear power\n" );
        s := [ q ];
    else
        s := R.operations.FactorsSquarefree( R, q );
    fi;

    # find factors of <g>
    for r  in s  do
        if 0 < Degree(g) and Degree(g) mod Degree(r) = 0  then
            q := Quotient( R, g, r );
            while 0 < Degree(g) and q <> false  do
                Add( s, r );
                g := q;
                if Degree(g) mod Degree(r) = 0  then
                    q := Quotient( R, g, r );
                else
                    q := false;
                fi;
            od;
        fi;
    od;

    # sort the factors
    Append( s, List( [ 1 .. v ],  f -> x ) );
    Sort(s);

    # return the (primitive) factors
    return s;

end;


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

    # handle trivial case
    if Degree(f) < 2  then
    f.factors := [ f ];
    elif Length(f.coefficients) = 1  then
    r := List( [ 1 .. f.valuation ], x -> Indeterminate(f.baseRing) );
    r[1] := r[1] * f.coefficients[1];
    f.factors := r;
    fi;

    # do we know the factors
    if IsBound(f.factors) and R.baseRing = f.baseRing  then
        return f.factors;
    fi;

    # compute the integer factors
    r := R.operations.IFactors( R, f );

    # convert into standard associates and sort
    r := List( r, x -> StandardAssociate(R,x) );
    Sort(r);

    # correct leading term
    r[1] := LeadingCoefficient(f) * r[1];

    # and return
    if f.baseRing = R.baseRing  then
    f.factors := r;
    fi;
    return r;

end;


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

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