////////////////////////////////////////////////////////////////
version="version olga.lib 4.1.2.0 Feb_2019 "; // $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
                if (defined(basering)) {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 (defined(basering))
                      {
                        if (size(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 (!defined(basering))
                {
                    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 (!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);
        }
        int n = nvars(basering);
        if (size(locData) < n) { // there are variables to eliminate
            ideal I = p;
            intvec modLocData = intvecComplement(locData, 1..nvars(basering));
            I = eliminateNC(I, modLocData);
            if (size(I)) {
                return(1);
            }
        } else {
            return(p != 0);
        }
    }
    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 transferred 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");
    }
    intvec w = 4,2,3,4,1;
    if (!isInS(x*y*Dx*Dy,2,w)) {
        ERROR("Weyl total rational isInS positive failed");
    }
    if (isInS(0, 2, w)) {
        ERROR("Weyl total 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");
}
//////////////////////////////////////////////////////////////////////