//////////////////////////////////////////////////////////////// version="version olga.lib 1.2.1.0 Jul_2017"; category="Noncommutative"; info=" LIBRARY: olga.lib Ore-localization in G-Algebras AUTHOR: Johannes Hoffmann, email: johannes.hoffmann at math.rwth-aachen.de Localization types: Type 0: Monoidal: represented by a list of polynomials g1...gk implemented and tested only for Weyl algebras of type K S = {g1^n1*...*gk^nk | n1...nk in N_0} Type 1: Geometrical: represented by a prime ideal p (Weyl-algebra only) implemented and tested only for Weyl algebras of type K S = K[x1,...,xN]\p Type 2: Rational: represented by an intvec v = [i_1,...,i_k] in the range 1..nvars(basering) some procedures require that i_1,...,i_k = 1..k (rationalForm,fractionLeadcoef,cancelLeftFraction,reduceLeftFraction) S = K[x_i_1,...,x_i_k]\{0} The type is represented by an int, whereas the localization data is a def to accommodate the different cases A fraction is represented as a vector with four entries: [s,r,p,t] Here, s^{-1}r is the left fraction representation while pt^{-1} is the right fraction representation. If s or t is zero, it means that the corresponding representation has not been evaluated. If both are zero, the fraction is not valid. PROCEDURES: checkFraction(vector, int, def); check if the given vector is a representation of a fraction in the specified localization checkLocData(int, def); check if the given data specifies a denominator set ore(poly, poly, int, def, int); compute Ore data for a given one-sided fraction 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 cancelLeftFraction(vector, int, def); cancels a given fraction in the specified localization reduceLeftFraction(vector, vector, int, def); completes one step of reducing a fraction in the specified localization isInS(poly, int, def); determine if a polynomial is in a denominator set isInvertibleLeftFraction(vector, int, def); check if a fraction is invertible in the specified localization (NOTE: check description) invertLeftFraction(vector, int, def); invert a fraction in the specified localization (NOTE: check description) fractionLeadcoef(vector, int, def); determine the leading coefficient of a fraction in the specified rational localization fractionLeadmonom(vector, int, def); determine the leading monomial of a fraction in the specified rational localization fractionLeadexp(vector, int, def); determine the leading exponent of a fraction in the specified rational localization rationalForm(def, int, def); factor out the monomials in non-invertable variables in the specified rational or geometrical localization "; LIB "dmodloc.lib"; LIB "ncpreim.lib"; LIB "elim.lib"; LIB "ncalg.lib"; //////////////////////////////////////////////////////////////////// proc testOlga() { example isInS; example ore; example convertRightToLeftFraction; example convertLeftToRightFraction; example addLeftFractions; example multiplyLeftFractions; example areEqualLeftFractions; example rationalForm; example fractionLeadcoef; example fractionLeadmonom; example fractionLeadexp; example reduceLeftFraction; example cancelLeftFraction; } //////////////////////////////////////////////////////////////////// proc checkFraction(vector a, int locType, def locData) "USAGE: checkFraction(a, locType, locData), vector a, int locType, list/intvec/vector locData PURPOSE: check if the given vector is a representation of a fraction in the specified localization ASSUME: RETURN: nothing NOTE: - throws error if checks were not successful, - checks if at least one denominator is specified, - checks if the specified denominators are in the denominator set specified by locType and locData, - if both denominators are specified, checks if left fraction and right fraction are equal. EXAMPLE: " { checkLocData(locType, locData); if( (a[1] == 0) && (a[4] == 0) ) { ERROR("vector is not a valid fraction " + print(a)); } if( (a[1] != 0) && (!isInS(a[1], locType, locData)) ) { ERROR("the left denominator " + print(a[1]) + " of fraction " + print(a) + " is not in the denominator set"); } if( (a[4] != 0) && (!isInS(a[4], locType, locData)) ) { ERROR("the right denominator " + print(a[4]) + " of fraction " + print(a) + " is not in the denominator set"); } if( (a[1] != 0) && (a[4] != 0) ) { if(a[2]*a[4]-a[1]*a[3] != 0) { ERROR("left and right fraction don't match in " + print(a)); } } } //////////////////////////////////////////////////////////////////// proc checkLocData(int locType, def locData) "USAGE: checkLocData(locType, locData), int locType, list/vector/intvec locData PURPOSE: check if the given data specifies a denominator set ASSUME: RETURN: nothing NOTE: - throws error if checks were not successful, - locType can only be 0, 1, or 2. EXAMPLE: " { if( (locType < 0) || (locType > 2) ) { ERROR("localization type invalid! valid types are: (0) for a monoidal localization | (1) for a geometrical localization | (2) for a rational localization"); } int i, j; int n = nvars(basering); string t = typeof(locData); if(locType == 0) { if(t != "list") { ERROR("for a monoidal localization, locData has to be a list, but is of type " + t); } if(size(locData) < 1) { ERROR("for a monoidal localization, locData has to be a non-empty list"); } intvec occurringVars = 0:n; //n = nvars(basering) div 2; //intvec comVars = 1..n; for(i = 1; i <= size(locData); i++) { if(typeof(locData[i]) != "poly") { ERROR("for a monoidal localization, locData has to be a list of polys, but entry " + print(i) + " is " + string(locData[i]) + ", which is of type " + t); } if(leadcoef(locData[i]) != 1) { ERROR("for a monoidal localization, locData has to be a list of monic polys, but entry " + print(i) + " has leading coefficient " + leadcoef(locData[i])); } for (j = 1; j <= n; j++) { if (subst(locData[i], var(j), 0) != locData[i]) { occurringVars[j] = 1; } } } for (i = 1; i <= n; i++) { for (j = i + 1; j <= n; j++) { if (occurringVars[i] == 1 && occurringVars[j] == 1) { if (var(i) * var(j) - var(j) * var(i) != 0) { ERROR("for a monoidal localization, locData has to be a list of polys contained in a commutative polynomial subring of basering"); } } } } } if(locType == 1) { // TODO: check order of variables if(!isWeyl()) { //ERROR("for a geometric localization, basering has to be a Weyl algebra"); } n = nvars(basering) div 2; for(i = 1; i <= n; i = i + 1) { if(var(n+i)*var(i)-var(i)*var(n+i) != 1) { //ERROR("for a geometrical localization, the Weyl algebra has to have the form K"); } } if(t != "ideal") { ERROR("for a geometrical localization, locData has to be of type ideal, but is of type " + t); } } if(locType == 2) { if(t != "intvec") { ERROR("for a rational localization, locData has to be an intvec, but is of type " + t); } if(!admissibleSub(locData)) { ERROR("provided localization data doesn't match: the indexed variables don't generate a sub-g-algebra"); } } } //////////////////////////////////////////////////////////////////// static proc checkMonDivisibility(poly p, list L) "USAGE: PURPOSE: ASSUME: RETURN: NOTE: this procedure will be replaced in the next version EXAMPLE: " { int n = nvars(basering) div 2; int i, j, k, l, bool; def bsRing = basering; list RL = ringlist(basering); list fact, result; RL = delete(RL,6); RL = delete(RL,5); ideal J = L[1]; for(i = 2; i <= size(L); i = i + 1) { J = J,L[i]; } def R = ring(RL); setring R; poly f = imap(bsRing, p); ideal I = imap(bsRing, J); fact = factorize(I[1],2); ideal irredFactors = fact[1]; for(i = 2; i <= size(I); i = i + 1) { // for every g_i fact = factorize(I[i],2); // list[ideal,intvec] l = size(irredFactors); for(j = 1; j <= size(fact[1]); j = j + 1) { // for every factor g_ij of g_i bool = 0; for(k = 1; k <= l; k = k + 1) { if(irredFactors[k] == fact[1][j]) { bool = 1; } } if(bool == 0) { irredFactors = irredFactors, fact[1][j]; } } } intvec exponent = 0:size(irredFactors); list expList; for(i = size(I); i >= 1; i = i - 1) { // for every g_i exponent = 0:size(irredFactors); fact = factorize(I[i],2); for(j = 1; j <= size(fact[1]); j = j + 1) { // for every factor g_ij of g_i for(k = 1; k <= size(irredFactors); k = k + 1) { if(irredFactors[k] == fact[1][j]) { exponent[k] = fact[2][j]; break; } } } expList = insert(expList, exponent); } exponent = 0:size(irredFactors); fact = factorize(f, 2); for(i = 1; i <= size(fact[1]); i = i + 1) { // for every factor of f bool = 0; for(j = 1; j <= size(irredFactors); j = j + 1) { if(irredFactors[j] == fact[1][i]) { bool = 1; exponent[j] = fact[2][i]; } } if(bool == 0) { result = 0; return (result); } } setring bsRing; ideal I = imap(R, irredFactors); list factors; for(i = 1; i <= size(I); i = i + 1) { factors[i] = I[i]; } result = 1, factors, exponent, expList; return (result); } 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; checkMonDivisibility(g, L); } //////////////////////////////////////////////////////////////////// proc isInS(poly p, int locType, def locData) "USAGE: isInS(p, locType, locData), poly p, int locType, list/vector/intvec locData PURPOSE: determine if a polynomial is in a denominator set ASSUME: RETURN: int NOTE: - returns 1, if p is in the denominator set specified by locType and locData, - returns 0, otherwise. EXAMPLE: example isInS; shows examples" { checkLocData(locType,locData); if (p == 0) { return (0); // the zero polynomial is never in S } if(locType == 0) { int n = nvars(basering) div 2; intvec comVars = 1..n; if(polyVars(p,comVars)) { // p is a poly in the first n variables int result = checkMonDivisibility(p, locData)[1]; return (result); } } if(locType == 1) { int n = nvars(basering) div 2; intvec comVars = 1..n; if (polyVars(p,comVars)) { ideal I = std(locData); if (NF(p,I) != 0) { return(1); } } } if(locType == 2) { if(polyVars(p, locData)) { return (1); // p is in the sub-algebra generated by the indexed variables } } 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); isInS(f,0,L); // geometrical localization ideal p = x-1,y-3; g = x^2+y-3; f = (x-1)*g; isInS(g,1,p); isInS(f,1,p); // rational localization intvec v = 2; g = y^5+17*y^2-4; f = x*y; isInS(g,2,v); isInS(f,2,v); } //////////////////////////////////////////////////////////////////// proc ore(poly g, poly f, int locType, def locData, int rightOre) "USAGE: ore(g, f, locType, locData, rightOre), poly g, f, int locType, list/vector/intvec locData, int rightOre PURPOSE: compute Ore data for a given one-sided fraction ASSUME: g is in the denominator set, determined via locType and locData RETURN: vector NOTE: - vector [s,r] will be returned with s in the denominator set, - if rightOre=0, sf=gr holds, - if rightOre=1, fs=rg holds. EXAMPLE: example ore; shows examples" { checkLocData(locType, locData); if(!isInS(g, locType, locData)) { ERROR("cannot find Ore-parameter as poly " + print(g) + " is not in the denominator set"); } if(rightOre == 0) { return (oreHelp(g, f, locType, locData, 0)); } if(rightOre == 1) { def bsRing = basering; if(locType == 0) { vector modLocData; for(int i = 1; i <= size(locData); i = i + 1) { modLocData = modLocData + locData[i] * gen(i); } } //if(locType == 1) { // vector modLocData = locData; //} def oppRing = opposite(bsRing); setring oppRing; if(locType == 0) { vector oppModLocData = oppose(bsRing, modLocData); list oppLocData; for(int j = size(oppModLocData); i >= 1; i = i - 1) { oppLocData = insert(oppLocData, oppModLocData[j]); } } if(locType == 1) { ideal oppLocData = oppose(bsRing, locData); } if(locType == 2) { intvec oppLocData = locData; } list oppResult = oreHelp(oppose(bsRing, g), oppose(bsRing, f), locType, oppLocData, 1); vector oppOreParas = oppResult[1]; ideal oppJ = oppResult[2]; setring bsRing; vector oreParas = oppose(oppRing, oppOreParas); ideal J = oppose(oppRing, oppJ); list result = oreParas, J; return (result); } } 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 = ore(g, f, 0, L, 0); print(rm[1]); rm[2]; rm[1][2]*g-rm[1][1]*f; // geometrical localization ideal p = x-1,y-3; f = dx; g = x^2+y; list rg = ore(g, f, 1, p, 0); 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 = ore(g, f, 2, rat, 0); print(rr[1]); rr[2]; rr[1][2]*g-rr[1][1]*f; // right Ore // monoidal localization g = x; f = dx; L = g; rm = ore(g, f, 0, L, 1); print(rm[1]); rm[2]; f*rm[1][1]-g*rm[1][2]; // geometrical localization g = x+y; f = dx+dy; p = x-1,y-3; rg = ore(g, f, 1, p, 1); print(rg[1]); rg[2]; f*rg[1][1]-g*rg[1][2]; // rational localization rat = 1; f = dx+dy; g = x; rr = ore(g, f, 2, rat, 1); print(rr[1]); rr[2]; f*rr[1][1]-g*rr[1][2]; } //////////////////////////////////////////////////////////////////// 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 representation of frac as a left fraction, - if the left representation of frac is already specified, frac will be returned. EXAMPLE: example convertRightToLeftFraction; shows examples" { checkLocData(locType,locData); checkFraction(frac, locType, locData); if(frac[1] != 0) { return (frac); } vector oreParas = ore(frac[4], frac[3], locType, locData, 0)[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; // geometrical 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 representation of frac as a right fraction, - if the right representation of frac is already specified, frac will be returned. EXAMPLE: example convertLeftToRightFraction; shows examples" { checkLocData(locType,locData); checkFraction(frac, locType, locData); if(frac[4] != 0) { return (frac); } vector oreParas = ore(frac[1], frac[2], locType, locData, 1)[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]; // geometrical 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]; } //////////////////////////////////////////////////////////////////// static proc oreHelp(poly g, poly f, int locType, def locData, int rightOre) "USAGE: ASSUME: RETURN: NOTE: EXAMPLE: " { // TODO: once ringlist-bug [#577] is fixed, use eliminateNC instead of eliminate in rightOre-rational and rightOre-geometrical poly s, r, h; poly gtemp = g; int i, j, n, nze, dif; // modification for monoidal localization if(locType == 0) { if(size(locData) > 0) { // bring the denominator on an equal exponent-level w.r.t. the g_i in locData list divGList = checkMonDivisibility(g, locData); list irredFactors = divGList[2]; intvec expG = divGList[3]; if(divGList[1] != 1) { ERROR("g is not in S"); } else { //"g has divList"; //divGList; } int maxPow = 0; for(i = 1; i <= size(irredFactors); i = i + 1) { if(expG[i] > maxPow) { maxPow = expG[i]; } } //vector testFrac = [g,f]; for(i = 1; i <= size(expG); i = i + 1) { dif = maxPow - expG[i]; g = g * irredFactors[i]^dif; // as for all t in S f = f * irredFactors[i]^dif; // s*f=r*g <=> s*f*t=r*g*t } } } // compute kernel def ker = modulo(f,g); ideal J = ker; ideal I = std(J); // calculate the poly in S if(locType == 0) { // type 0: monoidal localization intvec gle = leadexp(g); intvec fle; int m = deg(f) + 1; // upper bound for(i = 1; i <= size(I); i = i + 1) { fle = leadexp(I[i]); nze = 0; for(j = 1; j <= size(fle); j = j+1) { if(!gle[j] && fle[j]) { nze = 1; break; } } if(!nze) { int mf = 1; while(1) { if(mf * gle - fle >= 0) { break; } mf = mf + 1; } if(mf < m) { m = mf; } } } poly nf = NF(g^m, I); while(m <= deg(f) + 1) { if(nf == 0) { s = g^m; ideal Js = g^m; break; } nf = NF(g*nf, I); m = m + 1; } } if(locType == 1) { // type 1: geometrical localization n = nvars(basering) div 2; if(rightOre == 1) { poly elimVars = 1; for(i = 1; i <= n; i = i + 1) { elimVars = elimVars * var(i); } J = eliminate(I,elimVars); } else { intvec comVars = (n+1)..(2*n); J = eliminateNC(I,comVars); } J = std(J); // compute s, the generator of J with the smallest degree that is not contained in locData ideal K = intersect(J,locData); K = std(K); poly cand; ideal Js; ideal Q = std(locData); for (i = 1; i <= size(J); i++) { if (NF(J[i],Q) != 0) { // J[i] is not contained in locData cand = NF(J[i],K); Js = Js + cand; if (s == 0 || deg(cand) < deg(s)) { s = cand; } } } } if(locType == 2) { // type 2: rational localization n = nvars(basering); if(size(locData) < n) { // there are variables to eliminate intvec modLocData; int check; for(i = 1; i <= n; i = i + 1) { // calculate modLocData = {1...n}\locData check = 0; for(j = 1; j <= size(locData); j = j + 1) { if(i == locData[j]) { check = 1; break; } } if(check == 0) { if(modLocData[1] == 0) { modLocData = i; } else { modLocData = modLocData,i; } } } if(rightOre == 1) { poly elimVars = 1; for(i = 1; i <= size(locData); i = i + 1) { 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); s = J[1]; for(i = 2; i <= size(J); i = i + 1) { if(deg(J[i]) < deg(s)) { s = J[i]; // choose generator with lowest total degree } } ideal Js = J; } // calculate the other poly list divList = division(s*f,g); r = divList[1][1,1]; if(s == 0 || r*g-s*f != 0) { ERROR("no Ore data could be found"); } vector oreParas = [s,r]; list result = oreParas, Js; return (result); } //////////////////////////////////////////////////////////////////// proc addLeftFractions(vector a, vector b, int locType, def locData) "USAGE: addLeftFractions(a, b, locType, locData), vector a, b, int locType, list/vector/intvec locData 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" { checkLocData(locType,locData); checkFraction(a, locType, locData); checkFraction(b, locType, locData); if(a[1] == 0) { a = convertRightToLeftFraction(a, locType, locData); } if(b[1] == 0) { b = convertRightToLeftFraction(b, locType, locData); } if(a[2] == 0) { return (b); } if(b[2] == 0) { return (a); } if(a[1] == b[1]) { return ([a[1],a[2] + b[2],0,0]); } vector oreParas = ore(b[1], a[1], locType, locData, 0)[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); // geometrical 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) "USAGE: multiplyLeftFractions(a, b, locType, locData), vector a, b, int locType, list/vector/intvec locData 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" { checkLocData(locType,locData); checkFraction(a, locType, locData); checkFraction(b, locType, locData); if(a[1] == 0) { a = convertRightToLeftFraction(a, locType, locData); } if(b[1] == 0) { b = convertRightToLeftFraction(b, locType, locData); } if(a[1] == a[2]) { return (b); } if(b[1] == b[2]) { return (a); } if( (a[2] == 0) || (b[2] == 0) ) { return ([1,0,0,1]); } 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); // geometrical 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, if a=b as fractions in the localization specified by locType and locData, - returns 0, otherwise. EXAMPLE: example areEqualLeftFractions; shows examples" { checkLocData(locType, locData); checkFraction(a, locType, locData); checkFraction(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]]; checkFraction(negB, locType, locData); vector result = addLeftFractions(a, negB, locType, locData); if(result[2] == 0) { return (1); } else { return (0); } } 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 rationalForm(def input, int locType, def locData) "USAGE: rationalForm(input, locType, locData), vector/poly input, int locType, list/vector/intvec locData PURPOSE: factor out the monomials in non-invertible variables in the specified rational or geometrical localization ASSUME: - locType is 1 or 2 (geometrical or rational), - input is not the zero polynomial and not the zero vector, - if locType is 2 (rational), locData has to be 1..n for n<=nvars(basering). RETURN: list NOTE: - each entry of the returned list is a list with three entries: intvec alpha, poly d, poly c, - alpha is the exponent vector of the monomial d, - each of these lists represents a term of the input via c*d, - the returned list is ordered by the d-values via the ordering on the basering. EXAMPLE: example rationalForm; shows examples" { if(locType == 0) { ERROR("rationalForm not implemented for monoidal localizations"); } checkLocData(locType,locData); poly nom; if(typeof(input) == "vector") { checkFraction(input, locType, locData); if(input[1] == 0) { input = convertRightToLeftFraction(input, locType, locData); } nom = input[2]; } else { if(typeof(input) == "poly") { nom = input; } else { ERROR("unsupported type: rationalForm expects (poly,int,def) or (vector,int,def), but first type was: " + typeof(input)); } } if(nom == 0) { ERROR("the polynom " + print(nom) + " cannot be written in rational form"); } int n,m; if(locType == 1) { n = nvars(basering) div 2; m = n; } if(locType == 2) { n = size(locData); m = nvars(basering) - n; } int i; list L; poly lc, lt, p, d; intvec le, alpha, beta; list result; while(nom != 0) { lt = lead(nom); le = leadexp(lt); p = 1; d = 1; for(i = 1; i <= n; i = i + 1) { alpha[i] = le[i]; p = p * ((var(i))^le[i]); } for(i = 1; i <= m; i = i + 1) { beta[i] = le[n + i]; d = d * (var(n+i)^le[n+i]); } lc = leadcoef(lt); p = lc * p; L = beta, d, p; if(size(result) == 0) { //"size is zero"; result = insert(result, L); } else { //"size is not zero"; for(i = 1; i <= size(result); i = i + 1) { if(result[i][1] == beta) { //"found equal"; result[i][3] = result[i][3] + p; break; } if(result[i][2] < d) { //"found lower"; result = insert(result, L, i-1); break; } if(i == size(result)) { result = insert(result, L, size(result)); //"didn't find anything, inserted last"; break; } } } nom = nom - lt; } return (result); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,dx,dy),dp; def S = Weyl(); setring S; S; // monoidal localization not supported // geometrical localization ideal p = x-7,x+1; poly g = (4*x+17*y^3) * dx^2*dy + (-3*x+4) * dx + (4*y+9) * 1; rationalForm(g, 1, p); // rational localization intvec rLoc = 1,2; poly denom = x^2+y; poly nom = (4*x+17*y^3) * dx^2*dy + (-3*x+4) * dx + (4*y+9) * 1; vector frac = [denom,nom,0,0]; rationalForm(frac, 2, rLoc); } //////////////////////////////////////////////////////////////////// proc fractionLeadexp(vector frac, int locType, def locData) "USAGE: fractionLeadexp(frac, locType, locData), vector frac, int locType, list/vector/intvec locData PURPOSE: determine the leading exponent of a fraction in a rational localization specified by locData ASSUME: frac is not zero, locType is 2 (rational) RETURN: intvec NOTE: returns the leading exponent of frac EXAMPLE: example fractionLeadexp; shows examples" { if(locType != 2) { ERROR("fractionLeadexp only implemented for rational localizations"); } checkLocData(locType,locData); checkFraction(frac, locType, locData); list L = rationalForm(frac, locType, locData); return (L[1][1]); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,dx,dy),dp; def S = Weyl(); setring S; S; intvec locData = 1,2; vector frac = [x, 4*x*y+2*x*dx-6*y*dx, 0, 0]; fractionLeadexp(frac, 2, locData); } //////////////////////////////////////////////////////////////////// proc fractionLeadmonom(vector frac, int locType, def locData) "USAGE: fractionLeadmonom(frac, locType, locData), vector frac, int locType, list/vector/intvec locData PURPOSE: determine the leading monomial of a fraction in a rational localization specified by locData ASSUME: frac is not zero, locType is 2 (rational) RETURN: poly NOTE: returns the leading monomial of frac EXAMPLE: example fractionLeadmonom; shows examples" { if(locType != 2) { ERROR("fractionLeadmonom only implemented for rational localizations"); } checkLocData(locType,locData); checkFraction(frac, locType, locData); list L = rationalForm(frac, locType, locData); return (L[1][2]); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,dx,dy),dp; def S = Weyl(); setring S; S; intvec locData = 1,2; vector frac = [x, 4*x*y+2*x*dx-6*y*dx, 0, 0]; poly rr = fractionLeadmonom(frac, 2, locData); print(rr); } //////////////////////////////////////////////////////////////////// proc fractionLeadcoef(vector frac, int locType, def locData) "USAGE: fractionLeadcoef(frac, locType, locData), vector frac, int locType, list/vector/intvec locData PURPOSE: determine the leading coefficient of a fraction in a rational localization specified by locData ASSUME: frac is not zero, locType is 2 (rational) RETURN: vector NOTE: returns the leading coefficient of frac EXAMPLE: example fractionLeadcoef; shows examples" { if(locType != 2) { ERROR("fractionLeadcoef only implemented for rational localizations"); } checkLocData(locType,locData); checkFraction(frac, locType, locData); list L = rationalForm(frac, locType, locData); vector result = [frac[1],L[1][3],0,0]; return (result); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,dx,dy),dp; def S = Weyl(); setring S; S; intvec locData = 1,2; vector frac = [x, 4*x*y+2*x*dx-6*y*dx, 0, 0]; vector rr = fractionLeadcoef(frac, 2, locData); print(rr); } //////////////////////////////////////////////////////////////////// proc reduceLeftFraction(vector f, vector g, int locType, def locData) "USAGE: reduceLeftFraction(f, g, locType, locData), vector f, g, int locType, list/vector/intvec locData PURPOSE: completes one step of reducing a fraction in the specified localization ASSUME: locType is 2 (rational), lm(g) | lm(f) RETURN: vector NOTE: returns f-c*m*g, where m satisfies lm(f)=lm(lm(m)*lm(g)) and c satisfies lt(f)=c*lt(m*g) EXAMPLE: example reduceLeftFraction; shows examples" { if(locType != 2) { ERROR("reduceLeftFraction only implemented for rational localizations"); } checkLocData(locType,locData); checkFraction(f, locType, locData); checkFraction(g, locType, locData); if(f[1] == 0) { f = convertRightToLeftFraction(f, locType, locData); } if(f[2] == 0) { ERROR("cannot reduce f = " + print(f) + " because it is zero"); } if(g[1] == 0) { g = convertRightToLeftFraction(g, locType, locData); } if(g[2] == 0) { ERROR("cannot reduce with g = " + print(f) + " because it is zero"); } intvec fExp = fractionLeadexp(f, locType, locData); intvec gExp = fractionLeadexp(g, locType, locData); intvec expDiff = fExp - gExp; int bool = 1; for(int i = 1; i <= size(expDiff); i = i + 1) { if(expDiff[i] < 0) { bool = 0; break; } } if(bool == 0) { ERROR("cannot reduce f = " + print(f) + " with g = " + print(g) + " as lm(g) does not divide lm(f)"); } poly nom = 1; int n = nvars(basering) - size(expDiff); for(i = 1; i <= size(expDiff); i = i + 1) { nom = nom * var(n + i) ^ expDiff[i]; } vector m = [1,nom,0,0]; vector p = multiplyLeftFractions(m, g, locType, locData); vector lcf = fractionLeadcoef(f, locType, locData); vector lcp = fractionLeadcoef(p, locType, locData); vector c = [lcp[2]*lcf[1], -lcp[1]*lcf[2],0,0]; vector t = multiplyLeftFractions(c, p, locType, locData); vector result = addLeftFractions(f, t, locType, locData); return (result); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,dx,dy),dp; def S = Weyl(); setring S; S; intvec locData = 1,2; vector f = [x, 4*x*y+2*x*dx*dy-6*y*dx*dy, 0, 0]; vector g = [y, dx, 0, 0]; vector rr = reduceLeftFraction(f, g, 2, locData); print(rr); } //////////////////////////////////////////////////////////////////// proc cancelLeftFraction(vector frac, int locType, def locData) "USAGE: cancelLeftFraction(frac, locType, locData), vector frac, int locType, list/vector/intvec locData PURPOSE: cancel a given fraction in the specified localization ASSUME: RETURN: vector NOTE: the returned vector is a canceled representation of frac in the localization specified by locType and locData EXAMPLE: example cancelLeftFraction; shows examples " { checkLocData(locType,locData); checkFraction(frac, locType, locData); int i; poly lt; number currentCoef; if(frac[1] == 0) { frac = convertRightToLeftFraction(frac, locType, locData); } // simple cases if(frac[1] == 1) { return (frac); } if(frac[2] == 0) { return ([1,0,0,1]); } if(frac[1] == frac[2]) { return ([1,1,1,1]); } vector result = frac; poly denom = frac[1]; poly nom = frac[2]; // coefficient cancelling number cancelCoef = 0; while(denom != 0) { lt = lead(denom); currentCoef = leadcoef(lt); cancelCoef = gcd(cancelCoef, currentCoef); if(cancelCoef == 1) { break; } denom = denom - lt; } while(nom != 0) { lt = lead(nom); currentCoef = leadcoef(lt); cancelCoef = gcd(cancelCoef, currentCoef); if(cancelCoef == 1) { break; } nom = nom - lt; } denom = frac[1] / cancelCoef; nom = frac[2] / cancelCoef; result = [denom,nom,0,0]; // rational and geometrical localization if(locType >= 1) { list L = rationalForm(result, locType, locData); vector compolys; for(i = 1; i <= size(L); i = i + 1) { compolys = compolys + L[i][3]*gen(i); } def bsRing = basering; list RL = ringlist(basering); RL = delete(RL, 6); RL = delete(RL, 5); def R = ring(RL); setring R; poly gcdR = imap(bsRing, denom); vector compolysR = imap(bsRing, compolys); for(i = 1; i <= size(compolysR); i = i + 1) { gcdR = gcd(gcdR, compolysR[i]); } compolysR = compolysR/gcdR; setring bsRing; poly cancel = imap(R, gcdR); if(cancel == 1) { return (frac); } compolys = imap(R, compolysR); poly newDenom = denom/cancel; poly newNom = 0; for(i = 1; i <= size(L); i = i + 1) { newNom = newNom + compolys[i]*L[i][2]; } result = [newDenom, newNom, 0, 0]; } // monoidal localization if(locType == 0) { list divList = checkMonDivisibility(denom, locData); list irredFactors = divList[2]; intvec expDenom = divList[3]; list locDataExp = divList[4]; poly g = 1; int maxPow = 0; int dif; for(i = 1; i <= size(irredFactors); i = i + 1) { g = g * irredFactors[i]; if(expDenom[i] > maxPow) { maxPow = expDenom[i]; } } for(i = 1; i <= size(expDenom); i = i + 1) { dif = maxPow - expDenom[i]; nom = irredFactors[i]^dif * nom; } def bsRing = basering; list RL = ringlist(basering); list ordList = "c",0; RL[3] = insert(RL[3],ordList); def R = ring(RL); setring R; poly r = imap(bsRing, nom); poly gOrd = imap(bsRing, g); vector gVec = [gOrd,-1]; module gMod = gVec; vector red, rVec; vector resultOrd; for(i = 1; i <= maxPow; i = i + 1) { rVec = [r,0]; red = rightNF(rVec, gMod); if(red[1] != 0) { resultOrd = [gOrd^(maxPow - i + 1), r]; break; } else { r = red[2]; } } setring bsRing; result = imap(R, resultOrd); } // check if cancelled result represents the same fraction as the input if (areEqualLeftFractions(frac, result, locType, locData)) { return (result); } else { return (frac); } } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y,dx,dy),dp; def S = Weyl(); setring S; S; poly cancel, denom, nom; vector frac; // monoidal localization poly g1 = x+3; poly g2 = x*y; list L = g1,g2; poly g = g1^2*g2; poly f = dx; frac = [g^2, g*f, 0, 0]; print(frac); vector rm = cancelLeftFraction(frac, 0, L); print(rm); areEqualLeftFractions(frac, rm, 0, L); // geometrical localization ideal p = x,y-1; cancel = -3*(x^2 + 1); denom = cancel * (y^2 - 3); nom = cancel * x * (3*y*dx*dy + 7*dx^2 - 19); frac = [denom, nom, 0, 0]; vector rg = cancelLeftFraction(frac, 1, p); print(rg); areEqualLeftFractions(frac, rg, 1, p); // rational localization intvec rLoc = 1,2; cancel = 3 * x^2 * (x+y); denom = cancel * (x-y); nom = cancel * (4*x + 3*y*dx^2 - 42*x^3*y*dx^5*dy^2); frac = [denom,nom,0,0]; vector rr = cancelLeftFraction(frac, 2, rLoc); print(rr); areEqualLeftFractions(frac, rr, 2, rLoc); } //////////////////////////////////////////////////////////////////// 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: " { checkLocData(locType,locData); checkFraction(frac, locType, locData); if(frac[1] == 0) { frac = convertRightToLeftFraction(frac); } if(isInS(frac[2], locType, locData)) { return (1); } return (0); } 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); vector p = [1,0]; 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 localization 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: " { checkLocData(locType,locData); checkFraction(frac, locType, locData); if(frac[1] == 0) { frac = convertRightToLeftFraction(frac); } if (!isInvertibleLeftFraction(frac, locType, locData)) { return([1,0,0,1]); } else { return([frac[2],frac[1],frac[4],frac[2]]); } } 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)); vector p = [1,0]; 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)); } ////////////////////////////////////////////////////////////////////