#############################################################################
##
## HAPPRIME - polynomials.gi
## Functions, Operations and Methods to extend GAP's polynomials
## Paul Smith
##
## Copyright (C) 2008
## Paul Smith
## National University of Ireland Galway
##
## This file is part of HAPprime.
##
## HAPprime is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## HAPprime is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <
https://www.gnu.org/licenses/>.
##
##
#############################################################################
#####################################################################
## <#GAPDoc Label="TermsOfPolynomial_manPolynomial">
## <ManSection>
## <Attr Name="TermsOfPolynomial" Arg="poly"/>
##
## <Returns>
## List of pairs
## </Returns>
## <Description>
## Returns a list of the terms in the polynomial.
## This list is a list of pairs of the form <C>[mon, coeff]</C> where
## <C>mon</C> is a monomial and <C>coeff</C> is the coefficient of that
## monomial in the polynomial. The monomials are sorted according to
## the total degree/lexicographic order (the same as the in
## <Ref Oper="MonomialGrLexOrdering" BookName="ref"/>).
## </Description>
## </ManSection>
## <#/GAPDoc>
#####################################################################
InstallMethod(TermsOfPolynomial,
[IsPolynomial],
function(poly)
local terms, extrep, fam, one, i, term, mon, coeff;
terms := [];
extrep := ExtRepPolynomialRatFun(poly);
if IsEmpty(extrep) then
return [];
fi;
fam := FamilyObj(poly);
i := 1;
repeat
mon := extrep[i];
coeff := extrep[i+1];
term := [PolynomialByExtRep(fam, [mon, One(coeff)]), coeff];
Add(terms, term);
i := i + 2;
until i > Length(extrep);
return terms;
end
);
#####################################################################
#####################################################################
## <#GAPDoc Label="IsMonomial_manPolynomial">
## <ManSection>
## <Attr Name="IsMonomial" Arg="poly" Label="for polynomial"/>
##
## <Returns>
## Boolean
## </Returns>
## <Description>
## Returns <K>true</K> if <A>poly</A> is a monomial, i.e. the polynomial
## contains only one term.
## </Description>
## </ManSection>
## <#/GAPDoc>
#####################################################################
InstallMethod(IsMonomial,
"for polynomial",
[IsPolynomial],
function(poly)
return Length(ExtRepPolynomialRatFun(poly)) = 2;
end
);
#####################################################################
#####################################################################
## <#GAPDoc Label="UnivariateMonomialsOfMonomial_manPolynomial">
## <ManSection>
## <Attr Name="UnivariateMonomialsOfMonomial" Arg="mon"/>
##
## <Returns>
## List
## </Returns>
## <Description>
## Returns a list of the univariate monomials of the largest order
## whose product equals <A>mon</A>.
## The univariate monomials are sorted according to
## their indeterminate number.
## </Description>
## </ManSection>
## <#/GAPDoc>
#####################################################################
InstallMethod(UnivariateMonomialsOfMonomial,
[IsPolynomial],
function(poly)
local extrep, fam, coeffring, one, mon, i, unimons;
if IsOne(poly) or IsZero(poly) then
return [poly];
fi;
extrep := ExtRepPolynomialRatFun(poly);
fam := FamilyObj(poly);
coeffring := Ring(extrep[2]);
one := One(coeffring);
if Length(extrep) > 2 or extrep[2] <> one then
Error("<mon> must be a monomial");
fi;
mon := extrep[1];
i := 1;
unimons := [];
repeat
Add(unimons, PolynomialByExtRep(fam, [[mon[i], mon[i+1]], one]));
i := i + 2;
until i > Length(mon);
return unimons;
end
);
#####################################################################
#####################################################################
## <#GAPDoc Label="IndeterminateAndExponentOfUnivariateMonomial_manPolynomial">
## <ManSection>
## <Attr Name="IndeterminateAndExponentOfUnivariateMonomial" Arg="mon"/>
##
## <Returns>
## List
## </Returns>
## <Description>
## Returns a list <C>[indet, exp]</C> where <C>indet</C> is the indeterminate
## of the univariate monomial <A>mon</A> and <C>exp</C> is the exponent
## of that indeterminate in the monomial. If <A>mon</A> is an element in the
## coefficient ring (i.e. the monomial contains no indeterminates) then the
## first element will be <A>mon</A> with an exponent of zero.
## If <A>mon</A> is not a univariate monomial, then <K>fail</K> is returned.
## </Description>
## </ManSection>
## <#/GAPDoc>
#####################################################################
InstallMethod(IndeterminateAndExponentOfUnivariateMonomial,
[IsPolynomial],
function(poly)
local extrep, one, fam;
extrep := ExtRepPolynomialRatFun(poly);
if IsEmpty(extrep[1]) then
return [extrep[2], 0];
fi;
if Length(extrep) <> 2 then
return fail;
fi;
one := One(extrep[2]);
fam := FamilyObj(poly);
return [PolynomialByExtRep(fam, [[extrep[1][1], 1], one]), extrep[1][2]];
end
);
#####################################################################
#####################################################################
## <#GAPDoc Label="IndeterminatesOfPolynomial_manPolynomial">
## <ManSection>
## <Attr Name="IndeterminatesOfPolynomial" Arg="poly"/>
##
## <Returns>
## List
## </Returns>
## <Description>
## Returns a list of the indeterminates used in the polynomial <A>poly</A>.
## </Description>
## </ManSection>
## <#/GAPDoc>
#####################################################################
InstallMethod(IndeterminatesOfPolynomial,
[IsPolynomial],
function(poly)
# Is it a constant?
if IsEmpty(ExtRepPolynomialRatFun(poly)) or
IsEmpty(ExtRepPolynomialRatFun(poly)[1]) then
return [];
else
return IndeterminatesOfPolynomialRing(DefaultRing(poly));
fi;
end
);
#####################################################################
#####################################################################
## <#GAPDoc Label="ReduceIdeal_manPolynomial">
## <ManSection>
## <Heading>ReduceIdeal</Heading>
## <Oper Name="ReduceIdeal" Arg="I, O" Label="for Ideal"/>
## <Oper Name="ReduceIdeal" Arg="rels, O" Label="for list of relations"/>
##
## <Returns>
## Ideal or list
## </Returns>
## <Description>
## For an ideal <A>I</A> returns an ideal containing a reduced generating set
## for the ideal, i.e. one in which no monomial in a relation in <A>I</A> is
## divisible by the leading term of another polynomial in <A>I</A>.
## The monomial ordering to be used is
## specified by <A>O</A> (see <Ref Sect="Monomial Orderings" BookName="ref"/>).
## The ideal can instead be specified by a list of relations <A>rels</A>,
## in which case a reduced list of relations is returned.
## </Description>
## </ManSection>
## <#/GAPDoc>
#####################################################################
InstallOtherMethod(ReduceIdeal,
"for empty ideal",
[IsEmpty, IsMonomialOrdering],
function(ideal, order)
return [];
end
);
#####################################################################
InstallOtherMethod(ReduceIdeal,
"for Ideal",
[IsPolynomialRingIdeal, IsMonomialOrdering],
function(ideal, order)
return Ideal(
LeftActingRingOfIdeal(ideal),
ReduceIdeal(GeneratorsOfIdeal(ideal), order));
end
);
#####################################################################
InstallMethod(ReduceIdeal,
"for list of relations",
[IsHomogeneousList and IsRationalFunctionCollection, IsMonomialOrdering],
function(I, order)
local ideal, i, len, changed, poly;
ideal := ShallowCopy(I);
len := Length(ideal);
repeat
i := 1;
changed := false;
while i <= len do
poly := PolynomialReducedRemainder(
ideal[i], ideal{Difference([1..len],[i])}, order);
if poly <> ideal[i] then
changed := true;
fi;
if IsZero(poly) then
ideal[i] := ideal[len];
Unbind(ideal[len]);
len := len - 1;
else
ideal[i] := poly;
i := i + 1;
fi;
od;
until not changed;
for i in [1..Length(ideal)] do
ideal[i] := ideal[i] / LeadingCoefficientOfPolynomial(ideal[i], order);
od;
return ideal;
end
);
#####################################################################
#####################################################################
## <#GAPDoc Label="ReducedPolynomialRingPresentation_manPolynomial">
## <ManSection>
## <Oper Name="ReducedPolynomialRingPresentation" Arg="R, I[, avoid]"/>
## <Oper Name="ReducedPolynomialRingPresentationMap" Arg="R, I[, avoid]"/>
##
## <Returns>
## List
## </Returns>
## <Description>
## For a polynomial ring <A>R</A> and a list of relations <A>I</A> in that
## ring, returns a list <C>[S, J]</C> representing a polynomial quotient ring
## <M>S/J</M> which is isomorphic to the ring <M>R/I</M>, but which involves
## the minimal number of ring indeterminates. The indeterminates in <C>S</C>
## will be distinct from thise in <A>R</A>, and an optional argument
## <A>avoid</A> can be used to give a list of further indeterminates to avoid
## when creating the ring <C>S</C>.
## <P/>
## The extended version of this function,
## <Ref Oper="ReducedPolynomialRingPresentationMap"/>, returns an additional
## third element to the list, which contains two lists giving the mapping
## between the new ring indeterminates and the old ring indeterminates. The
## first list is of polynomials in the original indeterminates, the
## second the equivalent polynomials in the new ring indeterminates.
## </Description>
## </ManSection>
## <#/GAPDoc>
#####################################################################
InstallMethod(ReducedPolynomialRingPresentation,
"with nothing to avoid",
[IsPolynomialRing, IsHomogeneousList and IsRationalFunctionCollection, IsHomogeneous
List],
function(ring, I, avoid)
return ReducedPolynomialRingPresentationMap(ring, I, avoid){[1..2]};
end
);
#####################################################################
InstallOtherMethod(ReducedPolynomialRingPresentation,
"with nothing to avoid",
[IsPolynomialRing, IsHomogeneousList and IsRationalFunctionCollection],
function(ring, I)
return ReducedPolynomialRingPresentationMap(ring, I, []){[1..2]};
end
);
#####################################################################
InstallOtherMethod(ReducedPolynomialRingPresentation,
"for empty relations",
[IsPolynomialRing, IsEmpty, IsHomogeneousList],
function(ring, I, avoid)
local newring, newrelations, map;
# Just copy the ring
newring := PolynomialRing(
CoefficientsRing(ring), Length(IndeterminatesOfPolynomialRing(ring)),
Concatenation(avoid, IndeterminatesOfPolynomialRing(ring)));
newrelations := HAPPRIME_SwitchPolynomialIndeterminates(
ring, newring, I);
return [newring, newrelations ];
end
);
#####################################################################
InstallOtherMethod(ReducedPolynomialRingPresentationMap,
"with nothing to avoid",
[IsPolynomialRing, IsHomogeneousList and IsRationalFunctionCollection],
function(ring, I)
return ReducedPolynomialRingPresentationMap(ring, I, []);
end
);
#####################################################################
InstallOtherMethod(ReducedPolynomialRingPresentationMap,
"for empty relations and nothing to avoid",
[IsPolynomialRing, IsEmpty],
function(ring, I)
return ReducedPolynomialRingPresentationMap(ring, I, []);
end
);
#####################################################################
InstallOtherMethod(ReducedPolynomialRingPresentationMap,
"for empty relations",
[IsPolynomialRing, IsEmpty, IsHomogeneousList],
function(ring, I, avoid)
local newring, newrelations, map;
# Just copy the ring
newring := PolynomialRing(
CoefficientsRing(ring), Length(IndeterminatesOfPolynomialRing(ring)),
Concatenation(avoid, IndeterminatesOfPolynomialRing(ring)));
newrelations := HAPPRIME_SwitchPolynomialIndeterminates(
ring, newring, I);
return [newring, newrelations,
[IndeterminatesOfPolynomialRing(ring),
IndeterminatesOfPolynomialRing(newring)] ];
end
);
#####################################################################
InstallMethod(ReducedPolynomialRingPresentationMap,
[IsPolynomialRing, IsHomogeneousList and IsRationalFunctionCollection, IsHomogeneousList],
function(ring, I, avoid)
local ideal, indetorder, order, indet, i, t, t2, indets, allindets, ord,
removedindets, removedrelations, removedring, unimons, relation, poly, len,
nomod, newring, newrelations, map;
# Check for unit in the ideal - if so, we return a trivial ring
if Length(I) = 1 and IsOne(I[1]) then
return [Ring(Zero(I[1])), [ ] ];
fi;
ideal := ShallowCopy(I);
# indeterminates will be removed from indets and added to removedindets
# as we do our reduction
indets := ShallowCopy(IndeterminatesOfPolynomialRing(ring));
removedindets := [];
removedrelations := [];
# Are any of the leading terms of the ideal single indeterminates?
# If so then we can (after reduction) remove those indeterminates and those
# terms of the ideal
repeat
indet := false;
# Find a relation in the ideal that involves a solitary indeterminate
for i in [1..Length(ideal)] do
for t in TermsOfPolynomial(ideal[i]) do
unimons := UnivariateMonomialsOfMonomial(t[1]);
if Length(unimons) = 1 and unimons[1] = IndeterminateOfUnivariateRationalFunction(unimons[1]) then
# This term involves a solitary indeterminate.
# Check that no other term in this relation also involves this
# indeterminate
indet := unimons[1];
for t2 in TermsOfPolynomial(ideal[i] - t[1]*t[2]) do
if IsOne(DenominatorOfRationalFunction(t2[1] / indet)) then
indet := false;
break;
fi;
od;
if indet <> false then
# We're OK - no other term in this relation involves this
# indeterminate, so we can go on to remove it
break;
fi;
fi;
od;
# Have we found a solitary indeterminate?
if indet <> false then
relation := Remove(ideal, i);
Add(removedrelations, relation);
break;
fi;
od;
if indet <> false then
# Remove this indeterminate from the indets list and put it
# onto the removedindets list instead
Remove(indets, Position(indets, indet));
Add(removedindets, indet);
# And create a new (temporary) order which has indet as the most important
order := MonomialLexOrdering(Concatenation([indet], indets));
# And make sure that our relation has a unit coefficient
relation := relation / LeadingCoefficientOfPolynomial(relation, order);
# Now reduce all the other relations in the ideal with this one, which
# will get rid of this indeterminate
i := 1;
len := Length(ideal);
while i <= len do
poly := PolynomialReducedRemainder(ideal[i], [relation], order);
if IsZero(poly) then
ideal[i] := ideal[len];
Unbind(ideal[len]);
len := len - 1;
else
ideal[i] := poly;
i := i + 1;
fi;
od;
fi;
until indet = false;
# Tidy up the ideal
ideal := ReduceIdeal(ideal, MonomialLexOrdering());
# Now make the new ring and convert the indeterminates
newring := PolynomialRing(CoefficientsRing(ring), Length(indets),
Concatenation(avoid, IndeterminatesOfPolynomialRing(ring)));
newrelations := HAPPRIME_SwitchPolynomialIndeterminates(
PolynomialRing(CoefficientsRing(ring), indets), newring, ideal);
# And remember the mapping
map := [indets, ShallowCopy(IndeterminatesOfPolynomialRing(newring))];
if not IsEmpty(removedrelations) then
# Finally, sort out the map between the old and new indeterminates
# for the removed ones
# Start off by adding relations that tell us what the new indeterminates are
ideal := IndeterminatesOfPolynomialRing(newring) - indets;
# Now add the relations that we have removed
Append(ideal, removedrelations);
# Create an ordering that puts the removed relations largest so that
# they will be removed as much as possible, and the new relations smallest
# so they will be kept, then reduce this set of relations
ord := MonomialLexOrdering(Concatenation(removedindets, indets,
IndeterminatesOfPolynomialRing(newring)));
ideal := ReduceIdeal(ideal, ord);
# Now add to the map. The polynomials in the ideal that have
# leading monomials which still involve the removedindets are the ones
# we want
removedring := PolynomialRing(CoefficientsRing(ring), removedindets);
for i in ideal do
if LeadingMonomialOfPolynomial(i, ord) in removedring then
# is this a single indeterminate (i.e. a relation of the form x_i)?
if Length(TermsOfPolynomial(i)) = 1 then
poly := [[i], [Zero(i)]];
else
poly := [[], []];
for t in TermsOfPolynomial(i) do
if t[1] in ring then
Add(poly[1], t[1]*t[2]);
elif t[1] in newring then
Add(poly[2], t[1]*t[2]);
else
Error("Relation is not seperable when creating map. Please consult the package maintainer.");
fi;
od;
if IsEmpty(poly[2]) then
Error("Unexpected relation with no new indeterminates. Please consult the package maintainer.");
fi;
fi;
Add(map[1], Sum(poly[1]));
Add(map[2], Sum(poly[2]));
fi;
od;
fi;
return [newring, newrelations, map];
end
);
#####################################################################
#################################
## This functions new to GAP 4.4.10
## so needs defining if using an earlier version
if not IsBound(EmptyPlist) then
EmptyPlist := function(n)
return [];
end;
fi;
#################################
#####################################################################
## <#GAPDoc Label="HAPPRIME_SwitchPolynomialIndeterminates_manDTPolynomialInt">
## <ManSection>
## <Func Name="HAPPRIME_SwitchPolynomialIndeterminates" Arg="R, S, poly"/>
##
## <Returns>
## Polynomial or List of Polynomials
## </Returns>
## <Description>
## Changes the indeterminates in <A>poly</A>, which should be a polynomial or
## a list of polynomials, substituting the indeterminates of the polynomial
## ring <A>S</A> one-for-one for those in <A>R</A> (from which all polynomials
## in <A>poly</A> must come). The returned object is either a polynomial or a
## list of polynomials in the new indeterminates, depending on the input object.
## <P/>
## See <Ref Func="HAPPRIME_MapPolynomialIndeterminates"/> for a function that
## can work with more general indeterminate maps.
## </Description>
## </ManSection>
## <#/GAPDoc>
#####################################################################
InstallGlobalFunction(HAPPRIME_SwitchPolynomialIndeterminates,
function(R, S, polys)
local onepoly, Rindetnums, Sindetnums, newpolys, fam, p, extrep, extrepnew,
i, j, e;
if CoefficientsRing(R) <> CoefficientsRing(S) then
Error("<R> and <S> must have the same coefficient ring");
fi;
onepoly := false;
if Length(IndeterminatesOfPolynomialRing(R)) <>
Length(IndeterminatesOfPolynomialRing(S)) then
Error("<R> and <S> must have the same number of indeterminates");
fi;
if IsEmpty(polys) then
return polys;
fi;
if IsPolynomial(polys) then
polys := [polys];
onepoly := true;
elif not IsHomogeneousList(polys) then
Error("<polys> must be a polynomial or a list of polynomials in <R>");
fi;
Rindetnums := List(IndeterminatesOfPolynomialRing(R),
IndeterminateNumberOfUnivariateRationalFunction);
Sindetnums := List(IndeterminatesOfPolynomialRing(S),
IndeterminateNumberOfUnivariateRationalFunction);
newpolys := EmptyPlist(Length(polys));
fam := FamilyObj(polys[1]);
for p in polys do
if not p in R then
Error("<polys> must be a polynomial or a list of polynomials in <R>");
fi;
extrep := ExtRepPolynomialRatFun(p);
extrepnew := EmptyPlist(Length(extrep));
i := 1;
repeat
e := ShallowCopy(extrep[i]);
if not IsEmpty(e) then
j := 1;
repeat
# Swap the indeterminate number
e[j] := Sindetnums[Position(Rindetnums, e[j])];
j := j+2;
until j > Length(e);
fi;
Add(extrepnew, e);
Add(extrepnew, extrep[i+1]);
i := i+2;
until i > Length(extrep);
Add(newpolys, PolynomialByExtRep(fam, extrepnew));
od;
# And sort out the return type
if onepoly then
return newpolys[1];
else
return newpolys;
fi;
end);
######################################################
#####################################################################
## <#GAPDoc Label="HAPPRIME_MapPolynomialIndeterminates_manDTPolynomialInt">
## <ManSection>
## <Func Name="HAPPRIME_MapPolynomialIndeterminates" Arg="old, new, poly"/>
##
## <Returns>
## Polynomial or List of Polynomials
## </Returns>
## <Description>
## Changes the indeterminates in <A>poly</A>, which can be a polynomial or a
## list of polynomials, substituting the polynomials in <A>old</A> for those
## in <A>new</A>. The returned object is either a polynomial or a list of
## polynomials in the new indeterminates, depending on the input object.
## The change of variable arguments, <A>old</A> and <A>new</A>, do not
## have to be simply indeterminates: they can be can be lists of polynomials
## which are equivalent in the two different sets of indeterminates.
## If a polynomial cannot be converted (i.e. if it cannot be generated from the
## polynomials in <A>old</A>) then <K>fail</K> is returned for that polynomial.
## </Description>
## </ManSection>
## <#/GAPDoc>
#####################################################################
InstallGlobalFunction(HAPPRIME_MapPolynomialIndeterminates,
function(old, new, polys)
local old2, new2, oldindets, newindets, ord, newpolys, p, newp, onepoly, I,
lead, newI, ord2;
onepoly := false;
if not IsHomogeneousList(old) then
Error("<old> must a list of polynomials");
fi;
if not IsHomogeneousList(new) then
Error("<new> must a list of polynomials");
fi;
if Length(old) <> Length(new) then
Error("<old> and <new> must both be the same length");
fi;
if IsPolynomial(polys) then
polys := [polys];
onepoly := true;
elif not IsHomogeneousList(polys) then
Error("<polys> must be a polynomial or a list of polynomials in the old indeterminates");
fi;
oldindets := [];
for p in old do
UniteSet(oldindets, IndeterminatesOfPolynomial(p));
od;
newindets := [];
for p in new do
UniteSet(newindets, IndeterminatesOfPolynomial(p));
od;
if not IsEmpty(Intersection(oldindets, newindets)) then
Error("the indeterminates in <old> and <new> must be distinct sets.");
fi;
# Order the indeterminates so that the old ones are larger, so will
# be replaced
ord := MonomialLexOrdering(Concatenation(oldindets, newindets));
# Are there any zeros in the old list? If so, remove them and the
# corresponding element of the new list
old2 := ShallowCopy(old);
new2 := ShallowCopy(new);
repeat
p := Position(old2, Zero(old2[1]));
if p <> fail then
Remove(old2, p);
Remove(new2, p);
fi;
until p = fail;
# The conversion is relations in an ideal. Make sure it is a GroebnerBasis
I := HAPPRIME_SingularGroebnerBasis(old2 - new2, ord);
# Extract out the relations in the new ring
newI := Filtered(I, i->IsSubset(newindets, IndeterminatesOfPolynomial(i)));
# and turn this into a Groebner Basis with grevlex ordering
## TODO I would like to make it MonomialGrevlexOrdering here,
## but a bug in that function (reported by me on 27/6/08) means
## that that ordering doesn't work!
#ord2 := MonomialGrevlexOrdering();
ord2 := MonomialGrlexOrdering();
newI := HAPPRIME_SingularGroebnerBasis(newI, ord2);
newpolys := [];
for p in polys do
# is this just a constant?
lead := LeadingMonomial(p);
if IsEmpty(lead) or lead[2] = infinity or lead[2] = 0 then
Add(newpolys, p);
else
# Reduce this. Given the ordering, this will substitute the
# new in place of the old
# We further reduce this using the grevlex ordering to get the
# simplest (smallest degree) form in the new indeterminates
newp := PolynomialReducedRemainder(
PolynomialReducedRemainder(p, I, ord), newI, ord2);
# Check that this really is in the new indeterminates - it might not
# be if the map is not complete
if not IsSubset(newindets, IndeterminatesOfPolynomial(newp)) then
Add(newpolys, fail);
else
Add(newpolys, newp);
fi;
fi;
od;
# And sort out the return type
if onepoly then
return newpolys[1];
else
return newpolys;
fi;
end);
######################################################
#####################################################################
## <#GAPDoc Label="HAPPRIME_CombineIndeterminateMaps_manDTPolynomialInt">
## <ManSection>
## <Func Name="HAPPRIME_CombineIndeterminateMaps" Arg="coeff, M, N"/>
##
## <Returns>
## List
## </Returns>
## <Description>
## Returns the indeterminate map that results from applying map <A>M</A>
## followed by map <A>N</A>. An indeterminate map is a list containing two
## lists, the first of which is a list of polynomials in the original indeterminates,
## the second the equivalent polynomials in the new ring indeterminates.
## </Description>
## </ManSection>
## <#/GAPDoc>
#####################################################################
InstallGlobalFunction(HAPPRIME_CombineIndeterminateMaps,
function(M, N)
local indetsA, indetsB, indetsC, i, ord, I, J, relations, map, poly, t;
# Recreate the three rings
indetsA := [];
indetsB := [];
indetsC := [];
for i in M[1] do
Append(indetsA, IndeterminatesOfPolynomial(i));
od;
indetsA := AsSet(indetsA);
for i in M[2] do
Append(indetsB, IndeterminatesOfPolynomial(i));
od;
indetsB := AsSet(indetsB);
for i in N[1] do
Append(indetsC, IndeterminatesOfPolynomial(i));
od;
indetsC := AsSet(indetsC);
# Check that the second ring in M and the first in N are compatible
if not IsSubset(indetsC, indetsB) and not IsSubset(indetsB, indetsC) then
Error("the maps <M> and <N> are incompatible");
fi;
indetsC := [];
for i in N[2] do
Append(indetsC, IndeterminatesOfPolynomial(i));
od;
indetsC := AsSet(indetsC);
# And check that the three rings are distinct
if not IsEmpty(Intersection(indetsA, indetsB)) then
Error("the source and target rings of <M> are not different");
fi;
if not IsEmpty(Intersection(indetsA, indetsC)) then
Error("the source ring of <M> and target ring of <N> are not different");
fi;
if not IsEmpty(Intersection(indetsB, indetsC)) then
Error("the source and target rings of <N> are not different");
fi;
# Order the middle indeterminates first so that they are removed
ord := MonomialLexOrdering(Concatenation(indetsB, indetsA, indetsC));
# And get the two set of relations
I := M[1] - M[2];
J := N[1] - N[2];
# Now reduce each with the other
relations := [];
for i in I do
Add(relations, PolynomialReducedRemainder(i, J, ord));
od;
for i in J do
Add(relations, PolynomialReducedRemainder(i, I, ord));
od;
relations := AsSet(relations);
# Now create the new map. The polynomials in the ideal have
# leading monomials which involve the source ring
map := [[], []];
for i in relations do
poly := [[Zero(i)], [Zero(i)]];
for t in TermsOfPolynomial(i) do
if IsSubset(indetsA, IndeterminatesOfPolynomial(t[1])) then
Add(poly[1], t[1]*t[2]);
elif IsSubset(indetsC, IndeterminatesOfPolynomial(t[1])) then
Add(poly[2], t[1]*t[2]);
else
# This relation still involves an old indeterminate, so we
# can't do anything with it - ignore it
poly := "ignore";
break;
fi;
od;
if poly <> "ignore" then
Add(map[1], Sum(poly[1]));
Add(map[2], Sum(poly[2]));
fi;
od;
return map;
end);
######################################################
#####################################################################
## <#GAPDoc Label="HAPPRIME_SingularGroebnerBasis_manDTGroebnerInt">
## <ManSection>
## <Func Name="HAPPRIME_SingularGroebnerBasis" Arg="pols, O"/>
## <Func Name="HAPPRIME_SingularReducedGroebnerBasis" Arg="pols, O"/>
##
## <Returns>
## List
## </Returns>
## <Description>
## Returns the Gröbner basis (or reduced Gröbner basis) with respect to the
## ordering <A>O</A> for the ideal generated by the polynomials <A>pols</A>.
## This function uses the Gröbner basis implementation
## from &singular;, for preference, if available
## (<Ref Func="GroebnerBasis" BookName="singular"/>), and if so it also
## manuipulates the result to fix a bug in &singular; where the returned
## polynomials are not necessarily returned with a value external
## representation (see
## <Ref Sect="The Defining Attributes of Rational Functions" BookName="ref"/>).
## <P/>
## If the option <C>obeyGBASIS</C> is <K>true</K>, then this function will use
## whichever algorithm is specified by the <K>GBASIS</K> global variable
## (see <Ref Var="SINGULARGBASIS" BookName="singular"/>).
## </Description>
## </ManSection>
## <#/GAPDoc>
#####################################################################
#if LoadPackage("singular") = true then
if IsPackageMarkedForLoading("singular","0") then
InstallGlobalFunction(HAPPRIME_SingularGroebnerBasis,
function(pols, O)
local gbasis, ideal, fam;
# If empty, do nothing
if IsEmpty(pols) then
return pols;
fi;
# Make sure we use Singular
gbasis := GBASIS;
if not ValueOption("obeyGBASIS") = true then
GBASIS := SINGULARGBASIS;
fi;
ideal := GroebnerBasis(pols, O);
GBASIS := gbasis;
# Hack to make sure that the polynomials returned by Singular are correct
fam := FamilyObj(pols[1]);
return List(ideal, x->PolynomialByExtRep(fam, ExtRepPolynomialRatFun(x)));;
end
);
#####################################################################
InstallGlobalFunction(HAPPRIME_SingularReducedGroebnerBasis,
function(pols , O)
local ipr, mcf, R, I, input, out, fam;
# If empty, do nothing
if IsEmpty(pols) then
return pols;
fi;
if ValueOption("obeyGBASIS") = true and GBASIS = GAPGBASIS then
return ReducedGroebnerBasis(pols, O);
fi;
if IsPolynomialRingIdeal(pols) then
R := LeftActingRingOfIdeal(pols);
pols := GeneratorsOfTwoSidedIdeal(pols);
else
R := DefaultRing(pols);
fi;
if IsMonomialOrdering(O) then
ipr := ShallowCopy(IndeterminatesOfPolynomialRing(R));
mcf := MonomialComparisonFunction(O);
Sort(ipr, mcf);
ipr := Reversed(ipr);
R := PolynomialRing(LeftActingDomain(R), ipr);
fi;
if not(HasTermOrdering(R) and
IsIdenticalObj(TermOrdering(R), O)) then
SetTermOrdering(R, O);
SingularSetBaseRing(R);
fi;
I := Ideal(R, pols);
# Now use Singular to find the Groebner Basis
Info( InfoSingular, 2, "running GroebnerBasis..." );
SingularCommand("", "option(redSB)");
# preparing the input for Singular
input := "";
Append( input, "ideal GAP_groebner = simplify( groebner( " );
Append( input, ParseGapIdealToSingIdeal( I ) );
Append( input, " ), 1 );\n" );
out := SingularCommand( input, "string (GAP_groebner)" );
SingularCommand("", "option(noredSB)");
Info( InfoSingular, 2, "done GroebnerBasis." );
I := List( SplitString( out, ',' ), ParseSingPolyToGapPoly );
# Hack to make sure that the polynomials returned by Singular are correct
fam := FamilyObj(pols[1]);
return List(I, x->PolynomialByExtRep(fam, ExtRepPolynomialRatFun(x)));;
end );
else
InstallGlobalFunction(HAPPRIME_SingularGroebnerBasis,
function(I, O)
return GroebnerBasis(I, O);
end
);
InstallGlobalFunction(HAPPRIME_SingularReducedGroebnerBasis,
function(I, O)
return ReducedGroebnerBasis(I, O);
end
);
fi;
#####################################################################
