(**
 * Interface to Singular 
 * ---------------------
 * by Manuel Kauers and Viktor Levandovskyy
 *
 * TODO:
 * -- triangMH
 * -- timeout, full files (windows?)
 * -- timing
 **)

Singular`Private`SingularInterfaceVersion = "0.11 (2008-04-18)";
(* switched to MMA 6 from version 0.11 on *)

BeginPackage["Singular`"]

(* 0. Exported functions; declare public functions *)

Unprotect[SingularStd, SingularSlimgb, SingularGroebner, SingularEliminate, SingularReduce, SingularNF,
          SingularSyz, SingularQuotient, SingularSat, SingularIntersect, SingularDim, SingularVdim, 
          SingularLift, SingularFacstd, SingularMinAssGTZ, SingularPrimdecGTZ, SingularVersion, 
          SingularInterpolate, SingularNothing];

SingularStd::usage = "SingularStd[{f1,f2,...}, {x1,x2,...}..] computes a Groebner basis of the ideal (or module) generated by f1,f2... using Singular's std command";

SingularSlimgb::usage = "SingularSlimgb[{f1,f2,...}, {x1,x2,...}..] computes a Groebner basis of the ideal (or module) generated by f1,f2... using Singular's slimgb command";

SingularGroebner::usage = "SingularGroebner[{f1,f2,...}, {x1,x2,...}..] computes a Groebner basis of the ideal (or module) generated by f1,f2... using Singular's groebner command";

SingularEliminate::usage = "SingularEliminate[{f1,f2,...}, {x1,x2,...}, {y1,y2,...}] computes generators of the elimination ideal (or module) generated by f1,f2,... with the variables x1,x2,... being eliminated, using Singular's eliminate command.";

SingularReduce::usage = "SingularReduce[p, {f1,f2,...}, {x1,x2,...}..] computes the reductum of p with respect to the ideal generated by f1,f2,... using Singular's reduce command. If the f1,f2,... do not form a Groebner basis, the result may not be unique.";

SingularNF::usage = "SingularNF[p, {f1,f2,...}, {x1,x2,...}..] computes the reductum of p with respect to the ideal (or module) generated by f1,f2,... using Singular's NF command. If the f1,f2,... do not form a Groebner basis, the result may not be unique.";

SingularSyz::usage = "SingularSyz[{f1,f2,...}, {x1,x2,...}] computes a basis of the syzygy module of f1,f2,... using Singular's syz command.";

SingularPlus::usage = "SingularPlus[{f1,f2,...}, {g1,g2,...}, {x1,x2,...}] computes the sum of the two ideals (modules) generated by {f1,f2,..} and {g1,g2,...} respectively.";

SingularTimes::usage = "SingularTimes[{f1,f2,...}, {g1,g2,...}, {x1,x2,...}] computes the product of the two ideals generated by {f1,f2,..}, and {g1,g2,...} respectively.";

SingularQuotient::usage = "SingularQuotient[{f1,f2,...}, {g1,g2,...}, {x1,x2,...}] computes the quotient of the ideal (or module) generated by the f1,f2,.. and the ideal generated by g1,g2,...";

SingularSat::usage = "SingularSat[{f1,f2,...}, {g1,g2,...}, {x1,x2,...}] computes the saturation ideal (module) of the ideal (module) generated by f1,f2,... wrt. the ideal generated by g1,g2,... . Returns a pair {sat, exp} where sat is the result of the saturation and exp is the saturation exponent.";

SingularIntersect::usage = "SingularIntersect[{f1,f2,...}, {g1,g2,...}, ..., {x1,x2,...}] computes the intersection of the ideals (or modules) generated by f1,f2,... and g1,g2,... respectively.";

SingularMinAssGTZ::usage = "SingularMinAssGTZ[{f1,f2,...}, {x1,x2,...}] computes the minimal associated prime ideals of the ideal generated by f1,f2,...";

SingularPrimdecGTZ::usage = "SingularPrimdecGTZ[{f1,f2,...}, {x1,x2,...}] computes a primary decomposition of the ideal generated by f1,f2,...";

SingularDim::usage = "SingularDim[{f1,f2,...}, {x1,x2,...}] computes the Krull dimension of the ideal generated by f1,f2,... (of the annihilator of the module generated by f1,f2,...)"; 

SingularVdim::usage = "SingularVdim[{f1,f2,...}, {x1,x2,...}] computes the vector space dimension of the coordinate ring wrt. the ideal generated by f1,f2,...";

SingularFacstd::usage = "SingularFacstd[{f1,f2,...}, {x1,x2,...}] computes a factorizing Groebner basis";

