(* written by Manuel Kauers and Burkhard Zimmermann *)

(* Research Institute for Symbolic Computation (RISC) *)
(* J. Kepler University Linz, Austria *)

(* http://www.risc.uni-linz.ac.at/research/combinat/risc/software/ *)
(* email: manuel.kauers@risc.uni-linz.ac.at *)

(* ************************************************************************* *)
(*                                     Gb.m                                  *)
(*               An Interface to Faugere's amazing C++ library Gb            *)
(* ************************************************************************* *)

(* INSTALLATION INSTRUCTIONS:                                                *)
(*                                                                           *)
(*  (1)  Install Faugere's Gb library on your computer system                *)
(*  (2)  Put this file somewhere where Mathematica finds it                  *)
(*  (3)  Modify the last line of this file according to your needs           *)
(*                                                                           *)
(*                                                                           *)
(* NOTE:                                                                     *)
(*                                                                           *)
(*  This package has been developed for Mathematica 5.                       *)
(*  It should also run in the later versions of Mathematica 4, e.g., 4.2.    *)
(*  It definitely does not run under Mathematica 3 and older versions.       *)
(*                                                                           *)
(*                                                                           *)
(* BUGS:                                                                     *)
(*                                                                           *)
(*  Please report bugs concerning the interface package to Manuel Kauers     *)
(*  (manuel.kauers@risc.uni-linz.ac.at)                                      *)
(*                                                                           *)
(*                                                                           *)
(* LICENCE:                                                                  *)
(*                                                                           *)
(*  This software package is free; 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. The package 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.                             *)
(*                                                                           *)

(* ------------------------------------------------------------------------- *)

Gb`Private`Version = "1.0 (04-02-13)"; (* Version 1.0 *)

BeginPackage["Gb`"];

(* avoid problems on reload of the package *)
Unprotect[`GBasis, `Ideal, `GbCharacteristic, `GbDegree,
    `GbDimension, `GbEliminate, `GbFGLM, `GbHilbertFunction, `GbHilbertPolynomial,
    `GbIdeal, `GbMonomialOrder, `GbNormalForm, `GbPretend, `GbTime, `GbTriangular,
    `GbVariables, `GbVersion, `GbCofactors, `GbExtendedGBasis,
    MonomialOrder, MonomialOrdering, DRL, Lex, DegreeReverseLexicographic,
    Lexicographic, Modulus, Characteristic, Eliminate, Algorithm, Serveur,
    Variable, InFileName, OutFileName, CoefficientDomain, FiniteField
];

(************************)
(* Declaration of names *)
(************************)

(* organisatorical items *)
GbHome::usage = 
  "Directory of where to locate the serveur binaries of Gb";

$PrettyPrintIdeals::usage = 
  "Indicates if ideals should be prettyprinted or not.";
$PrettyPrintIdeals = True;

(* functionality *)

GbIdeal::usage =
  "Compute an Ideal object representing the polynomial ideal defined by the argument.";

GbBash::usage =
  "Command for executing a bash script. By default \"bash\"";

GbPretend::usage =
  "Compute an Ideal object of the given set of polynomials, assuming that these build a Groebner basis. A monomial order has to be specified.";

Ideal::usage =
  "Head of expressions of type Ideal. See GbIdeal.";

GBasis::usage = 
  "GBasis[{poly1, poly2, ...}, {x1, x2, ...}] gives a list of polynomials that form a Groebner basis for the set of polynoials polyi.";

GbMonomialOrder::usage = 
  "Returns the monomial order of a given ideal";

GbCharacteristic::usage =
  "Returns the characteristic of a given ideal";

GbHilbertPolynomial::usage = 
  "GbHilbertPolynomial[{poly1, poly2, ...}, {var1, var2, ...}] computes the Hilbert polynomial for the given polynomial ideal.";

GbHilbertFunction::usage =
  "GbHilbertFunction[{poly1, poly2, ...}, {var1, var2, ...}] computes the Hilbert function for the given polynomial ideal.";

GbDimension::usage = 
  "GbDimension[{poly1, poly2, ...}] computes the dimension of the given polynomial ideal.";

GbDegree::usage = 
  "GbDegree[{poly1, poly2, ...}] computes the degree of the given polynomial ideal.";

GbFGLM::usage = 
  "Transforms a Groebner basis in DRL order into Lex order";

GbTriangular::usage =
  "Compute a triangular system of a zero dimensional ideal.";

GbNormalForm::usage = 
  "Computes the normal form of a polynomial wrt. a given ideal.";

GbVersion::usage = 
  "Prints out the version of the underlying Gb and exists";

GbEliminate::usage = 
 "GbEliminate[{poly1, poly2, ...}, {x1, x2, ...}, {y1, y2, ...}] computes the elimination ideal of the given ideal in which x1, x2, ... are eliminated.";

GbTime::usage = 
 "Returns the time which was needed to construct the given ideal object.";

GbMonomialOrder::usage = 
 "Returns the monomial order of a given ideal";

GbCharacteristic::usage =
 "Returns the characteristic of a given ideal";

GbVariables::usage =
 "Returns the variables of a given ideal";

(* options and their values *)
MonomialOrder::usage = 
  "Option of GBasis, might be DRL or Lex. See also option Variables.";
MonomialOrdering::usage = 
  "Synonym for MonomialOrder";
DRL::usage = 
  "Possible value to the MonomialOrder option of GBasis. Specifies the degree reverse lexicographic ordering. See also Lex.";
Lex::usage = 
  "Possible value to the MonomialOrder option of GBasis. Specifies the lexicographic ordering. See also DRL.";
DegreeReverseLexicographic::usage =
  "Synonym for DRL";
Lexicographic::usage = 
  "Synonym for Lex";
Modulus::usage = 
  "Option of GBasis, might be any prime number <65522, or zero";
Characteristic::usage = "Synonym for Modulus.";
Eliminate::usage = 
  "Option of GBasis, might be True or False. This option is only applicable
if the MonomialOrder is an elimination ordering. See option Variables for
further information.";
Algorithm::usage = "Option to GBasis. Specifies the Gb serveur to be used.";
Serveur::usage = "Option to GbVersion. Specifies the program used to extract the version info. See also: Algorithm.";
Variable::usage = "Option to GbHilbertPolynomial. Specifies the variable in which to express the polynomial.";
InFileName::usage = "Option to ComputeIdealInformation. Specifies the file name of the input file to be created for Gb. (Defaults to a tmp file.)"; 
OutFileName::usage = "Option to ComputeIdealInformation. Specifies the file name of the output file to which Gb should write its results. (Defaults to a tmp file.)"; 
CoefficientDomain::usage = "Option to GBasis. Might be one of Rationals or FiniteField or Automatic. The latter determines the coefficient field by the given characteristics.";
FiniteField::usage = "Possible value for option CoefficientDomain in GBasis.";

(***************************)
(* Declaration of warnings *)
(***************************)

Gb::newvars = "New variables detected: `1`";
On[Gb::newvars];

Gb::noorder = "No monomial ordering specified. Taking `1`."; 
On[Gb::noorder];

Gb::novars = "No variable order specified.";
Off[Gb::novars]; (* because Gb::newvars is On *)

Gb::nochar = "No characteristic specified. Taking `1`.";
Off[Gb::nochar];

Gb::noalgo = "No algorithm specified. Taking `1`.";
Off[Gb::noalgo];

Gb::bogusopt = "Unknown Option(s) encountered: `1`";
On[Gb::bogusopt];

(*********************************)
(* Declaration of error messages *)
(*********************************)

Gb::boguselim = 
  "Eliminate -> True  makes only sense on elimination orderings.";
On[Gb::boguselim];

Gb::bogusfglm = "You cannot specify an elimination ordering for FGLM.";
On[Gb::bogusfglm];

Gb::boguschar = "Modulus `1` illegal. Prime <65522 or 0 expected.";
On[Gb::boguschar];

Gb::bogusord = "Termorder unknown or no termorder specified.";
On[Gb::bogusord];

Gb::timeout = "Computation aborted.";
On[Gb::timeout];

Gb::bogusalgo = "Algorithm `1` unknown or binary file not found. (Expecting, e.g., \"serveur\", \"serveur_t\", or \"accel_serveur\")";
On[Gb::bogusalgo];

Gb::bogusfield = "Field `1` not supported or incompatible to specified modulus `2`.";
On[Gb::bogusfield];

Gb::nyi = "Not yet implemented: `1`";
On[Gb::nyi];

Gb::bogustiming = "Timing specification `1` not understood. Expecting a positive integer or Infinity.";
On[Gb::bogustiming];

(**********************)
(* Start private part *)
(**********************)

Begin["`Private`"]; 

(* Definition of private identifiers. *)
(* Further private identifiers may be declared private elsewhere. *)

{`Pretend, `Variables, `DRL2Lex, `FGLM};

DRL2Lex::usage = 
  "Possible value to the MonomialOrder option of GBasis. Requests the conversion from a DRL Groebner basis to a Lex Groebner basis. The ideal has to be zero dimensional.";
FGLM::usage =
  "Synonym for DRL2Lex. See also: GbFGLM.";
Pretend::usage = "True or False. (private) Option to ComputeIdealInformation.";
Variables::usage = 
  "(Private) Option of GBasis, specifies the (order of the) variables.\nThe variables may be specified either as a list {x1,x2,...,xn} or as a list of (exactly) two lists {{x1,...xk},{y1,...,yn}}. The latter defines an elimination ordering with the first list of variables to be eliminated. The option Eliminate is applicable iff the variables are specified in the second form.\nFor reasons of compatibility to Mathematica's builtin GroebnerBasis, we have:\n (1) GBasis[basis, Variables -> list, ...] is equivalent to GBasis[basis, list, ...]\n (2) GBasis[basis, Variables -> {list1, list2}, ...] is equivalent to GBasis[basis, list1, list2, ...] See also option MonomialOrder.";

(****************************** Tools *****************************)

(**
 * ExtractVariables[expr]
 *
 * Returns a list of all variables in the expression. 
 * Every subexpression which is not a polynomial, is treated as a variable.
 * This is a private version of MMA's Variables. 
 *
 * Example:  ExtractVariables[ f[x]^2 + y*z - 17 ] == { f[x], y, z }
 * See also: IsPolynomial
 **)
ExtractVariables[expr_] := 
   Which[
     MemberQ[{Integer, Rational}, Head[expr]], (* numbers *)
       {},
     MemberQ[{Times, Plus}, Head[expr]], (* Plus and Times *)
       Apply[Union, Map[ExtractVariables, Level[expr, 1]]],
     MatchQ[expr, _^_Integer] && expr[[2]] > 1, (* poly ^ Integer *)
       ExtractVariables[expr[[1]]],
     True, {expr} (* not a polynomial: a variable *)
   ];

(**
 * IsPolynomial[expr]
 *
 * Checks if an expression is a multivariate polynomial. 
 * Plain symbols are polynomials, sums, products, rational multiples, 
 * and positive integer powers of polynomials are polynomials. 
 * Nothing else is a polynomial. 
 *
 * Example: IsPolynomial[ x*((y+1)^3)^3 ] == True
 *          IsPolynomial[ x*Sqrt[3] ] == False
 * Options: 
 *       -- Variables -> f
 *          f[expr] evaluates to True if expr should be considered
 *          a variable
 *          Default: Head[#]===Symbol&
 **)
Options[IsPolynomial] = {Variables -> (Head[#]===Symbol&)};
IsPolynomial[expr_, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {check},

   check = Variables /. {opts} 
                      /. Options[IsPolynomial];

      check[expr] (* base I: variables *)
   || MemberQ[{Integer, Rational}, Head[expr]] (* base II: constants *)
   || (MemberQ[{Times, Plus}, Head[expr]] (* poly + poly, poly * poly *)
       && Apply[And, Map[IsPolynomial[#, opts]&, Level[expr, 1]]])
   || ((Head[expr] === Power) (* poly ^ integer *)
       && (Head[expr[[2]]] === Integer) && (expr[[2]] > 0)
       && IsPolynomial[expr[[1]], opts])
];

(**
 * CleanVariables[basis]
 *
 * Computes two lists of transformation rules for translating the variables in 
 * the given basis into the variables x1, x2, x3, ... and back. 
 * ExtractVariables is used for extracting the variables from the basis.
 * The variables x1, x2, ... belong to the context Gb`Private`.
 * 
 * Example: CleanVariables[ {x*y, x+4/3} ] == { {x->x1, y->x2}, {x1->x, x2->y} }
 * Options: 
 *       -- Variables -> f
 *          A warning is printed for all detected variables x in basis for which 
 *          f[x] does not evaluate to True
 *          Default: True&
 **)
Options[CleanVariables] = {Variables -> (True&)};
CleanVariables[basis_, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {i, s, vars, back, forth, check, illvars},

   (* extract variables *)
   vars = Apply[Union, Map[ExtractVariables, basis]];

   (* complain about undeclared variables *)
   check = Variables /. {opts} /. Options[CleanVariables];
   illvars = Select[vars, !check[#]&];
   If [Length[illvars] > 0, Message[Gb::newvars, illvars]];

   (* construct transformation rules *)
   s[i_] := Symbol["Gb`Private`x" <> ToString[i]];
   forth = Table[vars[[i]] -> s[i], {i, 1, Length[vars]}];
   back = Table[s[i] -> vars[[i]], {i, 1, Length[vars]}];

   (* done *)
   {forth, back}   
];

(**
 * FileExists[name]
 *
 * Checks whether a given string corresponds to an existing file.
 * Returns True or False, accordingly.
 **)
FileExistsQ[name_String] := FileInformation[name] =!= {};

(**
 * MakeString[ex1, ex2, ex3,...]
 *
 * Produces the string obtained by concatenation of the arguments. 
 * Arguments which are not strings are converted into strings.
 * The substrings "Gb`Private`" and "Gb`" will be deleted from the result.
 **)
MakeString[args___] := 
   Module[ {l, ex}, 

   StringReplace[
     StringReplace[
       Apply[StringJoin,
         Map[ If[Head[#]=!=String, ToString[#, FormatType -> InputForm], #]&,
         {args}
         ] (* end map *)
       ], (* end apply *)
       "Gb`Private`" -> ""],
     "Gb`" -> ""]
];

(**
 * WriteInputFile[{data}]
 *
 * Writes data into an input file in the Gb syntax.
 * The option InFileName may specify a file name (as String). 
 * If the option does not point to a string, a fresh tmp file will be used.
 * 
 * The function returns the full name of the file that has been written.
 *
 * The data of the file is specified as a list of transformation rules whos
 * left hand sides correspond to the "#..." tags in Faugere's syntax as follows
 *
 *   Variables -> {x1, x2, ...}          ==>    #vars [x1, x2, ...]
 *   Characteristic -> 0                 ==>    #coeffring rational
 *   Characteristic -> prime             ==>    #coeffring modulo prime
 *   Properties -> {"p1", "p2", ...}     ==>    #properties [p1, p2, ...]
 *   List -> {{p1,p2,...}, list2, ...}   ==>    #list [p1, p2,...] list2...
 *   Error -> msg_String                 ==>    #error msg
 *
 * a single list "List -> {{p1...}}" may also be written simply "List -> {p1...}"
 *
 * Options:
 *      -- InFileName -> String specifying the desired filename, or any non-string
 *                 expression for using a tmp file. 
 * 
 * Unknown option specifications are ignored. 
 **)
Options[`WriteInputFile] = {InFileName -> Automatic};
WriteInputFile[{data:((_Rule|_RuleDelayed)...)}, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {a, in, inFileName, vars, char, props, list, i, error},

   inFileName = InFileName /. {opts} /. Options[WriteInputFile];
   vars = Variables /. {data};
   char = Characteristic /. {data};
   props = Properties /. {data};
   list = List /. {data};

   (* wrap a list around the list if there is only one *)
   If[ MatchQ[list, {___}] && !MatchQ[list, {{___}..}], 
       list = {list}
   ];

   (* create a tmp file if no file name was specified *)
   in = If[ Head[inFileName] === String, OpenWrite[inFileName], OpenTemporary[] ];

   (* write variable list *)
   If[ vars =!= Variables, 
       WriteString[in, MakeString["#", Apply[Symbol["Gb`Private`vars"], vars], "\n"]]
   ];

   (* write characteristics *)
   Which [ 
      char === 0,
         WriteString[in, "#coeffring rational\n"],
      IntegerQ[char] && char > 0 && char < 65522,
         WriteString[in, "#coeffring modulo " <> ToString[char] <> "\n"]
   ];

   (* write properties *)
   If[ props =!= Properties,
      WriteString[in, 
         MakeString["#properties[", 
            Apply[MakeString, 
              Drop[props /. a_String :> Sequence[a, ", "], -1]
            ],
            "]\n"
         ]
      ]
   ];

   (* write the data list *)
   If[ list =!= List, 
      WriteString[in, "#list\n"];
      Do[
         WriteString[in, 
            StringReplace[MakeString[Apply[`tmp, list[[i]] ], "\n"], "tmp[" -> "["]
         ],
         {i, 1, Length[list]}
      ]
   ];

   Close[in];

   in[[1]] (* return the file name *)
];

