source: git/Singular/LIB/divisors.lib

spielwiese
Last change on this file was 1243bc, checked in by Hans Schoenemann <hannes@…>, 3 years ago
more polymake dep.
  • Property mode set to 100644
File size: 19.3 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="version divisors.lib 4.2.0.1 Jan_2021 "; // $Id$
3category = "Commutative Algebra";
4info="
5LIBRARY:  divisors.lib      Divisors and P-Divisors
6
7AUTHORS:  Janko Boehm       boehm@mathematik.uni-kl.de
8@*        Lars Kastner      kastner@math.fu-berlin.de
9@*        Benjamin Lorenz   blorenz@math.uni-frankfurt.de
10@*        Hans Schoenemann  hannes@mathematik.uni-kl.de
11@*        Yue Ren           ren@mathematik.uni-kl.de
12
13OVERVIEW:
14
15We implement a class divisor on an algebraic variety and methods
16for computing with them. Divisors are represented by tuples of ideals
17defining the positive and the negative part. In particular, we implement the group
18structure on divisors, computing global sections and testing linear
19equivalence.
20
21In addition to this we provide a class formaldivisor which implements
22integer formal sums of divisors (not necessarily prime). A formal divisor
23can be evaluated to a divisor, and a divisor can be decomposed into
24a formal sum.
25
26Finally we provide a class pdivisor which implements polyhedral formal sums of
27divisors (P-divisors) where the coefficients are assumed to be polyhedra with fixed tail cone.
28There is a function to evaluate a P-divisor on a vector in the dual of the tail cone. The
29result will be a formal divisor.
30
31
32REFERENCES:
33
34For the class divisor we closely follow Macaulay2's tutorial on divisors.
35
36
37PROCEDURES:
38
39makeDivisor(ideal,ideal)                        create a divisor
40divisorplus(divisor,divisor)                    add two divisors
41multdivisor(int,divisor)                        multiply a divisor by an integer
42negativedivisor(divisor)                        compute the negative of the divisor
43normalForm(divisor)                             normal form of a divisor
44isEqualDivisor(divisor,divisor)                 test whether two divisors are equal
45globalSections(divisor)                         compute the global sections of a divisor
46degreeDivisor(divisor)                          degree of a divisor
47linearlyEquivalent(divisor,divisor)             test whether two divisors a linearly equivalent
48effective(divisor)                              compute an effective divisor
49                                                linearly equivalent to a divisor
50
51makeFormalDivisor(list)                         make a formal integer sum of divisors
52evaluateFormalDivisor(formaldivisor)            evaluate a formal sum of divisors to a divisor
53formaldivisorplus(formaldivisor,formaldivisor)  add two formal divisors
54negativeformaldivisor(formaldivisor)            compute the negative of the formal divisor
55multformaldivisor(int,formaldivisor)            multiply a formal divisor by an integer
56degreeFormalDivisor(formaldivisor)              degree of a formal divisor
57
58makePDivisor(List)                              make a formal polyhedral sum of divisors
59
60
61";
62
63
64////////////////////////////////////////////////////////////////////////////////
65
66LIB "primdec.lib";
67
68static proc mod_init()
69{
70LIB"customstd.so";
71LIB"gfanlib.so";
72newstruct("divisor","ideal den,ideal num");
73newstruct("formaldivisor","list summands");
74newstruct("pdivisor","list summands, cone tail");
75
76system("install","divisor","print",divisor_print,1);
77system("install","divisor","+",divisorplus,2);
78system("install","divisor","*",proxymultdivisor,2);
79system("install","formaldivisor","print",formaldivisor_print,1);
80system("install","formaldivisor","+",formaldivisorplus,2);
81system("install","formaldivisor","*",proxymultformaldivisor,2);
82}
83
84
85
86proc divisor_print(divisor D)
87"USAGE:  divisor_print(D); D; D = divisor. @*
88ASSUME:  D is a divisor.
89RETURN:  Will print D.
90KEYWORDS: divisors
91EXAMPLE: example divisor_print; shows an example
92"
93{
94"("+string(D.num)+") - ("+string(D.den)+")";
95}
96example
97{ "EXAMPLE:";
98ring r=31991,(x,y,z),dp;
99ideal I = y^2*z - x*(x-z)*(x+3*z);
100qring Q = std(I);
101divisor P = makeDivisor(ideal(x,z),ideal(1));
102P;
103}
104
105
106proc formaldivisor_print(formaldivisor fD)
107"USAGE:  formaldivisor_print(D); D; D = formaldivisor. @*
108ASSUME:  fD is a formaldivisor.
109RETURN:  Will print fD.
110KEYWORDS: divisors
111EXAMPLE: example formaldivisor_print; shows an example
112"
113{
114  int i; string s; list L=fD.summands;
115  list cd; int c; divisor d;
116  string linebreak =
117"
118";
119  for (i=1; i<=size(L); i++)
120  {
121    cd=L[i]; c=cd[1]; d=cd[2];
122    if (i>1 && c>=0) { s = s + "+"; }
123    s = s + string(c)+"*( ("+string(d.num)+") - ("+string(d.den)+") )";
124    if (i!=size(L)) { s = s + linebreak; }
125  }
126  s;
127}
128example
129{ "EXAMPLE:";
130ring r=31991,(x,y,z),dp;
131ideal I = y^2*z - x*(x-z)*(x+3*z);
132qring Q = std(I);
133divisor P = makeDivisor(ideal(x,z),ideal(1));
134P;
135}
136
137
138////////////////////////////////////////////////////////////////////////////////
139// divisors as numerator and denomiator ideals
140
141proc makeDivisor(ideal I, ideal J)
142"USAGE:  makeDivisor(I ,J); I = ideal, J = ideal. @*
143ASSUME:  I and J are ideals in a qring Q of a smooth irreducible variety X
144         such that any ideal in Q satisfies the S2 condition.
145RETURN:  a divisor on X
146THEORY:  The procedure will eliminate all components which are not of codimension 1.
147         The S2 condition requires that every proper nonzero principal ideal
148         has pure codimension 1.
149KEYWORDS: divisors
150EXAMPLE: example makeDivisor; shows an example
151"
152{
153divisor C;
154C.num = purify1(I);
155C.den = purify1(J);
156return(C);
157}
158
159example
160{ "EXAMPLE:";
161ring r=31991,(x,y,z),dp;
162ideal I = y^2*z - x*(x-z)*(x+3*z);
163qring Q = std(I);
164divisor P = makeDivisor(ideal(x,z),ideal(1));
165}
166
167
168proc divisorplus(divisor A, divisor B)
169"USAGE:  divisorplus(A ,B); A + B; A = divisor, B = divisor. @*
170ASSUME:  A and B are divisors on X.
171RETURN:  a divisor on X
172THEORY:  The procedure will compute the product of the numerator
173         and denominator ideals, respectively.
174KEYWORDS: divisors
175EXAMPLE: example divisorplus; shows an example
176"
177{
178return(makeDivisor(interred(A.num*B.num),interred(A.den*B.den)));
179}
180example
181{ "EXAMPLE:";
182ring r=31991,(x,y,z),dp;
183ideal I = y^2*z - x*(x-z)*(x+3*z);
184qring Q = std(I);
185divisor A = makeDivisor(ideal(x,z),ideal(1));
186divisor B = makeDivisor(ideal(x,y),ideal(1));
187A+B;
188}
189
190
191proc multdivisor(int n, divisor A)
192"USAGE:  multdivisor(n ,A); A*n; n = integer, A = divisor.@*
193ASSUME:  n is an integer and A is a divisor on X.
194RETURN:  a divisor on X
195THEORY:  The procedure will compute the n-th power of the numerator
196         and denominator ideals, respectively.
197KEYWORDS: divisors
198EXAMPLE: example multdivisor; shows an example
199"
200{
201if (n<0) {return(multdivisor(-n,negativedivisor(A)));}
202return(makeDivisor(interred((A.num)^n),interred((A.den)^n)));
203}
204example
205{ "EXAMPLE:";
206ring r=31991,(x,y,z),dp;
207ideal I = y^2*z - x*(x-z)*(x+3*z);
208qring Q = std(I);
209divisor A = makeDivisor(ideal(x,z),ideal(1));
210A;
211divisor D = multdivisor(4,A);
212D;
213A*4;
214}
215
216
217/***
218 * For operator overloading, which needs a procedure which takes a divisor first
219 * and integer second.
220 **/
221proc proxymultdivisor(divisor A, int n)
222{
223if (n<0) {return(multdivisor(-n,negativedivisor(A)));}
224return(makeDivisor(interred((A.num)^n),interred((A.den)^n)));
225}
226
227
228proc negativedivisor(divisor A)
229"USAGE:  negativedivisor(A); A*(-1); A = divisor.@*
230ASSUME:  A is a divisor on X.
231RETURN:  a divisor on X
232THEORY:  The procedure will interchange the numerator and denominator ideals.
233KEYWORDS: divisors
234EXAMPLE: example negativedivisor; shows an example
235"
236{
237divisor B;
238B.num=A.den;
239B.den=A.num;
240return(B);
241}
242example
243{ "EXAMPLE:";
244ring r=31991,(x,y,z),dp;
245ideal I = y^2*z - x*(x-z)*(x+3*z);
246qring Q = std(I);
247divisor A = makeDivisor(ideal(x,z),ideal(1));
248A;
249divisor D = negativedivisor(A);
250D;
251}
252
253
254proc normalForm(divisor A)
255"USAGE:  normalForm(A); A = divisor.@*
256ASSUME:  A is a divisor on X.
257RETURN:  different representative of the same divisor on X
258THEORY:  The procedure will cancel common components of numerator and denominator.
259KEYWORDS: divisors
260EXAMPLE: example normalForm; shows an example
261"
262{
263divisor B;
264B.num=quotient(A.num,A.den);
265B.den=quotient(A.den,A.num);
266return(B);
267}
268example
269{ "EXAMPLE:";
270ring r=31991,(x,y,z),dp;
271ideal I = y^2*z - x*(x-z)*(x+3*z);
272qring Q = std(I);
273divisor A = makeDivisor(ideal(x,z),ideal(1));
274divisor B = makeDivisor(ideal(x,y),ideal(1));
275divisor D = (A+B)+multdivisor(-1,B);
276D;
277normalForm(D);
278}
279
280
281
282
283static proc isEqualIdeal(ideal A,ideal B){
284return((size(NF(A,std(B)))==0) && (size(NF(B,std(A)))==0));
285}
286
287
288proc isEqualDivisor(divisor A,divisor B)
289"USAGE:  isEqualDivisor(A,B); A = divisor, B = divisor.@*
290ASSUME:  A and B are divisors on X.
291RETURN:  int 0 or 1, checks equality of A and B.
292THEORY:  The procedure will compute the normal forms of A and B and compare.
293KEYWORDS: divisors
294EXAMPLE: example isEqualDivisor; shows an example
295"
296{
297A=normalForm(A);
298B=normalForm(B);
299return((isEqualIdeal(A.num,B.num)) && (isEqualIdeal(A.den,B.den)));
300}
301example
302{ "EXAMPLE:";
303ring r=31991,(x,y,z),dp;
304ideal I = y^2*z - x*(x-z)*(x+3*z);
305qring Q = std(I);
306divisor A = makeDivisor(ideal(x,z),ideal(1));
307divisor B = makeDivisor(ideal(x,y),ideal(1));
308divisor D = (A+B)+multdivisor(-1,B);
309isEqualDivisor(A,D);
310}
311
312
313static proc purify1(ideal I)
314{
315I = simplify(I,2);
316if (I[1]==0){ERROR("expected a non-zero ideal");}
317ideal f = I[1];
318return(minbase(quotient(f,quotient(f,I))));
319}
320
321
322static proc basis(ideal I,int d)
323{
324I=simplify(jet(intersect(I,maxideal(d)),d),2);
325return(I)}
326
327//basis(ideal(x,y^3),2);
328
329
330proc globalSections(divisor D)
331"USAGE:  globalSections(A); A = divisor.@*
332ASSUME:  A is a divisor on X.
333RETURN:  a list with a basis of the space of global sections of D.
334THEORY:  We assume that the qring of X satisfies the S2-condition and that X is smooth. We compute sat((f*J) : I) /f
335         where D = (I)-(J).
336KEYWORDS: divisors
337EXAMPLE: example globalSections; shows an example
338"
339{
340poly f =(simplify(D.num,2))[1];
341ideal LD = basis(purify1(quotient(f*D.den,D.num)),deg(f));
342list L = simplify(LD,2),f;
343return(L);
344}
345example
346{ "EXAMPLE:";
347ring r=31991,(x,y,z),dp;
348ideal I = y^2*z - x*(x-z)*(x+3*z);
349qring Q = std(I);
350divisor P = makeDivisor(ideal(x,z),ideal(1));
351divisor D = multdivisor(4,P);
352globalSections(D);
353}
354
355
356static proc sectionIdeal(poly f, poly g, divisor D){
357return(purify1(quotient(quotient(f*D.num,g), D.den)));
358}
359
360proc degreeDivisor(divisor A)
361"USAGE:  degreeDivisor(A); A = divisor.@*
362ASSUME:  A is a divisor on X.
363RETURN:  The degree of A.
364THEORY:  We compute difference of the degrees of the numerator and denominator ideals.
365KEYWORDS: divisors
366EXAMPLE: example degreeDivisor; shows an example
367"
368{
369  return( mult(std(A.num))-mult(std(A.den)));
370}
371example
372{ "EXAMPLE:";
373ring r=31991,(x,y,z),dp;
374ideal I = y^2*z - x*(x-z)*(x+3*z);
375qring Q = std(I);
376divisor P = makeDivisor(ideal(x,z),ideal(1));
377degreeDivisor(P);
378}
379
380proc linearlyEquivalent(divisor A, divisor B)
381"USAGE:  linearlyEquivalent(A,B); A = divisor, B = divisor.@*
382ASSUME:  A and B are divisors on X.
383RETURN:  list if A and B a linearly equivalent and int 0 otherwise.
384THEORY:  Checks whether A-B is principle. If yes, returns a list L=(f,g) where
385         A - B = (f/g).
386KEYWORDS: divisors
387EXAMPLE: example linearlyEquivalent; shows an example
388"
389{
390divisor F = normalForm(divisorplus(A,negativedivisor(B)));
391list LB = globalSections(F);
392if (size(LB[1])!=1) {return(0);}
393ideal V = sectionIdeal(LB[1][1,1],LB[2],F);
394if (isEqualIdeal(V,ideal(1))==1) {return(list(LB[1][1,1],LB[2]));}
395return(0);
396}
397example
398{ "EXAMPLE:";
399ring r=31991,(x,y,z),dp;
400ideal I = y^2*z - x*(x-z)*(x+3*z);
401qring Q = std(I);
402divisor A = makeDivisor(ideal(x,z),ideal(1));
403divisor B = makeDivisor(ideal(x,y),ideal(1));
404linearlyEquivalent(A,B);
405linearlyEquivalent(multdivisor(2,A),multdivisor(2,B));
406}
407
408
409proc effective(divisor A)
410"USAGE:  effective(A); A = divisor.@*
411ASSUME:  A is a divisor on X which is linearly equivalent to an effective divisor.
412RETURN:  divisor on X.
413THEORY:  We compute an effective divisor linearly equivalent to A.
414KEYWORDS: divisors
415EXAMPLE: example effective; shows an example
416"
417{
418list LB = globalSections(A);
419if (size(LB[1])==0) {ERROR("the divisor is not linearly equivalent to an effective divisor");}
420ideal V = sectionIdeal(LB[1][1,1],LB[2],A);
421return(makeDivisor(V,ideal(1)));
422}
423example
424{ "EXAMPLE:";
425ring r=31991,(x,y,z),dp;
426ideal I = y^2*z - x*(x-z)*(x+3*z);
427qring Q = std(I);
428divisor A = makeDivisor(ideal(x,z),ideal(1));
429divisor B = makeDivisor(ideal(x,y),ideal(1));
430divisor D = divisorplus(multdivisor(2,B),negativedivisor(A));
431effective(D);
432}
433
434
435
436////////////////////////////////////////////////////////////////////////////////
437// formal sums of divisors
438
439proc makeFormalDivisor(list L)
440"USAGE:  makeFormalDivisor(L); L = list.@*
441ASSUME:  L is a list of tuples of an integer and a divisor.
442RETURN:  a formal divisor on X
443THEORY:  Represents an integer formal sum of divisors.
444KEYWORDS: divisors
445EXAMPLE: example makeFormalDivisor; shows an example
446"
447{
448formaldivisor C;
449C.summands = L;
450return(C);
451}
452example
453{ "EXAMPLE:";
454ring r=31991,(x,y,z),dp;
455ideal I = y^2*z - x*(x-z)*(x+3*z);
456qring Q = std(I);
457divisor A = makeDivisor(ideal(x,z),ideal(1));
458divisor B = makeDivisor(ideal(x,y),ideal(1));
459makeFormalDivisor(list(list(-5,A),list(2,B)));
460}
461
462
463proc evaluateFormalDivisor(formaldivisor D)
464"USAGE:  evaluateFormalDivisor(D); D = formal divisor.@*
465ASSUME:  D is a formal divisor on X.
466RETURN:  a divisor on X
467THEORY:  Will evaluate the formal sum.
468KEYWORDS: divisors
469EXAMPLE: example evaluateFormalDivisor; shows an example
470"
471{
472list L = D.summands;
473if (size(L)==0) {return(makeDivisor(ideal(1),ideal(1)));}
474int i;
475divisor E = multdivisor(L[1][1],L[1][2]);
476  for ( i=2; i <= size(L); i++ )
477  {
478    E = divisorplus(E, multdivisor(L[i][1],L[i][2]));
479  }
480return(E);
481}
482example
483{ "EXAMPLE:";
484ring r=31991,(x,y,z),dp;
485ideal I = y^2*z - x*(x-z)*(x+3*z);
486qring Q = std(I);
487divisor A = makeDivisor(ideal(x,z),ideal(1));
488divisor B = makeDivisor(ideal(x,y),ideal(1));
489formaldivisor fE= makeFormalDivisor(list(list(-5,A),list(2,B)));
490evaluateFormalDivisor(fE);
491}
492
493
494
495static proc position(divisor I,list L){
496int i;
497for (i=1; i <=size(L); i++){
498   if (isEqualDivisor(I,L[i][2])==1) {return(i);}
499}
500return(0);}
501
502
503proc formaldivisorplus(formaldivisor A, formaldivisor B)
504"USAGE:  formaldivisorplus(A ,B); A + B; A = formaldivisor, B = formaldivisor. @*
505ASSUME:  A and B are formal divisors on X.
506RETURN:  a formal divisor on X
507THEORY:  The procedure will add the formal sums.
508KEYWORDS: divisors
509EXAMPLE: example formaldivisorplus; shows an example
510"
511{
512formaldivisor C;
513int i,p;
514int j=1;
515list L;
516list LA=A.summands;
517list LB=B.summands;
518for (i=1; i<=size(LA);i++){
519   p=position(LA[i][2],L);
520   if (p==0) {
521        L[j]=list();
522        L[j][2]=LA[i][2];
523        L[j][1]=LA[i][1];
524        j=j+1;
525   } else {
526        L[p][1]=L[p][1]+LA[i][1];
527   };
528}
529for (i=1; i<=size(LB);i++){
530   p=position(LB[i][2],L);
531   if (p==0) {
532        L[j]=list();
533        L[j][2]=LB[i][2];
534        L[j][1]=LB[i][1];
535        j=j+1;
536   } else {
537        L[p][1]=L[p][1]+LB[i][1];
538   };
539}
540//C.summands = (A.summands)+(B.summands);
541return(L);
542}
543example
544{ "EXAMPLE:";
545ring r=31991,(x,y,z),dp;
546ideal I = y^2*z - x*(x-z)*(x+3*z);
547qring Q = std(I);
548divisor A = makeDivisor(ideal(x,z),ideal(1));
549divisor B = makeDivisor(ideal(x,y),ideal(1));
550divisor C = makeDivisor(ideal(x-z,y),ideal(1));
551formaldivisor fE= makeFormalDivisor(list(list(-5,A),list(2,B)));
552formaldivisor fE2= makeFormalDivisor(list(list(-5,A),list(2,C)));
553formaldivisorplus(fE,fE2);
554}
555
556
557proc degreeFormalDivisor(formaldivisor A)
558"USAGE:  degreeFormalDivisor(A); A = formaldivisor.@*
559ASSUME:  A is a formaldivisor on X.
560RETURN:  The degree of A.
561THEORY:  We compute degrees of the summands and return the weighted sum.
562KEYWORDS: divisors
563EXAMPLE: example degreeFormalDivisor; shows an example
564"
565{
566int i,s;
567list L = A.summands;
568for (i=1;i<=size(L);i++){
569    s=s+L[i][1]*degreeDivisor(L[i][2]);
570}
571return(s);
572}
573example
574{ "EXAMPLE:";
575ring r=31991,(x,y,z),dp;
576ideal I = y^2*z - x*(x-z)*(x+3*z);
577qring Q = std(I);
578divisor A = makeDivisor(ideal(x,z),ideal(1));
579divisor B = makeDivisor(ideal(x,y),ideal(1));
580formaldivisor fE= makeFormalDivisor(list(list(-5,A),list(2,B)));
581degreeFormalDivisor(fE);
582}
583
584
585proc multformaldivisor(int n,formaldivisor A)
586"USAGE:  multformaldivisor(n ,A); A*n; n = integer, A = formaldivisor.@*
587ASSUME:  n is an integer and A is a formal divisor on X.
588RETURN:  a formal divisor on X
589THEORY:  The procedure will multiply the formal sum with n.
590KEYWORDS: divisors
591EXAMPLE: example multformaldivisor; shows an example
592"
593{
594formaldivisor B;
595list L=A.summands;
596int i;
597for (i=1;i<=size(L);i++){
598  L[i][1]=n*L[i][1];
599}
600B.summands=L;
601return(B);
602}
603example
604{ "EXAMPLE:";
605ring r=31991,(x,y,z),dp;
606ideal I = y^2*z - x*(x-z)*(x+3*z);
607qring Q = std(I);
608divisor A = makeDivisor(ideal(x,z),ideal(1));
609divisor B = makeDivisor(ideal(x,y),ideal(1));
610formaldivisor fE= makeFormalDivisor(list(list(-5,A),list(2,B)));
611fE*2;
612}
613
614
615/***
616 * For operator overloading, which needs a procedure which takes a divisor first
617 * and integer second.
618 **/
619proc proxymultformaldivisor(formaldivisor A, int n)
620{
621formaldivisor B;
622list L=A.summands;
623int i;
624for (i=1;i<=size(L);i++){
625  L[i][1]=n*L[i][1];
626}
627B.summands=L;
628return(B);
629}
630
631
632
633proc negativeformaldivisor(formaldivisor A)
634"USAGE:  negativeformaldivisor(A); A = formaldivisor.@*
635ASSUME:  A is a formaldivisor on X.
636RETURN:  a formal divisor on X
637THEORY:  The procedure will change the signs of the coefficients.
638KEYWORDS: divisors
639EXAMPLE: example negativeformaldivisor; shows an example
640"
641{
642formaldivisor B;
643list L=A.summands;
644int i;
645for (i=1;i<=size(L);i++){
646  L[i][1]=-L[i][1];
647}
648B.summands=L;
649return(B);
650}
651example
652{ "EXAMPLE:";
653ring r=31991,(x,y,z),dp;
654ideal I = y^2*z - x*(x-z)*(x+3*z);
655qring Q = std(I);
656divisor A = makeDivisor(ideal(x,z),ideal(1));
657divisor B = makeDivisor(ideal(x,y),ideal(1));
658formaldivisor fE= makeFormalDivisor(list(list(-5,A),list(2,B)));
659negativeformaldivisor(fE);
660}
661
662
663static proc primDecDivisor(divisor D)
664"decompose a divisor into a formal divisor of primary divisors"
665{
666formaldivisor E;
667ideal I = D.num;
668ideal J = D.den;
669list L;
670int i;
671int j = 1;
672list LI = primdecGTZ(I);
673for (i=1;i<=size(LI);i++){
674   LI[i][2];
675   L[j]=list(1,makeDivisor(LI[i][1],ideal(1)));
676   j=j+1;
677};
678list LJ = primdecGTZ(J);
679for (i=1;i<=size(LJ);i++){
680   LJ[i][2];
681   L[j]=list(-1,makeDivisor(LJ[i][1],ideal(1)));
682   j=j+1;
683};
684E.summands=L;
685return(E);}
686
687
688
689
690////////////////////////////////////////////////////////////////////////////////
691// P-divisors
692
693proc makePDivisor(list L)
694"USAGE:  makePDivisor(L); L = list.@*
695ASSUME:  L is a list of tuples of a integral polyhedron and a divisor such that
696         all polyhedra have the same tail cone.
697RETURN:  a pdivisor on X
698THEORY:  Represents an polyhedral formal sum of divisors.
699KEYWORDS: divisors, polyhedra
700EXAMPLE: example makePDivisor; shows an example
701"
702{
703pdivisor P;
704list CP = decomposePolyhedron(L[1][1]);
705P.tail = CP[1];
706list LP;
707LP[1]=list(CP[2],L[1][2]);
708int i;
709for (i=2; i<=size(L);i++){
710   CP = decomposePolyhedron(L[i][1]);
711   if (!(CP[1]==P.tail)) {ERROR("All P-coefficients should have the same tail cone");}
712   LP[i]=list(CP[2],L[i][2]);
713}
714P.summands = LP;
715return(P);
716}
717example
718{ "EXAMPLE:";
719ring r=31991,(x,y,z),dp;
720ideal I = y^2*z - x*(x-z)*(x+3*z);
721qring Q = std(I);
722divisor A = makeDivisor(ideal(x,z),ideal(1));
723divisor B = makeDivisor(ideal(x,y),ideal(1));
724intmat M[4][4]= 1,4,0,0,
725                1,0,3,0,
726                0,0,0,2,
727                1,1,1,1;
728polytope PP = polytopeViaPoints(M);
729makePDivisor(list(list(PP,A),list(PP,B)));
730}
731
732static proc decomposePolyhedron(polytope P){
733bigintmat rays = vertices(P);
734bigintmat rays2 = rays;
735int i,j;
736for (i=1; i<=nrows(rays);i++){
737   if (rays[i,1]==1) {
738      for (j=1; j<=ncols(rays);j++){
739         rays[i,j]=0;
740      }
741   } else {
742      for (j=1; j<=ncols(rays);j++){
743         rays2[i,j]=0;
744      }
745   }
746}
747cone C = coneViaPoints(rays);
748polytope C2 = polytopeViaPoints(rays2);
749return(list(C,C2));
750}
751
752
753
Note: See TracBrowser for help on using the repository browser.