Index: Singular/LIB/ncfrac.lib
===================================================================
--- Singular/LIB/ncfrac.lib	(revision 3f58e0f44e7e559abc73c266dcd710d43ccfccb3)
+++ Singular/LIB/ncfrac.lib	(revision 3f58e0f44e7e559abc73c266dcd710d43ccfccb3)
@@ -0,0 +1,1023 @@
+//////////////////////////////////////////////////////////////////////
+version="version ncfrac.lib 4.0.0.0 Dec_2017 "; //$Id$
+category="Noncommutative";
+info="
+LIBRARY:   ncfrac.lib  object-oriented interface for olga.lib
+AUTHOR:    Johannes Hoffmann, email: johannes.hoffmann at math.rwth-aachen.de
+
+OVERVIEW:
+This library introduces a new type: ncfrac.
+This type wraps the data defining a (non-commutative) fraction in an Ore
+localization of a G-algebra as in olga.lib.
+An element of type ncfrac has five members:
+  - polys lnum, lden, rnum, rden
+  - ncloc loc
+
+OPERATIONS:
+string(ncfrac);
+  give a string representation of the data describing the fraction
+print(ncfrac);
+  prints the string representation of the fraction
+status(ncfrac);
+  report on the status/validity of the fraction
+test(ncfrac);
+  check if the fraction is valid
+
+INFIX OPERATIONS:
+ncfrac == ncfrac;
+  compare two fractions
+ncfrac != ncfrac;
+  compare two fractions
+ncfrac + ncfrac;
+  add two fractions
+ncfrac - ncfrac
+  subtract two fractions
+ncfrac * ncfrac
+  multiply two fractions
+ncfrac / ncfrac
+  divide two fractions
+ncfrac = int/number/poly
+  create a fraction with:
+    - left and right denominator equal to 1
+    - left and right numerator determined by the input
+    - localization data describing the trivial monoidal localization at 1
+ncfrac = vector
+  create a fraction from a vector v with unspecified localization such that
+  lden,lnum,rnum,rden = v[1],v[2],v[3],v[4]
+  (note: without specifying a localization afterwards this results is an
+  invalid fraction)
+ncfrac = list
+  create a fraction from a list L as follows:
+    - try to create a fraction from L[1] as above
+    - if L[2] is of type ncloc set the localization of the fraction to L[2]
+
+PROCEDURES:
+hasLeftDenom(ncfrac); check if the given fraction has a left representation
+hasRightDenom(ncfrac); check if the given fraction has a right representation
+isZeroNcfrac(ncfrac); check if the given fraction is equal to zero
+isOneNcfrac(ncfrac); check if the given fraction is equal to one
+zeroNcfrac(ncloc loc); create the fraction equal to zero in the given localization
+oneNcfrac(ncloc loc); create the fraction equal to one in the given localization
+ensureLeftNcfrac(ncfrac); compute a left representation of the given fraction if it does not have one
+ensureRightNcfrac(ncfrac); compute a right representation of the given fraction if it does not have one
+negateNcfrac(ncfrac); compute the additive inverse of the given fraction
+isInvertibleNcfrac(ncfrac); check if the given fraction is invertible (note: see the description of
+  isInvertibleLeftFraction from olga.lib for specific behaviour)
+invertNcfrac(ncfrac); compute the multiplicative inverse of the given fraction  (note: see the
+  description of invertLeftFraction from olga.lib for specific behaviour)
+testNcfrac(); execute a series of internal testing procedures
+testNcfracExamples(); execute the examples of all procedures in this library
+";
+//////////////////////////////////////////////////////////////////////
+proc testNcfracExamples()
+"USAGE:   testNcfracExamples()
+PURPOSE: execute the examples of all procedures in this library
+RETURN:  nothing
+NOTE:
+EXAMPLE: "
+{
+    example hasLeftDenom;
+    example hasRightDenom;
+    example isZeroNcfrac;
+    example isOneNcfrac;
+    example zeroNcfrac;
+    example oneNcfrac;
+    example ensureLeftNcfrac;
+    example ensureRightNcfrac;
+}
+//////////////////////////////////////////////////////////////////////
+proc testNcfrac()
+"USAGE:   testNcfrac()
+PURPOSE: execute a series of internal testing procedures
+RETURN:  nothing
+NOTE:
+EXAMPLE: "
+{
+    print("testing ncfrac.lib...");
+    testHasLeftHasRight();
+    testIsZeroIsOne();
+    testNcFracCreation();
+    testNcfracComparison();
+    testNcfracAddition();
+    testNcfracSubtraction();
+    testNcfracMultiplication();
+    testNcfracDivision();
+    testNcfracInversion();
+    testEnsureRightNcfrac();
+    testEnsureLeftNcfrac();
+    print("testing complete - ncfrac.lib OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc mod_init() {
+    LIB "ncloc.lib";
+    /* new type: ncfrac (non-commutative fraction) */
+    newstruct("ncfrac", "ncloc loc, poly rden, poly rnum, poly lden, poly lnum");
+    system("install", "ncfrac", "print", printNcfrac, 1);
+    system("install", "ncfrac", "string", ncfracToString, 1);
+    system("install", "ncfrac", "=", createNcfrac, 1);
+    system("install", "ncfrac", "==", compareNcfracs, 2);
+    system("install", "ncfrac", "!=", invertedCompareNcfracs, 2);
+    system("install", "ncfrac", "+", addNcfracs, 2);
+    system("install", "ncfrac", "-", subtractNcfracs, 2);
+    system("install", "ncfrac", "*", multiplyNcfracs, 2);
+    system("install", "ncfrac", "/", divideNcfracs, 2);
+    system("install", "ncfrac", "test", testNcfrac, 4);
+    system("install", "ncfrac", "status", statusNcfrac, 4);
+}
+//////////////////////////////////////////////////////////////////////
+proc hasLeftDenom(ncfrac frac)
+"USAGE:   hasLeftDenom(frac), ncfrac frac
+PURPOSE: checks if frac has a left representation
+RETURN:  int, 1 if frac has a left representation, 0 otherwise
+EXAMPLE: example hasLeftDenom; shows example"
+{
+    return(frac.lden != 0);
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    ncloc loc = ideal(x-3,y+7);
+    ncfrac noLeft = list([0,0,3*y*Dx,x+2], loc);
+    hasLeftDenom(noLeft);
+    ncfrac left = list([1,Dx,Dx,1], loc);
+    hasLeftDenom(left);
+}
+//////////////////////////////////////////////////////////////////////
+proc hasRightDenom(ncfrac frac)
+"USAGE:   hasRightDenom(frac), ncfrac frac
+PURPOSE: checks if frac has a right representation
+RETURN:  int, 1 if frac has a right representation, 0 otherwise
+EXAMPLE: example hasRightDenom; shows example"
+{
+    return(frac.rden != 0);
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    ncloc loc = ideal(x-3,y+7);
+    ncfrac noRight = list([x+2,3*y*Dx,0,0], loc);
+    hasRightDenom(noRight);
+    ncfrac right = list([1,Dx,Dx,1], loc);
+    hasRightDenom(right);
+}
+////////// status and printing ///////////////////////////////////////
+static proc toVector(ncfrac frac) {
+    return([frac.lden, frac.lnum, frac.rnum, frac.rden]);
+}
+//////////////////////////////////////////////////////////////////////
+static proc statusNcfrac(ncfrac frac) {
+    return(fracStatus(toVector(frac), frac.loc.locType, frac.loc.locData));
+    // fracStatus from olga.lib
+}
+//////////////////////////////////////////////////////////////////////
+static proc testNcfrac(ncfrac frac) {
+    list stat = status(frac);
+    if(!stat[1]) {
+        ERROR(stat[2]);
+    } else {
+        return();
+    }
+}
+//////////////////////////////////////////////////////////////////////
+static proc ncfracToString(ncfrac frac) {
+    list stat = status(frac);
+    if (!stat[1]) {
+        return(stat[2]);
+    } else {
+        return("left repr.: (" + string(frac.lden) + "," + string(frac.lnum)
+          + ")" + newline + "right repr.: (" + string(frac.rnum) + ","
+          + string(frac.rden) + ")");
+    }
+}
+//////////////////////////////////////////////////////////////////////
+static proc printNcfrac(ncfrac frac) {
+    string(frac);
+}
+//////////////////////////////////////////////////////////////////////
+proc isZeroNcfrac(ncfrac frac)
+"USAGE:   isZeroNcfrac(frac), ncfrac frac
+PURPOSE: checks if frac is zero
+RETURN:  int, 1 if frac is zero, 0 otherwise
+EXAMPLE: example isZeroNcfrac; shows example"
+{
+    testNcfrac(frac);
+    if (hasLeftDenom(frac)) { // frac has left representation
+        return(frac.lnum == 0);
+    } else { // frac has right representation, but no left representation
+        return(frac.rnum == 0);
+    }
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring Q = (0,q),(x,y,Qx,Qy),dp;
+    matrix C[4][4] = UpOneMatrix(4);
+    C[1,3] = q;
+    C[2,4] = q;
+    def ncQ = nc_algebra(C,0);
+    setring ncQ;
+    ncloc loc = intvec(2);
+    ncfrac frac = list([y^2+7*y+1,0,0,0], loc);
+    isZeroNcfrac(frac);
+    frac.lnum = 42*y*Qy+7*Qx+3*x+7;
+    isZeroNcfrac(frac);
+}
+//////////////////////////////////////////////////////////////////////
+proc isOneNcfrac(ncfrac frac)
+"USAGE:   isOneNcfrac(frac), ncfrac frac
+PURPOSE: checks if frac is one
+RETURN:  int, 1 if frac is one, 0 otherwise
+EXAMPLE: example isOneNcfrac; shows example"
+{
+    testNcfrac(frac);
+    if (hasLeftDenom(frac)) { // frac has left representation
+        return(frac.lden == frac.lnum);
+    } else { // frac has right representation, but no left representation
+        return(frac.rden == frac.rnum);
+    }
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring Q = (0,q),(x,y,Qx,Qy),dp;
+    matrix C[4][4] = UpOneMatrix(4);
+    C[1,3] = q;
+    C[2,4] = q;
+    def ncQ = nc_algebra(C,0);
+    setring ncQ;
+    ncloc loc = intvec(2);
+    ncfrac frac = list([y^2+7*y+1,y^2+7*y+1,0,0], loc);
+    isOneNcfrac(frac);
+    frac.lnum = 42*y*Qy+7*Qx+3*x+7;
+    isOneNcfrac(frac);
+}
+////////// initialization, comparison and declaration ////////////////
+static proc createNcfrac(def input) {
+    string inputType = typeof(input);
+    ncfrac result;
+    if (inputType == "list") {
+        if (size(input) >= 2) {
+            result = createNcfrac(input[1]);
+            if (typeof(input[2]) == "ncloc") {
+                result.loc = input[2];
+            }
+        }
+    }
+    if (inputType == "int" || inputType == "number" || inputType == "poly") {
+        result.lden = poly(1);
+        result.lnum = poly(input);
+        result.rden = poly(1);
+        result.rnum = poly(input);
+        result.loc = ncloc(list(1));
+    }
+    if (inputType == "vector") {
+        result.lden = input[1];
+        result.lnum = input[2];
+        result.rnum = input[3];
+        result.rden = input[4];
+    }
+    return(result);
+}
+//////////////////////////////////////////////////////////////////////
+proc zeroNcfrac(ncloc loc)
+"USAGE:   zeroNcfrac(loc), ncloc loc
+PURPOSE: returns the zero fraction in the localization loc
+RETURN:  ncfrac
+EXAMPLE: example zeroNcfrac; shows example"
+{
+    return(ncfrac(list([1,0,0,1], loc)));
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    ncloc loc = ideal(x-53,y-7);
+    zeroNcfrac(loc);
+}
+//////////////////////////////////////////////////////////////////////
+proc oneNcfrac(ncloc loc)
+"USAGE:   oneNcfrac(loc), ncloc loc
+PURPOSE: returns the one fraction in the localization loc
+RETURN:  ncfrac
+EXAMPLE: "
+{
+    return(ncfrac(list([1,1,1,1], loc)));
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    ncloc loc = ideal(x-42,y-17);
+    oneNcfrac(loc);
+}
+//////////////////////////////////////////////////////////////////////
+static proc trimNcfrac(ncfrac frac)
+"USAGE:   trimNcfrac(frac), ncfrac frac
+PURPOSE: simplifies complicated representations of zero and one
+RETURN:  ncfrac
+NOTE:    - if frac is zero, returns the default representation of zero
+         - if frac is one, returns the default representation of one
+         - otherwise returns frac
+EXAMPLE: example trimNcfrac; shows example"
+{
+    testNcfrac(frac);
+    if (isZeroNcfrac(frac)) {
+        return(zeroNcfrac(frac.loc));
+    }
+    if (isOneNcfrac(frac)) {
+        return(oneNcfrac(frac.loc));
+    }
+    return(frac);
+} example {
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    ncloc loc = poly(x);
+    loc;
+    ncfrac frac1 = list([x^2,x*x,0,0], loc);
+    trimNcfrac(frac1);
+    ncfrac frac2 = list([0,0,0,x], loc);
+    trimNcfrac(frac2);
+}
+//////////////////////////////////////////////////////////////////////
+static proc compareNcfracs(ncfrac frac1, ncfrac frac2) {
+    testNcfrac(frac1);
+    testNcfrac(frac2);
+    if (isZeroNcfrac(frac1)) {
+        return(isZeroNcfrac(frac2));
+    }
+    if (isOneNcfrac(frac1)) {
+        return(isOneNcfrac(frac2));
+    }
+    if (isZeroNcfrac(frac1 - frac2)) {
+        // frac1 and frac2 are equal iff their difference is zero
+        return (1);
+    } else {
+        return (0);
+    }
+}
+//////////////////////////////////////////////////////////////////////
+static proc invertedCompareNcfracs(ncfrac frac1, ncfrac frac2) {
+    return(!compareNcfracs(frac1, frac2));
+}
+////////// compute left resp. right representations //////////////////
+proc ensureLeftNcfrac(ncfrac frac)
+"USAGE:   ensureLeftNcfrac(frac), ncfrac frac
+PURPOSE: ensures that frac has a left representation (by computing it if not
+         alreaDy known)
+RETURN:  ncfrac, a representation of frac which has a left representation
+EXAMPLE: "
+{
+    testNcfrac(frac);
+    if (hasLeftDenom(frac)) {
+        return(frac);
+    }
+    ncloc loc = frac.loc;
+    vector f = toVector(frac);
+    vector result = convertRightToLeftFraction(f, loc.locType, loc.locData);
+    return(ncfrac(list(result, loc)));
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    // monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y;
+    list L = g1,g2;
+    ncloc loc0 = L;
+    poly g = g1^2*g2;
+    poly f = Dx;
+    ncfrac frac0 = [0,0,f,g];
+    frac0.loc = loc0;
+    ncfrac rm = ensureLeftNcfrac(frac0);
+    print(rm);
+    rm.lnum*g-rm.lden*f;
+    // geometric localization
+    ncloc loc1 = ideal(x-1,y-3);
+    f = Dx;
+    g = x^2+y;
+    ncfrac frac1 = [0,0,f,g];
+    frac1.loc = loc1;
+    ncfrac rg = ensureLeftNcfrac(frac1);
+    print(rg);
+    rg.lnum*g-rg.lden*f;
+    // rational localization
+    intvec rat = 1;
+    ncloc loc2 = rat;
+    f = Dx+Dy;
+    g = x;
+    ncfrac frac2 = [0,0,f,g];
+    frac2.loc = loc2;
+    ncfrac rr = ensureLeftNcfrac(frac2);
+    print(rr);
+    rr.lnum*g-rr.lden*f;
+}
+//////////////////////////////////////////////////////////////////////
+proc ensureRightNcfrac(ncfrac frac)
+"USAGE:   ensureLeftNcfrac(frac), ncfrac frac
+PURPOSE: ensures that frac has a right representation (by computing it if not
+         alreaDy known)
+RETURN:  ncfrac, a representation of frac which has a right representation
+EXAMPLE: "
+{
+    testNcfrac(frac);
+    if (hasRightDenom(frac)) {
+        return (frac);
+    }
+    ncloc loc = frac.loc;
+    vector f = toVector(frac);
+    vector result = convertLeftToRightFraction(f, loc.locType, loc.locData);
+    return(ncfrac(list(result, loc)));
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    // monoidal localization
+    poly g = x;
+    poly f = Dx;
+    ncloc loc0 = g;
+    ncfrac frac0 = [g,f,0,0];
+    frac0.loc = loc0;
+    ncfrac rm = ensureRightNcfrac(frac0);
+    print(rm);
+    f*rm.rden-g*rm.rnum;
+    // geometric localization
+    g = x+y;
+    f = Dx+Dy;
+    ncloc loc1 = ideal(x-1,y-3);
+    ncfrac frac1 = [g,f,0,0];
+    frac1.loc = loc1;
+    ncfrac rg = ensureRightNcfrac(frac1);
+    print(rg);
+    f*rg.rden-g*rg.rnum;
+    // rational localization
+    intvec rat = 1;
+    f = Dx+Dy;
+    g = x;
+    ncloc loc2 = rat;
+    ncfrac frac2 = [g,f,0,0];
+    frac2.loc = loc2;
+    ncfrac rr = ensureRightNcfrac(frac2);
+    print(rr);
+    f*rr.rden-g*rr.rnum;
+}
+////////// arithmetic ////////////////////////////////////////////////
+static proc addNcfracs(ncfrac frac1, ncfrac frac2) {
+    testNcfrac(frac1);
+    testNcfrac(frac2);
+    ncloc loc = frac1.loc;
+    if (loc != frac2.loc) {
+        ERROR("cannot add fractions: incompatible localizations");
+    }
+    vector f1 = toVector(frac1);
+    vector f2 = toVector(frac2);
+    vector result = addLeftFractions(f1, f2, loc.locType, loc.locData, 1);
+    return(ncfrac(list(result, loc)));
+}
+//////////////////////////////////////////////////////////////////////
+proc negateNcfrac(ncfrac frac)
+"USAGE:   negateNcfrac(frac), ncfrac frac
+PURPOSE: compute the negative (i.e. additive inverse) of frac
+RETURN:  ncfrac
+NOTE:    returns (-1)*frac
+EXAMPLE: example negateNcfrac; shows example"
+{
+    ncfrac result = frac;
+    result.lnum = result.lnum * (-1);
+    result.rnum = result.rnum * (-1);
+    return(result);
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    poly g = x*y^2+4*x+7*y-98;
+    ncloc loc = g;
+    ncfrac frac = list([g, 13*x^2], loc);
+    frac;
+    ncfrac negFrac = negateNcfrac(frac);
+    negFrac;
+    frac + negFrac;
+}
+//////////////////////////////////////////////////////////////////////
+static proc subtractNcfracs(ncfrac frac1, ncfrac frac2) {
+    testNcfrac(frac1);
+    testNcfrac(frac2);
+    if (frac1.loc != frac2.loc) {
+        ERROR("cannot subtract fractions: incompatible localizations");
+    }
+    return(frac1 + negateNcfrac(frac2));
+}
+//////////////////////////////////////////////////////////////////////
+static proc multiplyNcfracs(ncfrac frac1, ncfrac frac2) {
+    testNcfrac(frac1);
+    testNcfrac(frac2);
+    ncloc loc = frac1.loc;
+    if (loc != frac2.loc) {
+        ERROR("cannot multiply fractions: incompatible localizations");
+    }
+    vector f1 = toVector(frac1);
+    vector f2 = toVector(frac2);
+    vector result = multiplyLeftFractions(f1, f2, loc.locType, loc.locData, 1);
+    return(ncfrac(list(result, loc)));
+}
+//////////////////////////////////////////////////////////////////////
+proc isInvertibleNcfrac(ncfrac frac)
+"USAGE:   isInvertibleNcfrac(frac), ncfrac frac
+PURPOSE: checks if frac is invertible
+RETURN:  int, 1 if frac is invertible, 0 otherwise
+EXAMPLE: example isInvertibleNcfrac; shows example"
+{
+    testNcfrac(frac);
+    poly num;
+    if (hasLeftDenom(frac)) { // frac has left representation
+        num = frac.lnum;
+    } else { // frac has right representation, but no left representation
+        num = frac.rnum;
+    }
+    return(isDenom(num, frac.loc));
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    ncloc loc = intvec(2);
+    ncfrac frac = list([y,y+1,0,0], loc);
+    isInvertibleNcfrac(frac);
+    frac = list([y,x+1,0,0], loc);
+    isInvertibleNcfrac(frac);
+}
+//////////////////////////////////////////////////////////////////////
+proc invertNcfrac(ncfrac frac)
+"USAGE:   invertNcfrac(frac), ncfrac frac
+PURPOSE: compute the inverse of frac
+RETURN:  ncfrac
+NOTE:    returns the zero fraction if frac is not invertible
+EXAMPLE: example invertNcfrac; shows example"
+{
+    ncfrac result;
+    if (!isInvertibleNcfrac(frac)) {
+        result = zeroNcfrac(frac.loc);
+    } else {
+        result = list([frac.lnum, frac.lden, frac.rden, frac.rnum], frac.loc);
+    }
+    return(result);
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    ncloc loc = intvec(2);
+    ncfrac frac1 = list([y,y+1,0,0], loc);
+    // frac1 is invertible
+    ncfrac inv = invertNcfrac(frac1);
+    inv;
+    ncfrac frac2 = list([y,x+1,0,0], loc);
+    // frac2 is not invertible
+    inv = invertNcfrac(frac2);
+    inv;
+}
+//////////////////////////////////////////////////////////////////////
+static proc divideNcfracs(ncfrac frac1, ncfrac frac2) {
+    testNcfrac(frac1);
+    testNcfrac(frac2);
+    if (frac1.loc != frac2.loc) {
+        ERROR("cannot multiply fractions: incompatible localizations");
+    }
+    if(!isInvertibleNcfrac(frac2)) {
+        ERROR("division by non-invertible fraction");
+    } else {
+        return(frac1 * invertNcfrac(frac2));
+    }
+}
+//////////////////////////////////////////////////////////////////////
+////////// internal testing procedures ///////////////////////////////
+static proc testHasLeftHasRight()
+{
+    print("  testing ncfrac hasLeftDenom/hasRightDenom...");
+    ring r = 0,(x,y,Dx,Dy),dp;
+    def R = Weyl();
+    setring R;
+    ncloc loc = ideal(x-3,y+7);
+    ncfrac noLeft = list([0,0,3*y*Dx,x+2], loc);
+    if (hasLeftDenom(noLeft)) {
+        ERROR("hasLeftDenom noLeft failed");
+    }
+    ncfrac left = list([1,Dx,Dx,1], loc);
+    if (!hasLeftDenom(left)) {
+        ERROR("hasLeftDenom left failed");
+    }
+    ncfrac noRight = list([x+2,3*y*Dx,0,0], loc);
+    if (hasRightDenom(noRight)) {
+        ERROR("hasRightDenom noRight failed");
+    }
+    ncfrac right = list([1,Dx,Dx,1], loc);
+    if (!hasRightDenom(right)) {
+        ERROR("hasRightDenom right failed");
+    }
+    print("    ncloc hasLeftDenom/hasRightDenom OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testIsZeroIsOne()
+{
+    print("  testing ncfrac isZeroNcfrac/isOneNcfrac...");
+    ring Q = (0,q),(x,y,Qx,Qy),dp;
+    matrix C[4][4] = UpOneMatrix(4);
+    C[1,3] = q;
+    C[2,4] = q;
+    def ncQ = nc_algebra(C,0);
+    setring ncQ;
+    ncloc loc = intvec(2);
+    ncfrac frac = list([y^2+7*y+1,0,0,0], loc);
+    if (!isZeroNcfrac(frac)) {
+        ERROR("isZeroNcfrac zero failed");
+    }
+    frac.lnum = 42*y*Qy+7*Qx+3*x+7;
+    if (isZeroNcfrac(frac)) {
+        ERROR("isZeroNcfrac non-zero failed");
+    }
+    frac.lnum = frac.lden;
+    //frac = list([y^2+7*y+1,y^2+7*y+1,0,0], loc);
+    if (!isOneNcfrac(frac)) {
+        ERROR("isOneNcfrac one failed");
+    }
+    frac.lnum = 42*y*Qy+7*Qx+3*x+7;
+    if (isOneNcfrac(frac)) {
+        ERROR("isOneNcfrac non-one failed");
+    }
+    print("    ncloc isZeroNcfrac/isOneNcfrac OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testNcFracCreation() {
+    print("  testing ncfrac creation...");
+    ring r = 0,(x,y,Dx,Dy),dp;
+    def R = Weyl();
+    setring R;
+    ncloc loc = list(x);
+    ncfrac frac1 = [x,Dx,0,0]; // create from vector
+    frac1.loc = loc;
+    test(frac1);
+    ncfrac frac2 = 7; // create from interface
+    frac2.loc = loc;
+    test(frac2);
+    ncfrac frac3 = 4*x*Dx*Dy*Dy; // create from poly
+    frac3.loc = loc;
+    test(frac3);
+    ncfrac frac4 = list([x^2,Dx], loc); // create from list with vector
+    test(frac4);
+    ncfrac frac5 = list(42, loc); // create from list with int
+    test(frac5);
+    print("    ncfrac creation OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testNcfracComparison() {
+    print("  testing ncfrac comparison...");
+    ring r = 0,(x,y,Dx,Dy),dp;
+    def R = Weyl();
+    setring R;
+    // monoidal
+    poly g1 = x*y+3;
+    poly g2 = y^3;
+    ncloc locm = list(g1, g2);
+    ncfrac fracm1 = list([g1,Dx,0,0], locm);
+    ncfrac fracm2 = list([g1*g2,g2*Dx,0,0], locm);
+    ncfrac fracm3 = list([g1*g2,g1*Dx+3,0,0], locm);
+    if (!(fracm1 == fracm2)) {
+        ERROR("Weyl monoidal positive basic comparison error");
+    }
+    if (fracm1 == fracm3) {
+        ERROR("Weyl monoidal first negative basic comparison error");
+    }
+    if (fracm2 == fracm3) {
+        ERROR("Weyl monoidal second negative basic comparison error");
+    }
+    // geometric
+    ideal p = x+5, y-2;
+    ncloc locg = p;
+    ncfrac fracg1 = list([g1,Dx,0,0], locg);
+    ncfrac fracg2 = list([g1*g2,g2*Dx,0,0], locg);
+    ncfrac fracg3 = list([g1*g2,g1*Dx+3,0,0], locg);
+    if (!(fracg1 == fracg2)) {
+        ERROR("Weyl geometric positive basic comparison error");
+    }
+    if (fracg1 == fracg3) {
+        ERROR("Weyl geometric first negative basic comparison error");
+    }
+    if (fracg2 == fracg3) {
+        ERROR("Weyl geometric second negative basic comparison error");
+    }
+    // rational
+    intvec rat = 1,4;
+    ncloc locr = rat;
+    ncfrac fracr1 = list([x+Dy,Dx,0,0], locr);
+    ncfrac fracr2 = list([x*Dy*(x+Dy),x*Dx*Dy,0,0], locr);
+    ncfrac fracr3 = list([Dy*x*(x+Dy),x*Dx*Dy+1,0,0], locr);
+    if (!(fracr1 == fracr2)) {
+        ERROR("Weyl rational positive basic comparison error");
+    }
+    if (fracr1 == fracr3) {
+        ERROR("Weyl rational first negative basic comparison error");
+    }
+    if (fracr2 == fracr3) {
+        ERROR("Weyl rational second negative basic comparison error");
+    }
+    print("    ncfrac comparison OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testNcfracAddition() {
+    print("  testing ncfrac addition...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    poly g1 = x+3;
+    poly g2 = x*y+y;
+    ncloc loc = list(g1,g2);
+    ncfrac frac1 = list([g1,Dx,0,0], loc);
+    ncfrac frac2 = list([g2,Dy,0,0], loc);
+    ncfrac resu = frac1 + frac2;
+    if (resu.lden != g1*g2 || resu.lnum != g2*Dx+g1*Dy) {
+        ERROR("Weyl monoidal addition error");
+    }
+    loc = ideal(x-1,y-3);
+    frac1.loc = loc;
+    frac2.loc = loc;
+    resu = frac1 + frac2;
+    if (resu.lden != g1*g2 || resu.lnum != g2*Dx+g1*Dy) {
+        ERROR("Weyl geometric maximal addition error");
+    }
+    loc = ideal(y+3);
+    frac1.loc = loc;
+    frac2.loc = loc;
+    resu = frac1 + frac2;
+    if (resu.lden != g1*g2 || resu.lnum != g2*Dx+g1*Dy) {
+        ERROR("Weyl geometric prime addition error");
+    }
+    loc = intvec(2);
+    frac1 = list([y^2+y+1,Dx,0,0], loc);
+    frac2 = list([y-2,Dy,0,0], loc);
+    resu = frac1 + frac2;
+    if (resu.lden != (y^2+y+1)*(y-2) || resu.lnum != (y-2)*Dx+(y^2+y+1)*Dy) {
+        ERROR("Weyl rational addition error");
+    }
+    print("    ncfrac addition OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testNcfracSubtraction() {
+    print("  testing ncfrac subtraction...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    poly g1 = x+3;
+    poly g2 = x*y+y;
+    ncloc loc = list(g1,g2);
+    ncfrac frac1 = list([g1,Dx,0,0], loc);
+    ncfrac frac2 = list([g2,-Dy,0,0], loc);
+    ncfrac resu = frac1 - frac2;
+    if (resu.lden != g1*g2 || resu.lnum != g2*Dx+g1*Dy) {
+        ERROR("Weyl monoidal subtraction error");
+    }
+    loc = ideal(x-1,y-3);
+    frac1.loc = loc;
+    frac2.loc = loc;
+    resu = frac1 - frac2;
+    if (resu.lden != g1*g2 || resu.lnum != g2*Dx+g1*Dy) {
+        ERROR("Weyl geometric subtraction error");
+    }
+    loc = intvec(2);
+    frac1 = list([y^2+y+1,Dx,0,0], loc);
+    frac2 = list([y-2,-Dy,0,0], loc);
+    resu = frac1 - frac2;
+    if (resu.lden != (y^2+y+1)*(y-2) || resu.lnum != (y-2)*Dx+(y^2+y+1)*Dy) {
+        ERROR("Weyl rational subtraction error");
+    }
+    print("    ncfrac subtraction OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testNcfracMultiplication() {
+    print("  testing ncfrac multiplication...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    // monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y+y;
+    ncloc loc = list(g1,g2);
+    ncfrac frac1 = list([g1,Dx,0,0], loc);
+    ncfrac frac2 = list([g2,Dy,0,0], loc);
+    ncfrac resu = frac1 * frac2;
+    if (resu.lden != g1*g2^2 || resu.lnum != x*y*Dx*Dy+y*Dx*Dy-y*Dy) {
+        ERROR("Weyl monoidal multiplication error");
+    }
+    // geometric localization
+    loc = ideal(x-1,y-3);
+    frac1.loc = loc;
+    frac2.loc = loc;
+    resu = frac1 * frac2;
+    if (resu.lden != g1*g2*(x+1) || resu.lnum != x*Dx*Dy+Dx*Dy-Dy) {
+        ERROR("Weyl geometric multiplication error");
+    }
+    // rational localization
+    loc = intvec(2);
+    frac1 = list([y^2+y+1,Dx,0,0], loc);
+    frac2 = list([y-2,Dy,0,0], loc);
+    resu = frac1 * frac2;
+    if (resu.lden != (y^2+y+1)*(y-2) || resu.lnum != Dx*Dy) {
+        ERROR("Weyl rational multiplication (1*2) error");
+    }
+    resu = frac2 * frac1;
+    if (resu.lden != (y^2+y+1)^2*(y-2) || resu.lnum != y^2*Dx*Dy+y*Dx*Dy-2*y*Dx+Dx*Dy-Dx) {
+        ERROR("Weyl rational multiplication (2*1) error");
+    }
+    print("    ncfrac multiplication OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testNcfracDivision() {
+    print("  testing ncfrac division...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    // monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y+y;
+    ncloc loc = list(g1,g2);
+    ncfrac frac1 = list([g1,Dx,0,0], loc);
+    ncfrac frac2 = list([g2,g1,0,0], loc);
+    ncfrac resu = frac1 / frac2;
+    if (resu.lden != g1^3 || resu.lnum != x^2*y*Dx+4*x*y*Dx+3*y*Dx+2*y) {
+        ERROR("Weyl monoidal division error");
+    }
+    // geometric localization
+    loc = ideal(x-1,y-3);
+    frac1.loc = loc;
+    frac2.loc = loc;
+    resu = frac1 / frac2;
+    if (resu.lden != g1^3 || resu.lnum != x^2*y*Dx+4*x*y*Dx+3*y*Dx+2*y) {
+        ERROR("Weyl geometric division error");
+    }
+    // rational localization
+    loc = intvec(2);
+    frac1 = list([y^2+y+1,Dx,0,0], loc);
+    frac2 = list([y-2,3*y^2-21y+2,0,0], loc);
+    resu = frac1 / frac2;
+    if (resu.lden != 3*y^4-18*y^3-16*y^2-19*y+2 || resu.lnum != (y-2)*Dx) {
+        ERROR("Weyl geometric division error");
+    }
+    print("    ncfrac division OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testNcfracInversion() {
+    print("  testing ncfrac invertNcfrac...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    poly g1 = x+3;
+    poly g2 = x*y;
+    ncfrac frac = list([g1*g2, 17, 0, 0], ncloc(list(g1,g2)));
+    frac = invertNcfrac(frac);
+    if (frac.lden != 17 || frac.lnum != g1*g2) {
+        ERROR("Weyl monoidal inversion error");
+    }
+    frac = list([g1, 3*x, 0, 0], ncloc(ideal(x-1,y)));
+    frac = invertNcfrac(frac);
+    if (frac.lden != 3*x || frac.lnum != x+3) {
+        ERROR("Weyl geometric inversion error");
+    }
+    frac = list([g1*g2, y, 0, 0], ncloc(intvec(1,2)));
+    frac = invertNcfrac(frac);
+    if (frac.lden != y || frac.lnum != g1*g2) {
+        ERROR("Weyl rational inversion error");
+    }
+    print("    ncfrac invertNcfrac OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testEnsureRightNcfrac() {
+    print("  testing ncfrac ensureRightNcfrac...");
+    // Weyl
+    ring W = 0,(x,y,Dx,Dy),dp;
+    def ncW = Weyl();
+    setring ncW;
+    //// monoidal localization
+    ncloc monloc = list(x);
+    ncfrac monfrac = list([x,Dx,0,0], monloc);
+    ncfrac monresult = ensureRightNcfrac(monfrac);
+    if (!status(monresult)[1]) {
+        ERROR("Weyl monoidal ensureRightNcfrac failed");
+    }
+    //// geometrical localization
+    ncloc geoloc = ideal(x-1,y-3);
+    ncfrac geofrac = list([x+y,Dx+Dy,0,0], geoloc);
+    ncfrac georesult = ensureRightNcfrac(geofrac);
+    if (!status(georesult)[1]) {
+        ERROR("Weyl geometric ensureRightNcfrac failed");
+    }
+    //// rational localization
+    ncloc ratloc = intvec(1);
+    ncfrac ratfrac = list([x,Dx+Dy], ratloc);
+    ncfrac ratresult = ensureRightNcfrac(ratfrac);
+    if (!status(ratresult)[1]) {
+        ERROR("Weyl rational ensureRightNcfrac failed");
+    }
+    // shift rational localization
+    ring S = 0,(x,y,Sx,Sy),dp;
+    matrix D[4][4];
+    D[1,3] = Sx;
+    D[2,4] = Sy;
+    def ncS = nc_algebra(1, D);
+    setring ncS;
+    ncfrac shiftfrac = list([x,Sx+Sy], ratloc);
+    ncfrac shiftresult = ensureRightNcfrac(shiftfrac);
+    if (!status(shiftresult)[1]) {
+        ERROR("Shift rational ensureRightNcfrac failed");
+    }
+    // q-shift rational localization
+    ring Q = (0,q),(x,y,Qx,Qy),dp;
+    matrix C[4][4] = UpOneMatrix(4);
+    C[1,3] = q;
+    C[2,4] = q;
+    def ncQ = nc_algebra(C, 0);
+    setring ncQ;
+    ncfrac qshiftfrac = list([x,Qx+Qy], ratloc);
+    ncfrac qshiftresult = ensureRightNcfrac(qshiftfrac);
+    if (!status(qshiftresult)[1]) {
+        ERROR("q-shift rational ensureRightNcfrac failed");
+    }
+    print("    ncfrac ensureRightNcfrac OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testEnsureLeftNcfrac() {
+    print("  testing ncfrac ensureLeftNcfrac...");
+    // Weyl
+    ring W = 0,(x,y,Dx,Dy),dp;
+    def ncW = Weyl();
+    setring ncW;
+    //// monoidal localization
+    ncloc monloc = list(x+3, x*y);
+    ncfrac monfrac = list([0,0,Dx,(x+3)^2*x], monloc);
+    ncfrac monresult = ensureLeftNcfrac(monfrac);
+    if (!status(monresult)[1]) {
+        ERROR("Weyl monoidal ensureLeftNcfrac failed");
+    }
+    //// geometric localization
+    ncloc geoloc = ideal(x-1,y-3);
+    ncfrac geofrac = list([0,0,Dx,x^2+y], geoloc);
+    ncfrac georesult = ensureLeftNcfrac(geofrac);
+    if (!status(georesult)[1]) {
+        ERROR("Weyl monoidal ensureLeftNcfrac failed");
+    }
+    //// rational localization
+    ncloc ratloc = intvec(1);
+    ncfrac ratfrac = list([0,0,Dx+Dy,x], ratloc);
+    ncfrac ratresult = ensureLeftNcfrac(ratfrac);
+    if (!status(ratresult)[1]) {
+        ERROR("Weyl monoidal ensureLeftNcfrac failed");
+    }
+    // shift rational localization
+    ring S = 0,(x,y,Sx,Sy),dp;
+    matrix D[4][4];
+    D[1,3] = Sx;
+    D[2,4] = Sy;
+    def ncS = nc_algebra(1, D);
+    setring ncS;
+    ncfrac shiftfrac = list([0,0,Sx+Sy,x], ratloc);
+    ncfrac shiftresult = ensureLeftNcfrac(shiftfrac);
+    if (!status(shiftresult)[1]) {
+        ERROR("Shift rational ensureLeftNcfrac failed");
+    }
+    // q-shift rational localization
+    ring Q = (0,q),(x,y,Qx,Qy),dp;
+    matrix C[4][4] = UpOneMatrix(4);
+    C[1,3] = q;
+    C[2,4] = q;
+    def ncQ = nc_algebra(C, 0);
+    setring ncQ;
+    ncfrac qshiftfrac = list([0,0,Qx+Qy,x], ratloc);
+    ncfrac qshiftresult = ensureLeftNcfrac(qshiftfrac);
+    if (!status(qshiftresult)[1]) {
+        ERROR("q-shift rational ensureLeftNcfrac failed");
+    }
+    print("    ncfrac ensureLeftNcfrac OK");
+}
+//////////////////////////////////////////////////////////////////////
Index: Singular/LIB/ncloc.lib
===================================================================
--- Singular/LIB/ncloc.lib	(revision 3f58e0f44e7e559abc73c266dcd710d43ccfccb3)
+++ Singular/LIB/ncloc.lib	(revision 3f58e0f44e7e559abc73c266dcd710d43ccfccb3)
@@ -0,0 +1,361 @@
+//////////////////////////////////////////////////////////////////////
+version="version ncloc.lib 4.0.0.0 Dec_2017 "; // $Id$
+category="Noncommutative";
+info="
+LIBRARY:   ncloc.lib  Ore-localization in G-Algebras
+AUTHOR:    Johannes Hoffmann, email: johannes.hoffmann at math.rwth-aachen.de
+
+OVERVIEW:
+This library introduces a new type: ncloc.
+This type wraps the localization data defined as in olga.lib.
+An element of type ncloc has two members:
+  - int locType
+  - def locData
+
+OPERATIONS:
+string(ncloc);
+  give a string representation of the data describing the localization
+print(ncloc);
+  prints the string representation of the localization
+status(ncloc);
+  report on the status/validity of the localization
+test(ncloc);
+  check if the localization is valid
+
+INFIX OPERATIONS:
+ncloc == ncloc;
+  compare two localizations
+ncloc != ncloc;
+  compare two localizations
+ncloc = list/poly
+  create a monoidal localization from the given data
+ncloc = ideal
+  create a geometric localization from the given data
+ncloc = intvec
+  create a rational localization from the given data
+
+PROCEDURES:
+isDenom(poly, ncloc);
+  determine if a polynomial is in the denominator set of the localization
+testNcloc();
+  execute a series of internal testing procedures
+testNclocExamples();
+  execute the examples of all procedures in this library
+";
+//////////////////////////////////////////////////////////////////////
+proc testNclocExamples()
+"USAGE:   testNclocExamples()
+PURPOSE: execute the examples of all procedures in this library
+RETURN:  nothing
+NOTE:
+EXAMPLE: "
+{
+    example isDenom;
+}
+//////////////////////////////////////////////////////////////////////
+proc testNcloc()
+"USAGE:   testNcloc()
+PURPOSE: execute a series of internal testing procedures
+RETURN:  nothing
+NOTE:
+EXAMPLE: "
+{
+    print("testing ncloc.lib...");
+    testNclocCreation();
+    testNclocComparison();
+    testDenominatorDetection();
+    print("testing complete - ncloc.lib OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc mod_init() {
+    LIB "olga.lib";
+    /* new type: ncloc (non-commutative localization) */
+    newstruct("ncloc", "def locData, int locType");
+    // types: 0 - monoidal, 1 - geometric, 2 - rational
+    system("install", "ncloc", "string", nclocToString, 1);
+    system("install", "ncloc", "print", printNcloc, 1);
+    system("install", "ncloc", "=", createNcloc, 1);
+    system("install", "ncloc", "==", compareNclocs, 2);
+    system("install", "ncloc", "!=", invertedCompareNclocs, 2);
+    system("install", "ncloc", "test", testNcloc, 4);
+    system("install", "ncloc", "status", statusNcloc, 4);
+}
+////////// status and printing ///////////////////////////////////////
+static proc statusNcloc(ncloc loc) {
+    return(locStatus(loc.locType, loc.locData)); // locStatus from olga.lib
+}
+//////////////////////////////////////////////////////////////////////
+static proc testNcloc(ncloc loc) {
+    list stat = status(loc);
+    if(!stat[1]) {
+        ERROR(stat[2]);
+    } else {
+        return();
+    }
+}
+//////////////////////////////////////////////////////////////////////
+static proc nclocToString(ncloc loc) {
+    list stat = status(loc);
+    if (!stat[1]) {
+        return(stat[2]);
+    }
+    string result;
+    int locT = loc.locType;
+    int i;
+    if (locT == 0) {
+        list data = loc.locData;
+        result = result + "monoidal localization generated by:" + newline;
+        result = result + string(data[1]);
+        for (i = 2; i <= size(data); i++) {
+            result = result + ", " + string(data[i]);
+        }
+    }
+    if (locT == 1) {
+        ideal data = loc.locData;
+        result = result + "geometric localization at prime ideal:" + newline;
+        result = result + data;
+    }
+    if (locT == 2) {
+        intvec data = loc.locData;
+        result = result + "rational localization generated by:" + newline;
+        result = result + string(var(data[1]));
+        for (i = 2; i <= size(data); i++) {
+            result = result + ", " + string(var(data[i]));
+        }
+    }
+    if (result == "") {
+        result = "invalid localization";
+    }
+    return(result);
+}
+//////////////////////////////////////////////////////////////////////
+static proc printNcloc(ncloc loc) {
+    print(string(loc));
+}
+////////// initialization, comparison and declaration ////////////////
+static proc createNcloc(def input) {
+    string inputType = typeof(input);
+    ncloc result;
+    if (inputType == "list") {
+        result.locType = 0;
+        result.locData = normalizeMonoidal(input);
+    }
+    if (inputType == "poly") {
+        result.locType = 0;
+        result.locData = normalizeMonoidal(list(input));
+    }
+    if (inputType == "ideal") {
+        result.locType = 1;
+        result.locData = std(input);
+    }
+    if (inputType == "intvec") {
+        result.locType = 2;
+        result.locData = normalizeRational(input);
+    }
+    test(result);
+    return(result);
+}
+//////////////////////////////////////////////////////////////////////
+static proc compareNclocs(ncloc loc1, ncloc loc2) {
+    // only detects "naturally" isomorphic localizations
+    test(loc1);
+    test(loc2);
+    int locT = loc1.locType;
+    if (locT != loc2.locType) {
+        return(0); // different types of localizations
+    }
+    if (locT == 0) {
+        int i, j, found;
+        list L1 = loc1.locData;
+        list L2 = loc2.locData;
+        for (i = 1; i <= size(L1); i++) { // test if L1 subset of L2
+            found = 0;
+            for (j = 1; j <= size(L2); j++) {
+                if (L1[i] == L2[j]) {
+                    found = 1;
+                    break;
+                }
+            }
+            if (!found) {
+                return(0);
+            }
+        }
+        for (j = 1; j <= size(L2); j++) { // test if L2 subset of L1
+            found = 0;
+            for (i = 1; i <= size(L1); i++) {
+                if (L1[i] == L2[j]) {
+                    found = 1;
+                    break;
+                }
+            }
+            if (!found) {
+                return(0);
+            }
+        }
+        return(1);
+    }
+    if (locT == 1) {
+        ideal I1 = std(loc1.locData);
+        ideal I2 = std(loc2.locData);
+        if (size(NF(I1,I2)) != 0) { // test if I1 subset of I2
+            return(0);
+        }
+        if (size(NF(I2,I1)) != 0) { // test if I2 subset of I1
+            return(0);
+        }
+        return(1);
+    }
+    if (locT == 2) {
+        return(loc1.locData == loc2.locData); // direct comparison of intvec
+    }
+    ERROR("unreachable state reached"); // enforced by testNcloc
+}
+//////////////////////////////////////////////////////////////////////
+static proc invertedCompareNclocs(ncloc loc1, ncloc loc2) {
+    return(!(compareNclocs(loc1, loc2)));
+}
+//////////////////////////////////////////////////////////////////////
+proc isDenom(poly p, ncloc loc)
+"USAGE:   isDenom(p, loc), poly a, ncloc loc
+PURPOSE: check if p is a valid denominator in the localization loc
+RETURN:  int
+NOTE:    returns 1 or 0, depending whether p is a valid denominator
+EXAMPLE: example isDenom; shows examples"
+{
+    return(isInS(p, loc.locType, loc.locData, 1)); // isInS from olga.lib
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    // monoidal localization
+    ncloc loc;
+    poly g1 = x^2*y+x+2;
+    poly g2 = y^3+x*y;
+    list L = g1,g2;
+    loc = L;
+    poly g = g1^2*g2;
+    poly f = g - 1;
+    isDenom(g, loc);
+    isDenom(f, loc);
+    // geometrical localization
+    loc = ideal(x-1,y-3);
+    g = x^2+y-3;
+    f = (x-1)*g;
+    isDenom(g, loc);
+    isDenom(f, loc);
+    // rational localization
+    intvec v = 2;
+    loc = v;
+    g = y^5+17*y^2-4;
+    f = x*y;
+    isDenom(g, loc);
+    isDenom(f, loc);
+}
+//////////////////////////////////////////////////////////////////////
+////////// internal testing procedures ///////////////////////////////
+static proc testNclocCreation()
+{
+    print("  testing ncloc creation...");
+    ring r = 0,(x,y,Dx,Dy),dp;
+    def R = Weyl();
+    setring R;
+    ncloc loc1 = list(x*y,4*x); // create from list
+    ncloc loc2 = 3*Dx*Dy+2; // create from poly
+    ncloc loc3 = ideal(x-3,y2+y+1); // create from ideal
+    ncloc loc4 = intvec(1,2); // create from intvec
+    ncloc loc5 = intvec(3);
+    print("    ncloc creation OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testNclocComparison()
+{
+    print("  testing ncloc comparison...");
+    ring r = 0,(x,y,Dx,Dy),dp;
+    def R = Weyl();
+    setring R;
+    // monoidal localization
+    ncloc monloc1 = list(x,y-3);
+    ncloc monloc2 = list(x*(y-3)^2);
+    ncloc monloc3 = list(y-3);
+    if (!(monloc1 == monloc2)) {
+        ERROR("monoidal comparison failed");
+    }
+    if (!(monloc1 != monloc3)) {
+        ERROR("monoidal inverted comparison failed");
+    }
+    // geometric localization
+    ncloc geoloc1 = ideal(x*y-y+2, y^2-x-2);
+    ncloc geoloc2 = ideal(x*y-y+2, y^2-x-2, x^2+x+2*y-2);
+    ncloc geoloc3 = ideal(x*y-y+1, y^2-x-2);
+    if (!(geoloc1 == geoloc2)) {
+        ERROR("geometric comparison failed");
+    }
+    if (!(geoloc1 != geoloc3)) {
+        ERROR("geometric inverted comparison failed");
+    }
+    // rational localization
+    ncloc ratloc1 = intvec(1);
+    ncloc ratloc2 = intvec(1,1);
+    ncloc ratloc3 = intvec(2);
+    if (!(ratloc1 == ratloc2)) {
+        ERROR("rational comparison failed");
+    }
+    if (!(ratloc1 != ratloc3)) {
+        ERROR("rational inverted comparison failed");
+    }
+    print("    ncloc comparison OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testDenominatorDetection() {
+    print("  testing ncloc denominator detection...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    // monoidal localization
+    ncloc loc;
+    poly g1 = x^2*y+x+2;
+    poly g2 = y^3+x*y;
+    loc = list(g1,g2);
+    poly g = g1^2*g2;
+    poly f = g - 1;
+    if (!isDenom(g, loc)) {
+        ERROR("monoidal positive denominator detection error");
+    }
+    if (isDenom(f, loc)) {
+        ERROR("monoidal negative denominator detection error");
+    }
+    // geometric localization
+    loc = ideal(x-1,y-3);
+    g = x^2+y-3;
+    f = (x-1)*g;
+    if (!isDenom(g, loc)) {
+        ERROR("geometric maximal positive denominator detection error");
+    }
+    if (isDenom(f, loc)) {
+        ERROR("geometric maximal negative denominator detection error");
+    }
+    loc = ideal(x);
+    g = x*y^2-7*y+3;
+    f = x^2-3*x;
+    if (!isDenom(g, loc)) {
+        ERROR("geometric prime positive denominator detection error");
+    }
+    if (isDenom(f, loc)) {
+        ERROR("geometric prime negative denominator detection error");
+    }
+    // rational localization
+    loc = intvec(2);
+    g = y^5+17*y^2-4;
+    f = x*y;
+    if (!isDenom(g, loc)) {
+        ERROR("rational positive denominator detection error");
+    }
+    if (isDenom(f, loc)) {
+        ERROR("rational negative denominator detection error");
+    }
+    print("    ncloc denominator detection OK");
+}
+//////////////////////////////////////////////////////////////////////
Index: Singular/LIB/olga.lib
===================================================================
--- Singular/LIB/olga.lib	(revision 3f58e0f44e7e559abc73c266dcd710d43ccfccb3)
+++ Singular/LIB/olga.lib	(revision 3f58e0f44e7e559abc73c266dcd710d43ccfccb3)
@@ -0,0 +1,1879 @@
+////////////////////////////////////////////////////////////////
+version="version olga.lib 4.0.0.0 Dec_2017 "; // $Id$
+category="Noncommutative";
+info="
+LIBRARY:   olga.lib  Ore-localization in G-Algebras
+AUTHOR:    Johannes Hoffmann, email: johannes.hoffmann at math.rwth-aachen.de
+
+OVERVIEW:
+Let A be a G-algebra.
+
+Current localization types:
+Type 0: monoidal
+  - represented by a list of polys g_1,...,g_k that have to be contained in a
+    commutative polynomial subring of A generated by a subset of the variables
+    of A
+Type 1: geometric
+  - only for algebras with an even number of variables where the first half
+    induces a commutative polynomial subring B of A
+  - represented by an ideal p, which has to be a prime ideal in B
+Type 2: rational
+  - represented by an intvec v = [i_1,...,i_k] in the range 1..nvars(basering)
+
+Localization data is an int specifying the type and a def with the
+corresponding information.
+
+A fraction is represented as a vector with four entries: [s,r,p,t]
+Here, s^{-1}r is the left fraction representation, pt^{-1} is the right one.
+If s or t is zero, it means that the corresponding representation is not set.
+If both are zero, the fraction is not valid.
+
+A detailed description along with further examples can be found in our paper:
+Johannes Hoffmann, Viktor Levandovskyy:
+    Constructive Arithmetics in Ore Localizations of Domains
+    https://arxiv.org/abs/1712.01773
+
+PROCEDURES:
+locStatus(int, def);
+  report on the status/validity of the given localization data
+testLocData(int, def);
+  check if the given data specifies a denominator set
+isInS(poly, int, def);
+  determine if a polynomial is in a denominator set
+fracStatus(vector frac, int locType, def locData);
+  report on the status/validity of the given fraction wrt. to the given
+  localization data
+testFraction(vector, int, def);
+  check if the given vector is a representation of a fraction in the specified
+  localization
+leftOre(poly, poly, int, def)
+  compute left Ore data
+rightOre(poly, poly, int, def)
+  compute right Ore data
+convertRightToLeftFraction(vector, int, def);
+  determine a left fraction representation of a given fraction
+convertLeftToRightFraction(vector, int, def);
+  determine a right fraction representation of a given fraction
+addLeftFractions(vector, vector, int, def);
+  add two left fractions in the specified localization
+multiplyLeftFractions(vector, vector, int, def);
+  multiply two left fractions in the specified localization
+areEqualLeftFractions(vector, vector, int, def);
+  check if two given fractions are equal
+isInvertibleLeftFraction(vector, int, def);
+  check if a fraction is invertible in the specified localization
+  (NOTE: check description for specific behaviour)
+invertLeftFraction(vector, int, def);
+  invert a fraction in the specified localization
+  (NOTE: check description for specific behaviour)
+isZeroFraction(vector);
+  determine if the given fraction is equal to zero
+isOneFraction(vector);
+  determine if the given fraction is equal to one
+normalizeMonoidal(list);
+  determine a normal form for monoidal localization data
+normalizeRational(intvec);
+  determine a normal form for rational localization data
+testOlga();
+  execute a series of internal testing procedures
+testOlgaExamples();
+  execute the examples of all procedures in this library
+";
+
+LIB "dmodloc.lib"; // for polyVars
+LIB "ncpreim.lib";
+LIB "elim.lib";
+LIB "ncalg.lib";
+//////////////////////////////////////////////////////////////////////
+proc testOlgaExamples()
+"USAGE:   testOlgaExamples()
+PURPOSE: execute the examples of all procedures in this library
+RETURN:  nothing
+NOTE:
+EXAMPLE: "
+{
+    example isInS;
+    example leftOre;
+    example rightOre;
+    example convertRightToLeftFraction;
+    example convertLeftToRightFraction;
+    example addLeftFractions;
+    example multiplyLeftFractions;
+    example areEqualLeftFractions;
+    example isInvertibleLeftFraction;
+    example invertLeftFraction;
+}
+//////////////////////////////////////////////////////////////////////
+proc testOlga()
+"USAGE:   testOlga()
+PURPOSE: execute a series of internal testing procedures
+RETURN:  nothing
+NOTE:
+EXAMPLE: "
+{
+    print("testing olga.lib...");
+    testIsInS();
+    testLeftOre();
+    testRightOre();
+    testAddLeftFractions();
+    testMultiplyLeftFractions();
+    testAreEqualLeftFractions();
+    testConvertLeftToRightFraction();
+    testConvertRightToLeftFraction();
+    testIsInvertibleLeftFraction();
+    testInvertLeftFraction();
+    print("testing complete - olga.lib OK");
+}
+//////////////////////////////////////////////////////////////////////
+proc locStatus(int locType, def locData)
+"USAGE:   locStatus(locType, locData), int locType, list/vector/intvec locData
+PURPOSE: determine the status of a set of localization data
+ASSUME:
+RETURN:  list
+NOTE:    - the first entry is 0 or 1, depending whether the input represents
+           a valid localization
+         - the second entry is a string with a status/error message
+EXAMPLE: example locStatus; shows example"
+{
+    int i;
+    if (locType < 0 || locType > 2) {
+        string invalidTypeString = "invalid localization: type is "
+          + string(locType) + ", valid types are:";
+        invalidTypeString = invalidTypeString + newline
+          + "0 for a monoidal localization";
+        invalidTypeString = invalidTypeString + newline
+          + "1 for a geometric localization";
+        invalidTypeString = invalidTypeString + newline
+          + "2 for a rational localization";
+        return(list(0, invalidTypeString));
+    }
+    string t = typeof(locData);
+    if (t == "none") {
+        return(list(0, "uninitialized or invalid localization:"
+          + " locData has to be defined"));
+    }
+    if (locType == 0) { // monoidal localizations
+        if (t != "list") {
+            return(list(0, "for a monoidal localization, locData has to be of"
+              + " type list, but is of type " + t));
+        } else { // locData is of type list
+            if (size(locData) == 0) {
+                return(list(0, "for a monoidal localization, locData has to be"
+                  + " a non-empty list"));
+            } else { // locData is of type list and has at least one entry
+                ideal listEntries;
+                for (i = 1; i <= size(locData); i++) {
+                    t = typeof(locData[i]);
+                    if (t != "poly" && t != "int" && t != "number") {
+                        return(list(0, "for a monoidal localization, locData"
+                          + " has to be a list of polys, ints or numbers, but"
+                          + " entry " + string(i) + " is " + string(locData[i])
+                          + ", which is of type " + t));
+                    } else {
+                        if (listEntries == 0) {
+                            listEntries = locData[i];
+                        } else {
+                            listEntries = listEntries, locData[i];
+                        }
+                    }
+                }
+                // locData is of type list, has at least one entry and all
+                //   entries are polys
+                if (!inducesCommutativeSubring(listEntries)) {
+                    return(list(0, "for a monoidal localization, the variables"
+                      + " occurring in the polys in locData have to induce a"
+                      + " commutative polynomial subring of basering"));
+                }
+            }
+        }
+    }
+    if (locType == 1) { // geometric localizations
+        int n = nvars(basering) div 2;
+        if (2*n != nvars(basering)) {
+            return(list(0, "for a geometric localization, basering has to have"
+              + " an even number of variables"));
+        } else {
+            int j;
+            for (i = 1; i <= n; i++) {
+                for (j = i + 1; j <= n; j++) {
+                    if (var(i)*var(j) != var(j)*var(i)) {
+                        return(list(0, "for a geometric localization, the"
+                          + " first half of the variables of basering has to"
+                          + " induce a commutative polynomial subring of"
+                          + " basering"));
+                    }
+                }
+            }
+        }
+        if (t != "ideal") {
+            return(list(0, "for a geometric localization, locData has to be of"
+              + " type ideal, but is of type " + t));
+        }
+        for (i = 1; i <= size(locData); i++) {
+            if (!polyVars(locData[i],1..n)) {
+                return(list(0, "for a geometric localization, locData has to"
+                + " be an ideal generated by polynomials containing only"
+                + " variables from the first half of the variables"));
+            }
+        }
+    }
+    if (locType == 2) { // rational localizations
+        if (t != "intvec") {
+            return(list(0, "for a rational localization, locData has to be of"
+              + " type intvec, but is of type " + t));
+        } else { // locData is of type intvec
+            if(locData == 0) {
+                return(list(0, "for a rational localization, locData has to be"
+                  + " a non-zero intvec"));
+            } else { // locData is of type intvec and not zero
+                if (!admissibleSub(locData)) {
+                    return(list(0, "for a rational localization, the variables"
+                      + " indexed by locData have to generate a sub-G-algebra"
+                      + " of the basering"));
+                }
+            }
+        }
+    }
+    return(list(1, "valid localization"));
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    locStatus(42, list(1));
+    def undef;
+    locStatus(0, undef);
+    string s;
+    locStatus(0, s);
+    list L;
+    locStatus(0, L);
+    L = s;
+    print(L);
+    locStatus(0, L);
+    ring w = 0,(x,Dx,y,Dy),dp;
+    def W = Weyl(1);
+    setring W;
+    W;
+    locStatus(0, list(x, Dx));
+    ring R;
+    setring R;
+    R;
+    locStatus(1, s);
+    setring W;
+    locStatus(1, s);
+    ring t = 0,(x,y,Dx,Dy),dp;
+    def T = Weyl();
+    setring T;
+    T;
+    locStatus(1, s);
+    locStatus(1, ideal(Dx));
+    locStatus(2, s);
+    intvec v;
+    locStatus(2, v);
+    locStatus(2, intvec(1,2));
+}
+//////////////////////////////////////////////////////////////////////
+proc testLocData(int locType, def locData)
+"USAGE:   testLocData(locType, locData), int locType,
+         list/vector/intvec locData
+PURPOSE: test if the given data specifies a denominator set wrt. the checks
+         from locStatus
+ASSUME:
+RETURN:  nothing
+NOTE:    throws error if checks were not successful
+EXAMPLE: example testLocData; shows examples"
+{
+    list stat = locStatus(locType, locData);
+    if (!stat[1]) {
+        ERROR(stat[2]);
+    } else {
+        return();
+    }
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R; setring R;
+    testLocData(0, list(1)); // correct localization, no error
+    testLocData(42, list(1)); // incorrect localization, results in error
+}
+//////////////////////////////////////////////////////////////////////
+proc isInS(poly p, int locType, def locData, list #)
+"USAGE:   isInS(p, locType, locData(, override)), poly p, int locType,
+         list/vector/intvec locData(, int override)
+PURPOSE: determine if a polynomial is in a denominator set
+ASSUME:
+RETURN:  int
+NOTE:    - returns 0 or 1, depending whether p is in the denominator set
+           specified by locType and locData
+         - if override is set, will not normalize locData (use with care)
+EXAMPLE: example isInS; shows examples"
+{
+    testLocData(locType, locData);
+    if (p == 0) { // the zero polynomial is never a valid denominator
+        return(0);
+    }
+    if (number(p) != 0) {
+        // elements of the coefficient field are always valid denominators
+        return(1);
+    }
+    int override;
+    if (size(#) > 0) {
+        if(typeof(#[1]) == "int") {
+            override = #[1];
+        }
+    }
+    if (locType == 0) {
+        ideal pFactors = commutativeFactorization(p, 1);
+        list locFactors;
+        if (override) {
+            locFactors = locData;
+        } else {
+            locFactors = normalizeMonoidal(locData);
+        }
+        int i, j, foundFactor;
+        for (i = 1; i <= size(pFactors); i++) {
+            foundFactor = 0;
+            for (j = 1; j <= size(locFactors); j++) {
+                if (pFactors[i] == locFactors[j]) {
+                    foundFactor = 1;
+                    break;
+                }
+            }
+            if (!foundFactor) {
+                return(0);
+            }
+        }
+        return(1);
+    }
+    if (locType == 1) {
+        int n = nvars(basering) div 2;
+        ideal I = p;
+        ideal J = std(eliminateNC(I, (n+1)..(2*n)));
+        ideal K = std(intersect(J, locData));
+        int i;
+        for (i = 1; i <= size(J); i++) {
+            if (NF(J[i], K) != 0) {
+                return(1);
+            }
+        }
+    }
+    if (locType == 2) {
+        if (!override) {
+            locData = normalizeRational(locData);
+        }
+        ideal I = p;
+        intvec modLocData = intvecComplement(locData, 1..nvars(basering));
+        I = eliminateNC(I, modLocData);
+        if (size(I)) {
+            return(1);
+        }
+    }
+    return(0);
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    // monoidal localization
+    poly g1 = x^2*y+x+2;
+    poly g2 = y^3+x*y;
+    list L = g1,g2;
+    poly g = g1^2*g2;
+    poly f = g-1;
+    isInS(g, 0, L); // g is in the denominator set
+    isInS(f, 0, L); // f is NOT in the denominator set
+    // geometric localization
+    ideal p = x-1, y-3;
+    g = x^2+y-3;
+    f = (x-1)*g;
+    isInS(g, 1, p); // g is in the denominator set
+    isInS(f, 1, p); // f is NOT in the denominator set
+    // rational localization
+    intvec v = 2;
+    g = y^5+17*y^2-4;
+    f = x*y;
+    isInS(g, 2, v); // g is in the denominator set
+    isInS(f, 2, v); // f is NOT in the denominator set
+}
+//////////////////////////////////////////////////////////////////////
+proc fracStatus(vector frac, int locType, def locData)
+"USAGE:   fracStatus(frac, locType, locData), vector frac, int locType,
+         list/intvec/vector locData
+PURPOSE: determine if the given vector is a representation of a fraction in the
+         specified localization
+ASSUME:
+RETURN:  list
+NOTE:    - the first entry is 0 or 1, depending whether the input is valid
+         - the second entry is a string with a status message
+EXAMPLE: example fracStatus; shows examples"
+{
+    list locStat = locStatus(locType, locData);
+    if (!locStat[1]) { // there is a problem with the localization data
+        return(list(0, "invalid localization in fraction: "+ string(frac)
+          + newline + "      " + locStat[2]));
+    } else { // the specified localization is valid
+        if ((frac[1] == 0) && (frac[4] == 0)) {
+            return(list(0, "vector is not a valid fraction: no denominator"
+              + " specified in " + string(frac)));
+        }
+        if (frac[1] != 0) { // frac has a left representation
+            if (!isInS(frac[1], locType, locData)) {
+                return(list(0, "the left denominator " + string(frac[1])
+                  + " of fraction " + string(frac) + " is not in the"
+                  + " denominator set of type " + string(locType) + " given by "
+                  + string(locData)));
+            }
+        }
+        if (frac[4] != 0) { // frac has a right representation
+            if (!isInS(frac[4], locType, locData)) {
+                return(list(0, "the right denominator " + string(frac[4])
+                  + " of fraction " + string(frac) + " is not in the"
+                  + " denominator set of type " + string(locType) + " given by "
+                  + string(locData)));
+            }
+        }
+        if ((frac[1] != 0) && (frac[4] != 0)) {
+            // frac has left and right representations
+            if (frac[2]*frac[4] != frac[1]*frac[3]) {
+                // the representations are not equal
+                return(list(0, "left and right representation are not equal in:"
+                  + string(frac)));
+            }
+        }
+    }
+    return(list(1, "valid fraction"));
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring r = QQ[x,y,Dx,Dy];
+    def R = Weyl();
+    setring R;
+    fracStatus([1,0,0,0], 42, list(1));
+    list L = x;
+    fracStatus([0,7,x,0], 0, L);
+    fracStatus([Dx,Dy,0,0], 0, L);
+    fracStatus([0,0,Dx,Dy], 0, L);
+    fracStatus([x,Dx,Dy,x], 0, L);
+    fracStatus([x,Dx,x*Dx+2,x^2], 0, L);
+}
+//////////////////////////////////////////////////////////////////////
+proc testFraction(vector frac, int locType, def locData)
+"USAGE:   testFraction(frac, locType, locData), vector frac, int locType,
+         list/intvec/vector locData
+PURPOSE: test if the given vector is a representation of a fraction in the
+         specified localization wrt. the checks from fracStatus
+ASSUME:
+RETURN:  nothing
+NOTE:    throws error if checks were not successful
+EXAMPLE: example testFraction; shows examples"
+{
+    list stat = fracStatus(frac, locType, locData);
+    if (!stat[1]) {
+        ERROR(stat[2]);
+    } else {
+        return();
+    }
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring r = QQ[x,y,Dx,Dy];
+    def R = Weyl();
+    setring R;
+    list L = x;
+    vector frac = [x,Dx,x*Dx+2,x^2];
+    testFraction(frac, 0, L); // correct localization, no error
+    frac = [x,Dx,x*Dx,x^2];
+    testFraction(frac, 0, L); // incorrect localization, results in error
+}
+//////////////////////////////////////////////////////////////////////
+proc leftOre(poly s, poly r, int locType, def locData)
+"USAGE:   leftOre(s, r, locType, locData), poly s, r, int locType,
+         list/vector/intvec locData
+PURPOSE: compute left Ore data for a given tuple (s,r)
+ASSUME:  s is in the denominator set determined via locType and locData
+RETURN:  list
+NOTE:  - the first entry of the list is a vector [ts,tr] such that ts*r=tr*s
+       - the second entry of the list is a description of all choices for ts
+EXAMPLE: example leftOre; shows examples"
+{
+    testLocData(locType,locData);
+    locData = normalizeLocalization(locType, locData);
+    if(!isInS(s, locType, locData)) {
+        ERROR("cannot find Ore-parameter since poly " + string(s)
+          + " is not in the denominator set");
+    }
+    return(ore(s, r, locType, locData, 0));
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    // left Ore
+    // monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y;
+    list L = g1,g2;
+    poly g = g1^2*g2;
+    poly f = Dx;
+    list rm = leftOre(g, f, 0, L);
+    print(rm[1]);
+    rm[2];
+    rm[1][2]*g-rm[1][1]*f;
+    // geometric localization
+    ideal p = x-1, y-3;
+    f = Dx;
+    g = x^2+y;
+    list rg = leftOre(g, f, 1, p);
+    print(rg[1]);
+    rg[2];
+    rg[1][2]*g-rg[1][1]*f;
+    // rational localization
+    intvec rat = 1;
+    f = Dx+Dy;
+    g = x;
+    list rr = leftOre(g, f, 2, rat);
+    print(rr[1]);
+    rr[2];
+    rr[1][2]*g-rr[1][1]*f;
+}
+//////////////////////////////////////////////////////////////////////
+proc rightOre(poly s, poly r, int locType, def locData)
+"USAGE:   rightOre(s, r, locType, locData), poly s, r, int locType,
+         list/vector/intvec locData
+PURPOSE: compute right Ore data for a given tuple (s,r)
+ASSUME:  s is in the denominator set determined via locType and locData
+RETURN:  list
+NOTE:  - the first entry of the list is a vector [ts,tr] such that r*ts=s*tr
+       - the second entry of the list is a description of all choices for ts
+EXAMPLE: example rightOre; shows examples"
+{
+    testLocData(locType,locData);
+    locData = normalizeLocalization(locType, locData);
+    if(!isInS(s, locType, locData)) {
+        ERROR("cannot find Ore-parameter since poly " + string(s)
+          + " is not in the denominator set");
+    }
+    def bsRing = basering;
+    if (locType == 0) {
+        ideal modLocData = locData[1..size(locData)];
+    }
+    def oppRing = opposite(bsRing);
+    setring oppRing;
+    if (locType == 0) {
+        ideal oppModLocData = oppose(bsRing, modLocData);
+        list oppLocData = oppModLocData[1..size(oppModLocData)];
+    }
+    if (locType == 1) {
+        ideal oppLocData = oppose(bsRing,locData);
+    }
+    if (locType == 2) {
+        intvec oppLocData = locData;
+    }
+    poly oppS = oppose(bsRing, s);
+    poly oppR = oppose(bsRing, r);
+    list oppResult = ore(oppS, oppR, locType, oppLocData, 1);
+    vector oppOreParas = oppResult[1];
+    ideal oppJ = oppResult[2];
+    setring bsRing;
+    return(list(oppose(oppRing, oppOreParas), oppose(oppRing, oppJ)));
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    // monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y;
+    list L = g1,g2;
+    poly g = g1^2*g2;
+    poly f = Dx;
+    list rm = rightOre(g, f, 0, L);
+    print(rm[1]);
+    rm[2];
+    g*rm[1][2]-f*rm[1][1];
+    // geometric localization
+    ideal p = x-1, y-3;
+    f = Dx;
+    g = x^2+y;
+    list rg = rightOre(g, f, 1, p);
+    print(rg[1]);
+    rg[2];
+    g*rg[1][2]-f*rg[1][1];
+    // rational localization
+    intvec rat = 1;
+    f = Dx+Dy;
+    g = x;
+    list rr = rightOre(g, f, 2, rat);
+    print(rr[1]);
+    rr[2];
+    g*rr[1][2]-f*rr[1][1];
+}
+//////////////////////////////////////////////////////////////////////
+static proc ore(poly s, poly r, int locType, def locData, int rightOre) {
+    // TODO: once ringlist-bug [#577] is PROPERLY fixed (probability: low),
+    //   use eliminateNC instead of eliminate in rightOre-rational and
+    //   rightOre-geometric
+    // problem: the ordering on the opposite algebra is not recognized as
+    //   valid for eliminateNC, even though it is
+    if (r == 0) {
+        return([1,0], ideal(1));
+    }
+    int i, j;
+    // modification for monoidal localization
+    if (locType == 0) {
+        // computations will be carried out in the localization at the product
+        //   of all irreducible factors of s, thus all factors have to be
+        //   raised to the highest power occurring in the factorization of s
+        list factorList = commutativeFactorization(s);
+        ideal factors = factorList[1];
+        int maxPow = 0;
+        intvec exponents = factorList[2];
+        for (i = 1; i <= size(exponents); i++) {
+            if (exponents[i] > maxPow) {
+                maxPow = exponents[i];
+            }
+        }
+        int expDiff;
+        for (i = 1; i <= size(exponents); i++) {
+            expDiff = maxPow - exponents[i];
+            s = s * factors[i]^expDiff;
+            r = r * factors[i]^expDiff;
+        }
+    }
+    // compute kernel of x maps to x*r+R*s
+    ideal J = modulo(r, s);
+    ideal I = std(J);
+    // compute the poly in S
+    poly ts;
+    if (locType == 0) { // type 0: monoidal localization
+        intvec sle = leadexp(s);
+        intvec rle;
+        int m = deg(r) + 1; // upper bound
+        int dividesLeadingMonomial;
+        // find the minimal possible exponent m such that lm(f_i)|lm(s)^m
+        for (i = 1; i <= size(I); i++) {
+            rle = leadexp(I[i]);
+            dividesLeadingMonomial = 1;
+            for (j = 1; j <= size(rle); j++) {
+                if (sle[j] == 0 && rle[j] != 0) {
+                    // lm(f_i) divides no power of lm(s)
+                    dividesLeadingMonomial = 0;
+                    break;
+                }
+            }
+            if (dividesLeadingMonomial) {
+                int mf = 1;
+                while (1) {
+                    if (mf * sle - rle >= 0) {
+                        break; // now lm(r_i) divides lm(s)^mf
+                    } else {
+                        mf++;
+                    }
+                }
+                if (mf < m) {
+                    m = mf; // mf is lower than the current minimum m
+                }
+            }
+        }
+        poly nf = NF(s^m, I);
+        for (; m <= deg(r) + 1; m++) { // sic: m is initialized beforehand
+            if (nf == 0) { // s^m is in I
+                ts = s^m;
+                ideal Js = s^m;
+                break;
+            } else {
+                nf = NF(s*nf, I);
+            }
+        }
+    }
+    if (locType == 1) { // type 1: geometric localization
+        int n = nvars(basering) div 2;
+        if (rightOre == 1) {
+            // in the right Ore setting, the order of variables is inverted
+            poly elimVars = 1;
+            for (i = 1; i <= n; i++) {
+                elimVars = elimVars * var(i);
+            }
+            J = eliminate(I, elimVars);
+        } else {
+            J = eliminateNC(I, intvec((n+1)..(2*n)));
+        }
+        J = std(J);
+        // compute ts, the generator of J with the smallest degree that does
+        //   not vanish at locData
+        ideal K = intersect(J, locData);
+        K = std(K);
+        poly h, cand;
+        ideal Js;
+        ideal Q = std(locData);
+        for (i = 1; i <= size(J); i++) {
+            h = NF(J[i],Q);
+            if (h != 0) {
+                cand = NF(J[i],K);
+                Js = Js + cand;
+                if (ts == 0 || deg(cand) < deg(ts)) {
+                    ts = cand;
+                }
+            }
+        }
+    }
+    if (locType == 2) { // type 2: rational localization
+        int n = nvars(basering);
+        if (size(locData) < n) { // there are variables to eliminate
+            intvec modLocData;
+            int check;
+            modLocData = intvecComplement(locData, 1..n);
+            // calculate modLocData = {1...n}\locData
+            if (rightOre == 1) {
+                // in the right Ore setting, the order of variables is inverted
+                poly elimVars = 1;
+                for (i = 1; i <= size(locData); i++) {
+                    elimVars = elimVars * var(locData[i]);
+                }
+                J = eliminate(I, elimVars);
+            } else {
+                J = eliminateNC(I, modLocData);
+            }
+
+        } else { // no variables to eliminate (total localization)
+            J = I;
+        }
+        J = std(J);
+        ts = J[1];
+        for (i = 2; i <= size(J); i++) {
+            if (deg(J[i]) < deg(ts)) {
+                ts = J[i]; // choose generator with lowest total degree
+            }
+        }
+        ideal Js = J;
+    }
+    // calculate the other poly
+    poly tr = division(ts*r, s)[1][1,1];
+    if (ts == 0 || tr*s-ts*r != 0) {
+        string s;
+        if (rightOre) {
+            s = "right";
+        } else {
+            s = "left";
+        }
+        ERROR("no " + s + " Ore data could be found for s=" + string(s)
+          + " and r=" + string(r));
+    }
+    vector oreParas = [ts,tr];
+    list result = oreParas, Js;
+    return (result);
+}
+//////////////////////////////////////////////////////////////////////
+proc convertRightToLeftFraction(vector frac, int locType, def locData)
+"USAGE:   convertRightToLeftFraction(frac, locType, locData),
+         vector frac, int locType, list/vector/intvec locData
+PURPOSE: determine a left fraction representation of a given fraction
+ASSUME:
+RETURN:  vector
+NOTE:    - the returned vector contains a repr. of frac as a left fraction
+         - if the left representation of frac is already specified,
+           frac will be returned.
+EXAMPLE: example convertRightToLeftFraction; shows examples"
+{
+    testLocData(locType, locData);
+    testFraction(frac, locType, locData);
+    if (frac[1] != 0) { // frac already has a left representation
+        return (frac);
+    } else { // frac has no left representation, but a right one
+        vector oreParas = leftOre(frac[4], frac[3], locType, locData)[1];
+        vector result = [oreParas[1], oreParas[2], frac[3], frac[4]];
+        return (result);
+    }
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    // monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y;
+    list L = g1,g2;
+    poly g = g1^2*g2;
+    poly f = Dx;
+    vector fracm = [0,0,f,g];
+    vector rm = convertRightToLeftFraction(fracm, 0, L);
+    print(rm);
+    rm[2]*g-rm[1]*f;
+    // geometric localization
+    ideal p = x-1, y-3;
+    f = Dx;
+    g = x^2+y;
+    vector fracg = [0,0,f,g];
+    vector rg = convertRightToLeftFraction(fracg, 1, p);
+    print(rg);
+    rg[2]*g-rg[1]*f;
+    // rational localization
+    intvec rat = 1;
+    f = Dx+Dy;
+    g = x;
+    vector fracr = [0,0,f,g];
+    vector rr = convertRightToLeftFraction(fracr, 2, rat);
+    print(rr);
+    rr[2]*g-rr[1]*f;
+}
+//////////////////////////////////////////////////////////////////////
+proc convertLeftToRightFraction(vector frac, int locType, def locData)
+"USAGE:   convertLeftToRightFraction(frac, locType, locData), vector frac,
+         int locType, list/vector/intvec locData
+PURPOSE: determine a right fraction representation of a given fraction
+ASSUME:
+RETURN:  vector
+NOTE:    - the returned vector contains a repr. of frac as a right fraction,
+         - if the right representation of frac is already specified,
+           frac will be returned.
+EXAMPLE: example convertLeftToRightFraction; shows examples"
+{
+    testLocData(locType, locData);
+    testFraction(frac, locType, locData);
+    if (frac[4] != 0) { // frac already has a right representation
+        return (frac);
+    } else { // frac has no right representation, but a left one
+        vector oreParas = rightOre(frac[1], frac[2], locType, locData)[1];
+        vector result = [frac[1], frac[2], oreParas[2], oreParas[1]];
+        return (result);
+    }
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    // monoidal localization
+    poly g = x;
+    poly f = Dx;
+    vector fracm = [g,f,0,0];
+    list L = g;
+    vector rm = convertLeftToRightFraction(fracm, 0, L);
+    print(rm);
+    f*rm[4]-g*rm[3];
+    // geometric localization
+    g = x+y;
+    f = Dx+Dy;
+    vector fracg = [g,f,0,0];
+    ideal p = x-1, y-3;
+    vector rg = convertLeftToRightFraction(fracg, 1, p);
+    print(rg);
+    f*rg[4]-g*rg[3];
+    // rational localization
+    intvec rat = 1;
+    f = Dx+Dy;
+    g = x;
+    vector fracr = [g,f,0,0];
+    vector rr = convertLeftToRightFraction(fracr, 2, rat);
+    print(rr);
+    f*rr[4]-g*rr[3];
+}
+//////////////////////////////////////////////////////////////////////
+proc isZeroFraction(vector frac)
+"USAGE:   isZeroFraction(frac), vector frac
+PURPOSE: determine if the vector frac represents zero
+ASSUME:  frac is a valid fraction
+RETURN:  int
+NOTE:    returns 1, if frac == 0; 0 otherwise
+EXAMPLE: example isZeroFraction; shows examples"
+{
+    if (frac[1] == 0 && frac[4] == 0) {
+        return(0);
+    }
+    if (frac[1] == 0) { // frac has no left representation
+        if (frac[3] == 0) { // the right representation of frac is zero
+            return(1);
+        }
+    } else { // frac has a left representation
+        if (frac[2] == 0) { // the left representation of frac is zero
+            return(1);
+        }
+    }
+    return(0);
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    isZeroFraction([42,0,0,0]);
+    isZeroFraction([0,0,Dx,3]);
+    isZeroFraction([1,1,1,1]);
+}
+//////////////////////////////////////////////////////////////////////
+proc isOneFraction(vector frac)
+"USAGE:   isOneFraction(frac), vector frac
+PURPOSE: determine if the vector frac represents one
+ASSUME:  frac is a valid fraction
+RETURN:  int
+NOTE:    1, if frac == 1; 0 otherwise
+EXAMPLE: example isOneFraction; shows examples"
+{
+    if (frac[1] == 0 && frac[4] == 0) {
+        return(0);
+    }
+    if (frac[1] == 0) { // frac has no left representation
+        if (frac[3] == frac[4]) { // the right representation of frac is zero
+            return(1);
+        }
+    } else { // frac has a left representation
+        if (frac[2] == frac[1]) { // the left representation of frac is zero
+            return(1);
+        }
+    }
+    return(0);
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    isOneFraction([42,42,0,0]);
+    isOneFraction([0,0,Dx,3]);
+    isOneFraction([1,0,0,1]);
+}
+////////// arithmetic ////////////////////////////////////////////////
+proc addLeftFractions(vector a, vector b, int locType, def locData, list #)
+"USAGE:   addLeftFractions(a, b, locType, locData(, override)),
+         vector a, b, int locType, list/vector/intvec locData(, int override)
+PURPOSE: add two left fractions in the specified localization
+ASSUME:
+RETURN:  vector
+NOTE:    the returned vector is the sum of a and b as fractions in the
+         localization specified by locType and locData.
+EXAMPLE: example addLeftFractions; shows examples"
+{
+    int override = 0;
+    if (size(#) > 1) {
+        if(typeof(#[1]) == "int") {
+            override = #[1];
+        }
+    }
+    if (!override) {
+        testLocData(locType, locData);
+        testFraction(a, locType, locData);
+        testFraction(b, locType, locData);
+    }
+    // check for a shortcut
+    if (isZeroFraction(a)) {
+        return(b);
+    }
+    if (isZeroFraction(b)) {
+        return(a);
+    }
+    if (a[1] == 0) { // a has no left representation
+        a = convertRightToLeftFraction(a, locType, locData);
+    }
+    if (b[1] == 0) { // b has no left representation
+        b = convertRightToLeftFraction(b, locType, locData);
+    }
+    if (a[1] == b[1]) { // a and b have the same left denominator
+        return ([a[1],a[2] + b[2],0,0]);
+    }
+    // no shortcut found, use regular method
+    vector oreParas = leftOre(b[1], a[1], locType, locData)[1];
+    vector result = [oreParas[1]*a[1],oreParas[1]*a[2]+oreParas[2]*b[2],0,0];
+    return (result);
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    // monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y+y;
+    list L = g1,g2;
+    poly s1 = g1;
+    poly s2 = g2;
+    poly r1 = Dx;
+    poly r2 = Dy;
+    vector frac1 = [s1,r1,0,0];
+    vector frac2 = [s2,r2,0,0];
+    vector rm = addLeftFractions(frac1, frac2, 0, L);
+    print(rm);
+    // geometric localization
+    ideal p = x-1, y-3;
+    vector rg = addLeftFractions(frac1, frac2, 1, p);
+    print(rg);
+    // rational localization
+    intvec v = 2;
+    s1 = y^2+y+1;
+    s2 = y-2;
+    r1 = Dx;
+    r2 = Dy;
+    frac1 = [s1,r1,0,0];
+    frac2 = [s2,r2,0,0];
+    vector rr = addLeftFractions(frac1, frac2, 2, v);
+    print(rr);
+}
+//////////////////////////////////////////////////////////////////////
+proc multiplyLeftFractions(vector a, vector b, int locType, def locData, list #)
+"USAGE:   multiplyLeftFractions(a, b, locType, locData(, override)),
+         vector a, b, int locType, list/vector/intvec locData, int override
+PURPOSE: multiply two left fractions in the specified localization
+ASSUME:
+RETURN:  vector
+NOTE:    the returned vector is the product of a and b as fractions in the
+         localization specified by locType and locData.
+EXAMPLE: example multiplyLeftFractions; shows examples"
+{
+    int override = 0;
+    if (size(#) > 1) {
+        if(typeof(#[1]) == "int") {
+            override = #[1];
+        }
+    }
+    if (!override) {
+        testLocData(locType, locData);
+        testFraction(a, locType, locData);
+        testFraction(b, locType, locData);
+    }
+    // check for a shortcut
+    if (isZeroFraction(a) || isZeroFraction(b)) {
+        return([1,0,0,1]);
+    }
+    if (isOneFraction(a)) {
+        return(b);
+    }
+    if (isOneFraction(b)) {
+        return(a);
+    }
+    if(a[1] == 0) {
+        a = convertRightToLeftFraction(a, locType, locData);
+    }
+    if(b[1] == 0) {
+        b = convertRightToLeftFraction(b, locType, locData);
+    }
+    if( (a[2] == 0) || (b[2] == 0) ) {
+        return ([1,0,0,1]);
+    }
+    if (a[2]*b[1] == b[1]*a[2]) { // trivial solution of the Ore condition
+        return([b[1]*a[1],a[2]*b[2]]);
+    }
+    // no shortcut found, use regular method
+    vector oreParas = ore(b[1], a[2], locType, locData, 0)[1];
+    vector result = [oreParas[1]*a[1],oreParas[2]*b[2],0,0];
+    return (result);
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    // monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y+y;
+    list L = g1,g2;
+    poly s1 = g1;
+    poly s2 = g2;
+    poly r1 = Dx;
+    poly r2 = Dy;
+    vector frac1 = [s1,r1,0,0];
+    vector frac2 = [s2,r2,0,0];
+    vector rm = multiplyLeftFractions(frac1, frac2, 0, L);
+    print(rm);
+    // geometric localization
+    ideal p = x-1, y-3;
+    vector rg = multiplyLeftFractions(frac1, frac2, 1, p);
+    print(rg);
+    // rational localization
+    intvec v = 2;
+    s1 = y^2+y+1;
+    s2 = y-2;
+    r1 = Dx;
+    r2 = Dy;
+    frac1 = [s1,r1,0,0];
+    frac2 = [s2,r2,0,0];
+    vector rr1 = multiplyLeftFractions(frac1, frac2, 2, v);
+    print(rr1);
+    vector rr2 = multiplyLeftFractions(frac2, frac1, 2, v);
+    print(rr2);
+    areEqualLeftFractions(rr1, rr2, 2, v);
+}
+//////////////////////////////////////////////////////////////////////
+proc areEqualLeftFractions(vector a, vector b, int locType, def locData)
+"USAGE:   areEqualLeftFractions(a, b, locType, locData), vector a, b,
+         int locType, list/vector/intvec locData
+PURPOSE: check if two given fractions are equal
+ASSUME:
+RETURN:  int
+NOTE:    returns 1 or 0, depending whether a=b as fractions in the
+         localization specified by locType and locData
+EXAMPLE: example areEqualLeftFractions; shows examples"
+{
+    testLocData(locType, locData);
+    testFraction(a, locType, locData);
+    testFraction(b, locType, locData);
+    if(a[1] == 0) {
+        a = convertRightToLeftFraction(a, locType, locData);
+    }
+    if(b[1] == 0) {
+        b = convertRightToLeftFraction(b, locType, locData);
+    }
+    vector negB = [b[1], -b[2], -b[3], b[4]];
+    //testFraction(negB, locType, locData); // unnecessary check
+    vector result = addLeftFractions(a, negB, locType, locData);
+    return(isZeroFraction(result));
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    // monoidal
+    poly g1 = x*y+3;
+    poly g2 = y^3;
+    list L = g1,g2;
+    poly s1 = g1;
+    poly s2 = s1*g2;
+    poly s3 = s2;
+    poly r1 = Dx;
+    poly r2 = g2*r1;
+    poly r3 = s1*r1+3;
+    vector fracm1 = [s1,r1,0,0];
+    vector fracm2 = [s2,r2,0,0];
+    vector fracm3 = [s3,r3,0,0];
+    areEqualLeftFractions(fracm1, fracm2, 0, L);
+    areEqualLeftFractions(fracm1, fracm3, 0, L);
+    areEqualLeftFractions(fracm2, fracm3, 0, L);
+}
+//////////////////////////////////////////////////////////////////////
+proc isInvertibleLeftFraction(vector frac, int locType, def locData)
+"USAGE:   isInvertibleLeftFraction(frac, locType, locData), vector frac,
+         int locType, list/vector/intvec locData
+PURPOSE: check if a fraction is invertible in the specified localization
+ASSUME:
+RETURN:  int
+NOTE:    - returns 1, if the numerator of frac is in the denominator set,
+         - returns 0, otherwise (NOTE: this does NOT mean that the fraction is
+           not invertible, it just means it could not be determined by the
+           method above).
+EXAMPLE: example isInvertibleLeftFraction; shows examples"
+{
+    testLocData(locType, locData);
+    testFraction(frac, locType, locData);
+    locData = normalizeLocalization(locType, locData);
+    if(frac[1] != 0) { // frac has a left representation
+        return(isInS(frac[2], locType, locData));
+    } else { // frac has no left, but a right representation
+        return(isInS(frac[3], locType, locData));
+    }
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    poly g1 = x+3;
+    poly g2 = x*y;
+    list L = g1,g2;
+    vector frac = [g1*g2, 17, 0, 0];
+    isInvertibleLeftFraction(frac, 0, L);
+    ideal p = x-1, y;
+    frac = [g1, x, 0, 0];
+    isInvertibleLeftFraction(frac, 1, p);
+    intvec rat = 1,2;
+    frac = [g1*g2, Dx, 0, 0];
+    isInvertibleLeftFraction(frac, 2, rat);
+}
+//////////////////////////////////////////////////////////////////////
+proc invertLeftFraction(vector frac, int locType, def locData)
+"USAGE:   invertLeftFraction(frac, locType, locData), vector frac, int locType,
+         list/vector/intvec locData
+PURPOSE: invert a fraction in the specified localization
+ASSUME:  frac is invertible in the loc. specified by locType and locData
+RETURN:  vector
+NOTE:    - returns the multiplicative inverse of frac in the localization
+           specified by locType and locData,
+         - throws error if frac is not invertible (NOTE: this does NOT mean
+           that the fraction is not invertible, it just means it could not be
+           determined by the method listed above).
+EXAMPLE: example invertLeftFraction; shows examples"
+{
+    // standard tests are done by isInvertibleLeftFraction
+    if (isInvertibleLeftFraction(frac, locType, locData)) {
+        return([frac[2],frac[1],frac[4],frac[3]]);
+    } else {
+        return([1,0,0,1]);
+    }
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S; S;
+    poly g1 = x+3;
+    poly g2 = x*y;
+    list L = g1,g2;
+    vector frac = [g1*g2, 17, 0, 0];
+    print(invertLeftFraction(frac, 0, L));
+    ideal p = x-1, y;
+    frac = [g1, x, 0, 0];
+    print(invertLeftFraction(frac, 1, p));
+    intvec rat = 1,2;
+    frac = [g1*g2, y, 0, 0];
+    print(invertLeftFraction(frac, 2, rat));
+}
+//////////////////////////////////////////////////////////////////////
+proc normalizeMonoidal(list L)
+"USAGE:    normalizeMonoidal(L), list L
+PURPOSE:  compute a normal form of monoidal localization data
+RETURN:   list
+NOTE:     given a list of polys, returns a list of all unique factors appearing
+          in the given polys
+EXAMPLE:  example normalizeMonoidal; shows examples"
+{
+    ideal allFactors;
+    int i;
+    for (i = 1; i <= size(L); i++) {
+        allFactors = allFactors, commutativeFactorization(L[i],1);
+    }
+    allFactors = simplify(allFactors,1+2+4);
+    // simplify: divide by leading coefficients (1),
+    //   purge zero generators (2), purge double entries (4)
+    ideal rev = sort(allFactors)[1]; // sort sorts ascendingly
+    return(list(rev[size(rev)..1])); // reverse order and cast to list
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl(); setring S;
+    list L = x^2*y^3, (x+1)*(x*y-3*y^2+1);
+    L = normalizeMonoidal(L);
+    print(L);
+}
+//////////////////////////////////////////////////////////////////////
+proc normalizeRational(intvec v)
+"USAGE:    normalizeRational(v), intvec v
+PURPOSE:  compute a normal form of rational localization data
+RETURN:   intvec
+NOTE:     purges double entries and sorts ascendingly
+EXAMPLE:  example normalizeRational; shows examples"
+{
+    int n = nvars(basering);
+    int i;
+    intvec result;
+    intvec occurring = 0:n;
+    for (i = 1; i <= size(v); i++) {
+        occurring[v[i]] = 1;
+    }
+    for (i = 1; i <= size(occurring); i++) {
+        if (occurring[i]) {
+            if (result == 0) {
+                result = i;
+            } else {
+                result = result, i;
+            }
+        }
+    }
+    return(result);
+}
+example
+{
+    "EXAMPLE:"; echo = 2;
+    ring R; setring R;
+    intvec v = 9,5,9,3,1,5;
+    v = normalizeRational(v);
+    v;
+}
+////////// internal functions ////////////////////////////////////////
+static proc inducesCommutativeSubring(def input)
+{
+    ideal vars = variables(input);
+    int i, j;
+    for (i = 1; i <= size(vars); i++) {
+        for (j = i + 1; j <= size(vars); j++) {
+            if (vars[i]*vars[j] != vars[j]*vars[i]) {
+                return(0);
+            }
+        }
+    }
+    return(1);
+}
+//////////////////////////////////////////////////////////////////////
+static proc commutativeFactorization(poly p, list #)
+"USAGE:    commutativeFactorization(p[, #]), poly p[, list #]
+PURPOSE:  compute a factorization of p ignoring non-commutative relations
+RETURN:   list or ideal
+NOTE:     the optional parameter is passed to factorize after changing to a
+          commutative ring, the result of factorize is transfered back to
+          basering
+SEE ALSO: factorize
+EXAMPLE:  "
+{
+    int factorType = 0;
+    if (size(#)  > 0) {
+        if (typeof(#[1]) == "int") {
+            factorType = #[1];
+        }
+    }
+    list RL = ringlist(basering);
+    if (size(RL) > 4) {
+        def bsRing = basering;
+        RL = RL[1..4];
+        def commRing = ring(RL);
+        setring commRing;
+        poly commP = imap(bsRing, p);
+        if (factorType == 1) {
+            ideal commFactors = factorize(commP, 1);
+            setring bsRing;
+            return(imap(commRing, commFactors));
+        } else {
+            list commFac = factorize(commP, factorType);
+            intvec exponents = commFac[2];
+            ideal commFactors = commFac[1];
+            setring bsRing;
+            list result;
+            result[1] = imap(commRing, commFactors);
+            result[2] = exponents;
+            return(result);
+        }
+    } else {
+        return(factorize(p, factorType));
+    }
+}
+//////////////////////////////////////////////////////////////////////
+static proc normalizeLocalization(int locType, def locData) {
+    if (locType == 0) {
+        return(normalizeMonoidal(locData));
+    }
+    if (locType == 1) {
+        return(std(locData));
+    }
+    if (locType == 2) {
+        return(normalizeRational(locData));
+    }
+    return(locData);
+}
+//////////////////////////////////////////////////////////////////////
+static proc intvecComplement(intvec v, intvec w) { // complement of v in w as sets
+    intvec result;
+    int i;
+    int j;
+    int foundMatch;
+    for (i = 1; i <= size(w); i++) {
+        foundMatch = 0;
+        for (j = 1; j <= size(v); j++) {
+            if (w[i] == v[j]) {
+                foundMatch = 1;
+                break;
+            }
+        }
+        if (!foundMatch) { // v[i] is not in w
+            if (result == 0) {
+                result = w[i];
+            } else {
+                result = result, w[i];
+            }
+        }
+    }
+    return(result);
+}
+//////////////////////////////////////////////////////////////////////
+////////// internal testing procedures ///////////////////////////////
+static proc testIsInS()
+{
+    print("  testing isInS...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    // monoidal localization
+    poly g1 = x^2*y+x+2;
+    poly g2 = y^3+x*y;
+    list L = g1, g2;
+    poly g = g1^2*g2;
+    if (!isInS(g, 0, L)) {
+        ERROR("Weyl monoidal isInS direct positive failed");
+    }
+    if (!isInS(y^2+x, 0, L)) {
+        ERROR("Weyl monoidal isInS indirect positive failed");
+    }
+    if (isInS(g-1, 0, L)) {
+        ERROR("Weyl monoidal isInS negative failed");
+    }
+    // geometric localization
+    ideal p = x-1, y-3;
+    g = x^2+y-3;
+    if (!isInS(g, 1, p)) {
+        ERROR("Weyl geometric isInS positive failed");
+    }
+    if (isInS((x-1)*g, 1, p)) {
+        ERROR("Weyl geometric isInS negative failed");
+    }
+    // rational localization
+    intvec v = 2;
+    if (!isInS(y^5+17*y^2-4, 2, v)) {
+        ERROR("Weyl rational isInS positive failed");
+    }
+    if (isInS(x*y, 2, v)) {
+        ERROR("Weyl rational isInS negative failed");
+    }
+    print("    isInS OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testLeftOre() {
+    print("  testing leftOre...");
+    // Weyl
+    ring W = 0,(x,y,Dx,Dy),dp;
+    def ncW = Weyl();
+    setring ncW;
+    //// monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y;
+    list L = g1,g2;
+    poly g = g1^2*g2;
+    poly f = Dx;
+    vector rm = leftOre(g, f, 0, L)[1];
+    if (rm[1] == 0 || rm[2]*g-rm[1]*f != 0) {
+        ERROR("Weyl monoidal left Ore failed");
+    }
+    //// geometric localization
+    vector p1,p2,p3,p4 = [1,3],[0,0],[-1,-2],[10,-180];
+    vector rg;
+    ideal p;
+    list pVecs = p1,p2,p3,p4;
+    f = Dx;
+    g = x^2+y+3;
+    for(int i = 1; i <= 4; i = i + 1) {
+        p = x-pVecs[i][1], y-pVecs[i][2];
+        rg = leftOre(g, f, 1, p)[1];
+        if (rg[1] == 0 || rg[2]*g-rg[1]*f != 0) {
+            ERROR("Weyl geometric left Ore failed at maximal ideal"
+              + " induced by " + string(pVecs[i]));
+        }
+    }
+    //// rational localization
+    intvec rat = 1;
+    f = Dx+Dy;
+    g = x;
+    vector rr = leftOre(g, f, 2, rat)[1];
+    if (rr[1] == 0 || rr[2]*g-rr[1]*f != 0) {
+        ERROR("Weyl rational left Ore failed");
+    }
+    // shift rational localization
+    ring S = 0,(x,y,Sx,Sy),dp;
+    matrix D[4][4];
+    D[1,3] = Sx;
+    D[2,4] = Sy;
+    def ncS = nc_algebra(1, D);
+    setring ncS;
+    rat = 1;
+    poly f = Sx+Sy;
+    poly g = x;
+    vector rr = leftOre(g, f, 2, rat)[1];
+    if (rr[1] == 0 || rr[2]*g-rr[1]*f != 0) {
+        ERROR("shift rational left Ore failed");
+    }
+    // q-shift rational localization
+    ring Q = (0,q),(x,y,Qx,Qy),dp;
+    matrix C[4][4] = UpOneMatrix(4);
+    C[1,3] = q;
+    C[2,4] = q;
+    def ncQ = nc_algebra(C, 0);
+    setring ncQ;
+    rat = 1;
+    poly f = Qx+Qy;
+    poly g = x;
+    vector rr = leftOre(g, f, 2, rat)[1];
+    if (rr[1] == 0 || rr[2]*g-rr[1]*f != 0) {
+        ERROR("q-shift rational left Ore failed");
+    }
+    print("    leftOre OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testRightOre() {
+    print("  testing rightOre...");
+    // Weyl
+    ring W = 0,(x,y,Dx,Dy),dp;
+    def ncW = Weyl();
+    setring ncW;
+    //// monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y;
+    list L = g1,g2;
+    poly g = x;
+    poly f = Dx;
+    vector rm = rightOre(g, f, 0, L)[1];
+    if (rm[1] == 0 || f*rm[1]-g*rm[2] != 0) {
+        ERROR("Weyl monoidal right Ore failed");
+    }
+    //// geometric localization
+    g = x+y;
+    f = Dx+Dy;
+    ideal p = x-1,y-3;
+    vector rg = rightOre(g, f, 1, p)[1];
+    if (rg[1] == 0 || f*rg[1]-g*rg[2] != 0) {
+        ERROR("Weyl geometric right Ore failed");
+    }
+    //// rational localization
+    intvec rat = 1;
+    f = Dx+Dy;
+    g = x;
+    vector rr = rightOre(g, f, 2, rat)[1];
+    if (rr[1] == 0 || f*rr[1]-g*rr[2] != 0) {
+        ERROR("Weyl rational right Ore failed");
+    }
+    // shift rational localization
+    ring S = 0,(x,y,Sx,Sy),dp;
+    matrix D[4][4];
+    D[1,3] = Sx;
+    D[2,4] = Sy;
+    def ncS = nc_algebra(1, D);
+    setring ncS;
+    rat = 1;
+    poly f = Sx+Sy;
+    poly g = x;
+    vector rr = rightOre(g, f, 2, rat)[1];
+    if (rr[1] == 0 || f*rr[1]-g*rr[2] != 0) {
+        ERROR("shift rational right Ore failed");
+    }
+    ring Q = (0,q),(x,y,Qx,Qy),dp;
+    matrix C[4][4] = UpOneMatrix(4);
+    C[1,3] = q;
+    C[2,4] = q;
+    def ncQ = nc_algebra(C, 0);
+    setring ncQ;
+    rat = 1;
+    poly f = Qx+Qy;
+    poly g = x;
+    vector rr = rightOre(g, f, 2, rat)[1];
+    if (rr[1] == 0 || f*rr[1]-g*rr[2] != 0) {
+        ERROR("q-shift rational right Ore failed");
+    }
+    print("    rightOre OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testAddLeftFractions()
+{
+    print("  testing addLeftFractions...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    // monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y+y;
+    list L = g1,g2;
+    vector frac1 = [g1,Dx,0,0];
+    vector frac2 = [g2,Dy,0,0];
+    vector rm = addLeftFractions(frac1, frac2, 0, L);
+    if (rm[1] != x^2*y+4*x*y+3*y || rm[2] != x*y*Dx+y*Dx+x*Dy+3*Dy) {
+        ERROR("Weyl monoidal addition failed");
+    }
+    // geometric localization
+    ideal p = x-1,y-3;
+    vector rg = addLeftFractions(frac1, frac2, 1, p);
+    if (rg[1] != x^2*y+4*x*y+3*y || rg[2] != x*y*Dx+y*Dx+x*Dy+3*Dy) {
+        ERROR("Weyl geometric addition failed");
+    }
+    // rational localization
+    intvec v = 2;
+    frac1 = [y^2+y+1,Dx,0,0];
+    frac2 = [y-2,Dy,0,0];
+    vector rr = addLeftFractions(frac1, frac2, 2, v);
+    if (rr[1] != y^3-y^2-y-2 || rr[2] != y^2*Dy+y*Dx+y*Dy-2*Dx+Dy) {
+        ERROR("Weyl rational addition failed");
+    }
+    print("    addLeftFractions OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testMultiplyLeftFractions() {
+    print("  testing multiplyLeftFractions...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    // monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y+y;
+    list L = g1,g2;
+    vector frac1 = [g1,Dx,0,0];
+    vector frac2 = [g2,Dy,0,0];
+    vector rm = multiplyLeftFractions(frac1, frac2, 0, L);
+    if (rm[1] != g1*g2^2 || rm[2] != x*y*Dx*Dy+y*Dx*Dy-y*Dy) {
+        ERROR("Weyl monoidal multiplication error");
+    }
+    // geometric localization
+    ideal p = x-1,y-3;
+    vector rg = multiplyLeftFractions(frac1, frac2, 1, p);
+    if (rg[1] != g1*g2*(x+1) || rg[2] != x*Dx*Dy+Dx*Dy-Dy) {
+        ERROR("Weyl geometric multiplication error");
+    }
+    // rational localization
+    intvec v = 2;
+    frac1 = [y^2+y+1,Dx,0,0];
+    frac2 = [y-2,Dy,0,0];
+    vector rr = multiplyLeftFractions(frac1, frac2, 2, v);
+    if (rr[1] != (y^2+y+1)*(y-2) || rr[2] != Dx*Dy) {
+        ERROR("Weyl rational multiplication (1*2) error");
+    }
+    rr = multiplyLeftFractions(frac2, frac1, 2, v);
+    if (rr[1] != (y^2+y+1)^2*(y-2) || rr[2] != y^2*Dx*Dy+y*Dx*Dy-2*y*Dx+Dx*Dy-Dx) {
+        ERROR("Weyl rational multiplication (2*1) error");
+    }
+    print("    multiplyLeftFractions OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testAreEqualLeftFractions() {
+    print("  testing areEqualLeftFractions...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    // monoidal
+    poly g1 = x*y+3;
+    poly g2 = y^3;
+    list L = g1,g2;
+    vector fracm1 = [g1,Dx,0,0];
+    vector fracm2 = [g1*g2,g2*Dx,0,0];
+    vector fracm3 = [g1*g2,g1*Dx+3,0,0];
+    if (!areEqualLeftFractions(fracm1, fracm2, 0, L)) {
+        ERROR("Weyl monoidal positive basic comparison error");
+    }
+    if (areEqualLeftFractions(fracm1, fracm3, 0, L)) {
+        ERROR("Weyl monoidal first negative basic comparison error");
+    }
+    if (areEqualLeftFractions(fracm2, fracm3, 0, L)) {
+        ERROR("Weyl monoidal second negative basic comparison error");
+    }
+    // geometric
+    ideal p = x+5, y-2;
+    vector fracg1 = [g1,Dx,0,0];
+    vector fracg2 = [g1*g2,g2*Dx,0,0];
+    vector fracg3 = [g1*g2,g1*Dx+3,0,0];
+    if (!areEqualLeftFractions(fracg1, fracg2, 1, p)) {
+        ERROR("Weyl geometric positive basic comparison error");
+    }
+    if (areEqualLeftFractions(fracg1, fracg3, 1, p)) {
+        ERROR("Weyl geometric first negative basic comparison error");
+    }
+    if (areEqualLeftFractions(fracg2, fracg3, 1, p)) {
+        ERROR("Weyl geometric second negative basic comparison error");
+    }
+    // rational
+    intvec rat = 1,4;
+    vector fracr1 = [x+Dy,Dx,0,0];
+    vector fracr2 = [x*Dy*(x+Dy),x*Dx*Dy,0,0];
+    vector fracr3 = [Dy*x*(x+Dy),x*Dx*Dy+1,0,0];
+    if (!areEqualLeftFractions(fracr1, fracr2, 2, rat)) {
+        ERROR("Weyl rational positive basic comparison error");
+    }
+    if (areEqualLeftFractions(fracr1, fracr3, 2, rat)) {
+        ERROR("Weyl rational first negative basic comparison error");
+    }
+    if (areEqualLeftFractions(fracr2, fracr3, 2, rat)) {
+        ERROR("Weyl rational second negative basic comparison error");
+    }
+    print("    areEqualLeftFractions OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testConvertLeftToRightFraction() {
+    print("  testing convertLeftToRightFraction...");
+    // Weyl
+    ring W = 0,(x,y,Dx,Dy),dp;
+    def ncW = Weyl();
+    setring ncW;
+    //// monoidal localization
+    vector fracm = [x,Dx,0,0];
+    list L = x;
+    vector rm = convertLeftToRightFraction(fracm, 0, L);
+    if (!fracStatus(rm, 0, L)[1]) {
+        ERROR("Weyl monoidal convertLeftToRightFraction failed");
+    }
+    //// geometric localization
+    vector fracg = [x+y,Dx+Dy,0,0];
+    ideal p = x-1,y-3;
+    vector rg = convertLeftToRightFraction(fracg, 1, p);
+    if (!fracStatus(rg, 1, p)[1]) {
+        ERROR("Weyl geometric convertLeftToRightFraction failed");
+    }
+    //// rational localization
+    intvec rat = 1;
+    vector fracr = [x,Dx+Dy,0,0];
+    vector rr = convertLeftToRightFraction(fracr, 2, rat);
+    if (!fracStatus(rr, 2, rat)[1]) {
+        ERROR("Weyl rational convertLeftToRightFraction failed");
+    }
+    // shift rational localization
+    ring S = 0,(x,y,Sx,Sy),dp;
+    matrix D[4][4];
+    D[1,3] = Sx;
+    D[2,4] = Sy;
+    def ncS = nc_algebra(1, D);
+    setring ncS;
+    vector fracr = [x,Sx+Sy,0,0];
+    vector rr = convertLeftToRightFraction(fracr, 2, rat);
+    if (!fracStatus(rr, 2, rat)[1]) {
+        ERROR("Shift rational convertLeftToRightFraction failed");
+    }
+    // q-shift rational localization
+    ring Q = (0,q),(x,y,Qx,Qy),dp;
+    matrix C[4][4] = UpOneMatrix(4);
+    C[1,3] = q;
+    C[2,4] = q;
+    def ncQ = nc_algebra(C, 0);
+    setring ncQ;
+    vector fracr = [x,Qx+Qy,0,0];
+    vector rr = convertLeftToRightFraction(fracr, 2, rat);
+    if (!fracStatus(rr, 2, rat)[1]) {
+        ERROR("q-shift rational convertLeftToRightFraction failed");
+    }
+    print("    convertLeftToRightFraction OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testConvertRightToLeftFraction() {
+    print("  testing convertRightToLeftFraction...");
+    // Weyl
+    ring W = 0,(x,y,Dx,Dy),dp;
+    def ncW = Weyl();
+    setring ncW;
+    //// monoidal localization
+    poly g1 = x+3;
+    poly g2 = x*y;
+    list L = g1,g2;
+    vector fracm = [0,0,Dx,g1^2*g2];
+    vector rm = convertRightToLeftFraction(fracm, 0, L);
+    if (!fracStatus(rm, 0, L)[1]) {
+        ERROR("Weyl monoidal convertRightToLeftFraction failed");
+    }
+    //// geometric localization
+    ideal p = x-1,y-3;
+    vector fracg = [0,0,Dx,x^2+y];
+    vector rg = convertRightToLeftFraction(fracg, 1, p);
+    if (!fracStatus(rg, 1, p)[1]) {
+        ERROR("Weyl geometric convertRightToLeftFraction failed");
+    }
+    //// rational localization
+    intvec rat = 1;
+    vector fracr = [0,0,Dx+Dy,x];
+    vector rr = convertRightToLeftFraction(fracr, 2, rat);
+    if (!fracStatus(rr, 2, rat)[1]) {
+        ERROR("Weyl rational convertRightToLeftFraction failed");
+    }
+    // shift rational localization
+    ring S = 0,(x,y,Sx,Sy),dp;
+    matrix D[4][4];
+    D[1,3] = Sx;
+    D[2,4] = Sy;
+    def ncS = nc_algebra(1, D);
+    setring ncS;
+    vector fracr = [0,0,Sx+Sy,x];
+    vector rr = convertRightToLeftFraction(fracr, 2, rat);
+    if (!fracStatus(rr, 2, rat)[1]) {
+        ERROR("Shift rational convertRightToLeftFraction failed");
+    }
+    // q-shift rational localization
+    ring Q = (0,q),(x,y,Qx,Qy),dp;
+    matrix C[4][4] = UpOneMatrix(4);
+    C[1,3] = q;
+    C[2,4] = q;
+    def ncQ = nc_algebra(C, 0);
+    setring ncQ;
+    vector fracr = [0,0,Qx+Qy,x];
+    vector rr = convertRightToLeftFraction(fracr, 2, rat);
+    if (!fracStatus(rr, 2, rat)[1]) {
+        ERROR("q-shift rational convertRightToLeftFraction failed");
+    }
+    print("    convertRightToLeftFraction OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testIsInvertibleLeftFraction() {
+    print("  testing isInvertibleLeftFraction...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    poly g1 = x+3;
+    poly g2 = x*y;
+    // monoidal
+    list L = g1, g2;
+    if (!isInvertibleLeftFraction([g1*g2,17,0,0], 0, L)) {
+        ERROR("Weyl monoidal positive invertibility test error");
+    }
+    if (isInvertibleLeftFraction([g1*g2,Dx,0,0], 0, L)) {
+        ERROR("Weyl monoidal negative invertibility test error");
+    }
+    if (!isInvertibleLeftFraction([1,1,1,1], 0, L)) {
+        ERROR("Weyl monoidal one invertibility test error");
+    }
+    if (isInvertibleLeftFraction([1,0,0,1], 0, L)) {
+        ERROR("Weyl monoidal zero invertibility test error");
+    }
+    // geometric
+    ideal p = x-1, y;
+    if (!isInvertibleLeftFraction([g1,3*x,0,0], 1, p)) {
+        ERROR("Weyl geometric positive invertibility test error");
+    }
+    if (isInvertibleLeftFraction([g1,Dx,0,0], 0, L)) {
+        ERROR("Weyl geometric negative invertibility test error");
+    }
+    if (!isInvertibleLeftFraction([1,1,1,1], 1, p)) {
+        ERROR("Weyl geometric one invertibility test error");
+    }
+    if (isInvertibleLeftFraction([1,0,0,1], 1, p)) {
+        ERROR("Weyl geometric zero invertibility test error");
+    }
+    // rational
+    intvec rat = 1,2;
+    if (!isInvertibleLeftFraction([g1*g2,y,0,0], 2, rat)) {
+        ERROR("Weyl rational positive invertibility test error");
+    }
+    if (isInvertibleLeftFraction([g1*g2,Dx,0,0], 0, L)) {
+        ERROR("Weyl rational negative invertibility test error");
+    }
+    if (!isInvertibleLeftFraction([1,1,1,1], 2, rat)) {
+        ERROR("Weyl rational one invertibility test error");
+    }
+    if (isInvertibleLeftFraction([1,0,0,1], 2, rat)) {
+        ERROR("Weyl rational zero invertibility test error");
+    }
+    print("    isInvertibleLeftFraction OK");
+}
+//////////////////////////////////////////////////////////////////////
+static proc testInvertLeftFraction() {
+    print("  testing invertLeftFraction...");
+    ring R = 0,(x,y,Dx,Dy),dp;
+    def S = Weyl();
+    setring S;
+    poly g1 = x+3;
+    poly g2 = x*y;
+    // monoidal
+    list L = g1, g2;
+    vector rm = [g1*g2, 17, 0, 0];
+    vector rmInv = invertLeftFraction(rm, 0 , L);
+    if (!isOneFraction(multiplyLeftFractions(rm, rmInv, 0, L))) {
+        ERROR("Weyl monoidal inversion error");
+    }
+    // geometric
+    ideal p = x-1, y;
+    vector rg = [g1, 3*x, 0, 0];
+    vector rgInv = invertLeftFraction(rg, 1, p);
+    if (!isOneFraction(multiplyLeftFractions(rg, rgInv, 1, p))) {
+        ERROR("Weyl geometric inversion error");
+    }
+    // rational
+    intvec rat = 1,2;
+    vector rr = [g1*g2, y, 0, 0];
+    vector rrInv = invertLeftFraction(rr, 2, rat);
+    if (!isOneFraction(multiplyLeftFractions(rr, rrInv, 2, rat))) {
+        ERROR("Weyl rational inversion error");
+    }
+    print("    invertLeftFraction OK");
+}
+//////////////////////////////////////////////////////////////////////