(**
 * ServeurCall[serveur]
 * 
 * Calls the specified serveur. The serveur string should be such that
 * GbHome <> "/" <> serveur   yields an executable binary file.
 *
 * The function returns the file name of the file to which the output was
 * piped. If an error occurs, then the value $Failed is returned instead.
 * An output file may or may not exist in this case. 
 *
 * Options:
 *      -- InFileName -> name of the input file to parse to serveur with "-read"
 *                If this option is not a string, then no input will be read.
 *      -- OutFileName -> name of the file to which to pipe the output of the serveur
 *                If this option is not a string, then a tmp file will be used.
 *      -- TimeConstraint -> positive integer or Infinity for specifying timeout.
 *      -- Flags -> list of strings specifying the desired flags 
 *
 * Additional option specs are ignored. 
 **)
Options[`ServeurCall] = {
   InFileName -> Automatic,
   OutFileName -> Automatic,
   TimeConstraint -> Infinity,
   `Flags -> {}
};
ServeurCall[serveur_String, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {command, inFileName, outFileName, flags, timeout, out, tmp, result},

   (* extract options *)
   inFileName = InFileName /. {opts} /. Options[ServeurCall];
   outFileName = OutFileName /. {opts} /. Options[ServeurCall];
   flags = Flags /. {opts} /. Options[ServeurCall];
   timeout = TimeConstraint /. {opts} /. Options[ServeurCall];

   (* command = path + serveur, for the moment *)
   command = If[ Head[GbHome] === String, 
                 GbHome <> If[StringTake[GbHome, -1] == "/", 
                              "", 
                              "/"],
                 "./"] <> serveur;

   (* check the existence of the serveur *)
   If [ !FileExistsQ[command],

        (* if this doesn't exist, 
            try "__INT" instead of "__GINT"
            and/or "__BigLexp" instead of "__Lexp" *)

        Which[
          FileExistsQ[StringReplace[command, "__GINT" -> "__INT"]],
             command = StringReplace[command, "__GINT" -> "__INT"],
          FileExistsQ[StringReplace[command, "__Lexp" -> "__BigLexp"]],
             command = StringReplace[command, "__GINT" -> "__INT"],
          FileExistsQ[StringReplace[command, {"__Lexp" -> "__BigLexp", 
                                              "__GINT" -> "__INT"}]],
             command = StringReplace[command, 
                {"__Lexp" -> "__BigLexp", "__GINT" -> "__INT"}],
          True, (* otherwise give up: no serveur found *)
             Message[Gb::bogusalgo, serveur];
             Return[$Failed];
        ];
   ];

   (* prepare a tmp file for the output, if neccessary *)
   If[ Head[outFileName] =!= String,
       (* then: no outFileName given, use name of a tmp file *)
       out = OpenTemporary[];
       outFileName = out[[1]];
       Close[out]
   ];

   (* convert flag list into flag string *)
   flags = Apply[StringJoin, flags /. {a_String -> Sequence[a, " "]}];

   (* compose command *)
   command = command <> (* executable file name *)
             " " <> flags <>
             If[ Head[inFileName] === String, " -read " <> inFileName, "" ] <>
             " > " <> outFileName;

   (* run the command *)
   Which[ 
     timeout === Infinity, (* no time constraint *)
       Run[command],
     IntegerQ[timeout] && timeout > 0, (* time constraint *)
       (* execute command within wrapping watchdog bash script *)
       tmp = OpenTemporary[];
       WriteString[tmp, "(" <> command <> ") & cmdpid=$!\n"];
       WriteString[tmp, "(sleep " <> ToString[timeout] <> "; kill -9 $cmdpid ) &\n"];
       WriteString[tmp, "watchdogpid=$!\n"];
       WriteString[tmp, "wait $cmdpid 2>&1\n"];
       WriteString[tmp, "kill $watchdogpid > /dev/null 2>&1 \n"];
       Close[tmp];
       result = Run[GbBash <> " " <> tmp[[1]] ]; 
       (* if the result is different from zero, then we assume a timeout and 
          return $Failed *)
       If[ result =!= 0, 
          Message[Gb::timeout];
          Return[$Failed] 
       ],
     True, (* bogus time constraint spec. *)
       Message[Gb::bogustiming, timeout];
       Return[$Failed]
   ];

   outFileName (* return the file name of the output file *)
];