SingularLift::usage = "SingularLift[{f1,f2,...}, {g1,g2,...}, {x1,x2,...}] computes a representation of g1,g2,.. in terms of the f1,f2,...";

SingularInterpolation::usage = "SingularInterpolation[{p1,p2,...},{e1,e2,...}, {x1,x2,...}] computes a basis for the zero dimensional ideal having the solutions p1 with multiplicity e1, p2 with multiplicity e2, ... Here p1, p2, ... are ideals representing single points.";

SingularVersion::usage = "SingularVersion[] prints the version of the underlying Singular platform.";

SingularCommand::usage = "Raw command to launch Singular, including path if necessary. Default: \"Singular\"";

SingularNothing::usage = "SingularNothing[{f1,f2,..}, {x1,x2,...}...] returns the ideal (module) {f1,f2,...}. (debug procedure.)";

Singular::noord = "No monomial ordering specified. Taking `1`";
Singular::error = "`1`";

Unprotect[Lexicographic, LexicographicSeries, DegreeReverseLexicographic, DegreeReverseLexicographicSeries,
          DegreeLexicographic, DegreeLexicographicSeries, WeightedLexicographic, WeightedLexicographicSeries,
          WeightedReverseLexicographic, WeightedReverseLexicographicSeries, ModuleAscending, ModuleDescending];

Lexicographic::usage =
  "Lexicographic is a setting for the MonomialOrder option of GroebnerBasis and PolynomialReduce. Also used to encode Singular's \"lp\".";
LexicographicSeries::usage =
  "LexicographicSeries is a setting for the MonomialOrder option of Singular commands (corresponds to Singular's \"ls\").";
DegreeReverseLexicographic::usage = 
  "DegreeReverseLexicographic is a setting for the MonomialOrder option of GroebnerBasis and PolynomialReduce. Also used to encode Singular's \"dp\".";
DegreeReverseLexicographicSeries::usage = 
  "DegreeReverseLexicographicSeries is a setting for the MonomialOrder option of Singular commands (corresponds to Singular's \"ds\").";
DegreeLexicographic::usage = 
  "DegreeLexicographic is a setting for the MonomialOrder option of Singular commands (corresponds to Singular's \"Dp\").";
DegreeLexicographicSeries::usage = 
  "DegreeLexicographicSeries is a setting for the MonomialOrder option of Singular commands (corresponds to Singular's \"Ds\").";
WeightedLexicographic::usage = 
  "WeightedLexicographic is a setting for the MonomialOrder option of Singular commands (corresponds to Singular's \"Wp\").";
WeightedLexicographicSeries::usage = 
  "WeightedLexicographicSeries is a setting for the MonomialOrder option of Singular commands (corresponds to Singular's \"Ws\").";
WeightedReverseLexicographic::usage = 
  "WeightedReverseLexicographic is a setting for the MonomialOrder option of Singular commands (corresponds to Singular's \"wp\")";
WeightedReverseLexicographicSeries::usage = 
  "WeightedReverseLexicographicSeries is a setting for the MonomialOrder option of Singular commands (corresponds to Singular's \"ws\").";
ModuleAscending::usage = 
  "ModuleAscending is a setting for the MonomialOrder option of Singular commands (corresponds to Singular's \"C\").";
ModuleDescending::usage =
  "ModuleDescending is a setting for the MonomialOrder option of Singular commands (corresponds to Singular's \"c\").";

InFileName::usage = 
  "InFileName is an option for Singular commands. It specifies the file name into which the input for Singular should be written. If set to Automatic, a temporary file will be used.";
OutFileName::usage = 
  "OutFileName is an option for Singular commands. It specifies the file name into which the output for Singular should be written. If set to Automatic, a temporary file will be used.";

DegBound::usage = "DegBound is an option to SingularStd";

Protect[Lexicographic, LexicographicSeries, DegreeReverseLexicographic, DegreeReverseLexicographicSeries,
        DegreeLexicographic, DegreeLexicographicSeries, WeightedLexicographic, WeightedLexicographicSeries,
        WeightedReverseLexicographic, WeightedReverseLexicographicSeries, ModuleAscending, ModuleDescending,
        DegBound];

Begin["`Private`"]

(* I. public interface *)

SingularStd[mod_?MatrixQ, vars__List, opts:((_Rule|_RuleDelayed)...)] := (
  SingularSession["std", {mod}, Module, Dimension -> ModuleDim[mod], Variables -> {vars}, opts]
);

