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

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