(**
 * FreeOfErrorsQ[list]
 * 
 * Checks if any of the lines in the given output represents an error. 
 * In this case, it reports the error by Print and returns False. 
 * Otherwise it returns True.
 *
 * The list is assumed to contain the lines of the output file as strings, 
 * according to the result of Import[filename, "Lines"].
 **)
FreeOfErrorsQ[file_List] := 
   Module[{errors},

   (* select error messages *)
   errors = Select[file, (Head[#]===String && Length[#] > 5 && StringTake[#, 6] == "#error")&];
   
   If[ Length[errors] === 0,
      True, (* everything ok. *)
      Print["Gb reported the following error(s):"]; (* nothing ok. *)
      Print /@ errors; 
      False
   ]   
];

(**************************** Interface ***************************)

(**
 * GbHome
 *
 * Path where to find the serveur binaries. 
 * GbHome <> "/filename" should evaluate to an existing file.
 * Default: ".";
 **)
If[ Head[GbHome] =!= String, GbHome = "."]; 

(**
 * GbBash
 *
 * command line prefix for the execution of bash scripts, including the full path
 * of the bash interpreter if necessary.
 * Default: "bash", another useful setting might be "/usr/bin/bash".
 **)
If[ Head[GbBash] =!= String, GbBash = "bash"]; 

(**
 * ComputeIdealInformation[basis]
 * 
 * Computes certain informations about the polynomial ideal defined by the given basis.
 * The function outputs a list of transaction rules with the following left hand sides.
 * The transformation rules specified by (#) may not appear in every output.
 * This function generalizes GBasis, GbDimension, GbHilbertPolynomial, GbEliminate,
 * and GbDegree. 
 * 
 * Output:
 *       -- GBasis -> a Groebner basis for for the ideal
 *       -- MonomialOrder -> (copied from input, see below)
 *       -- Characteristic -> (copied from input, see below)
 *       -- Variables -> (copied from input, see below)
 *       -- Dimension (#) -> the dimension of the ideal 
 *       -- Degree (#) -> the degree of the ideal 
 *       -- Time (#) -> number of seconds spent by Gb 
 *       -- HilbertPolynomial[1] (#) -> the first Hilbert polynomial in T
 *       -- HilbertPolynomial[2] (#) -> the second Hilbert polynomial in T
 *
 * Input Options: 
 *       -- MonomialOrder -> either Lex or DRL or DRL2Lex
 *          o  Lex : lexicographic ordering
 *          o  DRL : degree reverse lexicographic ordering
 *          o  `DRL2Lex : applies the FGLM algorithm to the basis. The basis is
 *                assumed to be a GB wrt. DRL ordering. This "ordering" is incompatible
 *                with the {{},{}} variable specification (see below)
 *          o  Lexicographic : synonym for Lex
 *          o  DegreeReverseLexicographic : synonym for DRL
 *          o  `FGLM : synonym for DRL2Lex
 *       -- MonomialOrdering : synonym for MonomialOrder
 *       -- `Variables -> list or {list1, list2}
 *          o  list : list of expressions to be treated as variables
 *          o  {list1, list2} : Join[list1, list2] is the list of expressions
 *                to be treated as variables. This defines an elimination 
 *                ordering where the variables in list1 are to be eliminated.
 *       -- Modulus -> 1 < p < 65522, p prime, or 0
 *       -- Characteristic : synonym for Modulus
 *       -- CoefficientDomain -> either Automatic or FiniteField or Rationals
 *          o  Automatic: determine field from given characteristic
 *          o  Rationals: yields error if characteristic =!= 0
 *          o  FiniteField: yields error if characteristic === 0
 *          o  everything else: yields error.
 *       -- Eliminate -> either True or False
 *          Only applicable for elimination orderings.
 *          If True, polynomials containing variables from list1 (see Variables)
 *          will be omitted in the result. 
 *       -- Algorithm -> String prefix of the serveur to be used.
 *       -- InFileName -> a string specifying the file into which the input for Gb
 *          will be written. If this option is not given, a tmp file will be used
 *       -- OutFileName -> a string specifying the file into which Gb will be asked
 *          to write its output. If this option is not given, a tmp file will be used.
 *       -- TimeContraint -> a positive integer number (or Infinity) specifying
 *                the maximum number of seconds to be spent for the computation.
 **)

Options[Gb] = {

    (* Modulus and MonomialOrder intensionally not specified *)
    Characteristic :> (Message[Gb::nochar, 0]; 0),
    MonomialOrdering :> (Message[Gb::noorder, DRL]; DRL), 
    Variables :> (Message[Gb::novars]; {}),
    Algorithm :> (Message[Gb::noalgo, "serveur"]; "serveur"),
    Eliminate -> False,
    CoefficientDomain -> Automatic,
    TimeConstraint -> Infinity,
    Pretend -> False,
    InFileName -> Automatic,
    OutFileName -> Automatic
};

ComputeIdealInformation[basis_List, opts:((_Rule|_RuleDelayed)...)] := 
   Module[ {vars1, vars2, ord, eord, algo, char, field, elim, pretend,
            inFileName, outFileName, unknownOptions, variableTranslator, 
            timeout, change, serveur, in, out, result, i, h, p},

   (* extract information of options *)
   elim = Eliminate /. {opts} /. Options[Gb];
   vars1 = Variables /. {opts} /. Options[Gb];
   vars2 = {}; (* see below *)
   
   ord = MonomialOrder /. {opts} /. Options[Gb]
                       /. MonomialOrder -> MonomialOrdering 
                       /. {opts} /. Options[Gb]
                       /. {DegreeReverseLexicographic -> DRL,
                           Lexicographic -> Lex,
                           FGLM -> DRL2Lex};
   char = Modulus /. {opts} /. Options[Gb]
                  /. Modulus -> Characteristic
                  /. {opts} /. Options[Gb];
   algo = Algorithm /. {opts} /. Options[Gb];
   field = CoefficientDomain /. {opts} /. Options[Gb];
   timeout = TimeConstraint /. {opts} /. Options[Gb];
   pretend = Pretend /. {opts} /. Options[Gb];

   (* order change requested? *)
   If[ change = (ord === DRL2Lex),
       ord = DRL; (* because that's the serveur we want to call *)
   ];   

   eord = False; (* are we working with an elimination order? *)
   Which[ 
     MatchQ[vars1, {_List, _List}], (* if elimination order *)
       eord = True;
       If[ change, Message[Gb::bogusfglm]; Return[$Failed] ];
       vars2 = vars1[[2]];
       vars1 = vars1[[1]],
     elim === True,                 (* elif elimination switched on *)
       Message[Gb::boguselim];
       Return[$Failed];
   ];

   (* check options: unused/unknown lhs? *)
   unknownOptions = 
       Complement[Join[Cases[{opts}, (a_ -> b_) -> a],
                       Cases[{opts}, (a_ :> b_) -> a]],
                  {Eliminate, Variables, MonomialOrder, MonomialOrdering,
                   Modulus, Characteristic, Algorithm, CoefficientDomain,
                   Pretend, TimeConstraint, InFileName, OutFileName}];

   If[ unknownOptions =!= {}, 

       Message[Gb::bogusopt, unknownOptions];
   ];

   (* normalize variables *)
   variableTranslator = 
       CleanVariables[Join[basis, vars1, vars2], Variables -> (MemberQ[vars1,#] || MemberQ[vars2,#]&)];

   (* add unknown variables to vars2 *)
   vars2 = Join[vars2, 
                Complement[Flatten[Map[ExtractVariables, basis]], vars1, vars2]
               ];

   (* is the characteristic meaningful? *)
   If[ (char =!= 0 && !PrimeQ[char]) || char < 0 || char >= 65522, 

       Message[Gb::boguschar, char]; 
       Return[$Failed];
   ];

   (* check coefficient field option *)
   (* this option doesn't give additional information. *)
   (* its values are only checked for consistency *)
   If[
     (field === Rationals && char =!= 0) ||
     (field === FiniteField && char === 0) ||
     (!MemberQ[{Automatic, Rationals, FiniteField}, field]),
 
       Message[Gb::bogusfield, field, char];
       Return[$Failed];
   ];

   (* is the monomial ordering known? *)
   If[ !MemberQ[{Lex, DRL}, ord], 

       Message[Gb::bogusord]; 
       Return[$Failed];
   ];

   (* write the Gb input file *)
   inFileName = WriteInputFile[
      {Variables -> (Join[vars1, vars2] /. variableTranslator[[1]]),
       Characteristic -> char,
       Properties -> 
          If[change, 
             {"GB", "Dexp(" <> ToString[Length[vars1] + Length[vars2]] <> ")"}, 
             {"SYSTEM"}
          ],
       List -> Expand[basis /. variableTranslator[[1]] ]
      }, opts];       

   (* select the appropriate serveur *)
   serveur = algo <> "__DMP__" <> 
             Which[eord && ord === DRL, "DRLDRL",
                   ord === DRL, "Dexp",
                   ord === Lex, "Lexp"] <> "__" <>
             If[char === 0, "GINT", "GF"];

   (* call the serveur *)
   outFileName = ServeurCall[
      serveur,
      Flags -> {If[ eord, " -block " <> ToString[Length[vars1]], ""],
                If[ elim, " -elim ", ""],
                If[ pretend, " -pretend ", ""],
                If[ change, " -change", ""],
                If[ char > 0, " -prime " <> ToString[char], ""]},
      InFileName -> inFileName,
      opts (* possibly includes spec of outFileName and TimeConstraint *)
   ];

   If[ outFileName === $Failed, Return[$Failed] ];

   (* read the output file and extract the information *)
   out = Import[outFileName, "Lines"]; (* list of strings *)
   If[ !FreeOfErrorsQ[out], Return[$Failed] ];

   (* list of transformation rules *)
   result = {MonomialOrder -> If[change, Lex, ord],
             Characteristic -> char,
             Variables -> Which[
                            elim, vars2, 
                            vars2 === {}, vars1, 
                            True, {vars1, vars2}
                          ]};

   (* ------------------------------------------------------------------- *)
   (* We make the following assumptions upon the form of the output file: *)
   (* ------------------------------------------------------------------- *)
   (* (1) the first line consisting only of numbers and the symbols       *)
   (*     +,-,*,^,T declares the first Hilbert polynomial                 *)
   (* (2) the second line of this shape declares the second Hilbert       *)
   (*     polynomial.                                                     *)
   (* (3) A line starting with "Dim : " contains the dimen of the ideal   *)
   (* (4) A line starting with "Deg : " contains the degree of the ideal  *)
   (* (5) A line starting with "Output: " contains the polynomial list    *)
   (*     of the resulting Groebner basis.                                *)
   (* (6) A line starting with "#C Time" contains the number of seconds   *)
   (*     needed to complete the computation, and has the form            *)
   (*     #C Time ** sec **                                               *)
   (*     where the first "**" is a floating point number in MMA syntax   *)
   (*     and the second "**" is arbitrary.                               *)
   (* (7) There are no errors in the output file. (Checking for errors    *)
   (*     has been done some lines ago already.)                          *)
   (* (8) If there is no line starting with "Output:" but there is a line *)
   (*     starting with "[", then this line contains the GBasis result.   *)
   (* (9) If there is no line starting with "#C Time" but there is a line *)
   (*     starting with "End of FGLM" then this line contains the timing  *)
   (*     and is of the form                                              *)
   (*     "End of FGLM ... (... sec)!" where the ... are free of brackets *)
   (* ------------------------------------------------------------------- *)

   h = 1; (* Id of the next hilbert polynomial to be read *)
   Do[

     Which[

       (* is the line a polynomial in T with rational coefficients? *)
       Complement[ (* rule out "0 something" *)
         Characters[out[[i]]], 
         Join[CharacterRange["0", "9"], {"+", "*", "-", "/", "^", "T"}]
       ] === {} && 
       SyntaxQ[out[[i]]] && 
       IsPolynomial[ToExpression[out[[i]]], Variables -> (#===Symbol["T"]&)],
           (* if yes, then it is one of the Hilbert polynomials *)
           AppendTo[result, HilbertPolynomial[h++] -> ToExpression[out[[i]]]],

       StringMatchQ[out[[i]], "Dim :*"], 
           (* read the dimension *)
           AppendTo[result, Dimension -> ToExpression[StringDrop[out[[i]], 5]]],

       StringMatchQ[out[[i]], "Deg :*"],
           (* read the degree *)
           AppendTo[result, Degree -> ToExpression[StringDrop[out[[i]], 5]]],

       StringMatchQ[out[[i]], "Output:*"],
           (* read the Groebner basis *)
           AppendTo[result, 
              GBasis -> ((
                 StringReplace[ 
                    "List" <> StringDrop[out[[i]], 7], (* convert to MMA syntax *)
                    "x" -> "Gb`Private`x" (* localize variables *)
                 ]
                 // ToExpression) (* parse *)
                 /. variableTranslator[[2]]) (* translate to user's variables *)
            ], 

       StringMatchQ[out[[i]], "#C Time*"],
           (* read the time needed *)
           p = StringPosition[out[[i]], "sec"][[1, 1]]; (* end index *)
           AppendTo[result, 
                    Time -> ToExpression[StringTake[out[[i]], {8, p-1}]] * Second],

       StringMatchQ[out[[i]], "#error*"],
           (* Gb error, give up *)
           Print["Gb reported the following error(s)"];
           Print[StringDrop[out[[i]], 7]];
           Return[$Failed];
     ],
     {i, 1, Length[out]} (* for all lines *)
   ];

   (* if no GBasis result has been found, go through again and now look for lines
      starting on "[" *)
   If[ !MemberQ[Cases[result, (a_ -> _) -> a], GBasis],
     Do[
       If[StringMatchQ[out[[i]], "[*"] && SyntaxQ["List" <> out[[i]]],
         AppendTo[result,
            GBasis -> ((
              StringReplace[
                "List" <> out[[i]], "x" -> "Gb`Private`x"
              ] // ToExpression) /. variableTranslator[[2]])]
       ],
       {i, 1, Length[out]}
     ];
   ];

   (* if no Timing result has been found, go through again and now look for lines
      starting on "End of FGLM" *)
   If[ !MemberQ[Cases[result, (a_ -> _) -> a], Timing],
     Do[
       If[StringMatchQ[out[[i]], "End of FGLM *(* sec)*"],
         AppendTo[result,
            Time -> (
              StringReplace[
                StringTake[out[[i]], {StringPosition[out[[i]], "("][[1, 1]],
                                      StringPosition[out[[i]], ")"][[1, 1]]}],
                "sec" -> "*Second"]
              // ToExpression)
            ]
       ],
       {i, 1, Length[out]}
     ];
   ];

   (* if no Groebner basis has been found, return $Failed *)
   If[ Head[GBasis /. result] =!= List, Return[$Failed] ];

   (* return the result *)
   result
];

(********************************* The type Ideal **********************************)

(**
 * GbIdeal[basis, vars]
 *
 * Computes an object of type Ideal representing the ideal defined by the given 
 * basis.
 **)
GbIdeal[basis_List, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {m},

  m = ComputeIdealInformation[basis, opts];
  If[ m === $Failed, $Failed, Ideal[m]]
];
GbIdeal[basis_List, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {m},

  m = ComputeIdealInformation[basis, Variables -> vars, opts];
  If[ m === $Failed, $Failed, Ideal[m]]
];
GbIdeal[basis_List, vars1_List, vars2_List, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {m},

  m = ComputeIdealInformation[basis, Variables -> {vars1, vars2}, opts];
  If[ m === $Failed, $Failed, Ideal[m]]
];

(* An ideal looks as its basis, if $PrettyPrintIdeals is switched on *)
Ideal/: Format[Ideal[info_List]] := 
   If[ $PrettyPrintIdeals, 
       Apply[AngleBracket, ExtractInfoFromIdealObject[info, GBasis]],
       Ideal[Format[info]] (* inner "Format" avoids recursion *)
     ];

(**
 * ExtractInfoFromIdealObject[ideal, tag]
 *
 * Returns the value of <tag> in the <ideal>, or $Failed if no value for this tag
 * has been specified. 
 **)
ExtractInfoFromIdealObject[ideal_List, tag_] :=
   Module[ {result},

   result = tag /. ideal;

   If [ result === tag, $Failed, result]
];

GbMonomialOrder[Ideal[info_]] := ExtractInfoFromIdealObject[info, MonomialOrder];
GbCharacteristic[Ideal[info_]] := ExtractInfoFromIdealObject[info, Characteristic];
GbVariables[Ideal[info_]] := ExtractInfoFromIdealObject[info, Variables];
GbTime[Ideal[info_]] := ExtractInfoFromIdealObject[info, Time];

(**
 * Ideal[_] === Ideal[_]
 *
 * Two ideals are equal iff their orderings, their variable lists, their characteristic, 
 * and their Groebner bases agree literally. Otherwise they are not.
 **)
Ideal/: (Ideal[i1_] === Ideal[i2_]) := (   

   GbVariables[Ideal[i1]] === GbVariables[Ideal[i2]] &&
   GbMonomialOrder[Ideal[i1]] === GbMonomialOrder[Ideal[i2]] && 
   GbCharacteristic[Ideal[i1]] === GbCharacteristic[Ideal[i2]] &&
   Sort[GBasis[Ideal[i1]]] === Sort[GBasis[Ideal[i2]]]
);
Ideal/: (Ideal[i1_] =!= Ideal[i2_]) := !(Ideal[i1]===Ideal[i2])

(**
 * GbPretend[basis, vars, MonomialOrder -> ord]
 * 
 * Creates a raw ideal object, based on the given data.
 * It will be assumed without checking that the the given polynomial system is a 
 * Groebner basis wrt. the specified ordering.
 * 
 * The options MonomialOrder, MonomialOrdering, Modulus, Characteristic, CoefficientDomain
 * have the usual meaning. 
 **)
GbPretend[basis_, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {m},

   m = ComputeIdealInformation[basis, Pretend -> True, opts];
   If[ m === $Failed, $Failed, Ideal[m]]
];
GbPretend[basis_, vars_, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {m},

   m = ComputeIdealInformation[basis, Pretend -> True, Variables -> vars, opts];
   If[ m === $Failed, $Failed, Ideal[m]]
];
GbPretend[basis_, vars1_, vars2_, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {m},

   m = ComputeIdealInformation[basis, Pretend -> True, Variables -> {vars1,vars2}, opts];
   If[ m === $Failed, $Failed, Ideal[m]]
];

(******************************* Satellite Functions *******************************)

(**
 * GBasis[basis]
 *
 * Computes a Groebner basis for the given polynomial ideal.
 * See the documentation of ComputeIdealInformation for further information
 * (specification of MonomialOrderings etc.)
 *
 * The calls  GBasis[basis, vars]  and  GBasis[basis, vars1, vars2]
 * behave like the respective calls for MMA's builtin GroebnerBasis function. 
 * However, note that the default monomial ordering is DegreeReverseLexicographic,
 * as opposed to GroebnerBasis's default Lexicographic.
 **)
GBasis[Ideal[info_]] := ExtractInfoFromIdealObject[info, GBasis];
GBasis[basis_List, list1_List, list2_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {g},

   g = GBasis[GbIdeal[basis, list1, list2, opts]];
   If[ Head[g] === List, g, $Failed]
];
GBasis[basis_List, list_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {g},

   g = GBasis[GbIdeal[basis, list, opts]];
   If[ Head[g] === List, g, $Failed]
];

(**
 * GbEliminate[basis, away, stay]
 *
 * Computes a Groebnerbasis of <basis> in which the variables in <away> have been 
 * eliminated. 
 **)
GbEliminate[basis_, away_, stay_, opts:((_Rule|_RuleDelayed)...)] := 
   GBasis[basis, away, stay, Eliminate -> True, opts];
GbEliminate[Ideal[info_], away_List] :=
   Module[ {i, p, ord, vars, selector},

   (* If the monomial ordering of the ideal is already appropriate, we can do the
      elimination here, otherwise, we compute the elimination ideal by Gb. *)

   ord = GbMonomialOrder[Ideal[info]];
   vars = GbVariables[Ideal[info]];

   If[(MemberQ[{Lex, Lexicographic}, ord] &&
       Length[vars] >= away &&
       Table[vars[[i]], {i, 1, Length[away]}] === away) ||
      (MemberQ[{DRL, DegreeReverseLexicographic}, ord] &&
       Length[vars] === 2 &&
       vars[[1]] === away),

       (* then *)
       selector[p_] := Apply[And, Table[FreeQ[p, away[[i]]], {i, 1, Length[away]}]];

       Ideal[{MonomialOrder -> ord,
              GBasis -> Select[GBasis[Ideal[info]], selector],
              Variables -> If[ MatchQ[vars, {_List, _List}],
                               vars[[2]],
                               Table[vars[[i]], {i, Length[away] + 1, Length[vars]}]
                             ],
              Characteristic -> GbCharacteristic[Ideal[info]],
              Time -> 0. * Second
              }],

       (* else *)
       vars = Complement[Flatten[vars], away];
       GbIdeal[GBasis[Ideal[info]], 
               Variables -> {away, vars}, 
               MonomialOrder -> DRL, 
               Eliminate -> True,
               Characteristic -> GbCharacteristic[Ideal[info]]
              ]
   ];
];

(**
 * GbDimension[basis, vars]
 *
 * Returns the dimension of the ideal generated by the given basis.
 *
 * Example: GbDimension[{x*y}] == 1
 **)
GbDimension[Ideal[info_]] := ExtractInfoFromIdealObject[info, Dimension];
GbDimension[basis_List, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {d}, 

   d = GbDimension[GbIdeal[basis, vars, opts]];
   If[ Head[d] =!= Integer, $Failed, d]
];
GbDimension[basis_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {d},

   d = GbDimension[GbIdeal[basis, opts]];
   If[ Head[d] =!= Integer, $Failed, d]
];

(**
 * GbDegree[basis, vars]
 *
 * Returns the degree of the ideal given by the basis.
 **)
GbDegree[Ideal[info_]] := ExtractInfoFromIdealObject[info, Degree];
GbDegree[basis_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {d},

   d = GbDegree[GbIdeal[basis, opts]];
   If[ Head[d] =!= Integer, $Failed, d]
];
GbDegree[basis_List, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {d},

   d = GbDegree[GbIdeal[basis, vars, opts]];
   If[ Head[d] =!= Integer, $Failed, d]
];

(**
 * GbHilbertPolynomial[basis, vars]
 *
 * Returns the Hilbert polynomial of the ideal generated by the given basis.
 * 
 * Options:
 *     -- all the usual options (see docu of ComputeIdealInformation)
 *     -- Variable -> the variable in which the polynomial should be
 *        expressed. Default: T
 * See also: GbHilbertFunction
 **)
Options[GbHilbertPolynomial] = {Variable :> Symbol["T"]};
GbHilbertPolynomial[Ideal[info_], opts___] := 
   Module[ {h, d, v, var, a, k},

   var = Variable /. {opts} /. Options[GbHilbertFunction];
   
   h = ExtractInfoFromIdealObject[info, HilbertPolynomial[2]];
   d = ExtractInfoFromIdealObject[info, Dimension];

   a[k_] := If[k === 0, h /. Symbol["T"] -> 0, Coefficient[h, Symbol["T"]^k]];

   Sum[a[k]*Binomial[d + var - k - 1, d - 1], {k, 0, d}]      
];
GbHilbertPolynomial[basis_List, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {id, myOpts, var},

   myOpts = Apply[Sequence, DeleteCases[{opts}, Variable -> _]];

   id = GbIdeal[basis, vars, myOpts, MonomialOrdering -> DRL];

   myOpts = Apply[Sequence, Cases[{opts}, (Rule|RuleDelayed)[Variable, _]]];

   GbHilbertPolynomial[id, myOpts]
];

(**
 * GbHilbertFunction[basis, vars]
 *
 * Returns the Hilbert function of the ideal generated by the given basis.
 *
 * Options:
 *     -- all the usual options (see docu of ComputeIdealInformation)
 *     -- Variable -> the variable in which the polynomial should be
 *        expressed. Default: T
 * See also: GbHilbertPolynomial
 **)
Options[GbHilbertFunction] = {Variable :> Symbol["T"]};
GbHilbertFunction[Ideal[info_], opts___] :=
   Module[ {h, d, v, var},

   var = Variable /. {opts} /. Options[GbHilbertFunction];

   h = ExtractInfoFromIdealObject[info, HilbertPolynomial[2]];
   d = ExtractInfoFromIdealObject[info, Dimension];
   v = ExtractInfoFromIdealObject[info, Variables];

   h / (1 - Symbol["T"])^(Length[v] - d) /. Symbol["T"] -> var
];
GbHilbertFunction[basis_List, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {id, myOpts, var},

   myOpts = Apply[Sequence, DeleteCases[{opts}, Variable -> _]];

   id = GbIdeal[basis, vars, myOpts, MonomialOrdering -> DRL];

   myOpts = Apply[Sequence, Cases[{opts}, (Rule|RuleDelayed)[Variable, _]]];

   GbHilbertFunction[id, myOpts]
];

(**
 * GbFGLM[gbasis]
 * 
 * Converts a given Groebner basis for DRL into a Groebner basis for Lex.
 * The ideal has to be zero dimensional.
 **)
GbFGLM[basis_List, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
   GbFGLM[basis, Variables -> vars, opts];
GbFGLM[basis_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {info},
   
   info = ComputeIdealInformation[basis, MonomialOrder -> DRL2Lex, opts];

   If[!MatchQ[info, {___Rule}], Return[$Failed]];

   GBasis /. info
];
GbFGLM::bogusorder = "Illegal order: `1`. Expecting DegreeReverseLexicographic.";
GbFGLM::bogusdim = "Illegal dimension: `1`. Expecting 0.";
GbFGLM[Ideal[info_List], opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {ord, dim},

   If[ !MemberQ[{MonomialOrder, DRL, DegreeReverseLexicographic},
                 ord = MonomialOrder /. info],
       Message[GbFGLM::bogusorder, ord];
       Return[$Failed];
   ];

   If[ !MemberQ[{Dimension, 0}, dim = Dimension /. info],
       Message[GbFGLM::bogusdim, dim];
       Return[$Failed];
   ];

   Ideal[ComputeIdealInformation[
        GBasis[Ideal[info]],
        Variables -> GbVariables[Ideal[info]], 
        Characteristic -> GbCharacteristic[Ideal[info]],
        MonomialOrder -> DRL2Lex, 
        opts]]
];

(**
 * GbNormalForm[poly, ideal]
 *
 * Computes the normal form of the given <poly> wrt. the given <ideal>.
 *
 * Options:
 *        -- Algorithm -> specifies the prefix of the serveur to be used
 *        -- InFileName -> desired file name for input file of Gb
 *        -- OutFileName -> desired file name for output file of Gb
 **)
GbNormalForm[poly_, ideal_Ideal, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {ord, vars, basis, char, algo, variableTranslator,
            vars1, vars2, nf,
            in, out, inFileName, outFileName, command},

   ord = GbMonomialOrder[ideal];
   vars = GbVariables[ideal];
   basis = GBasis[ideal];
   char = GbCharacteristic[ideal];

   algo = Algorithm /. {opts} /. Options[Gb];
   inFileName = InFileName /. {opts} /. Options[Gb];
   outFileName = OutFileName /. {opts} /. Options[Gb];

   If[ MatchQ[vars, {_List, _List}],

      (* then *) vars1 = vars[[1]]; vars2 = vars[[2]],
      (* else *) vars1 = vars; vars2 = {}
   ];

   variableTranslator = 
       CleanVariables[Join[basis, vars1, vars2, {poly}], 
                      Variables -> (MemberQ[vars1, #] || MemberQ[vars2, #]&)];

   vars2 = Join[vars2, 
                Complement[Flatten[Map[ExtractVariables, Append[basis, poly]]], 
                           vars1, vars2]
               ];

   inFileName = WriteInputFile[{
      Variables -> (Join[vars1, vars2] /. variableTranslator[[1]]),
      List -> ({{Expand[poly]}, basis} /. variableTranslator[[1]])
   }, opts];

   outFileName = ServeurCall[
      algo <> "__DMP__" <> 
      Which[ ord === DRL && MatchQ[vars, {_List, _List}],
                "DRLDRL",
             ord === DRL,
                "Dexp",
             ord === Lex,
                "Lexp",
             True,
                Message[Gb::noorder, DRL]; "Dexp"] <> "__" <>
      If[ char === 0, "GINT", "GF" ], 
      Flags -> {"-mapnf", "-pipe", If[ MatchQ[vars, {_List, _List}], "-block " <> ToString[Length[vars[[1]]]], "" ] },
      InFileName -> inFileName,
      opts
   ];

   If[outFileName === $Failed, Return[$Failed]];

   (* read the input and return it. *)
   (* ASSUMPTION: The last line of the output file has the form *)
   (* [poly, number] *)
   (* where poly is the desired reduced polynomial *)

   out = Import[outFileName, "Lines"];
   If[ !FreeOfErrorsQ[out], Return[$Failed] ];
   nf = ToExpression[StringReplace["List" <> Last[out], "x" -> "Gb`Private`x"]];

   nf[[1]]/nf[[2]] /. variableTranslator[[2]]
];

(**
 * GbTriangular[ideal]
 *
 * computes a system of triangular sets for a given zero dimensional ideal. 
 **)
GbTri::bogusorder = "Illegal order: `1`. Expecting Lexicographic.";
GbTri::bogusdim = "Illegal dimension: `1`. Expecting 0.";
GbTriangular[Ideal[info_List], opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {algo, ord, dim, vars, inFileName, trans, out, i, result},

   algo = Algorithm /. {opts} /. Options[Gb];
   char = GbCharacteristic[Ideal[info]];
   vars = GbVariables[Ideal[info]];

   If[ !MemberQ[{MonomialOrder, Lex, Lexicographic},
                 ord = MonomialOrder /. info],
       Message[GbTri::bogusorder, ord];
       Return[$Failed];
   ];

   If[ !MemberQ[{Dimension, 0}, dim = Dimension /. info],
       Message[GbTri::bogusdim, dim];
       Return[$Failed];
   ];

   trans = CleanVariables[Join[GBasis[Ideal[info]], GbVariables[Ideal[info]]]];

   inFileName = WriteInputFile[{
           Variables -> (GbVariables[Ideal[info]] /. trans[[1]]),
           List -> (GBasis[Ideal[info]] /. trans[[1]]),
           Characteristic -> char, 
           Properties -> {"GB", "Lexp(" <> ToString[Length[vars]] <> ")"}
         }, opts];

   outFileName = ServeurCall[
      algo <> "__DMP__Lexp__" <> If[char === 0, "GINT", "GF"],
      Flags -> {"-lextri", "-pipe", If[char === 0, "", " -prime " <> ToString[char]]},
      InFileName -> inFileName,
      opts
   ];

   If[outFileName === $Failed, Return[$Failed]];

   out = Import[outFileName, "Lines"];
   If[ !FreeOfErrorsQ[out], Return[$Failed] ];

   (* Assumptions on format of the output file                  *)
   (* (1) there exists a line starting with #list               *)
   (* (2) the output data consists of all lines after that line *)
   (* (3) the concatenation of the output data lines gives a    *)
   (*     nested list [[...], ..., [...]] of polynomials        *)

   result = $Failed; 
   Do[
     If[ StringMatchQ[out[[i]], "#list*"], 
        (* Get the last (Length - i) lines. *)
        (* These are all lines behind the current line.*)
        result = Apply[StringJoin, Take[out, i - Length[out]]];
        (* Replace "[" by "List[" and try to parse. *)
        result = StringReplace[result, "[" -> "List["];
        result = ToExpression[result] /. trans[[2]];
        Break[];
     ],
     {i, 1, Length[out] - 1}
   ];

   result
];

(**
 * GbVersion[]
 *
 * Prints information about the version of the underlying Gb library.
 **)
Options[GbVersion] = {Serveur -> "serveur__DMP__Dexp__GF"};
GbVersion[opts:((_Rule|_RuleDelayed)...)] := 
   Module[ {tmp, tmpFile},

   tmp = OpenTemporary[];
   tmpFile = tmp[[1]];
   Close[tmp];

   Run[GbHome <> 
       "/" <> 
       (Serveur /. {opts} /. Options[GbVersion]) <> 
       " -version > " <>
       tmpFile];

   tmp = Import[tmpFile, "Lines"];
   If[ !FreeOfErrorsQ[tmp], Return[$Failed] ];
   Print[First[tmp]];
   Print["Gb Interface to Mathematica: Version " <> Version];
];

(**************************** Supplementary functions ***************************)

(**
 * RandomPolynomial[vars]
 *
 * Creates a random polynomial in the given variables. 
 * The polynomial is created by adding a certain number of random monomials with rational 
 * coefficients, each having a random total degree.
 * 
 * This function is intended for the generation of test cases. 
 *
 * Options: 
 *     -- Degree -> maximum total degree of the constructed polynomial
 *     -- Monomials -> number of random monomials to add up
 *     -- CoefficientBound -> maximum abs value for numerator and denominator of a single monomial. 
 **)
Options[RandomPolynomial] = {Degree -> 3, Monomials -> 5, CoefficientBound -> 10};
RandomPolynomial[vars_List, opts___] := 
   Module[ {d, m, c, i, j, p, deg, mon},

   d = Degree /. {opts} /. Options[RandomPolynomial];
   m = Monomials /. {opts} /. Options[RandomPolynomial];
   c = CoefficientBound /. {opts} /. Options[RandomPolynomial];

   Sum[
     deg = Random[Integer, d]; (* choose a total degree for the next monomial *)
     (* create a random monomial *)
     (2 * Random[Integer] - 1) (* random sign *)
       * Random[Integer, c]/Random[Integer, {1, 1}] (* random coeff *)
       * Product[vars[[Random[Integer, {1, Length[vars]}]]], {j, 1, deg}] (* random term *)
     , {i, 1, m}
   ]
];

(**
 * ExtendedGBasis[left, right][{b1, b2, ..., bn}, {x1, x2, ...}]
 *
 * Computes an extended Groebner basis in which the input polynomials are tagged by 
 * new variables right[1], right[2], ... which accumulate the actions that are taken
 * during the Groebner basis computation. 
 * 
 * It is required for the correctness of this function that none of the bi is the zero
 * polynomial.
 * 
 * The result of the computation is an ideal over the elimination order 
 * {left, right[1],...right[n]}, {x1, x2, ...} whose Groebner basis is a list of 
 * polynomials of the form 
 *
 *   left * p + right[1] * q1 + right[2] * q2 + ... + right[n] * qn
 *
 * such that b1*q1 + b2*q2 + ... + bn*qn == p and the list of all p builds a Groebner 
 * basis for {b1, b2, ..., bn}.
 *
 * (This is an internal function which is used by GbCofactors and GbExtendedGBasis.)
 * 
 * The usual options are applicable. Note that only one list of variables can be 
 * specified, i.e., no elimination orders are allowed for this function. 
 **)
ExtendedGBasis[left_, right_][basis_List, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {i},

   GbIdeal[
        Table[basis[[i]] * left - right[i], {i, 1, Length[basis]}], (* extended basis *)
        Join[{left}, Table[right[i], {i, 1, Length[basis]}]], (* markup variables are heavy *)
        vars, (* user's variables are not so heavy *)
        opts
      ]
];

(**
 * GbExtendedGBasis[{p1, p2, ...}, {x1, x2, ...}]
 *
 * Computes a Groebner basis for {p1, p2, ...} supplemented with the linear combination of each
 * element in terms of the original input basis. The result is a list of objects of the following
 * type
 *
 *   {b, {q1, q2, ...}}
 *
 * where b is the element of the Groebner basis and 
 *
 *   b == p1*q1 + p2*q2 + ...
 *
 * The usual options are applicable. 
 *
 * See also: GbCofactors
 **)
GbExtendedGBasis[basis_List, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {gb, b, f, e, i, b, j, eval},

   b = DeleteCases[basis, 0];

   gb = GBasis[ExtendedGBasis[f, e][b, vars, opts]];

   If[ gb === GBasis[$Failed], Return[$Failed] ];

   eval = Join[{f -> 0}, Table[e[i] -> 0, {i, 1, Length[b]}]];

   Table[
     {gb[[i]] /. (f -> 1) /. eval, Table[gb[[i]] /. (e[j] -> i) /. eval, {j, 1, Length[b]}]}
     , {i, 1, Length[gb]}
   ]
];

(**
 * GbCofactors[poly, {p1, p2, ...}, {x1, x2, ...}]
 * 
 * Computes cofactors of a given polynomial wrt. the given basis {p1, p2, ...}, and
 * returns $Failed if no linear combination exists. 
 * The computed cofactors are not necessarily unique.
 *
 * If the given poly is in the ideal generated by the given basis, then the function
 * will return a list {q1, q2, ...} of polynomials such that 
 *
 *   q1*p1 + q2*p2 + ... == poly
 * 
 * The usual options are applicable. Specified names for infiles and outfiles will be used twice. 
 **)
GbCofactors[poly_, basis_List, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {id, b, e, f, red, eval, i, lc, p}, 

   b = DeleteCases[basis, 0];

   id = ExtendedGBasis[f, e][b, vars, opts];
   red = GbNormalForm[f * poly, id, opts];
   eval = Join[{f -> 0}, Table[e[i] -> 0, {i, 1, Length[b]}]];

   If[ (red /. f -> 1 /. eval) =!= 0, 
       (* poly is not in the ideal *)
       $Failed
   ,
       (* the i-th cofactor is obtained by evaluating e[i] to 1 and the others to 0 *)
       lc = Table[red /. e[i] -> 1 /. eval, {i, 1, Length[b]}];

       (* insert the zeros which have been deleted before *)
       If[ Length[basis] > Length[b],
         p = Position[basis, 0, 1]; (* == {{firstindex}, {secondindex}, ...} *)
         Do[
           lc = Insert[lc, 0, p[[i]]];
         , {i, 1, Length[p]}
         ];
       ];

       lc
   ]
];

(******************************* End of private part ****************************)

End[];

(* protect defined names *)
Protect[GBasis, Ideal, GbCharacteristic, GbDegree,
    GbDimension, GbEliminate, GbFGLM, GbHilbertFunction, GbHilbertPolynomial,
    GbIdeal, GbMonomialOrder, GbNormalForm, GbPretend, GbTime, GbTriangular,
    GbVariables, GbVersion, GbCofactors, GbExtendedGBasis,
    MonomialOrder, MonomialOrdering, DRL, Lex, DegreeReverseLexicographic,
    Lexicographic, Modulus, Characteristic, Eliminate, Algorithm, Serveur,
    Variable, InFileName, OutFileName, CoefficientDomain, FiniteField
];

(* copyright note *)

(* Faugere header in green *)
If[$Notebooks,
   CellPrint[Cell[#, "Print", FontColor -> RGBColor[0, 0, 0], 
                              CellFrame -> 0.5, 
                              Background -> RGBColor[0.6, 0.8, 0.6]]]&,
   Print
  ]["Gb Library by Jean-Charles Faug\[EGrave]re \[LongDash] \[Copyright] 2003\n" <>
    "http://fgbrs.lip6.fr/jcf/"];

(* RISC header *)
If[$Notebooks,
   CellPrint[Cell[#, "Print", FontColor -> RGBColor[0, 0, 0], 
                              CellFrame -> 0.5, 
                              Background -> RGBColor[0.796887, 0.789075, 0.871107]]]&,
   Print
  ]["Gb Interface to Mathematica by Manuel Kauers and Burkhard Zimmermann " <> 
    "\[LongDash] \[Copyright] RISC Linz \[LongDash] V " <> Gb`Private`Version];


(********************************** End of package ******************************)

EndPackage[];

(* You may insert the the appropriate value here *)

GbHome = "/home/mkauers/scratch/systems/maple/faugere";