source: git/Singular/LIB/olga.lib @ 44168a

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