source: git/Singular/LIB/olga.lib @ a1b40a

spielwiese
Last change on this file since a1b40a was a1b40a, checked in by Hans Schoenemann <hannes@…>, 4 years ago
spelling p2
  • Property mode set to 100644
File size: 61.2 KB
Line 
1////////////////////////////////////////////////////////////////
2version="version olga.lib 4.1.2.0 Feb_2019 "; // $Id$
3category="Noncommutative";
4info="
5LIBRARY:   olga.lib  Ore-localization in G-Algebras
6AUTHOR:    Johannes Hoffmann, email: johannes.hoffmann at math.rwth-aachen.de
7
8OVERVIEW:
9Let A be a G-algebra.
10
11Current localization types:
12Type 0: monoidal
13  - represented by a list of polys g_1,...,g_k that have to be contained in a
14    commutative polynomial subring of A generated by a subset of the variables
15    of A
16Type 1: geometric
17  - only for algebras with an even number of variables where the first half
18    induces a commutative polynomial subring B of A
19  - represented by an ideal p, which has to be a prime ideal in B
20Type 2: rational
21  - represented by an intvec v = [i_1,...,i_k] in the range 1..nvars(basering)
22
23Localization data is an int specifying the type and a def with the
24corresponding information.
25
26A fraction is represented as a vector with four entries: [s,r,p,t]
27Here, s^{-1}r is the left fraction representation, pt^{-1} is the right one.
28If s or t is zero, it means that the corresponding representation is not set.
29If both are zero, the fraction is not valid.
30
31A detailed description along with further examples can be found in our paper:
32Johannes Hoffmann, Viktor Levandovskyy:
33    Constructive Arithmetics in Ore Localizations of Domains
34    https://arxiv.org/abs/1712.01773
35
36PROCEDURES:
37locStatus(int, def);
38  report on the status/validity of the given localization data
39testLocData(int, def);
40  check if the given data specifies a denominator set
41isInS(poly, int, def);
42  determine if a polynomial is in a denominator set
43fracStatus(vector frac, int locType, def locData);
44  report on the status/validity of the given fraction wrt. to the given
45  localization data
46testFraction(vector, int, def);
47  check if the given vector is a representation of a fraction in the specified
48  localization
49leftOre(poly, poly, int, def)
50  compute left Ore data
51rightOre(poly, poly, int, def)
52  compute right Ore data
53convertRightToLeftFraction(vector, int, def);
54  determine a left fraction representation of a given fraction
55convertLeftToRightFraction(vector, int, def);
56  determine a right fraction representation of a given fraction
57addLeftFractions(vector, vector, int, def);
58  add two left fractions in the specified localization
59multiplyLeftFractions(vector, vector, int, def);
60  multiply two left fractions in the specified localization
61areEqualLeftFractions(vector, vector, int, def);
62  check if two given fractions are equal
63isInvertibleLeftFraction(vector, int, def);
64  check if a fraction is invertible in the specified localization
65  (NOTE: check description for specific behaviour)
66invertLeftFraction(vector, int, def);
67  invert a fraction in the specified localization
68  (NOTE: check description for specific behaviour)
69isZeroFraction(vector);
70  determine if the given fraction is equal to zero
71isOneFraction(vector);
72  determine if the given fraction is equal to one
73normalizeMonoidal(list);
74  determine a normal form for monoidal localization data
75normalizeRational(intvec);
76  determine a normal form for rational localization data
77testOlga();
78  execute a series of internal testing procedures
79testOlgaExamples();
80  execute the examples of all procedures in this library
81";
82
83LIB "dmodloc.lib"; // for polyVars
84LIB "ncpreim.lib";
85LIB "elim.lib";
86LIB "ncalg.lib";
87//////////////////////////////////////////////////////////////////////
88proc testOlgaExamples()
89"USAGE:   testOlgaExamples()
90PURPOSE: execute the examples of all procedures in this library
91RETURN:  nothing
92NOTE:
93EXAMPLE: "
94{
95    example isInS;
96    example leftOre;
97    example rightOre;
98    example convertRightToLeftFraction;
99    example convertLeftToRightFraction;
100    example addLeftFractions;
101    example multiplyLeftFractions;
102    example areEqualLeftFractions;
103    example isInvertibleLeftFraction;
104    example invertLeftFraction;
105}
106//////////////////////////////////////////////////////////////////////
107proc testOlga()
108"USAGE:   testOlga()
109PURPOSE: execute a series of internal testing procedures
110RETURN:  nothing
111NOTE:
112EXAMPLE: "
113{
114    print("testing olga.lib...");
115    testIsInS();
116    testLeftOre();
117    testRightOre();
118    testAddLeftFractions();
119    testMultiplyLeftFractions();
120    testAreEqualLeftFractions();
121    testConvertLeftToRightFraction();
122    testConvertRightToLeftFraction();
123    testIsInvertibleLeftFraction();
124    testInvertLeftFraction();
125    print("testing complete - olga.lib OK");
126}
127//////////////////////////////////////////////////////////////////////
128proc locStatus(int locType, def locData)
129"USAGE:   locStatus(locType, locData), int locType, list/vector/intvec locData
130PURPOSE: determine the status of a set of localization data
131ASSUME:
132RETURN:  list
133NOTE:    - the first entry is 0 or 1, depending whether the input represents
134           a valid localization
135         - the second entry is a string with a status/error message
136EXAMPLE: example locStatus; shows example"
137{
138    int i;
139    if (locType < 0 || locType > 2) {
140        string invalidTypeString = "invalid localization: type is "
141          + string(locType) + ", valid types are:";
142        invalidTypeString = invalidTypeString + newline
143          + "0 for a monoidal localization";
144        invalidTypeString = invalidTypeString + newline
145          + "1 for a geometric localization";
146        invalidTypeString = invalidTypeString + newline
147          + "2 for a rational localization";
148        return(list(0, invalidTypeString));
149    }
150    string t = typeof(locData);
151    if (t == "none") {
152        return(list(0, "uninitialized or invalid localization:"
153          + " locData has to be defined"));
154    }
155    if (locType == 0)
156    { // monoidal localizations
157        if (t != "list")
158        {
159            return(list(0, "for a monoidal localization, locData has to be of"
160              + " type list, but is of type " + t));
161        }
162        else
163        { // locData is of type list
164            if (size(locData) == 0)
165            {
166                return(list(0, "for a monoidal localization, locData has to be"
167                  + " a non-empty list"));
168            }
169            else
170            { // locData is of type list and has at least one entry
171                if (defined(basering)) {ideal listEntries;}
172                for (i = 1; i <= size(locData); i++)
173                {
174                    t = typeof(locData[i]);
175                    if (t != "poly" && t != "int" && t != "number")
176                    {
177                        return(list(0, "for a monoidal localization, locData"
178                          + " has to be a list of polys, ints or numbers, but"
179                          + " entry " + string(i) + " is " + string(locData[i])
180                          + ", which is of type " + t));
181                    }
182                    else
183                    {
184                      if (defined(basering))
185                      {
186                        if (size(listEntries) == 0)
187                        {
188                            listEntries = locData[i];
189                        }
190                        else
191                        {
192                            listEntries = listEntries, locData[i];
193                        }
194                      }
195                    }
196                }
197                // locData is of type list, has at least one entry and all
198                //   entries are polys
199                if (!defined(basering))
200                {
201                    return(list(0, "for a monoidal localization, the variables"
202                      + " occurring in the polys in locData have to induce a"
203                      + " commutative polynomial subring of basering"));
204                }
205                if (!inducesCommutativeSubring(listEntries))
206                {
207                    return(list(0, "for a monoidal localization, the variables"
208                      + " occurring in the polys in locData have to induce a"
209                      + " commutative polynomial subring of basering"));
210                }
211            }
212        }
213    }
214    if (locType == 1)
215    { // geometric localizations
216        int n = nvars(basering) div 2;
217        if (2*n != nvars(basering))
218        {
219            return(list(0, "for a geometric localization, basering has to have"
220              + " an even number of variables"));
221        }
222        else
223        {
224            int j;
225            for (i = 1; i <= n; i++)
226            {
227                for (j = i + 1; j <= n; j++)
228                {
229                    if (var(i)*var(j) != var(j)*var(i))
230                    {
231                        return(list(0, "for a geometric localization, the"
232                          + " first half of the variables of basering has to"
233                          + " induce a commutative polynomial subring of"
234                          + " basering"));
235                    }
236                }
237            }
238        }
239        if (t != "ideal")
240        {
241            return(list(0, "for a geometric localization, locData has to be of"
242              + " type ideal, but is of type " + t));
243        }
244        for (i = 1; i <= size(locData); i++)
245        {
246            if (!polyVars(locData[i],1..n))
247            {
248                return(list(0, "for a geometric localization, locData has to"
249                + " be an ideal generated by polynomials containing only"
250                + " variables from the first half of the variables"));
251            }
252        }
253    }
254    if (locType == 2)
255    { // rational localizations
256        if (t != "intvec")
257        {
258            return(list(0, "for a rational localization, locData has to be of"
259              + " type intvec, but is of type " + t));
260        }
261        else
262        { // locData is of type intvec
263            if(locData == 0)
264            {
265                return(list(0, "for a rational localization, locData has to be"
266                  + " a non-zero intvec"));
267            }
268            else
269            { // locData is of type intvec and not zero
270                if (!admissibleSub(locData))
271                {
272                    return(list(0, "for a rational localization, the variables"
273                      + " indexed by locData have to generate a sub-G-algebra"
274                      + " of the basering"));
275                }
276            }
277        }
278    }
279    return(list(1, "valid localization"));
280}
281example
282{
283    "EXAMPLE:"; echo = 2;
284    locStatus(42, list(1));
285    def undef;
286    locStatus(0, undef);
287    string s;
288    locStatus(0, s);
289    list L;
290    locStatus(0, L);
291    L = s;
292    print(L);
293    locStatus(0, L);
294    ring w = 0,(x,Dx,y,Dy),dp;
295    def W = Weyl(1);
296    setring W;
297    W;
298    locStatus(0, list(x, Dx));
299    ring R;
300    setring R;
301    R;
302    locStatus(1, s);
303    setring W;
304    locStatus(1, s);
305    ring t = 0,(x,y,Dx,Dy),dp;
306    def T = Weyl();
307    setring T;
308    T;
309    locStatus(1, s);
310    locStatus(1, ideal(Dx));
311    locStatus(2, s);
312    intvec v;
313    locStatus(2, v);
314    locStatus(2, intvec(1,2));
315}
316//////////////////////////////////////////////////////////////////////
317proc testLocData(int locType, def locData)
318"USAGE:   testLocData(locType, locData), int locType,
319         list/vector/intvec locData
320PURPOSE: test if the given data specifies a denominator set wrt. the checks
321         from locStatus
322ASSUME:
323RETURN:  nothing
324NOTE:    throws error if checks were not successful
325EXAMPLE: example testLocData; shows examples"
326{
327    list stat = locStatus(locType, locData);
328    if (!stat[1]) {
329        ERROR(stat[2]);
330    } else {
331        return();
332    }
333}
334example
335{
336    "EXAMPLE:"; echo = 2;
337    ring R; setring R;
338    testLocData(0, list(1)); // correct localization, no error
339    testLocData(42, list(1)); // incorrect localization, results in error
340}
341//////////////////////////////////////////////////////////////////////
342proc isInS(poly p, int locType, def locData, list #)
343"USAGE:   isInS(p, locType, locData(, override)), poly p, int locType,
344         list/vector/intvec locData(, int override)
345PURPOSE: determine if a polynomial is in a denominator set
346ASSUME:
347RETURN:  int
348NOTE:    - returns 0 or 1, depending whether p is in the denominator set
349           specified by locType and locData
350         - if override is set, will not normalize locData (use with care)
351EXAMPLE: example isInS; shows examples"
352{
353    testLocData(locType, locData);
354    if (p == 0) { // the zero polynomial is never a valid denominator
355        return(0);
356    }
357    if (number(p) != 0) {
358        // elements of the coefficient field are always valid denominators
359        return(1);
360    }
361    int override;
362    if (size(#) > 0) {
363        if(typeof(#[1]) == "int") {
364            override = #[1];
365        }
366    }
367    if (locType == 0) {
368        ideal pFactors = commutativeFactorization(p, 1);
369        list locFactors;
370        if (override) {
371            locFactors = locData;
372        } else {
373            locFactors = normalizeMonoidal(locData);
374        }
375        int i, j, foundFactor;
376        for (i = 1; i <= size(pFactors); i++) {
377            foundFactor = 0;
378            for (j = 1; j <= size(locFactors); j++) {
379                if (pFactors[i] == locFactors[j]) {
380                    foundFactor = 1;
381                    break;
382                }
383            }
384            if (!foundFactor) {
385                return(0);
386            }
387        }
388        return(1);
389    }
390    if (locType == 1) {
391        int n = nvars(basering) div 2;
392        ideal I = p;
393        ideal J = std(eliminateNC(I, (n+1)..(2*n)));
394        ideal K = std(intersect(J, locData));
395        int i;
396        for (i = 1; i <= size(J); i++) {
397            if (NF(J[i], K) != 0) {
398                return(1);
399            }
400        }
401    }
402    if (locType == 2) {
403        if (!override) {
404            locData = normalizeRational(locData);
405        }
406        ideal I = p;
407        intvec modLocData = intvecComplement(locData, 1..nvars(basering));
408        I = eliminateNC(I, modLocData);
409        if (size(I)) {
410            return(1);
411        }
412    }
413    return(0);
414}
415example
416{
417    "EXAMPLE:"; echo = 2;
418    ring R = 0,(x,y,Dx,Dy),dp;
419    def S = Weyl();
420    setring S; S;
421    // monoidal localization
422    poly g1 = x^2*y+x+2;
423    poly g2 = y^3+x*y;
424    list L = g1,g2;
425    poly g = g1^2*g2;
426    poly f = g-1;
427    isInS(g, 0, L); // g is in the denominator set
428    isInS(f, 0, L); // f is NOT in the denominator set
429    // geometric localization
430    ideal p = x-1, y-3;
431    g = x^2+y-3;
432    f = (x-1)*g;
433    isInS(g, 1, p); // g is in the denominator set
434    isInS(f, 1, p); // f is NOT in the denominator set
435    // rational localization
436    intvec v = 2;
437    g = y^5+17*y^2-4;
438    f = x*y;
439    isInS(g, 2, v); // g is in the denominator set
440    isInS(f, 2, v); // f is NOT in the denominator set
441}
442//////////////////////////////////////////////////////////////////////
443proc fracStatus(vector frac, int locType, def locData)
444"USAGE:   fracStatus(frac, locType, locData), vector frac, int locType,
445         list/intvec/vector locData
446PURPOSE: determine if the given vector is a representation of a fraction in the
447         specified localization
448ASSUME:
449RETURN:  list
450NOTE:    - the first entry is 0 or 1, depending whether the input is valid
451         - the second entry is a string with a status message
452EXAMPLE: example fracStatus; shows examples"
453{
454    list locStat = locStatus(locType, locData);
455    if (!locStat[1]) { // there is a problem with the localization data
456        return(list(0, "invalid localization in fraction: "+ string(frac)
457          + newline + "      " + locStat[2]));
458    } else { // the specified localization is valid
459        if ((frac[1] == 0) && (frac[4] == 0)) {
460            return(list(0, "vector is not a valid fraction: no denominator"
461              + " specified in " + string(frac)));
462        }
463        if (frac[1] != 0) { // frac has a left representation
464            if (!isInS(frac[1], locType, locData)) {
465                return(list(0, "the left denominator " + string(frac[1])
466                  + " of fraction " + string(frac) + " is not in the"
467                  + " denominator set of type " + string(locType) + " given by "
468                  + string(locData)));
469            }
470        }
471        if (frac[4] != 0) { // frac has a right representation
472            if (!isInS(frac[4], locType, locData)) {
473                return(list(0, "the right denominator " + string(frac[4])
474                  + " of fraction " + string(frac) + " is not in the"
475                  + " denominator set of type " + string(locType) + " given by "
476                  + string(locData)));
477            }
478        }
479        if ((frac[1] != 0) && (frac[4] != 0)) {
480            // frac has left and right representations
481            if (frac[2]*frac[4] != frac[1]*frac[3]) {
482                // the representations are not equal
483                return(list(0, "left and right representation are not equal in:"
484                  + string(frac)));
485            }
486        }
487    }
488    return(list(1, "valid fraction"));
489}
490example
491{
492    "EXAMPLE:"; echo = 2;
493    ring r = QQ[x,y,Dx,Dy];
494    def R = Weyl();
495    setring R;
496    fracStatus([1,0,0,0], 42, list(1));
497    list L = x;
498    fracStatus([0,7,x,0], 0, L);
499    fracStatus([Dx,Dy,0,0], 0, L);
500    fracStatus([0,0,Dx,Dy], 0, L);
501    fracStatus([x,Dx,Dy,x], 0, L);
502    fracStatus([x,Dx,x*Dx+2,x^2], 0, L);
503}
504//////////////////////////////////////////////////////////////////////
505proc testFraction(vector frac, int locType, def locData)
506"USAGE:   testFraction(frac, locType, locData), vector frac, int locType,
507         list/intvec/vector locData
508PURPOSE: test if the given vector is a representation of a fraction in the
509         specified localization wrt. the checks from fracStatus
510ASSUME:
511RETURN:  nothing
512NOTE:    throws error if checks were not successful
513EXAMPLE: example testFraction; shows examples"
514{
515    list stat = fracStatus(frac, locType, locData);
516    if (!stat[1]) {
517        ERROR(stat[2]);
518    } else {
519        return();
520    }
521}
522example
523{
524    "EXAMPLE:"; echo = 2;
525    ring r = QQ[x,y,Dx,Dy];
526    def R = Weyl();
527    setring R;
528    list L = x;
529    vector frac = [x,Dx,x*Dx+2,x^2];
530    testFraction(frac, 0, L); // correct localization, no error
531    frac = [x,Dx,x*Dx,x^2];
532    testFraction(frac, 0, L); // incorrect localization, results in error
533}
534//////////////////////////////////////////////////////////////////////
535proc leftOre(poly s, poly r, int locType, def locData)
536"USAGE:   leftOre(s, r, locType, locData), poly s, r, int locType,
537         list/vector/intvec locData
538PURPOSE: compute left Ore data for a given tuple (s,r)
539ASSUME:  s is in the denominator set determined via locType and locData
540RETURN:  list
541NOTE:  - the first entry of the list is a vector [ts,tr] such that ts*r=tr*s
542       - the second entry of the list is a description of all choices for ts
543EXAMPLE: example leftOre; shows examples"
544{
545    testLocData(locType,locData);
546    locData = normalizeLocalization(locType, locData);
547    if(!isInS(s, locType, locData)) {
548        ERROR("cannot find Ore-parameter since poly " + string(s)
549          + " is not in the denominator set");
550    }
551    return(ore(s, r, locType, locData, 0));
552}
553example
554{
555    "EXAMPLE:"; echo = 2;
556    ring R = 0,(x,y,Dx,Dy),dp;
557    def S = Weyl();
558    setring S; S;
559    // left Ore
560    // monoidal localization
561    poly g1 = x+3;
562    poly g2 = x*y;
563    list L = g1,g2;
564    poly g = g1^2*g2;
565    poly f = Dx;
566    list rm = leftOre(g, f, 0, L);
567    print(rm[1]);
568    rm[2];
569    rm[1][2]*g-rm[1][1]*f;
570    // geometric localization
571    ideal p = x-1, y-3;
572    f = Dx;
573    g = x^2+y;
574    list rg = leftOre(g, f, 1, p);
575    print(rg[1]);
576    rg[2];
577    rg[1][2]*g-rg[1][1]*f;
578    // rational localization
579    intvec rat = 1;
580    f = Dx+Dy;
581    g = x;
582    list rr = leftOre(g, f, 2, rat);
583    print(rr[1]);
584    rr[2];
585    rr[1][2]*g-rr[1][1]*f;
586}
587//////////////////////////////////////////////////////////////////////
588proc rightOre(poly s, poly r, int locType, def locData)
589"USAGE:   rightOre(s, r, locType, locData), poly s, r, int locType,
590         list/vector/intvec locData
591PURPOSE: compute right Ore data for a given tuple (s,r)
592ASSUME:  s is in the denominator set determined via locType and locData
593RETURN:  list
594NOTE:  - the first entry of the list is a vector [ts,tr] such that r*ts=s*tr
595       - the second entry of the list is a description of all choices for ts
596EXAMPLE: example rightOre; shows examples"
597{
598    testLocData(locType,locData);
599    locData = normalizeLocalization(locType, locData);
600    if(!isInS(s, locType, locData)) {
601        ERROR("cannot find Ore-parameter since poly " + string(s)
602          + " is not in the denominator set");
603    }
604    def bsRing = basering;
605    if (locType == 0) {
606        ideal modLocData = locData[1..size(locData)];
607    }
608    def oppRing = opposite(bsRing);
609    setring oppRing;
610    if (locType == 0) {
611        ideal oppModLocData = oppose(bsRing, modLocData);
612        list oppLocData = oppModLocData[1..size(oppModLocData)];
613    }
614    if (locType == 1) {
615        ideal oppLocData = oppose(bsRing,locData);
616    }
617    if (locType == 2) {
618        intvec oppLocData = locData;
619    }
620    poly oppS = oppose(bsRing, s);
621    poly oppR = oppose(bsRing, r);
622    list oppResult = ore(oppS, oppR, locType, oppLocData, 1);
623    vector oppOreParas = oppResult[1];
624    ideal oppJ = oppResult[2];
625    setring bsRing;
626    return(list(oppose(oppRing, oppOreParas), oppose(oppRing, oppJ)));
627}
628example
629{
630    "EXAMPLE:"; echo = 2;
631    ring R = 0,(x,y,Dx,Dy),dp;
632    def S = Weyl();
633    setring S; S;
634    // monoidal localization
635    poly g1 = x+3;
636    poly g2 = x*y;
637    list L = g1,g2;
638    poly g = g1^2*g2;
639    poly f = Dx;
640    list rm = rightOre(g, f, 0, L);
641    print(rm[1]);
642    rm[2];
643    g*rm[1][2]-f*rm[1][1];
644    // geometric localization
645    ideal p = x-1, y-3;
646    f = Dx;
647    g = x^2+y;
648    list rg = rightOre(g, f, 1, p);
649    print(rg[1]);
650    rg[2];
651    g*rg[1][2]-f*rg[1][1];
652    // rational localization
653    intvec rat = 1;
654    f = Dx+Dy;
655    g = x;
656    list rr = rightOre(g, f, 2, rat);
657    print(rr[1]);
658    rr[2];
659    g*rr[1][2]-f*rr[1][1];
660}
661//////////////////////////////////////////////////////////////////////
662static proc ore(poly s, poly r, int locType, def locData, int rightOre) {
663    // TODO: once ringlist-bug [#577] is PROPERLY fixed (probability: low),
664    //   use eliminateNC instead of eliminate in rightOre-rational and
665    //   rightOre-geometric
666    // problem: the ordering on the opposite algebra is not recognized as
667    //   valid for eliminateNC, even though it is
668    if (r == 0) {
669        return([1,0], ideal(1));
670    }
671    int i, j;
672    // modification for monoidal localization
673    if (locType == 0) {
674        // computations will be carried out in the localization at the product
675        //   of all irreducible factors of s, thus all factors have to be
676        //   raised to the highest power occurring in the factorization of s
677        list factorList = commutativeFactorization(s);
678        ideal factors = factorList[1];
679        int maxPow = 0;
680        intvec exponents = factorList[2];
681        for (i = 1; i <= size(exponents); i++) {
682            if (exponents[i] > maxPow) {
683                maxPow = exponents[i];
684            }
685        }
686        int expDiff;
687        for (i = 1; i <= size(exponents); i++) {
688            expDiff = maxPow - exponents[i];
689            s = s * factors[i]^expDiff;
690            r = r * factors[i]^expDiff;
691        }
692    }
693    // compute kernel of x maps to x*r+R*s
694    ideal J = modulo(r, s);
695    ideal I = std(J);
696    // compute the poly in S
697    poly ts;
698    if (locType == 0) { // type 0: monoidal localization
699        intvec sle = leadexp(s);
700        intvec rle;
701        int m = deg(r) + 1; // upper bound
702        int dividesLeadingMonomial;
703        // find the minimal possible exponent m such that lm(f_i)|lm(s)^m
704        for (i = 1; i <= size(I); i++) {
705            rle = leadexp(I[i]);
706            dividesLeadingMonomial = 1;
707            for (j = 1; j <= size(rle); j++) {
708                if (sle[j] == 0 && rle[j] != 0) {
709                    // lm(f_i) divides no power of lm(s)
710                    dividesLeadingMonomial = 0;
711                    break;
712                }
713            }
714            if (dividesLeadingMonomial) {
715                int mf = 1;
716                while (1) {
717                    if (mf * sle - rle >= 0) {
718                        break; // now lm(r_i) divides lm(s)^mf
719                    } else {
720                        mf++;
721                    }
722                }
723                if (mf < m) {
724                    m = mf; // mf is lower than the current minimum m
725                }
726            }
727        }
728        poly nf = NF(s^m, I);
729        for (; m <= deg(r) + 1; m++) { // sic: m is initialized beforehand
730            if (nf == 0) { // s^m is in I
731                ts = s^m;
732                ideal Js = s^m;
733                break;
734            } else {
735                nf = NF(s*nf, I);
736            }
737        }
738    }
739    if (locType == 1) { // type 1: geometric localization
740        int n = nvars(basering) div 2;
741        if (rightOre == 1) {
742            // in the right Ore setting, the order of variables is inverted
743            poly elimVars = 1;
744            for (i = 1; i <= n; i++) {
745                elimVars = elimVars * var(i);
746            }
747            J = eliminate(I, elimVars);
748        } else {
749            J = eliminateNC(I, intvec((n+1)..(2*n)));
750        }
751        J = std(J);
752        // compute ts, the generator of J with the smallest degree that does
753        //   not vanish at locData
754        ideal K = intersect(J, locData);
755        K = std(K);
756        poly h, cand;
757        ideal Js;
758        ideal Q = std(locData);
759        for (i = 1; i <= size(J); i++) {
760            h = NF(J[i],Q);
761            if (h != 0) {
762                cand = NF(J[i],K);
763                Js = Js + cand;
764                if (ts == 0 || deg(cand) < deg(ts)) {
765                    ts = cand;
766                }
767            }
768        }
769    }
770    if (locType == 2) { // type 2: rational localization
771        int n = nvars(basering);
772        if (size(locData) < n) { // there are variables to eliminate
773            intvec modLocData;
774            int check;
775            modLocData = intvecComplement(locData, 1..n);
776            // calculate modLocData = {1...n}\locData
777            if (rightOre == 1) {
778                // in the right Ore setting, the order of variables is inverted
779                poly elimVars = 1;
780                for (i = 1; i <= size(locData); i++) {
781                    elimVars = elimVars * var(locData[i]);
782                }
783                J = eliminate(I, elimVars);
784            } else {
785                J = eliminateNC(I, modLocData);
786            }
787
788        } else { // no variables to eliminate (total localization)
789            J = I;
790        }
791        J = std(J);
792        ts = J[1];
793        for (i = 2; i <= size(J); i++) {
794            if (deg(J[i]) < deg(ts)) {
795                ts = J[i]; // choose generator with lowest total degree
796            }
797        }
798        ideal Js = J;
799    }
800    // calculate the other poly
801    poly tr = division(ts*r, s)[1][1,1];
802    if (ts == 0 || tr*s-ts*r != 0) {
803        string s;
804        if (rightOre) {
805            s = "right";
806        } else {
807            s = "left";
808        }
809        ERROR("no " + s + " Ore data could be found for s=" + string(s)
810          + " and r=" + string(r));
811    }
812    vector oreParas = [ts,tr];
813    list result = oreParas, Js;
814    return (result);
815}
816//////////////////////////////////////////////////////////////////////
817proc convertRightToLeftFraction(vector frac, int locType, def locData)
818"USAGE:   convertRightToLeftFraction(frac, locType, locData),
819         vector frac, int locType, list/vector/intvec locData
820PURPOSE: determine a left fraction representation of a given fraction
821ASSUME:
822RETURN:  vector
823NOTE:    - the returned vector contains a repr. of frac as a left fraction
824         - if the left representation of frac is already specified,
825           frac will be returned.
826EXAMPLE: example convertRightToLeftFraction; shows examples"
827{
828    testLocData(locType, locData);
829    testFraction(frac, locType, locData);
830    if (frac[1] != 0) { // frac already has a left representation
831        return (frac);
832    } else { // frac has no left representation, but a right one
833        vector oreParas = leftOre(frac[4], frac[3], locType, locData)[1];
834        vector result = [oreParas[1], oreParas[2], frac[3], frac[4]];
835        return (result);
836    }
837}
838example
839{
840    "EXAMPLE:"; echo = 2;
841    ring R = 0,(x,y,Dx,Dy),dp;
842    def S = Weyl();
843    setring S; S;
844    // monoidal localization
845    poly g1 = x+3;
846    poly g2 = x*y;
847    list L = g1,g2;
848    poly g = g1^2*g2;
849    poly f = Dx;
850    vector fracm = [0,0,f,g];
851    vector rm = convertRightToLeftFraction(fracm, 0, L);
852    print(rm);
853    rm[2]*g-rm[1]*f;
854    // geometric localization
855    ideal p = x-1, y-3;
856    f = Dx;
857    g = x^2+y;
858    vector fracg = [0,0,f,g];
859    vector rg = convertRightToLeftFraction(fracg, 1, p);
860    print(rg);
861    rg[2]*g-rg[1]*f;
862    // rational localization
863    intvec rat = 1;
864    f = Dx+Dy;
865    g = x;
866    vector fracr = [0,0,f,g];
867    vector rr = convertRightToLeftFraction(fracr, 2, rat);
868    print(rr);
869    rr[2]*g-rr[1]*f;
870}
871//////////////////////////////////////////////////////////////////////
872proc convertLeftToRightFraction(vector frac, int locType, def locData)
873"USAGE:   convertLeftToRightFraction(frac, locType, locData), vector frac,
874         int locType, list/vector/intvec locData
875PURPOSE: determine a right fraction representation of a given fraction
876ASSUME:
877RETURN:  vector
878NOTE:    - the returned vector contains a repr. of frac as a right fraction,
879         - if the right representation of frac is already specified,
880           frac will be returned.
881EXAMPLE: example convertLeftToRightFraction; shows examples"
882{
883    testLocData(locType, locData);
884    testFraction(frac, locType, locData);
885    if (frac[4] != 0) { // frac already has a right representation
886        return (frac);
887    } else { // frac has no right representation, but a left one
888        vector oreParas = rightOre(frac[1], frac[2], locType, locData)[1];
889        vector result = [frac[1], frac[2], oreParas[2], oreParas[1]];
890        return (result);
891    }
892}
893example
894{
895    "EXAMPLE:"; echo = 2;
896    ring R = 0,(x,y,Dx,Dy),dp;
897    def S = Weyl();
898    setring S; S;
899    // monoidal localization
900    poly g = x;
901    poly f = Dx;
902    vector fracm = [g,f,0,0];
903    list L = g;
904    vector rm = convertLeftToRightFraction(fracm, 0, L);
905    print(rm);
906    f*rm[4]-g*rm[3];
907    // geometric localization
908    g = x+y;
909    f = Dx+Dy;
910    vector fracg = [g,f,0,0];
911    ideal p = x-1, y-3;
912    vector rg = convertLeftToRightFraction(fracg, 1, p);
913    print(rg);
914    f*rg[4]-g*rg[3];
915    // rational localization
916    intvec rat = 1;
917    f = Dx+Dy;
918    g = x;
919    vector fracr = [g,f,0,0];
920    vector rr = convertLeftToRightFraction(fracr, 2, rat);
921    print(rr);
922    f*rr[4]-g*rr[3];
923}
924//////////////////////////////////////////////////////////////////////
925proc isZeroFraction(vector frac)
926"USAGE:   isZeroFraction(frac), vector frac
927PURPOSE: determine if the vector frac represents zero
928ASSUME:  frac is a valid fraction
929RETURN:  int
930NOTE:    returns 1, if frac == 0; 0 otherwise
931EXAMPLE: example isZeroFraction; shows examples"
932{
933    if (frac[1] == 0 && frac[4] == 0) {
934        return(0);
935    }
936    if (frac[1] == 0) { // frac has no left representation
937        if (frac[3] == 0) { // the right representation of frac is zero
938            return(1);
939        }
940    } else { // frac has a left representation
941        if (frac[2] == 0) { // the left representation of frac is zero
942            return(1);
943        }
944    }
945    return(0);
946}
947example
948{
949    "EXAMPLE:"; echo = 2;
950    ring R = 0,(x,y,Dx,Dy),dp;
951    def S = Weyl();
952    setring S; S;
953    isZeroFraction([42,0,0,0]);
954    isZeroFraction([0,0,Dx,3]);
955    isZeroFraction([1,1,1,1]);
956}
957//////////////////////////////////////////////////////////////////////
958proc isOneFraction(vector frac)
959"USAGE:   isOneFraction(frac), vector frac
960PURPOSE: determine if the vector frac represents one
961ASSUME:  frac is a valid fraction
962RETURN:  int
963NOTE:    1, if frac == 1; 0 otherwise
964EXAMPLE: example isOneFraction; shows examples"
965{
966    if (frac[1] == 0 && frac[4] == 0) {
967        return(0);
968    }
969    if (frac[1] == 0) { // frac has no left representation
970        if (frac[3] == frac[4]) { // the right representation of frac is zero
971            return(1);
972        }
973    } else { // frac has a left representation
974        if (frac[2] == frac[1]) { // the left representation of frac is zero
975            return(1);
976        }
977    }
978    return(0);
979}
980example
981{
982    "EXAMPLE:"; echo = 2;
983    ring R = 0,(x,y,Dx,Dy),dp;
984    def S = Weyl();
985    setring S; S;
986    isOneFraction([42,42,0,0]);
987    isOneFraction([0,0,Dx,3]);
988    isOneFraction([1,0,0,1]);
989}
990////////// arithmetic ////////////////////////////////////////////////
991proc addLeftFractions(vector a, vector b, int locType, def locData, list #)
992"USAGE:   addLeftFractions(a, b, locType, locData(, override)),
993         vector a, b, int locType, list/vector/intvec locData(, int override)
994PURPOSE: add two left fractions in the specified localization
995ASSUME:
996RETURN:  vector
997NOTE:    the returned vector is the sum of a and b as fractions in the
998         localization specified by locType and locData.
999EXAMPLE: example addLeftFractions; shows examples"
1000{
1001    int override = 0;
1002    if (size(#) > 1) {
1003        if(typeof(#[1]) == "int") {
1004            override = #[1];
1005        }
1006    }
1007    if (!override) {
1008        testLocData(locType, locData);
1009        testFraction(a, locType, locData);
1010        testFraction(b, locType, locData);
1011    }
1012    // check for a shortcut
1013    if (isZeroFraction(a)) {
1014        return(b);
1015    }
1016    if (isZeroFraction(b)) {
1017        return(a);
1018    }
1019    if (a[1] == 0) { // a has no left representation
1020        a = convertRightToLeftFraction(a, locType, locData);
1021    }
1022    if (b[1] == 0) { // b has no left representation
1023        b = convertRightToLeftFraction(b, locType, locData);
1024    }
1025    if (a[1] == b[1]) { // a and b have the same left denominator
1026        return ([a[1],a[2] + b[2],0,0]);
1027    }
1028    // no shortcut found, use regular method
1029    vector oreParas = leftOre(b[1], a[1], locType, locData)[1];
1030    vector result = [oreParas[1]*a[1],oreParas[1]*a[2]+oreParas[2]*b[2],0,0];
1031    return (result);
1032}
1033example
1034{
1035    "EXAMPLE:"; echo = 2;
1036    ring R = 0,(x,y,Dx,Dy),dp;
1037    def S = Weyl();
1038    setring S; S;
1039    // monoidal localization
1040    poly g1 = x+3;
1041    poly g2 = x*y+y;
1042    list L = g1,g2;
1043    poly s1 = g1;
1044    poly s2 = g2;
1045    poly r1 = Dx;
1046    poly r2 = Dy;
1047    vector frac1 = [s1,r1,0,0];
1048    vector frac2 = [s2,r2,0,0];
1049    vector rm = addLeftFractions(frac1, frac2, 0, L);
1050    print(rm);
1051    // geometric localization
1052    ideal p = x-1, y-3;
1053    vector rg = addLeftFractions(frac1, frac2, 1, p);
1054    print(rg);
1055    // rational localization
1056    intvec v = 2;
1057    s1 = y^2+y+1;
1058    s2 = y-2;
1059    r1 = Dx;
1060    r2 = Dy;
1061    frac1 = [s1,r1,0,0];
1062    frac2 = [s2,r2,0,0];
1063    vector rr = addLeftFractions(frac1, frac2, 2, v);
1064    print(rr);
1065}
1066//////////////////////////////////////////////////////////////////////
1067proc multiplyLeftFractions(vector a, vector b, int locType, def locData, list #)
1068"USAGE:   multiplyLeftFractions(a, b, locType, locData(, override)),
1069         vector a, b, int locType, list/vector/intvec locData, int override
1070PURPOSE: multiply two left fractions in the specified localization
1071ASSUME:
1072RETURN:  vector
1073NOTE:    the returned vector is the product of a and b as fractions in the
1074         localization specified by locType and locData.
1075EXAMPLE: example multiplyLeftFractions; shows examples"
1076{
1077    int override = 0;
1078    if (size(#) > 1) {
1079        if(typeof(#[1]) == "int") {
1080            override = #[1];
1081        }
1082    }
1083    if (!override) {
1084        testLocData(locType, locData);
1085        testFraction(a, locType, locData);
1086        testFraction(b, locType, locData);
1087    }
1088    // check for a shortcut
1089    if (isZeroFraction(a) || isZeroFraction(b)) {
1090        return([1,0,0,1]);
1091    }
1092    if (isOneFraction(a)) {
1093        return(b);
1094    }
1095    if (isOneFraction(b)) {
1096        return(a);
1097    }
1098    if(a[1] == 0) {
1099        a = convertRightToLeftFraction(a, locType, locData);
1100    }
1101    if(b[1] == 0) {
1102        b = convertRightToLeftFraction(b, locType, locData);
1103    }
1104    if( (a[2] == 0) || (b[2] == 0) ) {
1105        return ([1,0,0,1]);
1106    }
1107    if (a[2]*b[1] == b[1]*a[2]) { // trivial solution of the Ore condition
1108        return([b[1]*a[1],a[2]*b[2]]);
1109    }
1110    // no shortcut found, use regular method
1111    vector oreParas = ore(b[1], a[2], locType, locData, 0)[1];
1112    vector result = [oreParas[1]*a[1],oreParas[2]*b[2],0,0];
1113    return (result);
1114}
1115example
1116{
1117    "EXAMPLE:"; echo = 2;
1118    ring R = 0,(x,y,Dx,Dy),dp;
1119    def S = Weyl();
1120    setring S; S;
1121    // monoidal localization
1122    poly g1 = x+3;
1123    poly g2 = x*y+y;
1124    list L = g1,g2;
1125    poly s1 = g1;
1126    poly s2 = g2;
1127    poly r1 = Dx;
1128    poly r2 = Dy;
1129    vector frac1 = [s1,r1,0,0];
1130    vector frac2 = [s2,r2,0,0];
1131    vector rm = multiplyLeftFractions(frac1, frac2, 0, L);
1132    print(rm);
1133    // geometric localization
1134    ideal p = x-1, y-3;
1135    vector rg = multiplyLeftFractions(frac1, frac2, 1, p);
1136    print(rg);
1137    // rational localization
1138    intvec v = 2;
1139    s1 = y^2+y+1;
1140    s2 = y-2;
1141    r1 = Dx;
1142    r2 = Dy;
1143    frac1 = [s1,r1,0,0];
1144    frac2 = [s2,r2,0,0];
1145    vector rr1 = multiplyLeftFractions(frac1, frac2, 2, v);
1146    print(rr1);
1147    vector rr2 = multiplyLeftFractions(frac2, frac1, 2, v);
1148    print(rr2);
1149    areEqualLeftFractions(rr1, rr2, 2, v);
1150}
1151//////////////////////////////////////////////////////////////////////
1152proc areEqualLeftFractions(vector a, vector b, int locType, def locData)
1153"USAGE:   areEqualLeftFractions(a, b, locType, locData), vector a, b,
1154         int locType, list/vector/intvec locData
1155PURPOSE: check if two given fractions are equal
1156ASSUME:
1157RETURN:  int
1158NOTE:    returns 1 or 0, depending whether a=b as fractions in the
1159         localization specified by locType and locData
1160EXAMPLE: example areEqualLeftFractions; shows examples"
1161{
1162    testLocData(locType, locData);
1163    testFraction(a, locType, locData);
1164    testFraction(b, locType, locData);
1165    if(a[1] == 0) {
1166        a = convertRightToLeftFraction(a, locType, locData);
1167    }
1168    if(b[1] == 0) {
1169        b = convertRightToLeftFraction(b, locType, locData);
1170    }
1171    vector negB = [b[1], -b[2], -b[3], b[4]];
1172    //testFraction(negB, locType, locData); // unnecessary check
1173    vector result = addLeftFractions(a, negB, locType, locData);
1174    return(isZeroFraction(result));
1175}
1176example
1177{
1178    "EXAMPLE:"; echo = 2;
1179    ring R = 0,(x,y,Dx,Dy),dp;
1180    def S = Weyl();
1181    setring S; S;
1182    // monoidal
1183    poly g1 = x*y+3;
1184    poly g2 = y^3;
1185    list L = g1,g2;
1186    poly s1 = g1;
1187    poly s2 = s1*g2;
1188    poly s3 = s2;
1189    poly r1 = Dx;
1190    poly r2 = g2*r1;
1191    poly r3 = s1*r1+3;
1192    vector fracm1 = [s1,r1,0,0];
1193    vector fracm2 = [s2,r2,0,0];
1194    vector fracm3 = [s3,r3,0,0];
1195    areEqualLeftFractions(fracm1, fracm2, 0, L);
1196    areEqualLeftFractions(fracm1, fracm3, 0, L);
1197    areEqualLeftFractions(fracm2, fracm3, 0, L);
1198}
1199//////////////////////////////////////////////////////////////////////
1200proc isInvertibleLeftFraction(vector frac, int locType, def locData)
1201"USAGE:   isInvertibleLeftFraction(frac, locType, locData), vector frac,
1202         int locType, list/vector/intvec locData
1203PURPOSE: check if a fraction is invertible in the specified localization
1204ASSUME:
1205RETURN:  int
1206NOTE:    - returns 1, if the numerator of frac is in the denominator set,
1207         - returns 0, otherwise (NOTE: this does NOT mean that the fraction is
1208           not invertible, it just means it could not be determined by the
1209           method above).
1210EXAMPLE: example isInvertibleLeftFraction; shows examples"
1211{
1212    testLocData(locType, locData);
1213    testFraction(frac, locType, locData);
1214    locData = normalizeLocalization(locType, locData);
1215    if(frac[1] != 0) { // frac has a left representation
1216        return(isInS(frac[2], locType, locData));
1217    } else { // frac has no left, but a right representation
1218        return(isInS(frac[3], locType, locData));
1219    }
1220}
1221example
1222{
1223    "EXAMPLE:"; echo = 2;
1224    ring R = 0,(x,y,Dx,Dy),dp;
1225    def S = Weyl();
1226    setring S; S;
1227    poly g1 = x+3;
1228    poly g2 = x*y;
1229    list L = g1,g2;
1230    vector frac = [g1*g2, 17, 0, 0];
1231    isInvertibleLeftFraction(frac, 0, L);
1232    ideal p = x-1, y;
1233    frac = [g1, x, 0, 0];
1234    isInvertibleLeftFraction(frac, 1, p);
1235    intvec rat = 1,2;
1236    frac = [g1*g2, Dx, 0, 0];
1237    isInvertibleLeftFraction(frac, 2, rat);
1238}
1239//////////////////////////////////////////////////////////////////////
1240proc invertLeftFraction(vector frac, int locType, def locData)
1241"USAGE:   invertLeftFraction(frac, locType, locData), vector frac, int locType,
1242         list/vector/intvec locData
1243PURPOSE: invert a fraction in the specified localization
1244ASSUME:  frac is invertible in the loc. specified by locType and locData
1245RETURN:  vector
1246NOTE:    - returns the multiplicative inverse of frac in the localization
1247           specified by locType and locData,
1248         - throws error if frac is not invertible (NOTE: this does NOT mean
1249           that the fraction is not invertible, it just means it could not be
1250           determined by the method listed above).
1251EXAMPLE: example invertLeftFraction; shows examples"
1252{
1253    // standard tests are done by isInvertibleLeftFraction
1254    if (isInvertibleLeftFraction(frac, locType, locData)) {
1255        return([frac[2],frac[1],frac[4],frac[3]]);
1256    } else {
1257        return([1,0,0,1]);
1258    }
1259}
1260example
1261{
1262    "EXAMPLE:"; echo = 2;
1263    ring R = 0,(x,y,Dx,Dy),dp;
1264    def S = Weyl();
1265    setring S; S;
1266    poly g1 = x+3;
1267    poly g2 = x*y;
1268    list L = g1,g2;
1269    vector frac = [g1*g2, 17, 0, 0];
1270    print(invertLeftFraction(frac, 0, L));
1271    ideal p = x-1, y;
1272    frac = [g1, x, 0, 0];
1273    print(invertLeftFraction(frac, 1, p));
1274    intvec rat = 1,2;
1275    frac = [g1*g2, y, 0, 0];
1276    print(invertLeftFraction(frac, 2, rat));
1277}
1278//////////////////////////////////////////////////////////////////////
1279proc normalizeMonoidal(list L)
1280"USAGE:    normalizeMonoidal(L), list L
1281PURPOSE:  compute a normal form of monoidal localization data
1282RETURN:   list
1283NOTE:     given a list of polys, returns a list of all unique factors appearing
1284          in the given polys
1285EXAMPLE:  example normalizeMonoidal; shows examples"
1286{
1287    ideal allFactors;
1288    int i;
1289    for (i = 1; i <= size(L); i++) {
1290        allFactors = allFactors, commutativeFactorization(L[i],1);
1291    }
1292    allFactors = simplify(allFactors,1+2+4);
1293    // simplify: divide by leading coefficients (1),
1294    //   purge zero generators (2), purge double entries (4)
1295    ideal rev = sort(allFactors)[1]; // sort sorts ascendingly
1296    return(list(rev[size(rev)..1])); // reverse order and cast to list
1297}
1298example
1299{
1300    "EXAMPLE:"; echo = 2;
1301    ring R = 0,(x,y,Dx,Dy),dp;
1302    def S = Weyl(); setring S;
1303    list L = x^2*y^3, (x+1)*(x*y-3*y^2+1);
1304    L = normalizeMonoidal(L);
1305    print(L);
1306}
1307//////////////////////////////////////////////////////////////////////
1308proc normalizeRational(intvec v)
1309"USAGE:    normalizeRational(v), intvec v
1310PURPOSE:  compute a normal form of rational localization data
1311RETURN:   intvec
1312NOTE:     purges double entries and sorts ascendingly
1313EXAMPLE:  example normalizeRational; shows examples"
1314{
1315    int n = nvars(basering);
1316    int i;
1317    intvec result;
1318    intvec occurring = 0:n;
1319    for (i = 1; i <= size(v); i++) {
1320        occurring[v[i]] = 1;
1321    }
1322    for (i = 1; i <= size(occurring); i++) {
1323        if (occurring[i]) {
1324            if (result == 0) {
1325                result = i;
1326            } else {
1327                result = result, i;
1328            }
1329        }
1330    }
1331    return(result);
1332}
1333example
1334{
1335    "EXAMPLE:"; echo = 2;
1336    ring R; setring R;
1337    intvec v = 9,5,9,3,1,5;
1338    v = normalizeRational(v);
1339    v;
1340}
1341////////// internal functions ////////////////////////////////////////
1342static proc inducesCommutativeSubring(def input)
1343{
1344    ideal vars = variables(input);
1345    int i, j;
1346    for (i = 1; i <= size(vars); i++) {
1347        for (j = i + 1; j <= size(vars); j++) {
1348            if (vars[i]*vars[j] != vars[j]*vars[i]) {
1349                return(0);
1350            }
1351        }
1352    }
1353    return(1);
1354}
1355//////////////////////////////////////////////////////////////////////
1356static proc commutativeFactorization(poly p, list #)
1357"USAGE:    commutativeFactorization(p[, #]), poly p[, list #]
1358PURPOSE:  compute a factorization of p ignoring non-commutative relations
1359RETURN:   list or ideal
1360NOTE:     the optional parameter is passed to factorize after changing to a
1361          commutative ring, the result of factorize is transferred back to
1362          basering
1363SEE ALSO: factorize
1364EXAMPLE:  "
1365{
1366    int factorType = 0;
1367    if (size(#)  > 0) {
1368        if (typeof(#[1]) == "int") {
1369            factorType = #[1];
1370        }
1371    }
1372    list RL = ringlist(basering);
1373    if (size(RL) > 4) {
1374        def bsRing = basering;
1375        RL = RL[1..4];
1376        def commRing = ring(RL);
1377        setring commRing;
1378        poly commP = imap(bsRing, p);
1379        if (factorType == 1) {
1380            ideal commFactors = factorize(commP, 1);
1381            setring bsRing;
1382            return(imap(commRing, commFactors));
1383        } else {
1384            list commFac = factorize(commP, factorType);
1385            intvec exponents = commFac[2];
1386            ideal commFactors = commFac[1];
1387            setring bsRing;
1388            list result;
1389            result[1] = imap(commRing, commFactors);
1390            result[2] = exponents;
1391            return(result);
1392        }
1393    } else {
1394        return(factorize(p, factorType));
1395    }
1396}
1397//////////////////////////////////////////////////////////////////////
1398static proc normalizeLocalization(int locType, def locData) {
1399    if (locType == 0) {
1400        return(normalizeMonoidal(locData));
1401    }
1402    if (locType == 1) {
1403        return(std(locData));
1404    }
1405    if (locType == 2) {
1406        return(normalizeRational(locData));
1407    }
1408    return(locData);
1409}
1410//////////////////////////////////////////////////////////////////////
1411static proc intvecComplement(intvec v, intvec w) { // complement of v in w as sets
1412    intvec result;
1413    int i;
1414    int j;
1415    int foundMatch;
1416    for (i = 1; i <= size(w); i++) {
1417        foundMatch = 0;
1418        for (j = 1; j <= size(v); j++) {
1419            if (w[i] == v[j]) {
1420                foundMatch = 1;
1421                break;
1422            }
1423        }
1424        if (!foundMatch) { // v[i] is not in w
1425            if (result == 0) {
1426                result = w[i];
1427            } else {
1428                result = result, w[i];
1429            }
1430        }
1431    }
1432    return(result);
1433}
1434//////////////////////////////////////////////////////////////////////
1435////////// internal testing procedures ///////////////////////////////
1436static proc testIsInS()
1437{
1438    print("  testing isInS...");
1439    ring R = 0,(x,y,Dx,Dy),dp;
1440    def S = Weyl();
1441    setring S;
1442    // monoidal localization
1443    poly g1 = x^2*y+x+2;
1444    poly g2 = y^3+x*y;
1445    list L = g1, g2;
1446    poly g = g1^2*g2;
1447    if (!isInS(g, 0, L)) {
1448        ERROR("Weyl monoidal isInS direct positive failed");
1449    }
1450    if (!isInS(y^2+x, 0, L)) {
1451        ERROR("Weyl monoidal isInS indirect positive failed");
1452    }
1453    if (isInS(g-1, 0, L)) {
1454        ERROR("Weyl monoidal isInS negative failed");
1455    }
1456    // geometric localization
1457    ideal p = x-1, y-3;
1458    g = x^2+y-3;
1459    if (!isInS(g, 1, p)) {
1460        ERROR("Weyl geometric isInS positive failed");
1461    }
1462    if (isInS((x-1)*g, 1, p)) {
1463        ERROR("Weyl geometric isInS negative failed");
1464    }
1465    // rational localization
1466    intvec v = 2;
1467    if (!isInS(y^5+17*y^2-4, 2, v)) {
1468        ERROR("Weyl rational isInS positive failed");
1469    }
1470    if (isInS(x*y, 2, v)) {
1471        ERROR("Weyl rational isInS negative failed");
1472    }
1473    print("    isInS OK");
1474}
1475//////////////////////////////////////////////////////////////////////
1476static proc testLeftOre() {
1477    print("  testing leftOre...");
1478    // Weyl
1479    ring W = 0,(x,y,Dx,Dy),dp;
1480    def ncW = Weyl();
1481    setring ncW;
1482    //// monoidal localization
1483    poly g1 = x+3;
1484    poly g2 = x*y;
1485    list L = g1,g2;
1486    poly g = g1^2*g2;
1487    poly f = Dx;
1488    vector rm = leftOre(g, f, 0, L)[1];
1489    if (rm[1] == 0 || rm[2]*g-rm[1]*f != 0) {
1490        ERROR("Weyl monoidal left Ore failed");
1491    }
1492    //// geometric localization
1493    vector p1,p2,p3,p4 = [1,3],[0,0],[-1,-2],[10,-180];
1494    vector rg;
1495    ideal p;
1496    list pVecs = p1,p2,p3,p4;
1497    f = Dx;
1498    g = x^2+y+3;
1499    for(int i = 1; i <= 4; i = i + 1) {
1500        p = x-pVecs[i][1], y-pVecs[i][2];
1501        rg = leftOre(g, f, 1, p)[1];
1502        if (rg[1] == 0 || rg[2]*g-rg[1]*f != 0) {
1503            ERROR("Weyl geometric left Ore failed at maximal ideal"
1504              + " induced by " + string(pVecs[i]));
1505        }
1506    }
1507    //// rational localization
1508    intvec rat = 1;
1509    f = Dx+Dy;
1510    g = x;
1511    vector rr = leftOre(g, f, 2, rat)[1];
1512    if (rr[1] == 0 || rr[2]*g-rr[1]*f != 0) {
1513        ERROR("Weyl rational left Ore failed");
1514    }
1515    // shift rational localization
1516    ring S = 0,(x,y,Sx,Sy),dp;
1517    matrix D[4][4];
1518    D[1,3] = Sx;
1519    D[2,4] = Sy;
1520    def ncS = nc_algebra(1, D);
1521    setring ncS;
1522    rat = 1;
1523    poly f = Sx+Sy;
1524    poly g = x;
1525    vector rr = leftOre(g, f, 2, rat)[1];
1526    if (rr[1] == 0 || rr[2]*g-rr[1]*f != 0) {
1527        ERROR("shift rational left Ore failed");
1528    }
1529    // q-shift rational localization
1530    ring Q = (0,q),(x,y,Qx,Qy),dp;
1531    matrix C[4][4] = UpOneMatrix(4);
1532    C[1,3] = q;
1533    C[2,4] = q;
1534    def ncQ = nc_algebra(C, 0);
1535    setring ncQ;
1536    rat = 1;
1537    poly f = Qx+Qy;
1538    poly g = x;
1539    vector rr = leftOre(g, f, 2, rat)[1];
1540    if (rr[1] == 0 || rr[2]*g-rr[1]*f != 0) {
1541        ERROR("q-shift rational left Ore failed");
1542    }
1543    print("    leftOre OK");
1544}
1545//////////////////////////////////////////////////////////////////////
1546static proc testRightOre() {
1547    print("  testing rightOre...");
1548    // Weyl
1549    ring W = 0,(x,y,Dx,Dy),dp;
1550    def ncW = Weyl();
1551    setring ncW;
1552    //// monoidal localization
1553    poly g1 = x+3;
1554    poly g2 = x*y;
1555    list L = g1,g2;
1556    poly g = x;
1557    poly f = Dx;
1558    vector rm = rightOre(g, f, 0, L)[1];
1559    if (rm[1] == 0 || f*rm[1]-g*rm[2] != 0) {
1560        ERROR("Weyl monoidal right Ore failed");
1561    }
1562    //// geometric localization
1563    g = x+y;
1564    f = Dx+Dy;
1565    ideal p = x-1,y-3;
1566    vector rg = rightOre(g, f, 1, p)[1];
1567    if (rg[1] == 0 || f*rg[1]-g*rg[2] != 0) {
1568        ERROR("Weyl geometric right Ore failed");
1569    }
1570    //// rational localization
1571    intvec rat = 1;
1572    f = Dx+Dy;
1573    g = x;
1574    vector rr = rightOre(g, f, 2, rat)[1];
1575    if (rr[1] == 0 || f*rr[1]-g*rr[2] != 0) {
1576        ERROR("Weyl rational right Ore failed");
1577    }
1578    // shift rational localization
1579    ring S = 0,(x,y,Sx,Sy),dp;
1580    matrix D[4][4];
1581    D[1,3] = Sx;
1582    D[2,4] = Sy;
1583    def ncS = nc_algebra(1, D);
1584    setring ncS;
1585    rat = 1;
1586    poly f = Sx+Sy;
1587    poly g = x;
1588    vector rr = rightOre(g, f, 2, rat)[1];
1589    if (rr[1] == 0 || f*rr[1]-g*rr[2] != 0) {
1590        ERROR("shift rational right Ore failed");
1591    }
1592    ring Q = (0,q),(x,y,Qx,Qy),dp;
1593    matrix C[4][4] = UpOneMatrix(4);
1594    C[1,3] = q;
1595    C[2,4] = q;
1596    def ncQ = nc_algebra(C, 0);
1597    setring ncQ;
1598    rat = 1;
1599    poly f = Qx+Qy;
1600    poly g = x;
1601    vector rr = rightOre(g, f, 2, rat)[1];
1602    if (rr[1] == 0 || f*rr[1]-g*rr[2] != 0) {
1603        ERROR("q-shift rational right Ore failed");
1604    }
1605    print("    rightOre OK");
1606}
1607//////////////////////////////////////////////////////////////////////
1608static proc testAddLeftFractions()
1609{
1610    print("  testing addLeftFractions...");
1611    ring R = 0,(x,y,Dx,Dy),dp;
1612    def S = Weyl();
1613    setring S;
1614    // monoidal localization
1615    poly g1 = x+3;
1616    poly g2 = x*y+y;
1617    list L = g1,g2;
1618    vector frac1 = [g1,Dx,0,0];
1619    vector frac2 = [g2,Dy,0,0];
1620    vector rm = addLeftFractions(frac1, frac2, 0, L);
1621    if (rm[1] != x^2*y+4*x*y+3*y || rm[2] != x*y*Dx+y*Dx+x*Dy+3*Dy) {
1622        ERROR("Weyl monoidal addition failed");
1623    }
1624    // geometric localization
1625    ideal p = x-1,y-3;
1626    vector rg = addLeftFractions(frac1, frac2, 1, p);
1627    if (rg[1] != x^2*y+4*x*y+3*y || rg[2] != x*y*Dx+y*Dx+x*Dy+3*Dy) {
1628        ERROR("Weyl geometric addition failed");
1629    }
1630    // rational localization
1631    intvec v = 2;
1632    frac1 = [y^2+y+1,Dx,0,0];
1633    frac2 = [y-2,Dy,0,0];
1634    vector rr = addLeftFractions(frac1, frac2, 2, v);
1635    if (rr[1] != y^3-y^2-y-2 || rr[2] != y^2*Dy+y*Dx+y*Dy-2*Dx+Dy) {
1636        ERROR("Weyl rational addition failed");
1637    }
1638    print("    addLeftFractions OK");
1639}
1640//////////////////////////////////////////////////////////////////////
1641static proc testMultiplyLeftFractions() {
1642    print("  testing multiplyLeftFractions...");
1643    ring R = 0,(x,y,Dx,Dy),dp;
1644    def S = Weyl();
1645    setring S;
1646    // monoidal localization
1647    poly g1 = x+3;
1648    poly g2 = x*y+y;
1649    list L = g1,g2;
1650    vector frac1 = [g1,Dx,0,0];
1651    vector frac2 = [g2,Dy,0,0];
1652    vector rm = multiplyLeftFractions(frac1, frac2, 0, L);
1653    if (rm[1] != g1*g2^2 || rm[2] != x*y*Dx*Dy+y*Dx*Dy-y*Dy) {
1654        ERROR("Weyl monoidal multiplication error");
1655    }
1656    // geometric localization
1657    ideal p = x-1,y-3;
1658    vector rg = multiplyLeftFractions(frac1, frac2, 1, p);
1659    if (rg[1] != g1*g2*(x+1) || rg[2] != x*Dx*Dy+Dx*Dy-Dy) {
1660        ERROR("Weyl geometric multiplication error");
1661    }
1662    // rational localization
1663    intvec v = 2;
1664    frac1 = [y^2+y+1,Dx,0,0];
1665    frac2 = [y-2,Dy,0,0];
1666    vector rr = multiplyLeftFractions(frac1, frac2, 2, v);
1667    if (rr[1] != (y^2+y+1)*(y-2) || rr[2] != Dx*Dy) {
1668        ERROR("Weyl rational multiplication (1*2) error");
1669    }
1670    rr = multiplyLeftFractions(frac2, frac1, 2, v);
1671    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) {
1672        ERROR("Weyl rational multiplication (2*1) error");
1673    }
1674    print("    multiplyLeftFractions OK");
1675}
1676//////////////////////////////////////////////////////////////////////
1677static proc testAreEqualLeftFractions() {
1678    print("  testing areEqualLeftFractions...");
1679    ring R = 0,(x,y,Dx,Dy),dp;
1680    def S = Weyl();
1681    setring S;
1682    // monoidal
1683    poly g1 = x*y+3;
1684    poly g2 = y^3;
1685    list L = g1,g2;
1686    vector fracm1 = [g1,Dx,0,0];
1687    vector fracm2 = [g1*g2,g2*Dx,0,0];
1688    vector fracm3 = [g1*g2,g1*Dx+3,0,0];
1689    if (!areEqualLeftFractions(fracm1, fracm2, 0, L)) {
1690        ERROR("Weyl monoidal positive basic comparison error");
1691    }
1692    if (areEqualLeftFractions(fracm1, fracm3, 0, L)) {
1693        ERROR("Weyl monoidal first negative basic comparison error");
1694    }
1695    if (areEqualLeftFractions(fracm2, fracm3, 0, L)) {
1696        ERROR("Weyl monoidal second negative basic comparison error");
1697    }
1698    // geometric
1699    ideal p = x+5, y-2;
1700    vector fracg1 = [g1,Dx,0,0];
1701    vector fracg2 = [g1*g2,g2*Dx,0,0];
1702    vector fracg3 = [g1*g2,g1*Dx+3,0,0];
1703    if (!areEqualLeftFractions(fracg1, fracg2, 1, p)) {
1704        ERROR("Weyl geometric positive basic comparison error");
1705    }
1706    if (areEqualLeftFractions(fracg1, fracg3, 1, p)) {
1707        ERROR("Weyl geometric first negative basic comparison error");
1708    }
1709    if (areEqualLeftFractions(fracg2, fracg3, 1, p)) {
1710        ERROR("Weyl geometric second negative basic comparison error");
1711    }
1712    // rational
1713    intvec rat = 1,4;
1714    vector fracr1 = [x+Dy,Dx,0,0];
1715    vector fracr2 = [x*Dy*(x+Dy),x*Dx*Dy,0,0];
1716    vector fracr3 = [Dy*x*(x+Dy),x*Dx*Dy+1,0,0];
1717    if (!areEqualLeftFractions(fracr1, fracr2, 2, rat)) {
1718        ERROR("Weyl rational positive basic comparison error");
1719    }
1720    if (areEqualLeftFractions(fracr1, fracr3, 2, rat)) {
1721        ERROR("Weyl rational first negative basic comparison error");
1722    }
1723    if (areEqualLeftFractions(fracr2, fracr3, 2, rat)) {
1724        ERROR("Weyl rational second negative basic comparison error");
1725    }
1726    print("    areEqualLeftFractions OK");
1727}
1728//////////////////////////////////////////////////////////////////////
1729static proc testConvertLeftToRightFraction() {
1730    print("  testing convertLeftToRightFraction...");
1731    // Weyl
1732    ring W = 0,(x,y,Dx,Dy),dp;
1733    def ncW = Weyl();
1734    setring ncW;
1735    //// monoidal localization
1736    vector fracm = [x,Dx,0,0];
1737    list L = x;
1738    vector rm = convertLeftToRightFraction(fracm, 0, L);
1739    if (!fracStatus(rm, 0, L)[1]) {
1740        ERROR("Weyl monoidal convertLeftToRightFraction failed");
1741    }
1742    //// geometric localization
1743    vector fracg = [x+y,Dx+Dy,0,0];
1744    ideal p = x-1,y-3;
1745    vector rg = convertLeftToRightFraction(fracg, 1, p);
1746    if (!fracStatus(rg, 1, p)[1]) {
1747        ERROR("Weyl geometric convertLeftToRightFraction failed");
1748    }
1749    //// rational localization
1750    intvec rat = 1;
1751    vector fracr = [x,Dx+Dy,0,0];
1752    vector rr = convertLeftToRightFraction(fracr, 2, rat);
1753    if (!fracStatus(rr, 2, rat)[1]) {
1754        ERROR("Weyl rational convertLeftToRightFraction failed");
1755    }
1756    // shift rational localization
1757    ring S = 0,(x,y,Sx,Sy),dp;
1758    matrix D[4][4];
1759    D[1,3] = Sx;
1760    D[2,4] = Sy;
1761    def ncS = nc_algebra(1, D);
1762    setring ncS;
1763    vector fracr = [x,Sx+Sy,0,0];
1764    vector rr = convertLeftToRightFraction(fracr, 2, rat);
1765    if (!fracStatus(rr, 2, rat)[1]) {
1766        ERROR("Shift rational convertLeftToRightFraction failed");
1767    }
1768    // q-shift rational localization
1769    ring Q = (0,q),(x,y,Qx,Qy),dp;
1770    matrix C[4][4] = UpOneMatrix(4);
1771    C[1,3] = q;
1772    C[2,4] = q;
1773    def ncQ = nc_algebra(C, 0);
1774    setring ncQ;
1775    vector fracr = [x,Qx+Qy,0,0];
1776    vector rr = convertLeftToRightFraction(fracr, 2, rat);
1777    if (!fracStatus(rr, 2, rat)[1]) {
1778        ERROR("q-shift rational convertLeftToRightFraction failed");
1779    }
1780    print("    convertLeftToRightFraction OK");
1781}
1782//////////////////////////////////////////////////////////////////////
1783static proc testConvertRightToLeftFraction() {
1784    print("  testing convertRightToLeftFraction...");
1785    // Weyl
1786    ring W = 0,(x,y,Dx,Dy),dp;
1787    def ncW = Weyl();
1788    setring ncW;
1789    //// monoidal localization
1790    poly g1 = x+3;
1791    poly g2 = x*y;
1792    list L = g1,g2;
1793    vector fracm = [0,0,Dx,g1^2*g2];
1794    vector rm = convertRightToLeftFraction(fracm, 0, L);
1795    if (!fracStatus(rm, 0, L)[1]) {
1796        ERROR("Weyl monoidal convertRightToLeftFraction failed");
1797    }
1798    //// geometric localization
1799    ideal p = x-1,y-3;
1800    vector fracg = [0,0,Dx,x^2+y];
1801    vector rg = convertRightToLeftFraction(fracg, 1, p);
1802    if (!fracStatus(rg, 1, p)[1]) {
1803        ERROR("Weyl geometric convertRightToLeftFraction failed");
1804    }
1805    //// rational localization
1806    intvec rat = 1;
1807    vector fracr = [0,0,Dx+Dy,x];
1808    vector rr = convertRightToLeftFraction(fracr, 2, rat);
1809    if (!fracStatus(rr, 2, rat)[1]) {
1810        ERROR("Weyl rational convertRightToLeftFraction failed");
1811    }
1812    // shift rational localization
1813    ring S = 0,(x,y,Sx,Sy),dp;
1814    matrix D[4][4];
1815    D[1,3] = Sx;
1816    D[2,4] = Sy;
1817    def ncS = nc_algebra(1, D);
1818    setring ncS;
1819    vector fracr = [0,0,Sx+Sy,x];
1820    vector rr = convertRightToLeftFraction(fracr, 2, rat);
1821    if (!fracStatus(rr, 2, rat)[1]) {
1822        ERROR("Shift rational convertRightToLeftFraction failed");
1823    }
1824    // q-shift rational localization
1825    ring Q = (0,q),(x,y,Qx,Qy),dp;
1826    matrix C[4][4] = UpOneMatrix(4);
1827    C[1,3] = q;
1828    C[2,4] = q;
1829    def ncQ = nc_algebra(C, 0);
1830    setring ncQ;
1831    vector fracr = [0,0,Qx+Qy,x];
1832    vector rr = convertRightToLeftFraction(fracr, 2, rat);
1833    if (!fracStatus(rr, 2, rat)[1]) {
1834        ERROR("q-shift rational convertRightToLeftFraction failed");
1835    }
1836    print("    convertRightToLeftFraction OK");
1837}
1838//////////////////////////////////////////////////////////////////////
1839static proc testIsInvertibleLeftFraction() {
1840    print("  testing isInvertibleLeftFraction...");
1841    ring R = 0,(x,y,Dx,Dy),dp;
1842    def S = Weyl();
1843    setring S;
1844    poly g1 = x+3;
1845    poly g2 = x*y;
1846    // monoidal
1847    list L = g1, g2;
1848    if (!isInvertibleLeftFraction([g1*g2,17,0,0], 0, L)) {
1849        ERROR("Weyl monoidal positive invertibility test error");
1850    }
1851    if (isInvertibleLeftFraction([g1*g2,Dx,0,0], 0, L)) {
1852        ERROR("Weyl monoidal negative invertibility test error");
1853    }
1854    if (!isInvertibleLeftFraction([1,1,1,1], 0, L)) {
1855        ERROR("Weyl monoidal one invertibility test error");
1856    }
1857    if (isInvertibleLeftFraction([1,0,0,1], 0, L)) {
1858        ERROR("Weyl monoidal zero invertibility test error");
1859    }
1860    // geometric
1861    ideal p = x-1, y;
1862    if (!isInvertibleLeftFraction([g1,3*x,0,0], 1, p)) {
1863        ERROR("Weyl geometric positive invertibility test error");
1864    }
1865    if (isInvertibleLeftFraction([g1,Dx,0,0], 0, L)) {
1866        ERROR("Weyl geometric negative invertibility test error");
1867    }
1868    if (!isInvertibleLeftFraction([1,1,1,1], 1, p)) {
1869        ERROR("Weyl geometric one invertibility test error");
1870    }
1871    if (isInvertibleLeftFraction([1,0,0,1], 1, p)) {
1872        ERROR("Weyl geometric zero invertibility test error");
1873    }
1874    // rational
1875    intvec rat = 1,2;
1876    if (!isInvertibleLeftFraction([g1*g2,y,0,0], 2, rat)) {
1877        ERROR("Weyl rational positive invertibility test error");
1878    }
1879    if (isInvertibleLeftFraction([g1*g2,Dx,0,0], 0, L)) {
1880        ERROR("Weyl rational negative invertibility test error");
1881    }
1882    if (!isInvertibleLeftFraction([1,1,1,1], 2, rat)) {
1883        ERROR("Weyl rational one invertibility test error");
1884    }
1885    if (isInvertibleLeftFraction([1,0,0,1], 2, rat)) {
1886        ERROR("Weyl rational zero invertibility test error");
1887    }
1888    print("    isInvertibleLeftFraction OK");
1889}
1890//////////////////////////////////////////////////////////////////////
1891static proc testInvertLeftFraction() {
1892    print("  testing invertLeftFraction...");
1893    ring R = 0,(x,y,Dx,Dy),dp;
1894    def S = Weyl();
1895    setring S;
1896    poly g1 = x+3;
1897    poly g2 = x*y;
1898    // monoidal
1899    list L = g1, g2;
1900    vector rm = [g1*g2, 17, 0, 0];
1901    vector rmInv = invertLeftFraction(rm, 0 , L);
1902    if (!isOneFraction(multiplyLeftFractions(rm, rmInv, 0, L))) {
1903        ERROR("Weyl monoidal inversion error");
1904    }
1905    // geometric
1906    ideal p = x-1, y;
1907    vector rg = [g1, 3*x, 0, 0];
1908    vector rgInv = invertLeftFraction(rg, 1, p);
1909    if (!isOneFraction(multiplyLeftFractions(rg, rgInv, 1, p))) {
1910        ERROR("Weyl geometric inversion error");
1911    }
1912    // rational
1913    intvec rat = 1,2;
1914    vector rr = [g1*g2, y, 0, 0];
1915    vector rrInv = invertLeftFraction(rr, 2, rat);
1916    if (!isOneFraction(multiplyLeftFractions(rr, rrInv, 2, rat))) {
1917        ERROR("Weyl rational inversion error");
1918    }
1919    print("    invertLeftFraction OK");
1920}
1921//////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.