SingularStd[ideal_?VectorQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  Module2Ideal[SingularStd[Ideal2Module[ideal], vars, opts]];

SingularSlimgb[mod_?MatrixQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["slimgb", {mod}, Module, Dimension -> ModuleDim[mod], Variables -> {vars}, opts];

SingularSlimgb[ideal_?VectorQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  Module2Ideal[SingularSlimgb[Ideal2Module[ideal], vars, opts]];

SingularGroebner[mod_?MatrixQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["groebner", {mod}, Module, Dimension -> ModuleDim[mod], Variables -> {vars}, opts];

SingularGroebner[ideal_?VectorQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  Module2Ideal[SingularGroebner[Ideal2Module[ideal], vars, opts]];

SingularEliminate[mod_?MatrixQ, elimvars_List, othervars_List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["eliminate", {mod, {Times@@elimvars}}, Module, Dimension -> ModuleDim[mod], 
    Variables -> {elimvars, othervars}, opts];

SingularEliminate[ideal_?VectorQ, elimvars_List, othervars_List, opts:((_Rule|_RuleDelayed)...)] :=
  Module2Ideal[SingularEliminate[Ideal2Module[ideal], elimvars, othervars, opts]];

SingularReduce[vector_?VectorQ, mod_?MatrixQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["reduce", {vector, mod}, Vector, Dimension -> Length[vector], Variables -> {vars}, 
     opts];

SingularReduce[poly_, ideal_?VectorQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  Vector2Poly[SingularReduce[Poly2Vector[poly], Ideal2Module[ideal], vars, opts]];

SingularNF = SingularReduce;

SingularSyz[mod_?MatrixQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["syz", {mod}, Module, Dimension -> Length[mod], Variables -> {vars}, opts];

SingularSyz[ideal_?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSyz[Ideal2Module[ideal], vars, opts]; (* no Module2Ideal conversion in this case! *)

SingularPlus[mod1_?MatrixQ, mod2_?MatrixQ, vars_List, opts___] /; ModuleDim[mod1] === ModuleDim[mod2] :=
  Join[mod1, mod2];

SingularPlus[id1_?VectorQ, id2_?VectorQ, vars_List, opts___] := Join[id1, id2];

SingularTimes[id1_?VectorQ, id2_?VectorQ, vars_List, opts___] := Module[ {i,j},
    Flatten[Table[ id1[[i]] * id2[[j]], {i, 1, Length[id1]}, {j, 1, Length[id2]}]]
  ];

SingularQuotient[mod_?MatrixQ, ideal_?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["quotient", {mod, ideal}, Module, Dimension -> ModuleDim[mod], Variables -> {vars}, opts];

SingularQuotient[ideal_?VectorQ, ideal2_?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  Module2Ideal[SingularQuotient[Ideal2Module[ideal], ideal2, vars, opts]];

SingularSat[mod_?MatrixQ, ideal_?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["sat", {mod, Ideal2Module[ideal]}, {Module, Integer}, Dimension -> ModuleDim[mod], 
       Variables -> {vars}, Library -> {"elim.lib"}, opts];

SingularSat[ideal_?VectorQ, ideal2_?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  MapAt[Module2Ideal, SingularSat[Ideal2Module[ideal], ideal2, vars, opts], 1];

SingularIntersect[mod__?MatrixQ, vars_List, opts:((_Rule|_RuleDelayed)...)] /; 
    Equal@@(ModuleDim /@ {mod}) :=
  SingularSession["intersect", {mod}, Module, Dimension -> ModuleDim[First[{mod}]], Variables -> {vars}, opts];

SingularIntersect[id__?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  Module2Ideal[SingularIntersect[##, vars, opts]&@@(Ideal2Module /@ {id})];

SingularMinAssGTZ[id_?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  Module2Ideal /@
    SingularSession["minAssGTZ", {Ideal2Module[id]}, List[Module], 
      Dimension -> 1, Variables -> {vars}, Library -> {"primdec.lib"}, opts];

SingularPrimdecGTZ[id_?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  Map[Module2Ideal,
    SingularSession["primdecGTZ", {Ideal2Module[id]}, List[Module], 
      Dimension -> 1, Variables -> {vars}, Library -> {"primdec.lib"}, opts], {2}];

SingularDim[mod_?MatrixQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["dim", {mod}, Integer, Variables -> {vars}, opts];

SingularDim[id_?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularDim[Ideal2Module[id], vars, opts];

SingularVdim[mod_?MatrixQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["vdim", {mod}, Integer, Variables -> {vars}, opts] /. -1 -> Infinity;

SingularVdim[id_?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularVdim[Ideal2Module[id], vars, opts];

SingularFacstd[id_?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  Map[Module2Ideal,
    SingularSession["facstd", {Ideal2Module[id]}, List[Module], Dimension -> 1, 
       Variables -> {vars}, opts]
  ];

SingularFacstd[id_?VectorQ, id2_?VectorQ, vars_List, opts:((_Rule|_RuleDelayed)...)] :=
  Map[Module2Ideal,
    SingularSession["facstd", {Ideal2Module[id], Ideal2Module[id2]}, List[Module], 
       Variables -> {vars}, opts]
  ];

SingularLift[id_?VectorQ, id2_?VectorQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularLift[Ideal2Module[id], Ideal2Module[id2], vars, opts];

SingularLift[mod_?MatrixQ, mod2_?MatrixQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["lift", {mod, mod2}, Module, Dimension -> Length[mod],
       Variables -> {vars}, opts];

SingularInterpolation[ids_?MatrixQ, mult_?VectorQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  Module2Ideal[SingularSession["interpolation", {ids, mult}, Module, Dimension -> 1, 
       Variables -> {vars}, opts]];

SingularVersion[opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["system", {"version"}, Integer, Variables -> {{x}}, opts];

SingularNothing[mod_?MatrixQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  SingularSession["", {mod}, Module, Dimension -> ModuleDim[mod], Variables -> {vars}, opts];

SingularNothing[id_?VectorQ, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  Module2Ideal[SingularNothing[Ideal2Module[id], vars, opts]];

SingularNothing[poly_, vars__List, opts:((_Rule|_RuleDelayed)...)] :=
  First[SingularNothing[{poly}, vars, opts]];

SingularCommand = 
  Switch[$System,
    "Microsoft Windows",
       (* assuming presence of cygwin and Singular in cygwin's $PATH *) 
       "bash Singular", 
    _, "Singular" (* all other systems *)
  ];

(* II. private routines *)

Protect[`a, `gen];

(* 1. "type conversions" between module and ideal stuff *)

(* a) {{p1}, {p2}, {p3}, ...} --> {p1, p2, ...} *)
Module2Ideal[$Failed] = $Failed;
Module2Ideal[mod_?MatrixQ] := First[Transpose[mod]];
Module2Ideal[{}] = {};

(* b) {p1, p2, p3, ...} -> {{p1}, {p2}, ...} *)
Ideal2Module[$Failed] = $Failed;
Ideal2Module[ideal_?VectorQ] := Transpose[{ideal}];
Ideal2Module[{}] = {{0}}; (* trivial generator included to make module distinguishable from ideal *)

(* c) {p} -> p *)
Vector2Poly[$Failed] = $Failed;
Vector2Poly[{poly_}] := poly;

(* d) p -> {p} *)
Poly2Vector[$Failed] = $Failed;
Poly2Vector[poly_] := {poly};

(* Determines the dimension (= number of coordinates) of a module *)
ModuleDim[mod_?MatrixQ] := Length[First[mod]];

(**
 * MyVariables[expr]
 *
 * Returns the largest subexpressions of expr whose head is neither List nor 
 * an arithmetic operation. 
 **)
ClearAll[MyVariables];
MyVariables[{a__}] := Union@@(MyVariables /@ {a});
MyVariables[a_^_Integer] := MyVariables[a];
MyVariables[a_ * b__] := MyVariables[{a, b}];
MyVariables[a_ + b__] := MyVariables[{a, b}];
MyVariables[_?NumberQ] = {};
MyVariables[_?NumericQ] = {};
MyVariables[_String] = {};
MyVariables[a_] := {a}; (* variable found *)

(**
 * MakeString[ex1, ex2, ex3,...]
 *
 * Produces the string obtained by concatenation of the arguments. 
 * Arguments which are not strings are converted into strings.
 * The substrings "Singular`Private`" and "Singular`" will be deleted from the result.
 **)
MakeString[args___] := 
   StringReplace[
     Apply[StringJoin,
       Map[ If[Head[#]=!=String, ToString[#, FormatType -> InputForm], #]&,
         {args}
       ] (* end map *)
     ], (* end apply *)
     {"Singular`Private`" -> "", "Singular`" -> ""}
   ];

ClearAll[MyPoly2String];
Poly2String[p_] := 
   Block[ {$RecursionLimit = Infinity},
     StringReplace[StringReplace[MyPoly2String[p], "-1*" -> "-"], "+-" -> "-"]
   ];
MyPoly2String[a_+b_] := Poly2String[a] <> "+" <> Poly2String[b];
MyPoly2String[a_*b_] /; Complement[{Head[a], Head[b]}, {Times, Power Symbol, Integer, Rational, gen}] === {} := 
   Poly2String[a] <> "*" <> Poly2String[b];
MyPoly2String[a_*b_] /; MemberQ[{Times, Power, Symbol, Integer, Rational, gen}, Head[a]] :=
   Poly2String[a] <> "*(" <> Poly2String[b] <> ")";
MyPoly2String[a_*b_] := "(" <> Poly2String[a] <> ")*(" <> Poly2String[b] <> ")";
MyPoly2String[a_?NumberQ*b_] := ToString[a] <> "*" <> Poly2String[b];
MyPoly2String[a_Symbol^b_] /; b > 0 := Poly2String[a] <> "^" <> Poly2String[b];
MyPoly2String[a_Symbol^b_] /; b < 0 := Poly2String[a] <> "^(" <> Poly2String[b] <> ")";
MyPoly2String[a_^b_Integer] /; b > 0 := "(" <> Poly2String[a] <> ")^" <> Poly2String[b];
MyPoly2String[a_^b_Integer] /; b < 0 := "(" <> Poly2String[a] <> ")^(" <> Poly2String[b] <> ")";
MyPoly2String[a_Symbol] := MakeString[a];
MyPoly2String[a_?NumberQ] := MakeString[a];
MyPoly2String[gen[i_]] := "gen(" <> ToString[i] <> ")";
MyPoly2String[a_] := Throw["Unable to convert " <> ToString[a] <> " to Singular expression."];

(**
 * GetType[expr]
 *
 * Detects whether the given expression encodes an object of one of the following types:
 * "poly", "ideal", "vector", "module", "int", and returns the appropriate string.
 * An exception is thrown if no type matches.
 * It is assumed that polynomials and ideals are encoded in module form, i.e.,
 * the ideal generated by a,b,c is represented by {{a},{b},{c}}, while {a,b,c} is a vector.
 **)
ClearAll[GetType];
GetType[{{_}..}] = "ideal";
GetType[_?MatrixQ] = "module";
GetType[{_}] = "poly";
GetType[_?VectorQ] = "vector";
GetType[_Integer] = "int";
GetType[_String] = "string";

GetType[Module, dim_Integer] := If[ dim==1, "ideal", "module"];
GetType[Vector, dim_Integer] := If[ dim==1, "poly", "vector"];
GetType[Integer, _] = "int";
GetType[String, _] = "string";
GetType[{___}, _] = "list";

GetType[x_] := (Throw["Unknown type encountered."]);

(**
 * GetString[expr]
 * 
 * Converts an expression for a module or a vector into Singular synatx for an object of the
 * type that GetType returns for expr.
 **)
ClearAll[GetString];
GetString[id:{{_}..}] := StringTake[ToString[GetString /@ id], {2, -2}];
GetString[mod_?MatrixQ] := StringTake[ToString[GetString /@ mod], {2, -2}];
GetString[{poly_}] := Poly2String[poly];
GetString[vec_?VectorQ] := Module[ {i},
  Poly2String[Sum[vec[[i]]*gen[i],{i, 1, Length[vec]}]]
];
GetString[int_Integer] := MakeString[int];
GetString[str_String] := "\"" <> str <> "\"";
GetString[x_] := (Throw["Unknown type encountered."]);

(**
 * ParseOutput[type, listofstrings]
 *
 * Constructs an MMA expression which corresponds to the Singular expression encoded in the 
 * list of strings.
 *
 * Type might be Module, Vector, Integer, String.
 * Module also accepts ideals, Vector also accepts polynoials. 
 **)
ClearAll[ParseOutput];
ParseOutput[Integer, dim_Integer, {int_String}] := ToExpression[int];
ParseOutput[Vector, dim_Integer, {vec_String}] := Module[ {i, v},
  If[ StringMatchQ[vec, "[*"],
    ToExpression[StringReplace[
      vec,
      {"[" -> "{", "]" -> "}", "x" -> "Singular`Private`x", 
       "t" -> "Singular`Private`t", "a" -> "Singular`Private`a"}
    ]]
  ,
    v = ToExpression[StringReplace[vec, 
          Join[
            Table["gen(" <> ToString[i] <> ")" -> "Singular`Private`gen[" <> ToString[i] <> "]", {i, 1, dim}],
            {"x" -> "Singular`Private`x", "t" -> "Singular`Private`t", "a" -> "Singular`Private`a"}
          ]]
        ];
    If[ dim > 1, Table[Coefficient[v, gen[i]], {i, 1, dim}], {v} ]
  ]
];
ParseOutput[Module, dim_Integer, {mod_String}] := Module[ {s, i, p, gens},
  If[ StringMatchQ[mod, "[*"],
    Join[#, Table[0, {dim - Length[#]}]]& /@ 
      ToExpression[StringReplace[
        "{" <> mod <> "}",
        {"[" -> "{", "]" -> "}", "x" -> "Singular`Private`x", 
         "t" -> "Singular`Private`t", "a" -> "Singular`Private`a"}
      ]]
  ,
    s = "," <> mod <> ",";
    p = Union[Flatten[StringPosition[s, ","]]];
    gens = Table[ 
      ParseOutput[Vector, dim, {StringTake[s, {p[[i]] + 1, p[[i + 1]] - 1}]}],
      {i, 1, Length[p] - 1}
    ]
  ]
];
ParseOutput[List[type_], dim_Integer, out_List] := Module[ {i}, 
   Table[ParseOutput[type, dim, {out[[i]]}], {i, 1, Length[out]}]
];
ParseOutput[List[type_, types__], dim_Integer, out_List] := (
   Prepend[ParseOutput[{types}, dim, Rest[out]], ParseOutput[type, dim, {First[out]}]]
);

ParseOutput[__] := Throw["Unable to read Singular's output."];

(* 3. Singular session
      given a ring (characteristic, algebraic extension, list of parameters, list of variables, order spec.)
      a command, and a number of arguments to that command (might be modules or vectors),
      this command composes a suitable Singular session, executes it, and re-formulates the result
      in MMA language. The result type may be Integer, {Module, dim}, {Vector, dim}.
 *)
Options[SingularSession] = 
  {Variables :> Throw["No variables specified."], InFileName -> Automatic, OutFileName -> Automatic,
   Dimension -> 1, Modulus -> 0, MonomialOrder -> Automatic, Library -> {}, Throw -> False, DegBound -> 0};
SingularSession[command_String, args_List, resulttype_, opts___] :=
   Module[ {myargs, alg, trans, varBlocks, t, x, minpoly, char, order, variableConverter, i, j, vars,
           infile, outfile, in, out, keepInfile, keepOutfile, root, result, dim, lib, throw},

   result = Catch[

     (***** a. extract and check options *****)
     {varBlocks, infile, outfile, char, order, char, dim, lib, throw} =
     {Variables, InFileName, OutFileName, Modulus, MonomialOrder, Modulus, Dimension, Library, Throw} 
         /. {opts} /. Options[SingularSession];

     If[ !IntegerQ[char] || char < 0 || (char > 0 && !PrimeQ[char]), 
       Throw["Illegal modulus: " <> ToString[char]]
     ];

     If[ Complement[First /@ {opts}, First /@ Options[SingularSession]] =!= {},
       Throw["Illegal options encountered: " <> 
         ToString[Complement[First /@ {opts}, First /@ Options[SingularSession]]]
       ]
     ];

     If[ !FreeQ[ varBlocks, {} ],
       Throw["Illegal variable declaration: " <> 
         ToString[varBlocks]];
     ];

     (***** b. put algebraic numbers into a common field. *****)
     alg = Union[Cases[args, (AlgebraicNumber[_,_]|Root[_,__]|Complex[_,_]|Power[_,_Rational]|Sqrt[_]), Infinity]];
     myargs = args;

     If[ Length[alg] > 0,
       alg = Apply[Rule, Transpose[{alg, ToNumberField[alg]}], {1}];
       root = alg[[1,2,1]];
       minpoly = MinimalPolynomial[root];
       (* insert minimal polynomial to all coordinates of all modules appearing in the argument list *)
       Do[
         If[ GetType[myargs[[i]]] == "module" || GetType[myargs[[i]]] == "ideal",
           myargs[[i]] = Join[myargs[[i]], 
             minpoly[a] * IdentityMatrix[ModuleDim[myargs[[i]]]]
           ];
         ]
       , {i, 1, Length[myargs]}];
     ,
       alg = {};
       minpoly = 0&;
       root = 1;
     ];

     myargs = myargs /. alg /. AlgebraicNumber[_, x_List] :> Sum[x[[i]]*a^(i-1), {i,1,Length[x]}];

     (***** c. determine transcendental elements besides variables *****)
     vars = Flatten[varBlocks];
     Block[ {$RecursionLimit = 10^5},
       trans = Complement[MyVariables[myargs], vars, {a}]
     ];

     (***** d. find nicer variables: a for algebraic, t1,t2,... for parameters, x1,x2,... for variables *****)
     variableConverter = 
        Join[
          Table[trans[[i]] ->  Symbol["Singular`Private`t" <> ToString[i]], {i, 1, Length[trans]}],
          Table[vars[[i]] -> Symbol["Singular`Private`x" <> ToString[i]], {i, 1, Length[vars]}]
        ];

     {myargs, trans, vars} = {myargs, trans, vars} /. variableConverter;

     (***** e. choose term order *****)
     order = order /. 
         Automatic :> 
           (If[ MemberQ[{"groebner", "std", "slimgb", "NF", "reduce", "dim", "vdim"}, command],
              Message[Singular::noord, DegreeReverseLexicographic]
            ]; DegreeReverseLexicographic);

     order = order /. {DegreeReverseLexicographic -> "dp",   DegreeReverseLexicographicSeries -> "ds",
                       DegreeLexicographic -> "Dp",          DegreeLexicographicSeries -> "Ds",
                       Lexicographic -> "lp",                LexicographicSeries -> "ls",
                       WeightedReverseLexicographic -> "wp", WeightedReverseLexicographicSeries -> "ws",
                       WeightedLexicographic -> "Wp",        WeightedLexicographicSeries -> "Ws",
                       ModuleAscending -> "C",               ModuleDescending -> "c"
                      };

     If[ Length[varBlocks] > 1 && Head[order] =!= List,
       order = Table[order, {Length[varBlocks]}];
     ];

     If[ Length[alg] > 0, 
       AppendTo[varBlocks, {a}]; AppendTo[vars, a];
       Which[ 
         Head[order] === List && !MatrixQ[order],
           AppendTo[order, "lp"],
         True,
           order = {order, "lp"}
       ];
     ];

     order = Which[
       MemberQ[{"dp", "ds", "Dp", "Ds", "lp", "ls"}, order], 
         (* one block with standard order only *)
         order,
       MatchQ[order, ("wp"|"Wp"|"ws"|"Ws")[__Integer]] && Length[order] === Length[vars],
         (* weighted order *)
         MakeString[order],
       MatrixQ[order, NumberQ] && Length[order] === Length[vars], 
         (* explicit matrix order, one block only *)
         "M" <> MakeString[Flatten[order]],
       Head[order] === List,
         (* block order *)
         If[ FreeQ[order, "c"] && FreeQ[order, "C"], AppendTo[order, "C"] ];
         If[ Length[order] =!= Length[varBlocks] + 1, 
           Throw["Number of blocks must agree with number of orders"] 
         ];
         j = 0; (* index shift due to appearence of C in the order list *)
         ToString[
           Table[
             Which[
               MemberQ[{"dp", "ds", "Dp", "Ds", "lp", "ls"}, order[[i]]],
                 order[[i]] <> "(" <> ToString[Length[varBlocks[[i + j]]]] <> ")", 
               (order[[i]] === "C" || order[[i]] === "c") && j === 0,
                 j--; order[[i]],
               MatchQ[order[[i]], ("wp"|"Wp"|"ws"|"Ws")[__Integer]] && Length[order[[i]]] === Length[varBlocks[[i + j]]],
                 MakeString[order[[i]]],
               MatrixQ[order[[i]], NumberQ] && Length[order[[i]]] === Length[varBlocks[[i + j]]],
                 "M" <> MakeString[Flatten[order[[i]]]],
               True,
                 Throw["Unknown sub term order in block order, or block size mismatch: "]
             ],
           {i, 1, Length[varBlocks] + 1}]
         ],
       True,
         Throw["Unknown term order: " <> ToString[order]];
     ];

     order = StringReplace[order, {"{" -> "(", "}" -> ")", "\"" -> "", "[" -> "(", "]" -> ")"}];

     (***** f. write input file *****)
     keepInfile = False; keepOutfile = False;
     in = If[ Head[infile] === String, keepInfile = True; OpenWrite[infile], OpenTemporary[] ];
     infile = First[in];
     out = If[ Head[outfile] === String, keepOutfile = True; OpenWrite[outfile], OpenTemporary[] ];
     outfile = First[out];
     Close[out];

     (* --> load libraries *)
     Do[
       WriteString[in, "LIB \"" <> lib[[i]] <> "\";\n"]
     , {i, 1, Length[lib]}];     

     (* --> ring declaration *)
     WriteString[in, 
       StringReplace[
         MakeString["ring r = ", (Prepend[trans, char] /. {a_} -> a), ", ", vars, ", ", order, ";\n"],
         {"{" -> "(", "}" -> ")"}
       ]
     ];

     (* --> flags *)
     WriteString[in, "option(redSB); option(redTail);\n"];
     WriteString[in, "degBound = ", DegBound /. {opts} /. Options[SingularSession], ";\n"];

     (* --> arguments *)
     Which[ 
       command == "interpolation",
         If[ Length[myargs] != 2, Throw["illegal arguments to 'interpolation'"] ];
         Do[
           WriteString[in,
             MakeString["ideal p", i, " = ", GetString[{#}& /@ myargs[[1, i]]], ";\n"]]
         , {i, 1, Length[First[myargs]]}];
         WriteString[in, "list a1 = "];
         Do[ WriteString[in, "p" <> ToString[i] <> ", "], {i, 1, Length[First[myargs]] - 1} ];
         WriteString[in, MakeString["p", Length[First[myargs]], ";\n"]];
         WriteString[in, MakeString["intvec a2 = ", GetString[{#}& /@ myargs[[2]]], ";\n"]],
       True, (* default argument queue *)
         Do[
           WriteString[in, 
             MakeString[GetType[myargs[[i]]], " a", i, " = ", GetString[myargs[[i]]], ";\n"]
           ]
         , {i, 1, Length[myargs]}]
     ];

     (* --> command *)
     WriteString[in,
       StringReplace[
         MakeString[GetType[resulttype, dim], " result = ", command,
                      Table[Symbol["a" <> ToString[i]], {i, 1, Length[myargs]}], ";\n"],
         {"{" -> "(", "}" -> ")"}
       ]               
     ];

     (* --> export *)
     If[ Head[resulttype] === List, 
       WriteString[in, "for(int i=1; i <= size(result); i++) {\n"];
       If[ StringMatchQ[command, "primdec*"], 
          (* Singular is too stupid to write nested lists properly into a file. 
             Argh... We need a hack here. *)
         WriteString[in, "for(int j=1; j <= 2; j++) {\n"];
       ]
     ];
     WriteString[in, 
       MakeString[
         "write(\":a ", 
         StringReplace[outfile, "\\" -> "\\\\"], (* '\' must be escaped in windows *)
         "\", result",
         Which[
           StringMatchQ[command, "primdec*"], "[i][j]",
           Head[resulttype] === List, "[i]",
           True, ""
         ], ");\n"
       ]
     ];
     If[ Head[resulttype] === List,
       WriteString[in, "};\n"];
       If[ StringMatchQ[command, "primdec*"], 
         WriteString[in, "};\n"];
       ];
     ];

     (* --> quit *)
     WriteString[in, "quit;\n"];

     Close[in];

     (***** g. execute singular *****)
     (* option -b deleted, causes problems in windows *)
     Run[SingularCommand <> " -tq --no-warn ", First[in]]; 
     If[ !keepInfile, DeleteFile[infile] ];

     (***** h. read and parse output file *****)
     out = Import[outfile, "Lines"]; (* list of strings *)
     If[ !keepOutfile, DeleteFile[outfile] ];

     result = ParseOutput[resulttype, dim, out];

     (* postprocessing hack for primdecGTZ *)
     If[ StringMatchQ[command, "primdec*"], 
       result = Table[{result[[2i-1]], result[[2i]]}, {i, 1, Length[result]/2}];
     ];

     (***** i. switch to user's variables, and introduce algebraic numbers *****)
     result = result /. (Reverse /@ variableConverter);
     If[ Length[alg] > 0, 
       result = result /. a -> AlgebraicNumber[root, {0, 1}] // RootReduce
     ];
     result = DeleteCases[result, {0..}];

     (***** k. return result *****)
     Return[result]
   ];

   If[ Head[result] === String,
     If[ throw, 
       Throw["Singular exception: " <> result]
     ,
       Message[Singular::error, result]; Return[$Failed]
     ]
   ,
     Return[result]
   ]
];

End[]

Protect[SingularStd, SingularSlimgb, SingularGroebner, SingularEliminate, SingularReduce, SingularNF,
        SingularSyz, SingularQuotient, SingularSat, SingularIntersect, SingularDim, SingularVdim, 
        SingularMinAssGTZ, SingularPrimdecGTZ, SingularLift, SingularFacstd, SingularVersion, 
        SingularInterpolate, SingularNothing];

If[$Notebooks,
   CellPrint[Cell[#, "Print", FontColor -> RGBColor[0, 0, 0], 
                              CellFrame -> 0.5, 
                              Background -> RGBColor[0.796887, 0.789075, 0.871107]]]&,
   Print
  ]["Singular -- Interface to Mathematica   Package by Manuel Kauers (mkauers@risc.uni-linz.ac.at) and Viktor Levandovskyy (levandov@risc.uni-linz.ac.at)\nhttp://www.risc.uni-linz.ac.at/research/combinat/software/Singular/ \[LongDash] \[Copyright] RISC Linz \[LongDash] V " 
    <> Singular`Private`SingularInterfaceVersion]

EndPackage[]

(* uncomment the following line and modify it according to your needs, if appropriate *)
(* SingularCommand = "Singular" *)